Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision 3f832f95a1224aa53015664f7faacbe6d8a3ce2d)
+++ sasmodels/models/ellipsoid.c	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
@@ -7,14 +7,9 @@
 double _ellipsoid_kernel(double q, double rpolar, double requatorial, double sin_alpha)
 {
-    double sn, cn;
     double ratio = rpolar/requatorial;
     const double u = q*requatorial*sqrt(1.0
                    + sin_alpha*sin_alpha*(ratio*ratio - 1.0));
-    SINCOS(u, sn, cn);
-    //const double f = ( u==0.0 ? 1.0 : 3.0*(sn-u*cn)/(u*u*u) );
-    const double usq = u*u;
-    const double f = (u < 1.e-1)
-        ? 1.0 + usq*(-3./30. + usq*(3./840. + usq*(-3./45360.)))// + qrsq*(3./3991680.))))
-        : 3.0*(sn/u - cn)/usq;
+    const double f = J1c(u);
+
     return f*f;
 }
Index: sasmodels/models/ellipsoid.py
===================================================================
--- sasmodels/models/ellipsoid.py	(revision cade6209f5b85c2af69cfe7e1237f41ffe36e1dc)
+++ sasmodels/models/ellipsoid.py	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
@@ -152,5 +152,5 @@
              ]
 
-source = ["lib/J1.c", "lib/gauss76.c", "ellipsoid.c"]
+source = ["lib/J1.c", "lib/J1c.c", "lib/gauss76.c", "ellipsoid.c"]
 
 def ER(rpolar, requatorial):
Index: sasmodels/models/gel_fit.py
===================================================================
--- sasmodels/models/gel_fit.py	(revision 30b4ddfb88f4c0a299757def2ee628c2a3d50a42)
+++ sasmodels/models/gel_fit.py	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
@@ -18,5 +18,6 @@
 .. math::
 
