Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
+++ sasmodels/models/ellipsoid.c	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
@@ -10,5 +10,5 @@
     const double u = q*requatorial*sqrt(1.0
                    + sin_alpha*sin_alpha*(ratio*ratio - 1.0));
-    const double f = J1c(u);
+    const double f = sph_j1c(u);
 
     return f*f;
Index: sasmodels/models/ellipsoid.py
===================================================================
--- sasmodels/models/ellipsoid.py	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
+++ sasmodels/models/ellipsoid.py	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
@@ -152,5 +152,5 @@
              ]
 
-source = ["lib/J1.c", "lib/J1c.c", "lib/gauss76.c", "ellipsoid.c"]
+source = ["lib/J1.c", "lib/sph_j1c.c", "lib/gauss76.c", "ellipsoid.c"]
 
 def ER(rpolar, requatorial):
Index: smodels/models/lib/J1c.c
===================================================================
--- sasmodels/models/lib/J1c.c	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
+++ 	(revision )
@@ -1,21 +1,0 @@
-// Use taylor series 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.  Single precision only
-// requires the first 3 terms.  Double precision requires the 4th term.
-// The fifth term is not needed, and is commented out.
-double J1c(double qr);
-double J1c(double qr)
-{
-    const double qr2 = qr*qr;
-    double sn, cn;
-
-    SINCOS(qr, sn, cn);
-
-    const double bes = (qr < 1.e-1) ?
-        1.0 + qr2*(-3./30. + qr2*(3./840. + qr2*(-3./45360.)))
-        // + qr2*(3./3991680.)))
-        : 3.0*(sn/qr - cn)/qr2;
-
-    return bes;
-}
Index: sasmodels/models/lib/sph_j1c.c
===================================================================
--- sasmodels/models/lib/sph_j1c.c	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
+++ sasmodels/models/lib/sph_j1c.c	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
@@ -0,0 +1,25 @@
+/**
+* Spherical Bessel function 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.  Single precision only
+* requires the first 3 terms.  Double precision requires the 4th term.
+* The fifth term is not needed, and is commented out.
+*/
+
+double sph_j1c(double q);
+double sph_j1c(double q)
+{
+    const double q2 = q*q;
+    double sin_q, cos_q;
+
+    SINCOS(q, sin_q, cos_q);
+
+    const double bessel = (q < 1.e-1) ?
+        1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))) // + q2*(3./3991680.)))
+        : 3.0*(sin_q/q - cos_q)/q2;
+
+    return bessel;
+}
Index: sasmodels/models/mass_fractal.c
===================================================================
--- sasmodels/models/mass_fractal.c	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
+++ sasmodels/models/mass_fractal.c	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
@@ -17,10 +17,7 @@
           double cutoff_length)
 {
-
-    //calculate P(q) for the spherical subunits; not normalized
-    //double pq = pow((3.0*(sin(q*radius) - q*radius*cos(q*radius))/pow((q*radius),3)),2);
-    double pq = J1c(q*radius);
-
-    //pq = 1.0;
+    //calculate P(q)
+    double pq = sph_j1c(q*radius);
+    pq = pq*pq;
 
     //calculate S(q)
@@ -30,6 +27,4 @@
     sq /= pow((1.0 + (q*cutoff_length)*(q*cutoff_length)),(mmo/2.0));
     sq /= q;
-
-    //sq = 1.0;
 
     //combine and return
Index: sasmodels/models/mass_fractal.py
===================================================================
--- sasmodels/models/mass_fractal.py	(revision 9eaadb373d3543d4f7dda6dafe056914c2f74b6b)
+++ sasmodels/models/mass_fractal.py	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
@@ -80,14 +80,11 @@
 
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["radius",        "Ang",  10.0, [0.0, inf], "",
-               "Particle radius"],
-              ["mass_dim",      "",      1.9, [1.0, 6.0],   "",
-               "Mass fractal dimension"],
-              ["cutoff_length", "Ang", 100.0, [0.0, inf], "",
-               "Cut-off length"],
+parameters = [["radius",        "Ang",  10.0, [0.0, inf], "", "Particle radius"],
+              ["mass_dim",      "",      1.9, [1.0, 6.0], "", "Mass fractal dimension"],
+              ["cutoff_length", "Ang", 100.0, [0.0, inf], "", "Cut-off length"],
               ]
 
 
-source = ["lib/J1c.c", "lib/lanczos_gamma.c", "mass_fractal.c"]
+source = ["lib/sph_j1c.c", "lib/lanczos_gamma.c", "mass_fractal.c"]
 
 demo = dict(scale=1, background=0,
@@ -104,5 +101,5 @@
            'mass_dim':      3.3,
            'cutoff_length': 1.0,
-           }, 0.001, 2.68342996226],
+           }, 0.5, 1.29016774904],
 
          [{'radius':        1.0,
@@ -116,4 +113,4 @@
            'cutoff_length': 1.0,
            'scale':        10.0,
-           }, 0.051, 11.6258202095],
+           }, 0.051, 11.6227966145],
          ]
