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/bcc.py
===================================================================
--- sasmodels/models/bcc.py	(revision 2f0c07deb2cb7a8c0c93476fbf1e0ffbe6c278cb)
+++ sasmodels/models/bcc.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -124,5 +124,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/J1.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc.c"]
+source = ["lib/sph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc.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/core_shell_bicelle.c
===================================================================
--- sasmodels/models/core_shell_bicelle.c	(revision cf85329d4eadae1d4c59b4a1c0b6f34770d49de5)
+++ sasmodels/models/core_shell_bicelle.c	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -56,6 +56,6 @@
     double sinarg2 = qq*(length+facthick)*cn;
 
-    be1 = J1c(besarg1);
-    be2 = J1c(besarg2);
+    be1 = sas_J1c(besarg1);
+    be2 = sas_J1c(besarg2);
     si1 = sinc(sinarg1);
     si2 = sinc(sinarg2);
Index: sasmodels/models/core_shell_bicelle.py
===================================================================
--- sasmodels/models/core_shell_bicelle.py	(revision 2f0c07deb2cb7a8c0c93476fbf1e0ffbe6c278cb)
+++ sasmodels/models/core_shell_bicelle.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -87,5 +87,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/Si.c", "lib/J1.c", "lib/J1c.c", "lib/gauss76.c", "core_shell_bicelle.c"]
+source = ["lib/Si.c","lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "core_shell_bicelle.c"]
 
 demo = dict(scale=1, background=0,
Index: sasmodels/models/core_shell_cylinder.c
===================================================================
--- sasmodels/models/core_shell_cylinder.c	(revision f4878dc778e82ee49c11c59cc0517797917067cc)
+++ sasmodels/models/core_shell_cylinder.c	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -11,5 +11,5 @@
 double _cyl(double twovd, double besarg, double siarg)
 {
-    const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg);
+    const double bj = (besarg == 0.0 ? 0.5 : 0.5*sas_J1c(besarg));
     const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg);
     return twovd*si*bj;
Index: sasmodels/models/core_shell_cylinder.py
===================================================================
--- sasmodels/models/core_shell_cylinder.py	(revision 01eece695fca4715baea58b0e86d78de1ed9ca7c)
+++ sasmodels/models/core_shell_cylinder.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -130,5 +130,5 @@
              ]
 
-source = ["lib/J1.c", "lib/gauss76.c", "core_shell_cylinder.c"]
+source = ["lib/polevl.c","lib/sas_J1.c", "lib/gauss76.c", "core_shell_cylinder.c"]
 
 def ER(radius, thickness, length):
Index: sasmodels/models/core_shell_parallelepiped.py
===================================================================
--- sasmodels/models/core_shell_parallelepiped.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ sasmodels/models/core_shell_parallelepiped.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -141,5 +141,5 @@
              ]
 
-source = ["lib/J1.c", "lib/gauss76.c", "core_shell_parallelepiped.c"]
+source = ["lib/gauss76.c", "core_shell_parallelepiped.c"]
 
 
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/ellipsoid.py
===================================================================
--- sasmodels/models/ellipsoid.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ sasmodels/models/ellipsoid.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -132,5 +132,5 @@
              ]
 
-source = ["lib/J1.c", "lib/sph_j1c.c", "lib/gauss76.c", "ellipsoid.c"]
+source = ["lib/sph_j1c.c", "lib/gauss76.c", "ellipsoid.c"]
 
 def ER(rpolar, requatorial):
Index: sasmodels/models/elliptical_cylinder.c
===================================================================
--- sasmodels/models/elliptical_cylinder.c	(revision a8b3cdb60df4f2e22c9b3c5c3110bb51b7dba2c6)
+++ sasmodels/models/elliptical_cylinder.c	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -22,5 +22,6 @@
         retval = 1.0;
     }else{
-        retval = 2.0*NR_BessJ1(arg)/arg;
+        //retval = 2.0*NR_BessJ1(arg)/arg;
+        retval = sas_J1c(arg);
     }
     return retval*retval ;
@@ -155,5 +156,6 @@
       Be = 0.5;
     }else{
-      Be = NR_BessJ1(qr)/qr;
+      //Be = NR_BessJ1(qr)/qr;
+      Be = 0.5*sas_J1c(qr);
     }
 