-    I(Q) = I(0)_L \frac{1}{\left( 1+\left[ ((D+1/3)Q^2a_{1}^2 \right]\right)^{D/2}} + I(0)_G exp\left( -Q^2a_{2}^2\right) + B
+    I(Q) = I(0)_L \frac{1}{\left( 1+\left[ ((D+1/3)Q^2a_{1}^2
+    \right]\right)^{D/2}} + I(0)_G exp\left( -Q^2a_{2}^2\right) + B
 
 where
@@ -26,7 +27,8 @@
     a_{2}^2 \approx \frac{R_{g}^2}{3}
 
-Note that the first term reduces to the Ornstein-Zernicke equation when $D = 2$;
-ie, when the Flory exponent is 0.5 (theta conditions).
-In gels with significant hydrogen bonding $D$ has been reported to be ~2.6 to 2.8.
+Note that the first term reduces to the Ornstein-Zernicke equation
+when $D = 2$; ie, when the Flory exponent is 0.5 (theta conditions).
+In gels with significant hydrogen bonding $D$ has been reported to be
+~2.6 to 2.8.
 
 
@@ -38,7 +40,9 @@
 ---------
 
-Mitsuhiro Shibayama, Toyoichi Tanaka, Charles C Han, J. Chem. Phys. 1992, 97 (9), 6829-6841
+Mitsuhiro Shibayama, Toyoichi Tanaka, Charles C Han, J. Chem. Phys. 1992, 97 (9),
+6829-6841
 
-Simon Mallam, Ferenc Horkay, Anne-Marie Hecht, Adrian R Rennie, Erik Geissler, Macromolecules 1991, 24, 543-548
+Simon Mallam, Ferenc Horkay, Anne-Marie Hecht, Adrian R Rennie, Erik Geissler,
+Macromolecules 1991, 24, 543-548
 
 """
Index: sasmodels/models/lib/J1c.c
===================================================================
--- sasmodels/models/lib/J1c.c	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
+++ sasmodels/models/lib/J1c.c	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
@@ -0,0 +1,21 @@
+// 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/polymer_excl_volume.py
===================================================================
--- sasmodels/models/polymer_excl_volume.py	(revision 23d3b41b22d87ce9f10b06dde71e4b25ffbe7197)
+++ sasmodels/models/polymer_excl_volume.py	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
@@ -1,10 +1,11 @@
 r"""
-This model describes the scattering from polymer chains subject to excluded volume effects
-and has been used as a template for describing mass fractals.
+This model describes the scattering from polymer chains subject to excluded
+volume effects and has been used as a template for describing mass fractals.
 
 Definition
 ----------
 
-The form factor was originally presented in the following integral form (Benoit, 1957)
+The form factor was originally presented in the following integral form
+(Benoit, 1957)
 
 .. math::
@@ -16,5 +17,6 @@
 $a$ is the statistical segment length of the polymer chain,
 and $n$ is the degree of polymerization.
-This integral was later put into an almost analytical form as follows (Hammouda, 1993)
+This integral was later put into an almost analytical form as follows
+(Hammouda, 1993)
 
 .. math::
@@ -43,14 +45,16 @@
 Note that this model applies only in the mass fractal range (ie, $5/3<=m<=3$ )
 and **does not apply** to surface fractals ( $3<m<=4$ ).
-It also does not reproduce the rigid rod limit (m=1) because it assumes chain flexibility
-from the outset. It may cover a portion of the semi-flexible chain range ( $1<m<5/3$ ).
+It also does not reproduce the rigid rod limit (m=1) because it assumes chain
+flexibility from the outset. It may cover a portion of the semi-flexible chain
+range ( $1<m<5/3$ ).
 
-A low-Q expansion yields the Guinier form and a high-Q expansion yields the Porod form
-which is given by
+A low-Q expansion yields the Guinier form and a high-Q expansion yields the
+Porod form which is given by
 
 .. math::
 
-    P(Q\rightarrow \infty) = \frac{1}{\nu U^{1/2\nu}}\Gamma\left(\frac{1}{2\nu}\right) -
-    \frac{1}{\nu U^{1/\nu}}\Gamma\left(\frac{1}{\nu}\right)
+    P(Q\rightarrow \infty) = \frac{1}{\nu U^{1/2\nu}}\Gamma\left(
+    \frac{1}{2\nu}\right) - \frac{1}{\nu U^{1/\nu}}\Gamma\left(
+    \frac{1}{\nu}\right)
 
 Here $\Gamma(x) = \gamma(x,\infty)$ is the gamma function.
@@ -102,5 +106,6 @@
 name = "polymer_excl_volume"
 title = "Polymer Excluded Volume model"
-description = """Compute the scattering intensity from polymers with excluded volume effects.
+description = """Compute the scattering intensity from polymers with excluded
+                volume effects.
                 rg:         radius of gyration
                 porod_exp:  Porod exponent
@@ -127,6 +132,6 @@
 
     intensity = ((1.0/(nu*power(u, o2nu))) * (gamma(o2nu)*gammainc(o2nu, u) -
-                 1.0/power(u, o2nu) * gamma(porod_exp)*gammainc(porod_exp, u)))*(q > 0) + \
-                 1.0*(q <= 0)
+                  1.0/power(u, o2nu) * gamma(porod_exp) *
+                  gammainc(porod_exp, u))) * (q > 0) + 1.0*(q <= 0)
 
     return intensity
@@ -156,4 +161,5 @@
          [{'rg': 2.2, 'porod_exp': 22.0, 'background': 100.0}, 5.0, 100.0],
 
-         [{'rg': 1.1, 'porod_exp': 1, 'background': 10.0, 'scale': 1.25}, 20000., 10.0000712097]
+         [{'rg': 1.1, 'porod_exp': 1, 'background': 10.0, 'scale': 1.25},
+         20000., 10.0000712097]
          ]
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision d53d3cd773e86c4090a2a101b72ec7579ccbf132)
+++ sasmodels/models/sphere.py	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
@@ -78,4 +78,5 @@
              ]
 
+source = ["lib/J1c.c"]
 
 # No volume normalization despite having a volume parameter