Index: sasmodels/models/mass_surface_fractal.py
===================================================================
--- sasmodels/models/mass_surface_fractal.py	(revision 7ed702f2d0d4ecbcb0115970cdabf0d05ac08de0)
+++ sasmodels/models/mass_surface_fractal.py	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
@@ -91,5 +91,5 @@
 
 
-source = ["lib/J1c.c", "mass_surface_fractal.c"]
+source = ["mass_surface_fractal.c"]
 
 demo = dict(scale=1, background=0,
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
+++ sasmodels/models/sphere.py	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
@@ -78,5 +78,5 @@
              ]
 
-source = ["lib/J1c.c"]
+source = ["lib/sph_j1c.c"]
 
 # No volume normalization despite having a volume parameter
@@ -88,5 +88,5 @@
 Iq = """
     const double qr = q*radius;
-    const double bes = J1c(qr);
+    const double bes = sph_j1c(qr);
     const double fq = bes * (sld - solvent_sld) * form_volume(radius);
     return 1.0e-4*fq*fq;
Index: sasmodels/models/surface_fractal.c
===================================================================
--- sasmodels/models/surface_fractal.c	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
+++ sasmodels/models/surface_fractal.c	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
@@ -21,5 +21,5 @@
     //pq = pow((3.0*(sin(q*radius) - q*radius*cos(q*radius))/pow((q*radius),3)),2);
 
-    pq = J1c(q*radius);
+    pq = sph_j1c(q*radius);
     pq = pq*pq;
 
Index: sasmodels/models/surface_fractal.py
===================================================================
--- sasmodels/models/surface_fractal.py	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
+++ sasmodels/models/surface_fractal.py	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
@@ -90,5 +90,5 @@
 
 
-source = ["lib/J1c.c", "lib/lanczos_gamma.c", "surface_fractal.c"]
+source = ["lib/sph_j1c.c", "lib/lanczos_gamma.c", "surface_fractal.c"]
 
 demo = dict(scale=1, background=0,
Index: sasmodels/models/triaxial_ellipsoid.c
===================================================================
--- sasmodels/models/triaxial_ellipsoid.c	(revision 7ed702f2d0d4ecbcb0115970cdabf0d05ac08de0)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
@@ -36,5 +36,5 @@
             const double y = 0.5*(Gauss76Z[j] + 1.0);
             const double t = q*sqrt(acosx2 + bsinx2*(1.0-y*y) + c2*y*y);
-            const double fq = J1c(t);
+            const double fq = sph_j1c(t);
             inner += Gauss76Wt[j] * fq * fq ;
         }
Index: sasmodels/models/triaxial_ellipsoid.py
===================================================================
--- sasmodels/models/triaxial_ellipsoid.py	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
+++ sasmodels/models/triaxial_ellipsoid.py	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
@@ -117,5 +117,5 @@
              ]
 
-source = ["lib/J1.c", "lib/J1c.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
+source = ["lib/J1.c", "lib/sph_j1c.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
 
 def ER(req_minor, req_major, rpolar):
