Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision 7fcdc9f91a5fd5d9d1fa2a34356d607aca2be682)
+++ sasmodels/generate.py	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
@@ -139,5 +139,5 @@
     square(x) = x*x
     cube(x) = x*x*x
-    sinc(x) = sin(x)/x, with sin(0)/0 -> 1
+    sas_sinx_x(x) = sin(x)/x, with sin(0)/0 -> 1
     all double precision constants must include the decimal point
     all double declarations may be converted to half, float, or long double
Index: sasmodels/kernel_header.c
===================================================================
--- sasmodels/kernel_header.c	(revision b7e8b94b477e6d9a7e4e2666bb218d3568ee66d2)
+++ sasmodels/kernel_header.c	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
@@ -146,5 +146,5 @@
 inline double square(double x) { return x*x; }
 inline double cube(double x) { return x*x*x; }
-inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; }
+inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; }
 
 #if 1
Index: sasmodels/kernel_template.c
===================================================================
--- sasmodels/kernel_template.c	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
+++ sasmodels/kernel_template.c	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
@@ -133,5 +133,5 @@
 inline double square(double x) { return x*x; }
 inline double cube(double x) { return x*x*x; }
-inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; }
+inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; }
 
 
Index: sasmodels/models/barbell.c
===================================================================
--- sasmodels/models/barbell.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/barbell.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
@@ -33,5 +33,5 @@
         const double t = Gauss76Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
-        const double bj = sas_J1c(qrst*sqrt(radical));
+        const double bj = sas_2J1x_x(qrst*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
         total += Gauss76Wt[i] * Fq;
@@ -49,6 +49,6 @@
 {
     const double bell_fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha);
-    const double bj = sas_J1c(q*radius*sin_alpha);
-    const double si = sinc(q*half_length*cos_alpha);
+    const double bj = sas_2J1x_x(q*radius*sin_alpha);
+    const double si = sas_sinx_x(q*half_length*cos_alpha);
     const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si;
     const double Aq = bell_fq + cyl_fq;
Index: sasmodels/models/bcc_paracrystal.py
===================================================================
--- sasmodels/models/bcc_paracrystal.py	(revision b0c42719ac9096c35137bf644db59e1c5915b06b)
+++ sasmodels/models/bcc_paracrystal.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -128,5 +128,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc_paracrystal.c"]
+source = ["lib/sas_3j1x_x.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc_paracrystal.c"]
 
 # parameters for demo
Index: sasmodels/models/binary_hard_sphere.c
===================================================================
--- sasmodels/models/binary_hard_sphere.c	(revision 4f79d94ef5835e630e5531fb9a69ea5d02515e60)
+++ sasmodels/models/binary_hard_sphere.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -58,6 +58,6 @@
     qr2 = r2*q;
 
-    sc1 = sph_j1c(qr1);
-    sc2 = sph_j1c(qr2);
+    sc1 = sas_3j1x_x(qr1);
+    sc2 = sas_3j1x_x(qr2);
     b1 = r1*r1*r1*(rho1-rhos)*M_4PI_3*sc1;
     b2 = r2*r2*r2*(rho2-rhos)*M_4PI_3*sc2;
Index: sasmodels/models/binary_hard_sphere.py
===================================================================
--- sasmodels/models/binary_hard_sphere.py	(revision b0c42719ac9096c35137bf644db59e1c5915b06b)
+++ sasmodels/models/binary_hard_sphere.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -108,5 +108,5 @@
              ]
 
-source = ["lib/sph_j1c.c", "binary_hard_sphere.c"]
+source = ["lib/sas_3j1x_x.c", "binary_hard_sphere.c"]
 
 # parameters for demo and documentation
Index: sasmodels/models/capped_cylinder.c
===================================================================
--- sasmodels/models/capped_cylinder.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/capped_cylinder.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
@@ -39,5 +39,5 @@
         const double t = Gauss76Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
-        const double bj = sas_J1c(qrst*sqrt(radical));
+        const double bj = sas_2J1x_x(qrst*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
         total += Gauss76Wt[i] * Fq;
@@ -54,6 +54,6 @@
 {
     const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha);
-    const double bj = sas_J1c(q*radius*sin_alpha);
-    const double si = sinc(q*half_length*cos_alpha);
+    const double bj = sas_2J1x_x(q*radius*sin_alpha);
+    const double si = sas_sinx_x(q*half_length*cos_alpha);
     const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si;
     const double Aq = cap_Fq + cyl_Fq;
Index: sasmodels/models/core_multi_shell.c
===================================================================
--- sasmodels/models/core_multi_shell.c	(revision c5ac2b2355ae3705fd04e87490d99bfe4140d545)
+++ sasmodels/models/core_multi_shell.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -3,5 +3,5 @@
 f_constant(double q, double r, double sld)
 {
-  const double bes = sph_j1c(q * r);
+  const double bes = sas_3j1x_x(q * r);
   const double vol = M_4PI_3 * cube(r);
   return sld * vol * bes;
@@ -33,9 +33,9 @@
   f = 0.;
   for (int i=0; i<n; i++) {
-    f += M_4PI_3 * cube(r) * (sld[i] - last_sld) * sph_j1c(q*r);
+    f += M_4PI_3 * cube(r) * (sld[i] - last_sld) * sas_3j1x_x(q*r);
     last_sld = sld[i];
     r += thickness[i];
   }
-  f += M_4PI_3 * cube(r) * (solvent_sld - last_sld) * sph_j1c(q*r);
+  f += M_4PI_3 * cube(r) * (solvent_sld - last_sld) * sas_3j1x_x(q*r);
   return f * f * 1.0e-4;
 }
Index: sasmodels/models/core_multi_shell.py
===================================================================
--- sasmodels/models/core_multi_shell.py	(revision 2d73a53de4a4fb58419bdaa7cf5c5e8cddaa8c5c)
+++ sasmodels/models/core_multi_shell.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -101,5 +101,5 @@
              ]
 
-source = ["lib/sph_j1c.c", "core_multi_shell.c"]
+source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"]
 
 def profile(sld_core, radius, sld_solvent, n, sld, thickness):
Index: sasmodels/models/core_shell_bicelle.c
===================================================================
--- sasmodels/models/core_shell_bicelle.c	(revision 5bddd89a3ed5b121f80ca5d5663daeb28ce822d8)
+++ sasmodels/models/core_shell_bicelle.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
@@ -55,8 +55,8 @@
     double sinarg2 = qq*(length+facthick)*cos_alpha;
 