@@ -87,18 +88,5 @@
 Iq = """
     const double qr = q*radius;
-    const double qrsq = qr*qr;
-    double sn, cn;
-    SINCOS(qr, sn, cn);
-    // Use taylor series for low q to avoid cancellation error.  Tested against
-    // the following expression in quad precision:
-    //     3.0*(sn-qr*cn)/(qr*qr*qr);
-    // 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 below.
-    const double bes = (qr < 1.e-1)
-        ? 1.0 + qrsq*(-3./30. + qrsq*(3./840. + qrsq*(-3./45360.)))// + qrsq*(3./3991680.))))
-        : 3.0*(sn/qr - cn)/qrsq;
+    const double bes = 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 da84551fb2482445240877cb0f409b5810adaa3e)
+++ sasmodels/models/surface_fractal.c	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
@@ -1,6 +1,13 @@
 double form_volume(double radius);
 
-double Iq(double q, double radius, double surface_dim, double cutoff_length);
-double Iqxy(double qx, double qy, double radius, double surface_dim, double cutoff_length);
+double Iq(double q,
+          double radius,
+          double surface_dim,
+          double cutoff_length);
+
+double Iqxy(double qx, double qy,
+            double radius,
+            double surface_dim,
+            double cutoff_length);
 
 static double _gamln(double q)
@@ -26,5 +33,5 @@
 }
 
-static double surface_fractal_kernel(double q,
+static double _surface_fractal_kernel(double q,
     double radius,
     double surface_dim,
@@ -33,6 +40,9 @@
     double pq, sq, mmo, result;
 
-    //calculate P(q) for the spherical subunits; not normalized
-	pq = pow((3.0*(sin(q*radius) - q*radius*cos(q*radius))/pow((q*radius),3)),2);
+    //Replaced the original formula with Taylor expansion near zero.
+    //pq = pow((3.0*(sin(q*radius) - q*radius*cos(q*radius))/pow((q*radius),3)),2);
+
+    pq = J1c(q*radius);
+    pq = pq*pq;
 
     //calculate S(q)
@@ -59,5 +69,5 @@
     )
 {
-    return surface_fractal_kernel(q, radius, surface_dim, cutoff_length);
+    return _surface_fractal_kernel(q, radius, surface_dim, cutoff_length);
 }
 
@@ -69,5 +79,5 @@
 {
     double q = sqrt(qx*qx + qy*qy);
-    return surface_fractal_kernel(q, radius, surface_dim, cutoff_length);
+    return _surface_fractal_kernel(q, radius, surface_dim, cutoff_length);
 }
 
Index: sasmodels/models/surface_fractal.py
===================================================================
--- sasmodels/models/surface_fractal.py	(revision 30b4ddfb88f4c0a299757def2ee628c2a3d50a42)
+++ sasmodels/models/surface_fractal.py	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
@@ -22,5 +22,6 @@
 .. math::
 
-    S(q) = \frac{\Gamma(5-D_S)\zeta^{5-D_S}}{\left[1+(q\zeta)^2\right]^{(5-D_S)/2}}
+    S(q) = \frac{\Gamma(5-D_S)\zeta^{5-D_S}}{\left[1+(q\zeta)^2
+    \right]^{(5-D_S)/2}}
     \frac{sin\left[(D_S - 5) tan^{-1}(q\zeta) \right]}{q}
 
@@ -33,11 +34,13 @@
     V = \frac{4}{3}\pi R^3
 
-where $R$ is the radius of the building block, $D_S$ is the **surface** fractal dimension,
-$\zeta$ is the cut-off length, $\rho_{solvent}$ is the scattering length density of the solvent,
+where $R$ is the radius of the building block, $D_S$ is the **surface** fractal
+dimension,$\zeta$ is the cut-off length, $\rho_{solvent}$ is the scattering
+length density of the solvent,
 and $\rho_{particle}$ is the scattering length density of particles.
 
 .. note::
     The surface fractal dimension $D_s$ is only valid if $1<surface\_dim<3$.
-    It is also only valid over a limited $q$ range (see the reference for details)
+    It is also only valid over a limited $q$ range (see the reference for
+    details)
 
 
@@ -78,11 +81,14 @@
 
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["radius",        "Ang", 10.0, [0, inf],   "", "Particle radius"],
-              ["surface_dim",   "",    2.0,  [0, inf],   "", "Surface fractal dimension"],
-              ["cutoff_length", "Ang", 500., [0.0, inf], "",  "Cut-off Length"],
+parameters = [["radius",        "Ang", 10.0, [0, inf],   "",
+               "Particle radius"],
+              ["surface_dim",   "",    2.0,  [0, inf],   "",
+               "Surface fractal dimension"],
+              ["cutoff_length", "Ang", 500., [0.0, inf], "",
+               "Cut-off Length"],
               ]
 
 
-source = ["surface_fractal.c"]
+source = ["lib/J1c.c", "surface_fractal.c"]
 
 demo = dict(scale=1, background=0,
Index: sasmodels/models/triaxial_ellipsoid.c
===================================================================
--- sasmodels/models/triaxial_ellipsoid.c	(revision d53d3cd773e86c4090a2a101b72ec7579ccbf132)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
@@ -36,10 +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);
-            SINCOS(t, st, ct);
-            //const double fq = ( t==0.0 ? 1.0 : 3.0*(st-t*ct)/(t*t*t) );
-            const double tsq = t*t;
-            const double fq = (t < 1.e-1)
-                ? 1.0 + tsq*(-3./30. + tsq*(3./840. + tsq*(-3./45360.)))// + tsq*(3./3991680.))))
-                : 3.0*(st/t - ct)/tsq;
+            const double fq = J1c(t);
             inner += Gauss76Wt[j] * fq * fq ;
         }
Index: sasmodels/models/triaxial_ellipsoid.py
===================================================================
--- sasmodels/models/triaxial_ellipsoid.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/triaxial_ellipsoid.py	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
@@ -117,5 +117,5 @@
              ]
 
-source = ["lib/J1.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
+source = ["lib/J1.c", "lib/J1c.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
 
 def ER(req_minor, req_major, rpolar):
Index: sasmodels/models/two_lorentzian.py
===================================================================
--- sasmodels/models/two_lorentzian.py	(revision 421e55c8d6051fdab7b8bee21bd622416a29c8b4)
+++ sasmodels/models/two_lorentzian.py	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
@@ -1,4 +1,5 @@
 r"""
