Index: sasmodels/models/barbell.c
===================================================================
--- sasmodels/models/barbell.c	(revision 50e1e40cf14d4ded28bf96e61a3cfaf35003b899)
+++ sasmodels/models/barbell.c	(revision 26141cba42ea53a3ca3b54e60b579f5525dd583d)
@@ -31,5 +31,5 @@
         const double t = Gauss76Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
-        const double bj = J1c(qrst*sqrt(radical));
+        const double bj = sas_J1c(qrst*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
         total += Gauss76Wt[i] * Fq;
@@ -73,5 +73,5 @@
 
         const double bell_Fq = _bell_kernel(q, h, bell_radius, half_length, sin_alpha, cos_alpha);
-        const double bj = J1c(q*radius*sin_alpha);
+        const double bj = sas_J1c(q*radius*sin_alpha);
         const double si = sinc(q*half_length*cos_alpha);
         const double cyl_Fq = M_PI*radius*radius*length*bj*si;
@@ -108,5 +108,5 @@
     SINCOS(alpha, sin_alpha, cos_alpha);
     const double bell_Fq = _bell_kernel(q, h, bell_radius, half_length, sin_alpha, cos_alpha);
-    const double bj = J1c(q*radius*sin_alpha);
+    const double bj = sas_J1c(q*radius*sin_alpha);
     const double si = sinc(q*half_length*cos_alpha);
     const double cyl_Fq = M_PI*radius*radius*length*bj*si;
Index: sasmodels/models/barbell.py
===================================================================
--- sasmodels/models/barbell.py	(revision 2f0c07deb2cb7a8c0c93476fbf1e0ffbe6c278cb)
+++ sasmodels/models/barbell.py	(revision 26141cba42ea53a3ca3b54e60b579f5525dd583d)
@@ -111,5 +111,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/J1c.c", "lib/gauss76.c", "barbell.c"]
+source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"]
 
 # parameters for demo
Index: sasmodels/models/capped_cylinder.c
===================================================================
--- sasmodels/models/capped_cylinder.c	(revision 50e1e40cf14d4ded28bf96e61a3cfaf35003b899)
+++ sasmodels/models/capped_cylinder.c	(revision 26141cba42ea53a3ca3b54e60b579f5525dd583d)
@@ -37,5 +37,5 @@
         const double t = Gauss76Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
-        const double bj = J1c(qrst*sqrt(radical));
+        const double bj = sas_J1c(qrst*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
         total += Gauss76Wt[i] * Fq;
@@ -94,5 +94,5 @@
 
         const double cap_Fq = _cap_kernel(q, h, cap_radius, half_length, sin_alpha, cos_alpha);
-        const double bj = J1c(q*radius*sin_alpha);
+        const double bj = sas_J1c(q*radius*sin_alpha);
         const double si = sinc(q*half_length*cos_alpha);
         const double cyl_Fq = M_PI*radius*radius*length*bj*si;
@@ -129,5 +129,5 @@
     SINCOS(alpha, sin_alpha, cos_alpha);
     const double cap_Fq = _cap_kernel(q, h, cap_radius, half_length, sin_alpha, cos_alpha);
-    const double bj = J1c(q*radius*sin_alpha);
+    const double bj = sas_J1c(q*radius*sin_alpha);
     const double si = sinc(q*half_length*cos_alpha);
     const double cyl_Fq = M_PI*radius*radius*length*bj*si;
Index: sasmodels/models/capped_cylinder.py
===================================================================
--- sasmodels/models/capped_cylinder.py	(revision 2f0c07deb2cb7a8c0c93476fbf1e0ffbe6c278cb)
+++ sasmodels/models/capped_cylinder.py	(revision 26141cba42ea53a3ca3b54e60b579f5525dd583d)
@@ -129,5 +129,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/J1c.c", "lib/gauss76.c", "capped_cylinder.c"]
+source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"]
 
 demo = dict(scale=1, background=0,
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision 50e1e40cf14d4ded28bf96e61a3cfaf35003b899)
+++ sasmodels/models/cylinder.c	(revision 26141cba42ea53a3ca3b54e60b579f5525dd583d)
@@ -29,5 +29,5 @@
         double sn, cn;
         SINCOS(alpha, sn, cn);
-        const double fq = sinc(qh*cn) * J1c(qr*sn);
+        const double fq = sinc(qh*cn) * sas_J1c(qr*sn);
         total += Gauss76Wt[i] * fq*fq * sn;
     }
@@ -57,5 +57,5 @@
 
     SINCOS(alpha, sn, cn);
-    const double fq = sinc(q*0.5*length*cn) * J1c(q*radius*sn);
+    const double fq = sinc(q*0.5*length*cn) * sas_J1c(q*radius*sn);
     const double s = (sld-solvent_sld) * form_volume(radius, length);
     return 1.0e-4 * square(s * fq);
Index: sasmodels/models/cylinder.py
===================================================================
--- sasmodels/models/cylinder.py	(revision 01eece695fca4715baea58b0e86d78de1ed9ca7c)
+++ sasmodels/models/cylinder.py	(revision 26141cba42ea53a3ca3b54e60b579f5525dd583d)
@@ -119,5 +119,5 @@
              ]
 
-source = ["lib/J1c.c", "lib/gauss76.c", "cylinder.c"]
+source = ["lib/polevl.c","lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"]
 
 def ER(radius, length):