-    be1 = sas_J1c(besarg1);
-    be2 = sas_J1c(besarg2);
-    si1 = sinc(sinarg1);
-    si2 = sinc(sinarg2);
+    be1 = sas_2J1x_x(besarg1);
+    be2 = sas_2J1x_x(besarg2);
+    si1 = sas_sinx_x(sinarg1);
+    si2 = sas_sinx_x(sinarg2);
 
     const double t = vol1*dr1*si1*be1 +
Index: sasmodels/models/core_shell_bicelle.py
===================================================================
--- sasmodels/models/core_shell_bicelle.py	(revision fcb33e4672814b4a8fe1ee865d430090b5fc03e6)
+++ sasmodels/models/core_shell_bicelle.py	(revision 4b541acdec7c4de7cd7f18ff5911816123ce8333)
@@ -141,5 +141,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c",
+source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c",
           "core_shell_bicelle.c"]
 
Index: sasmodels/models/core_shell_bicelle_elliptical.c
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical.c	(revision fcb33e4672814b4a8fe1ee865d430090b5fc03e6)
+++ sasmodels/models/core_shell_bicelle_elliptical.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
@@ -76,6 +76,6 @@
         double sinarg1 = qq*halfheight*cos_alpha;
         double sinarg2 = qq*(halfheight+facthick)*cos_alpha;
-        si1 = sinc(sinarg1);
-        si2 = sinc(sinarg2);
+        si1 = sas_sinx_x(sinarg1);
+        si2 = sas_sinx_x(sinarg2);
         for(int j=0;j<76;j++) {
             //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe)
@@ -85,6 +85,6 @@
             double besarg1 = qq*rr*sin_alpha;
             double besarg2 = qq*(rr+radthick)*sin_alpha;
-            be1 = sas_J1c(besarg1);
-            be2 = sas_J1c(besarg2);
+            be1 = sas_2J1x_x(besarg1);
+            be2 = sas_2J1x_x(besarg2);
             inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 +
                                               dr2*si2*be2 +
@@ -129,8 +129,8 @@
     // ASSUME the sin_alpha is included in the separate integration over orientation of rod angle
     const double r = rad*sqrt(square(x_core*cos_nu) + cos_mu*cos_mu);
-    const double be1 = sas_J1c(qq*r);
-    const double be2 = sas_J1c( qq*(r + radthick ) );
-    const double si1 = sinc( qq*halfheight*cos_val );
-    const double si2 = sinc( qq*(halfheight + facthick)*cos_val );
+    const double be1 = sas_2J1x_x( qq*r );
+    const double be2 = sas_2J1x_x( qq*(r + radthick ) );
+    const double si1 = sas_sinx_x( qq*halfheight*cos_val );
+    const double si2 = sas_sinx_x( qq*(halfheight + facthick)*cos_val );
     const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1);
     //const double vol = form_volume(radius_minor, r_ratio, length);
Index: sasmodels/models/core_shell_bicelle_elliptical.py
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical.py	(revision fcb33e4672814b4a8fe1ee865d430090b5fc03e6)
+++ sasmodels/models/core_shell_bicelle_elliptical.py	(revision 8e68ea06a7246c4bd95ce82c755719ced6fa5a60)
@@ -132,5 +132,6 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c","core_shell_bicelle_elliptical.c"]
+source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c",
+          "core_shell_bicelle_elliptical.c"]
 
 demo = dict(scale=1, background=0,
Index: sasmodels/models/core_shell_cylinder.c
===================================================================
--- sasmodels/models/core_shell_cylinder.c	(revision 9aa48810160434e0c0c0901a1d80874ca75855c9)
+++ sasmodels/models/core_shell_cylinder.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
@@ -11,5 +11,5 @@
 double _cyl(double vd, double besarg, double siarg)
 {
-    return vd * sinc(siarg) * sas_J1c(besarg);
+    return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg);
 }
 
Index: sasmodels/models/core_shell_ellipsoid.py
===================================================================
--- sasmodels/models/core_shell_ellipsoid.py	(revision 7c2935c1e6bf1f935aa7dd83519aeb98216c717a)
+++ sasmodels/models/core_shell_ellipsoid.py	(revision 8e68ea06a7246c4bd95ce82c755719ced6fa5a60)
@@ -137,5 +137,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/gfn.c", "lib/gauss76.c",
+source = ["lib/sas_3j1x_x.c", "lib/gfn.c", "lib/gauss76.c",
           "core_shell_ellipsoid.c"]
 
Index: sasmodels/models/core_shell_parallelepiped.c
===================================================================
--- sasmodels/models/core_shell_parallelepiped.c	(revision 14838a3913650894395800a6e3abd9c585ae40d8)
+++ sasmodels/models/core_shell_parallelepiped.c	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
@@ -87,8 +87,8 @@
             double sin_uu, cos_uu;
             SINCOS(M_PI_2*uu, sin_uu, cos_uu);
-            const double si1 = sinc(mu_proj * sin_uu * a_scaled);
-            const double si2 = sinc(mu_proj * cos_uu);
-            const double si3 = sinc(mu_proj * sin_uu * ta);
-            const double si4 = sinc(mu_proj * cos_uu * tb);
+            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
+            const double si2 = sas_sinx_x(mu_proj * cos_uu);
+            const double si3 = sas_sinx_x(mu_proj * sin_uu * ta);
+            const double si4 = sas_sinx_x(mu_proj * cos_uu * tb);
 
             // Expression in libCylinder.c (neither drC nor Vot are used)
@@ -109,5 +109,5 @@
 
         // now sum up the outer integral
-        const double si = sinc(mu * c_scaled * sigma);
+        const double si = sas_sinx_x(mu * c_scaled * sigma);
         outer_total += Gauss76Wt[i] * inner_total * si * si;
     }
@@ -160,10 +160,10 @@
     double tc = length_a + 2.0*thick_rim_c;
     //handle arg=0 separately, as sin(t)/t -> 1 as t->0
-    double siA = sinc(0.5*q*length_a*cos_val_a);
-    double siB = sinc(0.5*q*length_b*cos_val_b);
-    double siC = sinc(0.5*q*length_c*cos_val_c);
-    double siAt = sinc(0.5*q*ta*cos_val_a);
-    double siBt = sinc(0.5*q*tb*cos_val_b);
-    double siCt = sinc(0.5*q*tc*cos_val_c);
+    double siA = sas_sinx_x(0.5*q*length_a*cos_val_a);
+    double siB = sas_sinx_x(0.5*q*length_b*cos_val_b);
+    double siC = sas_sinx_x(0.5*q*length_c*cos_val_c);
+    double siAt = sas_sinx_x(0.5*q*ta*cos_val_a);
+    double siBt = sas_sinx_x(0.5*q*tb*cos_val_b);
+    double siCt = sas_sinx_x(0.5*q*tc*cos_val_c);
     
 
Index: sasmodels/models/core_shell_sphere.py
===================================================================
--- sasmodels/models/core_shell_sphere.py	(revision 4962519cc1910422a467e4f643dba223623e0e09)
+++ sasmodels/models/core_shell_sphere.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -75,5 +75,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/core_shell.c", "core_shell_sphere.c"]
+source = ["lib/sas_3j1x_x.c", "lib/core_shell.c", "core_shell_sphere.c"]
 
 demo = dict(scale=1, background=0, radius=60, thickness=10,
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision b829b161dd452ced7aab4d44d03905daab7234ac)
+++ sasmodels/models/cylinder.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
@@ -18,5 +18,5 @@
     const double qr = q*radius;
     const double qh = q*0.5*length; 