Index: sasmodels/models/elliptical_cylinder.py
===================================================================
--- sasmodels/models/elliptical_cylinder.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ sasmodels/models/elliptical_cylinder.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -105,5 +105,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/nr_bess_j1.c", "lib/gauss76.c", "lib/gauss20.c", "elliptical_cylinder.c"]
+source = ["lib/polevl.c","lib/sas_J1.c", "lib/gauss76.c", "lib/gauss20.c", "elliptical_cylinder.c"]
 
 demo = dict(scale=1, background=0, r_minor=100, r_ratio=1.5, length=400.0,
Index: sasmodels/models/fcc.py
===================================================================
--- sasmodels/models/fcc.py	(revision 2f0c07deb2cb7a8c0c93476fbf1e0ffbe6c278cb)
+++ sasmodels/models/fcc.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -116,5 +116,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/J1.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc.c"]
+source = ["lib/sph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc.c"]
 
 # parameters for demo
Index: sasmodels/models/flexible_cylinder.c
===================================================================
--- sasmodels/models/flexible_cylinder.c	(revision e7678b2f993854334aca4c45e1c7022b2fdbb509)
+++ sasmodels/models/flexible_cylinder.c	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -25,5 +25,5 @@
     const double qr = q*radius;
     //const double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr);
-    const double crossSect = J1c(qr);
+    const double crossSect = sas_J1c(qr);
     double flex = Sk_WR(q,length,kuhn_length);
     flex *= crossSect*crossSect;
Index: sasmodels/models/flexible_cylinder.py
===================================================================
--- sasmodels/models/flexible_cylinder.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ sasmodels/models/flexible_cylinder.py	(revision 26141cba42ea53a3ca3b54e60b579f5525dd583d)
@@ -84,5 +84,5 @@
     ]
 # pylint: enable=bad-whitespace, line-too-long
-source = ["lib/J1.c", "lib/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: sasmodels/models/flexible_cylinder_ex.c
===================================================================
--- sasmodels/models/flexible_cylinder_ex.c	(revision e7678b2f993854334aca4c45e1c7022b2fdbb509)
+++ sasmodels/models/flexible_cylinder_ex.c	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -24,5 +24,5 @@
         SINCOS(zi, sn, cn);
         double arg = q*sqrt(a*a*sn*sn+b*b*cn*cn);
-        double yyy = pow((double)J1c(arg),2);
+        double yyy = pow((double)sas_J1c(arg),2);
         yyy *= Gauss76Wt[i];
         summ += yyy;
Index: sasmodels/models/flexible_cylinder_ex.py
===================================================================
--- sasmodels/models/flexible_cylinder_ex.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ sasmodels/models/flexible_cylinder_ex.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -109,5 +109,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/J1.c", "lib/J1c.c", "lib/gauss76.c", "lib/wrc_cyl.c", "flexible_cylinder_ex.c"]
+source = ["lib/polevl.c","lib/sas_J1.c", "lib/gauss76.c", "lib/wrc_cyl.c", "flexible_cylinder_ex.c"]
 
 demo = dict(scale=1.0, background=0.0001,
Index: sasmodels/models/hollow_cylinder.c
===================================================================
--- sasmodels/models/hollow_cylinder.c	(revision 0420af743cbd4738603a3f63dde3bffd78ebfdcd)
+++ sasmodels/models/hollow_cylinder.c	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -24,10 +24,10 @@
     	lam1 = 1.0;
     }else{
-    	lam1 = 2.0*J1(arg1)/arg1;
+    	lam1 = sas_J1c(arg1);
     }
     if (arg2 == 0.0){
     	lam2 = 1.0;
     }else{
-    	lam2 = 2.0*J1(arg2)/arg2;
+    	lam2 = sas_J1c(arg2);
     }
     //Todo: Need to check psi behavior as gamma goes to 1.
Index: sasmodels/models/hollow_cylinder.py
===================================================================
--- sasmodels/models/hollow_cylinder.py	(revision 01eece695fca4715baea58b0e86d78de1ed9ca7c)
+++ sasmodels/models/hollow_cylinder.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -79,5 +79,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/J1.c", "lib/gauss76.c", "hollow_cylinder.c"]
+source = ["lib/polevl.c","lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"]
 
 # pylint: disable=W0613
Index: sasmodels/models/hollow_rectangular_prism.py
===================================================================
--- sasmodels/models/hollow_rectangular_prism.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ sasmodels/models/hollow_rectangular_prism.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -115,5 +115,5 @@
              ]
 