-This model calculates an empirical functional form for SAS data characterized by two Lorentzian-type functions.
+This model calculates an empirical functional form for SAS data characterized
+by two Lorentzian-type functions.
 
 Definition
@@ -12,6 +13,7 @@
 
 where $A$ = Lorentzian scale factor #1, $C$ = Lorentzian scale #2,
-$\xi_1$ and $\xi_2$ are the corresponding correlation lengths, and $n$ and $m$ are the respective
-power law exponents (set $n = m = 2$ for Ornstein-Zernicke behaviour).
+$\xi_1$ and $\xi_2$ are the corresponding correlation lengths, and $n$ and
+$m$ are the respective power law exponents (set $n = m = 2$ for
+Ornstein-Zernicke behaviour).
 
 For 2D data the scattering intensity is calculated in the same way as 1D,
@@ -52,15 +54,27 @@
 category = "shape-independent"
 
-#             ["name", "units", default, [lower, upper], "type", "description"],
-parameters = [["lorentz_scale_1",  "",     10.0, [-inf, inf], "", "First power law scale factor"],
-              ["lorentz_length_1", "Ang", 100.0, [-inf, inf], "", "First Lorentzian screening length"],
-              ["lorentz_exp_1",    "",      3.0, [-inf, inf], "", "First exponent of power law"],
-              ["lorentz_scale_2",  "",      1.0, [-inf, inf], "", "Second scale factor for broad Lorentzian peak"],
-              ["lorentz_length_2", "Ang",  10.0, [-inf, inf], "", "Second Lorentzian screening length"],
-              ["lorentz_exp_2",    "",      2.0, [-inf, inf], "", "Second exponent of power law"],
+#            ["name", "units", default, [lower, upper], "type", "description"],
+parameters = [["lorentz_scale_1",  "",     10.0, [-inf, inf], "",
+               "First power law scale factor"],
+              ["lorentz_length_1", "Ang", 100.0, [-inf, inf], "",
+               "First Lorentzian screening length"],
+              ["lorentz_exp_1",    "",      3.0, [-inf, inf], "",
+               "First exponent of power law"],
+              ["lorentz_scale_2",  "",      1.0, [-inf, inf], "",
+               "Second scale factor for broad Lorentzian peak"],
+              ["lorentz_length_2", "Ang",  10.0, [-inf, inf], "",
+               "Second Lorentzian screening length"],
+              ["lorentz_exp_2",    "",      2.0, [-inf, inf], "",
+               "Second exponent of power law"],
               ]
 
 
-def Iq(q, lorentz_scale_1, lorentz_length_1, lorentz_exp_1, lorentz_scale_2, lorentz_length_2, lorentz_exp_2):
+def Iq(q,
+       lorentz_scale_1,
+       lorentz_length_1,
+       lorentz_exp_1,
+       lorentz_scale_2,
+       lorentz_length_2,
+       lorentz_exp_2):
 
     """
@@ -75,6 +89,8 @@
     """
 
-    intensity  = lorentz_scale_1/(1.0 + power(q*lorentz_length_1, lorentz_exp_1))
-    intensity += lorentz_scale_2/(1.0 + power(q*lorentz_length_2, lorentz_exp_2))
+    intensity  = lorentz_scale_1/(1.0 +
+                                  power(q*lorentz_length_1, lorentz_exp_1))
+    intensity += lorentz_scale_2/(1.0 +
+                                  power(q*lorentz_length_2, lorentz_exp_2))
 
     return intensity
@@ -101,15 +117,27 @@
                lorentz_exp_1='exponent_1',  lorentz_exp_2='exponent_2')
 