-    return sas_J1c(qr*sn) * sinc(qh*cn);
+    return sas_2J1x_x(qr*sn) * sas_sinx_x(qh*cn);
 }
 
Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision 73e08aee242fc777e9b6a88ee4ab9fea7732e685)
+++ sasmodels/models/ellipsoid.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -18,5 +18,5 @@
     const double r = radius_equatorial
                      * sqrt(1.0 + sin_alpha*sin_alpha*(ratio*ratio - 1.0));
-    const double f = sph_j1c(q*r);
+    const double f = sas_3j1x_x(q*r);
 
     return f*f;
Index: sasmodels/models/ellipsoid.py
===================================================================
--- sasmodels/models/ellipsoid.py	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
+++ sasmodels/models/ellipsoid.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -135,5 +135,5 @@
              ]
 
-source = ["lib/sph_j1c.c", "lib/gauss76.c", "ellipsoid.c"]
+source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"]
 
 def ER(radius_polar, radius_equatorial):
Index: sasmodels/models/elliptical_cylinder.c
===================================================================
--- sasmodels/models/elliptical_cylinder.c	(revision 251f54b51119923e12cbc8eae49fefaeb6158a47)
+++ sasmodels/models/elliptical_cylinder.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
@@ -39,5 +39,5 @@
             const double theta = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
             const double r = sin_val*sqrt(rA - rB*cos(theta));
-            const double be = sas_J1c(q*r);
+            const double be = sas_2J1x_x(q*r);
             inner_sum += Gauss20Wt[j] * be * be;
         }
@@ -46,5 +46,5 @@
 
         //now calculate outer integral
-        const double si = sinc(q*0.5*length*cos_val);
+        const double si = sas_sinx_x(q*0.5*length*cos_val);
         outer_sum += Gauss76Wt[i] * inner_sum * si * si;
     }
@@ -73,6 +73,6 @@
     // Given:    radius_major = r_ratio * radius_minor
     const double r = radius_minor*sqrt(square(r_ratio*cos_nu) + cos_mu*cos_mu);
-    const double be = sas_J1c(q*r);
-    const double si = sinc(q*0.5*length*cos_val);
+    const double be = sas_2J1x_x(q*r);
+    const double si = sas_sinx_x(q*0.5*length*cos_val);
     const double Aq = be * si;
     const double delrho = sld - solvent_sld;
Index: sasmodels/models/fcc_paracrystal.py
===================================================================
--- sasmodels/models/fcc_paracrystal.py	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
+++ sasmodels/models/fcc_paracrystal.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -116,5 +116,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc_paracrystal.c"]
+source = ["lib/sas_3j1x_x.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc_paracrystal.c"]
 
 # parameters for demo
Index: sasmodels/models/flexible_cylinder.c
===================================================================
--- sasmodels/models/flexible_cylinder.c	(revision e6408d049d3598925e345d8adec1cdb0bb0c6adb)
+++ sasmodels/models/flexible_cylinder.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
@@ -14,5 +14,5 @@
 {
     const double contrast = sld - solvent_sld;
-    const double cross_section = sas_J1c(q*radius);
+    const double cross_section = sas_2J1x_x(q*radius);
     const double volume = M_PI*radius*radius*length;
     const double flex = Sk_WR(q, length, kuhn_length);
Index: sasmodels/models/flexible_cylinder_elliptical.c
===================================================================
--- sasmodels/models/flexible_cylinder_elliptical.c	(revision 92ce163837b3c51ecc732713d4a5ea719c037b29)
+++ sasmodels/models/flexible_cylinder_elliptical.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
@@ -22,5 +22,5 @@
         SINCOS(zi, sn, cn);
         const double arg = q*sqrt(a*a*sn*sn + b*b*cn*cn);
-        const double yyy = sas_J1c(arg);
+        const double yyy = sas_2J1x_x(arg);
         sum += Gauss76Wt[i] * yyy * yyy;
     }
Index: sasmodels/models/fractal.c
===================================================================
--- sasmodels/models/fractal.c	(revision 217590bce5ffcf2c4de704a4770a9cab003dd287)
+++ sasmodels/models/fractal.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -14,5 +14,5 @@
     //calculate P(q) for the spherical subunits
     const double V = M_4PI_3*cube(radius);
-    const double pq = V * square((sld_block-sld_solvent)*sph_j1c(q*radius));
+    const double pq = V * square((sld_block-sld_solvent)*sas_3j1x_x(q*radius));
 
     // scale to units cm-1 sr-1 (assuming data on absolute scale)
Index: sasmodels/models/fractal.py
===================================================================
--- sasmodels/models/fractal.py	(revision d1cfa86b03d71f574114b242dc91eefade357ba2)
+++ sasmodels/models/fractal.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -97,5 +97,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "lib/fractal_sq.c", "fractal.c"]
+source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "lib/fractal_sq.c", "fractal.c"]
 
 demo = dict(volfraction=0.05,
Index: sasmodels/models/fractal_core_shell.py
===================================================================
--- sasmodels/models/fractal_core_shell.py	(revision d6f60c35e7286f24ca0fdd9d5f740e642d007b8a)
+++ sasmodels/models/fractal_core_shell.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -95,5 +95,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "lib/core_shell.c",
+source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "lib/core_shell.c",
           "lib/fractal_sq.c", "fractal_core_shell.c"]
 
Index: sasmodels/models/fuzzy_sphere.py
===================================================================
--- sasmodels/models/fuzzy_sphere.py	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/fuzzy_sphere.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -81,5 +81,5 @@
 # pylint: enable=bad-whitespace,line-too-long
 
-source = ["lib/sph_j1c.c"]
+source = ["lib/sas_3j1x_x.c"]
 
 # No volume normalization despite having a volume parameter