Index: sasmodels/models/flexible_cylinder.py
===================================================================
--- sasmodels/models/flexible_cylinder.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
+++ sasmodels/models/flexible_cylinder.py	(revision 26141cba42ea53a3ca3b54e60b579f5525dd583d)
@@ -84,5 +84,5 @@
     ]
 # pylint: enable=bad-whitespace, line-too-long
-source = ["lib/polevl.c","lib/sas_J1c.c", "lib/wrc_cyl.c", "flexible_cylinder.c"]
+source = ["lib/polevl.c", "lib/sas_J1.c", "lib/wrc_cyl.c", "flexible_cylinder.c"]
 
 demo = dict(scale=1.0, background=0.0001,
Index: smodels/models/lib/J1.c
===================================================================
--- sasmodels/models/lib/J1.c	(revision 994d77f9b713691e41106d1f9af9661dfb2dc419)
+++ 	(revision )
@@ -1,38 +1,0 @@
-double J1(double x);
-double J1(double x)
-{
-  const double ax = fabs(x);
-  if (ax < 8.0) {
-    const double y = x*x;
-    const double ans1 = x*(72362614232.0
-              + y*(-7895059235.0
-              + y*(242396853.1
-              + y*(-2972611.439
-              + y*(15704.48260
-              + y*(-30.16036606))))));
-    const double ans2 = 144725228442.0
-              + y*(2300535178.0
-              + y*(18583304.74
-              + y*(99447.43394
-              + y*(376.9991397
-              + y))));
-    return ans1/ans2;
-  } else {
-    const double y = 64.0/(ax*ax);
-    const double xx = ax - 2.356194491;
-    const double ans1 = 1.0
-              + y*(0.183105e-2
-              + y*(-0.3516396496e-4
-              + y*(0.2457520174e-5
-              + y*-0.240337019e-6)));
-    const double ans2 = 0.04687499995
-              + y*(-0.2002690873e-3
-              + y*(0.8449199096e-5
-              + y*(-0.88228987e-6
-              + y*0.105787412e-6)));
-    double sn,cn;
-    SINCOS(xx, sn, cn);
-    const double ans = sqrt(0.636619772/ax) * (cn*ans1 - (8.0/ax)*sn*ans2);
-    return (x < 0.0) ? -ans : ans;
-  }
-}
Index: smodels/models/lib/J1c.c
===================================================================
--- sasmodels/models/lib/J1c.c	(revision e7678b2f993854334aca4c45e1c7022b2fdbb509)
+++ 	(revision )
@@ -1,44 +1,0 @@
-/**
-* Bessel function 2*J1(x)/x
-*
-* Used for low q to avoid cancellation error.
-*/
-double J1c(double x);
-double J1c(double x)
-{
-  const double ax = fabs(x);
-  if (ax < 8.0) {
-    const double y = x*x;
-    const double ans1 = 2.0*(72362614232.0
-              + y*(-7895059235.0
-              + y*(242396853.1
-              + y*(-2972611.439
-              + y*(15704.48260
-              + y*(-30.16036606))))));
-    const double ans2 = 144725228442.0
-              + y*(2300535178.0
-              + y*(18583304.74
-              + y*(99447.43394
-              + y*(376.9991397
-              + y))));
-    return ans1/ans2;
-  } else {
-    const double y = 64.0/(ax*ax);
-    const double xx = ax - 2.356194491;
-    const double ans1 = 1.0
-              + y*(0.183105e-2
-              + y*(-0.3516396496e-4
-              + y*(0.2457520174e-5
-              + y*-0.240337019e-6)));
-    const double ans2 = 0.04687499995
-              + y*(-0.2002690873e-3
-              + y*(0.8449199096e-5
-              + y*(-0.88228987e-6
-              + y*0.105787412e-6)));
-    double sn,cn;
-    SINCOS(xx, sn, cn);
-    const double ans = sqrt(0.636619772/ax)
-        * (cn*ans1 - (8.0/ax)*sn*ans2)*2.0/x;
-    return (x < 0.0) ? -ans : ans;
-  }
-}
Index: smodels/models/lib/nr_bess_j1.c
===================================================================
--- sasmodels/models/lib/nr_bess_j1.c	(revision a8b3cdb60df4f2e22c9b3c5c3110bb51b7dba2c6)
+++ 	(revision )
@@ -1,28 +1,0 @@
-double NR_BessJ1(double x);
-double NR_BessJ1(double x)
-{
-    double ax,z;
-    double xx,y,ans,ans1,ans2;
-
-    if ((ax=fabs(x)) < 8.0) {
-        y=x*x;
-        ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
-                                                  +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
-        ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
-                                               +y*(99447.43394+y*(376.9991397+y*1.0))));
-        ans=ans1/ans2;
-    } else {
-        z=8.0/ax;
-        y=z*z;
-        xx=ax-2.356194491;
-        ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
-                                   +y*(0.2457520174e-5+y*(-0.240337019e-6))));
-        ans2=0.04687499995+y*(-0.2002690873e-3
-                              +y*(0.8449199096e-5+y*(-0.88228987e-6
-                                                     +y*0.105787412e-6)));
-        ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
-        if (x < 0.0) ans = -ans;
-    }
-
-    return(ans);
-}