-source = ["lib/J1.c", "lib/gauss76.c", "hollow_rectangular_prism.c"]
+source = [ "lib/gauss76.c", "hollow_rectangular_prism.c"]
 
 def ER(a_side, b2a_ratio, c2a_ratio, thickness):
Index: sasmodels/models/hollow_rectangular_prism_infinitely_thin_walls.py
===================================================================
--- sasmodels/models/hollow_rectangular_prism_infinitely_thin_walls.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ sasmodels/models/hollow_rectangular_prism_infinitely_thin_walls.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -99,5 +99,5 @@
              ]
 
-source = ["lib/J1.c", "lib/gauss76.c", "hollow_rectangular_prism_infinitely_thin_walls.c"]
+source = ["lib/gauss76.c", "hollow_rectangular_prism_infinitely_thin_walls.c"]
 
 def ER(a_side, b2a_ratio, c2a_ratio):
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);
-}
Index: sasmodels/models/parallelepiped.py
===================================================================
--- sasmodels/models/parallelepiped.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ sasmodels/models/parallelepiped.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -193,5 +193,5 @@
              ]
 
-source = ["lib/J1.c", "lib/gauss76.c", "parallelepiped.c"]
+source = ["lib/gauss76.c", "parallelepiped.c"]
 
 def ER(a_side, b_side, c_side):
Index: sasmodels/models/rectangular_prism.py
===================================================================
--- sasmodels/models/rectangular_prism.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ sasmodels/models/rectangular_prism.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -107,5 +107,5 @@
              ]
 
-source = ["lib/J1.c", "lib/gauss76.c", "rectangular_prism.c"]
+source = ["lib/gauss76.c", "rectangular_prism.c"]
 
 def ER(a_side, b2a_ratio, c2a_ratio):
Index: sasmodels/models/stacked_disks.c
===================================================================
--- sasmodels/models/stacked_disks.c	(revision d507c3a946bb1b2103ca169efe83263b3955b43d)
+++ sasmodels/models/stacked_disks.c	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -52,6 +52,6 @@
 	const double sinarg2 = qq*(halfheight+layer_thick)*cos(zi);
 
-    const double be1 = J1(besarg1)/besarg1;
-	const double be2 = J1(besarg2)/besarg2;
+    const double be1 = sas_J1c(besarg1);
+	const double be2 = sas_J1c(besarg2);
 	const double si1 = sin(sinarg1)/sinarg1;
 	const double si2 = sin(sinarg2)/sinarg2;
@@ -62,6 +62,6 @@
 	const double totald=2.0*(layer_thick+halfheight);
 
-	const double t1 = 2.0*area*(2.0*halfheight)*dr1*(si1)*(be1);
-	const double t2 = 2.0*area*dr2*(totald*si2-2.0*halfheight*si1)*(be2);
+	const double t1 = area*(2.0*halfheight)*dr1*(si1)*(be1);
+	const double t2 = area*dr2*(totald*si2-2.0*halfheight*si1)*(be2);
 
 
Index: sasmodels/models/stacked_disks.py
===================================================================
--- sasmodels/models/stacked_disks.py	(revision 5111921153c056b59ff03411e4b7b028d238ae26)
+++ sasmodels/models/stacked_disks.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -133,5 +133,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/gauss76.c", "lib/J1.c", "stacked_disks.c"]
+source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "stacked_disks.c"]
 
 demo = dict(background=0.001,
Index: sasmodels/models/triaxial_ellipsoid.py
===================================================================
--- sasmodels/models/triaxial_ellipsoid.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ sasmodels/models/triaxial_ellipsoid.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
@@ -104,5 +104,5 @@
              ]
 
-source = ["lib/J1.c", "lib/sph_j1c.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
+source = ["lib/sph_j1c.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
 
 def ER(req_minor, req_major, rpolar):