@@ -91,5 +91,5 @@
 Iq = """
     const double qr = q*radius;
-    const double bes = sph_j1c(qr);
+    const double bes = sas_3j1x_x(qr);
     const double qf = q*fuzziness;
     const double fq = bes * (sld - sld_solvent) * form_volume(radius) * exp(-0.5*qf*qf);
Index: sasmodels/models/hollow_cylinder.c
===================================================================
--- sasmodels/models/hollow_cylinder.c	(revision f8f0991d67f2db3a0efcd8211a67ea2f580210a2)
+++ sasmodels/models/hollow_cylinder.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
@@ -20,11 +20,11 @@
 {
     const double qs = q*sin_val;
-    const double lam1 = sas_J1c((radius+thickness)*qs);
-    const double lam2 = sas_J1c(radius*qs);
+    const double lam1 = sas_2J1x_x((radius+thickness)*qs);
+    const double lam2 = sas_2J1x_x(radius*qs);
     const double gamma_sq = square(radius/(radius+thickness));
-    //Note: lim_{thickness -> 0} psi = J0(radius*qs)
-    //Note: lim_{radius -> 0} psi = sas_J1c(thickness*qs)
+    //Note: lim_{thickness -> 0} psi = sas_J0(radius*qs)
+    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qs)
     const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq);	//SRK 10/19/00
-    const double t2 = sinc(0.5*q*length*cos_val);
+    const double t2 = sas_sinx_x(0.5*q*length*cos_val);
     return psi*t2;
 }
Index: sasmodels/models/hollow_rectangular_prism.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism.c	(revision 6f676fb91e0720edd32b0faa11126cc646a2cad7)
+++ sasmodels/models/hollow_rectangular_prism.c	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
@@ -45,6 +45,6 @@
         SINCOS(theta, sin_theta, cos_theta);
 
-        const double termC1 = sinc(q * c_half * cos(theta));
-        const double termC2 = sinc(q * (c_half-thickness)*cos(theta));
+        const double termC1 = sas_sinx_x(q * c_half * cos(theta));
+        const double termC2 = sas_sinx_x(q * (c_half-thickness)*cos(theta));
 
         double inner_sum = 0.0;
@@ -57,9 +57,9 @@
             // Amplitude AP from eqn. (13), rewritten to avoid round-off effects when arg=0
 
-            const double termA1 = sinc(q * a_half * sin_theta * sin_phi);
-            const double termA2 = sinc(q * (a_half-thickness) * sin_theta * sin_phi);
+            const double termA1 = sas_sinx_x(q * a_half * sin_theta * sin_phi);
+            const double termA2 = sas_sinx_x(q * (a_half-thickness) * sin_theta * sin_phi);
 
-            const double termB1 = sinc(q * b_half * sin_theta * cos_phi);
-            const double termB2 = sinc(q * (b_half-thickness) * sin_theta * cos_phi);
+            const double termB1 = sas_sinx_x(q * b_half * sin_theta * cos_phi);
+            const double termB2 = sas_sinx_x(q * (b_half-thickness) * sin_theta * cos_phi);
 
             const double AP1 = vol_total * termA1 * termB1 * termC1;
Index: smodels/models/lib/Si.c
===================================================================
--- sasmodels/models/lib/Si.c	(revision f719764b9af71d472c265bb3ac96d855dfc1bf06)
+++ 	(revision )
@@ -1,24 +1,0 @@
-// integral of sin(x)/x Taylor series approximated to w/i 0.1%
-double Si(double x);
-double Si(double x)
-{
-    if (x >= M_PI*6.2/4.0) {
-        const double xxinv = 1./(x*x);
-        // Explicitly writing factorial values triples the speed of the calculation
-        const double out_cos = (((-720.*xxinv + 24.)*xxinv - 2.)*xxinv + 1.)/x;
-        const double out_sin = (((-5040.*xxinv + 120.)*xxinv - 6.)*xxinv + 1)*xxinv;
-
-        double sin_x, cos_x;
-        SINCOS(x, sin_x, cos_x);
-        return M_PI_2 - cos_x*out_cos - sin_x*out_sin;
-    } else {
-        const double xx = x*x;
-        // Explicitly writing factorial values triples the speed of the calculation
-        return (((((-1./439084800.*xx
-            + 1./3265920.)*xx
-            - 1./35280.)*xx
-            + 1./600.)*xx
-            - 1./18.)*xx
-            + 1.)*x;
-    }
-}
Index: sasmodels/models/lib/core_shell.c
===================================================================
--- sasmodels/models/lib/core_shell.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/lib/core_shell.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -16,5 +16,5 @@
     const double core_qr = q * radius;
     const double core_contrast = core_sld - shell_sld;
-    const double core_bes = sph_j1c(core_qr);
+    const double core_bes = sas_3j1x_x(core_qr);
     const double core_volume = M_4PI_3 * cube(radius);
     double f = core_volume * core_bes * core_contrast;
@@ -23,5 +23,5 @@
     const double shell_qr = q * (radius + thickness);
     const double shell_contrast = shell_sld - solvent_sld;
-    const double shell_bes = sph_j1c(shell_qr);
+    const double shell_bes = sas_3j1x_x(shell_qr);
     const double shell_volume = M_4PI_3 * cube(radius + thickness);
     f += shell_volume * shell_bes * shell_contrast;
Index: sasmodels/models/lib/gfn.c
===================================================================
--- sasmodels/models/lib/gfn.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/lib/gfn.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -18,5 +18,5 @@
     // changing to more accurate sph_j1c since the following inexplicably fails on Radeon Nano.
     //const double siq = (uq == 0.0 ? 1.0 : 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq);
-    const double siq = sph_j1c(uq);
+    const double siq = sas_3j1x_x(uq);
     const double vc = M_4PI_3*aa*aa*bb;
     const double gfnc = siq*vc*delpc;
@@ -26,5 +26,5 @@
     const double vt = M_4PI_3*trmaj*trmaj*trmin;
     //const double sit = (ut == 0.0 ? 1.0 : 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut);
-    const double sit = sph_j1c(ut);
+    const double sit = sas_3j1x_x(ut);
     const double gfnt = sit*vt*delps;
 
Index: sasmodels/models/lib/sas_3j1x_x.c
===================================================================
--- sasmodels/models/lib/sas_3j1x_x.c	(revision 473a9f131b0c8e4ae5fb862b046d8e30d4eddedb)
+++ sasmodels/models/lib/sas_3j1x_x.c	(revision 473a9f131b0c8e4ae5fb862b046d8e30d4eddedb)
@@ -0,0 +1,56 @@
+/**
+* Spherical Bessel function 3*j1(x)/x
+*
+* Used for low q to avoid cancellation error.
+* Note that the values differ from sasview ~ 5e-12 rather than 5e-14, but
+* in this case it is likely cancellation errors in the original expression
+* using double precision that are the source.
+*/
+double sas_3j1x_x(double q);
+
+// The choice of the number of terms in the series and the cutoff value for
+// switching between series and direct calculation depends on the numeric
+// precision.
+//
+// Point where direct calculation reaches machine precision:
+//
+//   single machine precision eps 3e-8 at qr=1.1 **
+//   double machine precision eps 4e-16 at qr=1.1
+//
+// Point where Taylor series reaches machine precision (eps), where taylor
+// series matches direct calculation (cross) and the error at that point:
+//
+//   prec   n eps  cross  error
+//   single 3 0.28  0.4   6.2e-7
+//   single 4 0.68  0.7   2.3e-7
+//   single 5 1.18  1.2   7.5e-8
+//   double 3 0.01  0.03  2.3e-13
+//   double 4 0.06  0.1   3.1e-14
+//   double 5 0.16  0.2   5.0e-15
+//
+// ** Note: relative error on single precision starts increase on the direct
+// method at qr=1.1, rising from 3e-8 to 5e-5 by qr=1e3.  This should be
+// safe for the sans range, with objects of 100 nm supported to a q of 0.1
+// while maintaining 5 digits of precision.  For usans/sesans, the objects
+// are larger but the q is smaller, so again it should be fine.
+//
+// See explore/sph_j1c.py for code to explore these ranges.
+
+// Use 4th order series
+#if FLOAT_SIZE>4
+#define SPH_J1C_CUTOFF 0.1
+#else
+#define SPH_J1C_CUTOFF 0.7
+#endif
+
+double sas_3j1x_x(double q)
+{
+    if (q < SPH_J1C_CUTOFF) {
+        const double q2 = q*q;
+        return (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))));// + q2*(3./3991680.)))));
+    } else {
+        double sin_q, cos_q;
+        SINCOS(q, sin_q, cos_q);
+        return 3.0*(sin_q/q - cos_q)/(q*q);
+    }
+}
Index: sasmodels/models/lib/sas_J1.c
===================================================================
--- sasmodels/models/lib/sas_J1.c	(revision c8902ac8c2f740f2cd6d59d72c25f8957b7bdd8e)
+++ sasmodels/models/lib/sas_J1.c	(revision 473a9f131b0c8e4ae5fb862b046d8e30d4eddedb)
@@ -217,6 +217,6 @@
 
 //Finally J1c function that equals 2*J1(x)/x
-double sas_J1c(double x);
-double sas_J1c(double x)
+double sas_2J1x_x(double x);
+double sas_2J1x_x(double x)
 {
     return (x != 0.0 ) ? 2.0*sas_J1(x)/x : 1.0;
Index: sasmodels/models/lib/sas_Si.c
===================================================================
--- sasmodels/models/lib/sas_Si.c	(revision 473a9f131b0c8e4ae5fb862b046d8e30d4eddedb)
+++ sasmodels/models/lib/sas_Si.c	(revision 473a9f131b0c8e4ae5fb862b046d8e30d4eddedb)
@@ -0,0 +1,24 @@
+// integral of sin(x)/x Taylor series approximated to w/i 0.1%
+double sas_Si(double x);
+double sas_Si(double x)
+{
+    if (x >= M_PI*6.2/4.0) {
+        const double xxinv = 1./(x*x);
+        // Explicitly writing factorial values triples the speed of the calculation
+        const double out_cos = (((-720.*xxinv + 24.)*xxinv - 2.)*xxinv + 1.)/x;
+        const double out_sin = (((-5040.*xxinv + 120.)*xxinv - 6.)*xxinv + 1)*xxinv;
+
+        double sin_x, cos_x;
+        SINCOS(x, sin_x, cos_x);
+        return M_PI_2 - cos_x*out_cos - sin_x*out_sin;
+    } else {
+        const double xx = x*x;
+        // Explicitly writing factorial values triples the speed of the calculation
+        return (((((-1./439084800.*xx
+            + 1./3265920.)*xx
+            - 1./35280.)*xx
+            + 1./600.)*xx
+            - 1./18.)*xx
+            + 1.)*x;
+    }
+}
Index: smodels/models/lib/sph_j1c.c
===================================================================
--- sasmodels/models/lib/sph_j1c.c	(revision ba32cddebeff04b66d828ca58b1843244d4425fb)
+++ 	(revision )
@@ -1,56 +1,0 @@
-/**
-* Spherical Bessel function 3*j1(x)/x
-*
-* Used for low q to avoid cancellation error.
-* Note that the values differ from sasview ~ 5e-12 rather than 5e-14, but
-* in this case it is likely cancellation errors in the original expression
-* using double precision that are the source.
-*/
-double sph_j1c(double q);
-
-// The choice of the number of terms in the series and the cutoff value for
-// switching between series and direct calculation depends on the numeric
-// precision.
-//
-// Point where direct calculation reaches machine precision:
-//
-//   single machine precision eps 3e-8 at qr=1.1 **
-//   double machine precision eps 4e-16 at qr=1.1
-//
-// Point where Taylor series reaches machine precision (eps), where taylor
-// series matches direct calculation (cross) and the error at that point:
-//
-//   prec   n eps  cross  error
-//   single 3 0.28  0.4   6.2e-7
-//   single 4 0.68  0.7   2.3e-7
-//   single 5 1.18  1.2   7.5e-8
-//   double 3 0.01  0.03  2.3e-13
-//   double 4 0.06  0.1   3.1e-14
-//   double 5 0.16  0.2   5.0e-15
-//
-// ** Note: relative error on single precision starts increase on the direct
-// method at qr=1.1, rising from 3e-8 to 5e-5 by qr=1e3.  This should be
-// safe for the sans range, with objects of 100 nm supported to a q of 0.1
-// while maintaining 5 digits of precision.  For usans/sesans, the objects
-// are larger but the q is smaller, so again it should be fine.
-//
-// See explore/sph_j1c.py for code to explore these ranges.
-
-// Use 4th order series
-#if FLOAT_SIZE>4
-#define SPH_J1C_CUTOFF 0.1
-#else
-#define SPH_J1C_CUTOFF 0.7
-#endif
-
-double sph_j1c(double q)
-{
-    if (q < SPH_J1C_CUTOFF) {
-        const double q2 = q*q;
-        return (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))));// + q2*(3./3991680.)))));
-    } else {
-        double sin_q, cos_q;
-        SINCOS(q, sin_q, cos_q);
-        return 3.0*(sin_q/q - cos_q)/(q*q);
-    }
-}
Index: sasmodels/models/lib/sphere_form.c
===================================================================
--- sasmodels/models/lib/sphere_form.c	(revision ba32cddebeff04b66d828ca58b1843244d4425fb)
+++ sasmodels/models/lib/sphere_form.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -9,5 +9,5 @@
 double sphere_form(double q, double radius, double sld, double solvent_sld)
 {
-    const double fq = sphere_volume(radius) * sph_j1c(q*radius);
+    const double fq = sphere_volume(radius) * sas_3j1x_x(q*radius);
     const double contrast = (sld - solvent_sld);
     return 1.0e-4*square(contrast * fq);
Index: sasmodels/models/linear_pearls.c
===================================================================
--- sasmodels/models/linear_pearls.c	(revision 4962519cc1910422a467e4f643dba223623e0e09)
+++ sasmodels/models/linear_pearls.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -44,5 +44,5 @@
 
     //sine functions of a pearl
-    double psi = sph_j1c(q * radius);
+    double psi = sas_3j1x_x(q * radius);
 
     // N pearls contribution
@@ -50,5 +50,5 @@
     n_contrib = num_pearls;
     for(int num=1; num<=n_max; num++) {
-        n_contrib += (2.0*(num_pearls-num)*sinc(q*separation*num));
+        n_contrib += (2.0*(num_pearls-num)*sas_sinx_x(q*separation*num));
     }
     // form factor for num_pearls
Index: sasmodels/models/linear_pearls.py
===================================================================
--- sasmodels/models/linear_pearls.py	(revision 4962519cc1910422a467e4f643dba223623e0e09)
+++ sasmodels/models/linear_pearls.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -63,5 +63,5 @@
 single = False
 
-source = ["lib/sph_j1c.c", "linear_pearls.c"]
+source = ["lib/sas_3j1x_x.c", "linear_pearls.c"]
 
 demo = dict(scale=1.0, background=0.0,
Index: sasmodels/models/mass_fractal.c
===================================================================
--- sasmodels/models/mass_fractal.c	(revision 6d96b66e1e58166fbc135a8df1ddf8d37f64b54a)
+++ sasmodels/models/mass_fractal.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -5,5 +5,5 @@
 {
     //calculate P(q)
-    const double pq = square(sph_j1c(q*radius));
+    const double pq = square(sas_3j1x_x(q*radius));
 
     //calculate S(q)
Index: sasmodels/models/mass_fractal.py
===================================================================
--- sasmodels/models/mass_fractal.py	(revision 6d96b66e1e58166fbc135a8df1ddf8d37f64b54a)
+++ sasmodels/models/mass_fractal.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -86,5 +86,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "mass_fractal.c"]
+source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "mass_fractal.c"]
 
 demo = dict(scale=1, background=0,
Index: sasmodels/models/multilayer_vesicle.c
===================================================================
--- sasmodels/models/multilayer_vesicle.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/multilayer_vesicle.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -20,5 +20,5 @@
         // layer 1
         voli = M_4PI_3*ri*ri*ri;
-        fval += voli*sldi*sph_j1c(ri*q);
+        fval += voli*sldi*sas_3j1x_x(ri*q);
 
         ri += thick_shell;
@@ -26,5 +26,5 @@
         // layer 2
         voli = M_4PI_3*ri*ri*ri;
-        fval -= voli*sldi*sph_j1c(ri*q);
+        fval -= voli*sldi*sas_3j1x_x(ri*q);
 
         //do 2 layers at a time
Index: sasmodels/models/multilayer_vesicle.py
===================================================================
--- sasmodels/models/multilayer_vesicle.py	(revision 1a6cd57bd7b12db89d69b335e984c55413295845)
+++ sasmodels/models/multilayer_vesicle.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -107,5 +107,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "multilayer_vesicle.c"]
+source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"]
 
 polydispersity = ["radius", "n_pairs"]
Index: sasmodels/models/onion.c
===================================================================
--- sasmodels/models/onion.c	(revision 97623415b6465256171bc225ecfa9d65aac160ce)
+++ sasmodels/models/onion.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -6,5 +6,5 @@
   const double vol = M_4PI_3 * cube(r);
   const double qr = q * r;
-  const double bes = sph_j1c(qr);
+  const double bes = sas_3j1x_x(qr);
   const double alpha = A * r/thickness;
   double result;
Index: sasmodels/models/onion.py
===================================================================
--- sasmodels/models/onion.py	(revision 97623415b6465256171bc225ecfa9d65aac160ce)
+++ sasmodels/models/onion.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -314,5 +314,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "onion.c"]
+source = ["lib/sas_3j1x_x.c", "onion.c"]
 single = False
 
Index: sasmodels/models/parallelepiped.c
===================================================================
--- sasmodels/models/parallelepiped.c	(revision 14838a3913650894395800a6e3abd9c585ae40d8)
+++ sasmodels/models/parallelepiped.c	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
@@ -39,11 +39,11 @@
             double sin_uu, cos_uu;
             SINCOS(M_PI_2*uu, sin_uu, cos_uu);
-            const double si1 = sinc(mu_proj * sin_uu * a_scaled);
-            const double si2 = sinc(mu_proj * cos_uu);
+            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
+            const double si2 = sas_sinx_x(mu_proj * cos_uu);
             inner_total += Gauss76Wt[j] * square(si1 * si2);
         }
         inner_total *= 0.5;
 
-        const double si = sinc(mu * c_scaled * sigma);
+        const double si = sas_sinx_x(mu * c_scaled * sigma);
         outer_total += Gauss76Wt[i] * inner_total * si * si;
     }
@@ -70,7 +70,7 @@
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a);
 
-    const double siA = sinc(0.5*q*length_a*cos_val_a);
-    const double siB = sinc(0.5*q*length_b*cos_val_b);
-    const double siC = sinc(0.5*q*length_c*cos_val_c);
+    const double siA = sas_sinx_x(0.5*q*length_a*cos_val_a);
+    const double siB = sas_sinx_x(0.5*q*length_b*cos_val_b);
+    const double siC = sas_sinx_x(0.5*q*length_c*cos_val_c);
     const double V = form_volume(length_a, length_b, length_c);
     const double drho = (sld - solvent_sld);
Index: sasmodels/models/pearl_necklace.c
===================================================================
--- sasmodels/models/pearl_necklace.c	(revision 21261311ccad3df50fbfc3f21eb04942110f1c9d)
+++ sasmodels/models/pearl_necklace.c	(revision 4b541acdec7c4de7cd7f18ff5911816123ce8333)
@@ -34,10 +34,10 @@
     // But there is a 1/(1-sinc) term below which blows up so don't bother
     const double q_edge = q * edge_sep;
-    const double beta = (Si(q*(A_s-radius)) - Si(q*radius)) / q_edge;
-    const double gamma = Si(q_edge) / q_edge;
-    const double psi = sph_j1c(q*radius);
+    const double beta = (sas_Si(q*(A_s-radius)) - sas_Si(q*radius)) / q_edge;
+    const double gamma = sas_Si(q_edge) / q_edge;
+    const double psi = sas_3j1x_x(q*radius);
 
     // Precomputed sinc terms
-    const double si = sinc(q*A_s);
+    const double si = sas_sinx_x(q*A_s);
     const double omsi = 1.0 - si;
     const double pow_si = pow(si, num_pearls);
@@ -54,5 +54,5 @@
         - 2.0 * (1.0 - pow_si/si)*beta*beta / (omsi*omsi)
         + 2.0 * num_strings*beta*beta / omsi
-        + num_strings * (2.0*gamma - square(sinc(q_edge/2.0)))
+        + num_strings * (2.0*gamma - square(sas_sinx_x(q_edge/2.0)))
         );
 
Index: sasmodels/models/pearl_necklace.py
===================================================================
--- sasmodels/models/pearl_necklace.py	(revision 21261311ccad3df50fbfc3f21eb04942110f1c9d)
+++ sasmodels/models/pearl_necklace.py	(revision 4b541acdec7c4de7cd7f18ff5911816123ce8333)
@@ -92,5 +92,5 @@
              ]
 
-source = ["lib/Si.c", "lib/sph_j1c.c", "pearl_necklace.c"]
+source = ["lib/sas_Si.c", "lib/sas_3j1x_x.c", "pearl_necklace.c"]
 single = False  # use double precision unless told otherwise
 
Index: sasmodels/models/polymer_micelle.c
===================================================================
--- sasmodels/models/polymer_micelle.c	(revision c3ebc71c4319402dc4ddde47d1ee1264eb354d01)
+++ sasmodels/models/polymer_micelle.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -31,5 +31,5 @@
 
     // Self-correlation term of the core
-    const double bes_core = sph_j1c(q*radius_core);
+    const double bes_core = sas_3j1x_x(q*radius_core);
     const double term1 = square(n_aggreg*beta_core*bes_core);
 
@@ -41,5 +41,5 @@
     // Interference cross-term between core and chains
     const double chain_ampl = (qrg2 == 0.0) ? 1.0 : -expm1(-qrg2)/qrg2;
-    const double bes_corona = sinc(q*(radius_core + d_penetration * rg));
+    const double bes_corona = sas_sinx_x(q*(radius_core + d_penetration * rg));
     const double term3 = 2.0 * n_aggreg * n_aggreg * beta_core * beta_corona *
                  bes_core * chain_ampl * bes_corona;
Index: sasmodels/models/polymer_micelle.py
===================================================================
--- sasmodels/models/polymer_micelle.py	(revision bba9361c45987d76f3673d7d0a4eb16d28a8b701)
+++ sasmodels/models/polymer_micelle.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -104,5 +104,5 @@
 single = False
 
-source = ["lib/sph_j1c.c", "polymer_micelle.c"]
+source = ["lib/sas_3j1x_x.c", "polymer_micelle.c"]
 
 demo = dict(scale=1, background=0,
Index: sasmodels/models/pringle.c
===================================================================
--- sasmodels/models/pringle.c	(revision 30fbe2eda4d6f97acc840bd5e605de3829caf3c7)
+++ sasmodels/models/pringle.c	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
@@ -91,5 +91,5 @@
         SINCOS(psi, sin_psi, cos_psi);
         double bessel_term = _sum_bessel_orders(radius, alpha, beta, q*sin_psi, q*cos_psi);
-        double sinc_term = square(sinc(q * thickness * cos_psi / 2.0));
+        double sinc_term = square(sas_sinx_x(q * thickness * cos_psi / 2.0));
         double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term;
         sum += Gauss76Wt[i] * pringle_kernel;
Index: sasmodels/models/raspberry.c
===================================================================
--- sasmodels/models/raspberry.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
+++ sasmodels/models/raspberry.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -51,11 +51,11 @@
 
     //Form factors for each particle
-    psiL = sph_j1c(q*rL);
-    psiS = sph_j1c(q*rS);
+    psiL = sas_3j1x_x(q*rL);
+    psiS = sas_3j1x_x(q*rS);
 
     //Cross term between large and small particles
-    sfLS = psiL*psiS*sinc(q*(rL+deltaS*rS));
+    sfLS = psiL*psiS*sas_sinx_x(q*(rL+deltaS*rS));
     //Cross term between small particles at the surface
-    sfSS = psiS*psiS*sinc(q*(rL+deltaS*rS))*sinc(q*(rL+deltaS*rS));
+    sfSS = psiS*psiS*sas_sinx_x(q*(rL+deltaS*rS))*sas_sinx_x(q*(rL+deltaS*rS));
 
     //Large sphere form factor term
Index: sasmodels/models/raspberry.py
===================================================================
--- sasmodels/models/raspberry.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/models/raspberry.py	(revision 8e68ea06a7246c4bd95ce82c755719ced6fa5a60)
@@ -129,5 +129,5 @@
             Ref: J. coll. inter. sci. (2010) vol. 343 (1) pp. 36-41."""
 category = "shape:sphere"