-tests = [[{'lorentz_scale_1': 10, 'lorentz_length_1': 100.0, 'lorentz_exp_1' : 3.0,
-          'lorentz_scale_2': 1, 'lorentz_length_2': 10.0, 'lorentz_exp_2' : 2.0,
+tests = [[{'lorentz_scale_1':   10.0,
+           'lorentz_length_1': 100.0,
+           'lorentz_exp_1':      3.0,
+           'lorentz_scale_2':    1.0,
+           'lorentz_length_2':  10.0,
+           'lorentz_exp_2':      2.0,
            }, 0.000332070182643, 10.9996228107],
 
-         [{'lorentz_scale_1': 0, 'lorentz_length_1': 0.0, 'lorentz_exp_1' : 0.0,
-          'lorentz_scale_2': 0, 'lorentz_length_2': 0.0, 'lorentz_exp_2' : 0.0,
-           'background':100.
+         [{'lorentz_scale_1':  0.0,
+           'lorentz_length_1': 0.0,
+           'lorentz_exp_1':    0.0,
+           'lorentz_scale_2':  0.0,
+           'lorentz_length_2': 0.0,
+           'lorentz_exp_2':    0.0,
+           'background':     100.0
            }, 5.0, 100.0],
 
-         [{'lorentz_scale_1': 200, 'lorentz_length_1': 10.0, 'lorentz_exp_1' : 0.1,
-          'lorentz_scale_2': 0.1, 'lorentz_length_2': 5.0, 'lorentz_exp_2' : 2.0
+         [{'lorentz_scale_1': 200.0,
+           'lorentz_length_1': 10.0,
+           'lorentz_exp_1':     0.1,
+           'lorentz_scale_2':   0.1,
+           'lorentz_length_2':  5.0,
+           'lorentz_exp_2':     2.0
            }, 20000., 45.5659201896],
          ]