-#single = False
+
 
 #             [ "name", "units", default, [lower, upper], "type", "description"],
@@ -152,5 +152,5 @@
              ]
 
-source = ["lib/sph_j1c.c", "raspberry.c"]
+source = ["lib/sas_3j1x_x.c", "raspberry.c"]
 
 # parameters for demo
Index: sasmodels/models/rectangular_prism.c
===================================================================
--- sasmodels/models/rectangular_prism.c	(revision ab2aea81f7fb392e5f24f1c76dd989914d8f7022)
+++ sasmodels/models/rectangular_prism.c	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
@@ -33,5 +33,5 @@
         SINCOS(theta, sin_theta, cos_theta);
 
-        const double termC = sinc(q * c_half * cos_theta);
+        const double termC = sas_sinx_x(q * c_half * cos_theta);
 
         double inner_sum = 0.0;
@@ -42,6 +42,6 @@
 
             // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0
-            const double termA = sinc(q * a_half * sin_theta * sin_phi);
-            const double termB = sinc(q * b_half * sin_theta * cos_phi);
+            const double termA = sas_sinx_x(q * a_half * sin_theta * sin_phi);
+            const double termB = sas_sinx_x(q * b_half * sin_theta * cos_phi);
             const double AP = termA * termB * termC;
             inner_sum += Gauss76Wt[j] * AP * AP;
Index: sasmodels/models/sc_paracrystal.py
===================================================================
--- sasmodels/models/sc_paracrystal.py	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
+++ sasmodels/models/sc_paracrystal.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -133,5 +133,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/sphere_form.c", "lib/gauss150.c", "sc_paracrystal.c"]
+source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c", "lib/gauss150.c", "sc_paracrystal.c"]
 
 demo = dict(scale=1, background=0,
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision 7e6bea817616b3ff05f1c96c2a11a94853dc6276)
+++ sasmodels/models/sphere.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -66,5 +66,5 @@
              ]
 
-source = ["lib/sph_j1c.c", "lib/sphere_form.c"]
+source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c"]
 
 # No volume normalization despite having a volume parameter
Index: sasmodels/models/spherical_sld.c
===================================================================
--- sasmodels/models/spherical_sld.c	(revision 54bcd4abeaba4faf03e24319c050286fc10969da)
+++ sasmodels/models/spherical_sld.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -34,5 +34,5 @@
     const double qr = q * r;
     const double qrsq = qr * qr;
-    const double bes = sph_j1c(qr);
+    const double bes = sas_3j1x_x(qr);
     double sinqr, cosqr;
     SINCOS(qr, sinqr, cosqr);
@@ -60,7 +60,7 @@
 
         // uniform shell; r=0 => r^3=0 => f=0, so works for core as well.
-        f -= M_4PI_3 * cube(r) * sld_l * sph_j1c(q*r);
+        f -= M_4PI_3 * cube(r) * sld_l * sas_3j1x_x(q*r);
         r += thickness[shell];
-        f += M_4PI_3 * cube(r) * sld_l * sph_j1c(q*r);
+        f += M_4PI_3 * cube(r) * sld_l * sas_3j1x_x(q*r);
 
         // iterate over sub_shells in the interface
@@ -92,5 +92,5 @@
     }
     // add in solvent effect
-    f -= M_4PI_3 * cube(r) * sld_solvent * sph_j1c(q*r);
+    f -= M_4PI_3 * cube(r) * sld_solvent * sas_3j1x_x(q*r);
 
     const double f2 = f * f * 1.0e-4;
Index: sasmodels/models/spherical_sld.py
===================================================================
--- sasmodels/models/spherical_sld.py	(revision 2d65d51719530a31ea7870f34610b834f36a1319)
+++ sasmodels/models/spherical_sld.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -214,5 +214,5 @@
              ]
 # pylint: enable=bad-whitespace, line-too-long
-source = ["lib/polevl.c", "lib/sas_erf.c", "lib/sph_j1c.c", "spherical_sld.c"]
+source = ["lib/polevl.c", "lib/sas_erf.c", "lib/sas_3j1x_x.c", "spherical_sld.c"]
 single = False  # TODO: fix low q behaviour
 
Index: sasmodels/models/stacked_disks.c
===================================================================
--- sasmodels/models/stacked_disks.c	(revision 98ce141a6e091b02f8139ab3b235c16dafbce663)
+++ sasmodels/models/stacked_disks.c	(revision 6c3e2666e0741d50f118672e164dfa561dcc7765)
@@ -53,9 +53,9 @@
     const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha;
 
-    const double be1 = sas_J1c(besarg1);
-    //const double be2 = sas_J1c(besarg2);
+    const double be1 = sas_2J1x_x(besarg1);
+    //const double be2 = sas_2J1x_x(besarg2);
     const double be2 = be1;
-    const double si1 = sinc(sinarg1);
-    const double si2 = sinc(sinarg2);
+    const double si1 = sas_sinx_x(sinarg1);
+    const double si2 = sas_sinx_x(sinarg2);
 
     const double dr1 = core_sld - solvent_sld;
Index: sasmodels/models/surface_fractal.c
===================================================================
--- sasmodels/models/surface_fractal.c	(revision b716cc63882c24230b143debfbf45e545ba6bac2)
+++ sasmodels/models/surface_fractal.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -15,5 +15,5 @@
 {
     // calculate P(q)
-    const double pq = square(sph_j1c(q*radius));
+    const double pq = square(sas_3j1x_x(q*radius));
 
     // calculate S(q)
Index: sasmodels/models/surface_fractal.py
===================================================================
--- sasmodels/models/surface_fractal.py	(revision 5c94f414f50d8c0d9e3d33e06b3e085b7730049a)
+++ sasmodels/models/surface_fractal.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -72,5 +72,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "surface_fractal.c"]
+source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "surface_fractal.c"]
 
 demo = dict(scale=1, background=1e-5,
Index: sasmodels/models/triaxial_ellipsoid.c
===================================================================
--- sasmodels/models/triaxial_ellipsoid.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -37,5 +37,5 @@
             const double ysq = square(Gauss76Z[j]*zm + zb);
             const double t = q*sqrt(acosx2 + bsinx2*(1.0-ysq) + c2*ysq);
-            const double fq = sph_j1c(t);
+            const double fq = sas_3j1x_x(t);
             inner += Gauss76Wt[j] * fq * fq ;
         }
@@ -64,5 +64,5 @@
                           + radius_equat_major*radius_equat_major*cmu*cmu
                           + radius_polar*radius_polar*calpha*calpha);
-    const double fq = sph_j1c(t);
+    const double fq = sas_3j1x_x(t);
     const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar);
 
Index: sasmodels/models/triaxial_ellipsoid.py
===================================================================
--- sasmodels/models/triaxial_ellipsoid.py	(revision 416f5c7f5af38d5e5d926a5b6b140a94e7e60502)
+++ sasmodels/models/triaxial_ellipsoid.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -104,5 +104,5 @@
              ]
 
-source = ["lib/sph_j1c.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
+source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
 
 def ER(radius_equat_minor, radius_equat_major, radius_polar):
Index: sasmodels/models/vesicle.c
===================================================================
--- sasmodels/models/vesicle.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/vesicle.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -30,10 +30,10 @@
     contrast = sld_solvent-sld;
     vol = M_4PI_3*cube(radius);
-    f = vol * sph_j1c(q*radius) * contrast;
+    f = vol * sas_3j1x_x(q*radius) * contrast;
  
     //now the shell. No volume normalization as this is done by the caller
     contrast = sld-sld_solvent;
     vol = M_4PI_3*cube(radius+thickness);
-    f += vol * sph_j1c(q*(radius+thickness)) * contrast;
+    f += vol * sas_3j1x_x(q*(radius+thickness)) * contrast;
 
     //rescale to [cm-1]. 
Index: sasmodels/models/vesicle.py
===================================================================
--- sasmodels/models/vesicle.py	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/vesicle.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
@@ -94,5 +94,5 @@
              ]
 
-source = ["lib/sph_j1c.c", "vesicle.c"]
+source = ["lib/sas_3j1x_x.c", "vesicle.c"]
 
 def ER(radius, thickness):
