Index: sasmodels/models/bcc.c
===================================================================
--- sasmodels/models/bcc.c	(revision 82d239a5cb4d8005a4b04bf42b1425792a28383d)
+++ sasmodels/models/bcc.c	(revision 7ed702f2d0d4ecbcb0115970cdabf0d05ac08de0)
@@ -100,5 +100,5 @@
 
   double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
-  double q_z;
+  //double q_z;
   double cos_val_b3, cos_val_b2, cos_val_b1;
   double a1_dot_q, a2_dot_q,a3_dot_q;
@@ -123,5 +123,5 @@
   const double latticescale = 2.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1);
   // q vector
-  q_z = 0.0; // for SANS; assuming qz is negligible
+  //q_z = 0.0; // for SANS; assuming qz is negligible
   /// Angles here are respect to detector coordinate
   ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
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/hardsphere.py
===================================================================
--- sasmodels/models/hardsphere.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/hardsphere.py	(revision 7ed702f2d0d4ecbcb0115970cdabf0d05ac08de0)
@@ -23,5 +23,5 @@
 
 
-.. figure:: img/HardSphere_1d.jpg
+.. figure:: img/hardSphere_1d.jpg
 
     1D plot using the default values (in linear scale).
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/lib/lanczos_gamma.c
===================================================================
--- sasmodels/models/lib/lanczos_gamma.c	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
+++ sasmodels/models/lib/lanczos_gamma.c	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
@@ -0,0 +1,21 @@
+double lanczos_gamma(double q);
+double lanczos_gamma(double q)
+{
+    // Lanczos approximation to the Gamma function.
+
+    double x,y,tmp,ser;
+    double coeff[6]=
+        {76.18009172947146,     -86.50532032941677,
+         24.01409824083091,     -1.231739572450155,
+          0.1208650973866179e-2,-0.5395239384953e-5};
+
+    y=x=q;
+    tmp  = x+5.5;
+    tmp -= (x+0.5)*log(tmp);
+    ser  = 1.000000000190015;
+    for (int j=0; j<=5; j++) {
+        y+=1.0;
+        ser += coeff[j]/y;
+    }
+    return -tmp+log(2.5066282746310005*ser/x);
+}
Index: sasmodels/models/mass_fractal.c
===================================================================
--- sasmodels/models/mass_fractal.c	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
+++ sasmodels/models/mass_fractal.c	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
@@ -0,0 +1,68 @@
+double form_volume(double radius);
+
+double Iq(double q,
+          double radius,
+          double mass_dim,
+          double cutoff_length);
+
+double Iqxy(double qx, double qy,
+          double radius,
+          double mass_dim,
+          double cutoff_length);
+
+
+static double _mass_fractal_kernel(double q,
+          double radius,
+          double mass_dim,
+          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 S(q)
+    double mmo = mass_dim-1.0;
+    double sq = exp(lanczos_gamma(mmo))*sin((mmo)*atan(q*cutoff_length));
+    sq *= pow(cutoff_length, mmo);
+    sq /= pow((1.0 + (q*cutoff_length)*(q*cutoff_length)),(mmo/2.0));
+    sq /= q;
+
+    //sq = 1.0;
+
+    //combine and return
+    double result = pq * sq;
+
+    return result;
+}
+double form_volume(double radius){
+
+    return 1.333333333333333*M_PI*radius*radius*radius;
+}
+
+double Iq(double q,
+          double radius,
+          double mass_dim,
+          double cutoff_length)
+{
+    return _mass_fractal_kernel(q,
+           radius,
+           mass_dim,
+           cutoff_length);
+}
+
+// Iqxy is never called since no orientation or magnetic parameters.
+double Iqxy(double qx, double qy,
+          double radius,
+          double mass_dim,
+          double cutoff_length)
+{
+    double q = sqrt(qx*qx + qy*qy);
+    return _mass_fractal_kernel(q,
+           radius,
+           mass_dim,
+           cutoff_length);
+}
+
Index: sasmodels/models/mass_fractal.py
===================================================================
--- sasmodels/models/mass_fractal.py	(revision 9eaadb373d3543d4f7dda6dafe056914c2f74b6b)
+++ sasmodels/models/mass_fractal.py	(revision 9eaadb373d3543d4f7dda6dafe056914c2f74b6b)
@@ -0,0 +1,119 @@
+r"""
+Calculates the scattering from fractal-like aggregates based on
+the Mildner reference.
+
+Definition
+----------
+
+The scattering intensity $I(q)$ is calculated as
+
+.. math::
+
+    I(q) = scale \times P(q)S(q) + background
+
+.. math::
+
+    P(q) = F(qR)^2
+
+.. math::
+
+    F(x) = \frac{3\left[sin(x)-xcos(x)\right]}{x^3}
+
+.. math::
+
+    S(q) = \frac{\Gamma(D_m-1)\zeta^{D_m-1}}{\left[1+(q\zeta)^2
+    \right]^{(D_m-1)/2}}
+    \frac{sin\left[(D_m - 1) tan^{-1}(q\zeta) \right]}{q}
+
+.. math::
+
+    scale = scale\_factor \times NV^2(\rho_{particle} - \rho_{solvent})^2
+
+.. math::
+
+    V = \frac{4}{3}\pi R^3
+
+where $R$ is the radius of the building block, $D_m$ is the **mass** 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 mass fractal dimension ( $D_m$ ) is only
+    valid if $0 < mass_dim < 6$. It is also only valid over a limited
+    $q$ range (see the reference for details).
+
+.. figure:: img/mass_fractal_1d.jpg
+
+    1D plot using the default values.
+
+Reference
+---------
+
+D Mildner and P Hall, *J. Phys. D: Appl. Phys.*,  19 (1986) 1535-1545 Equation(9)
+
+
+"""
+
+from numpy import inf
+
+name = "mass_fractal"
+title = "Mass Fractal model"
+description = """
+        The scattering intensity  I(x) = scale*P(x)*S(x) + background, where
+        scale = scale_factor  * V * delta^(2)
+        p(x)=  F(x*radius)^(2)
+        F(x) = 3*[sin(x)-x cos(x)]/x**3
+        S(x) = [(gamma(Dm-1)*colength^(Dm-1)*[1+(x^2*colength^2)]^((1-Dm)/2)
+        * sin[(Dm-1)*arctan(x*colength)])/x]
+        where delta = sldParticle -sldSolv.
+        radius       =  Particle radius
+        mass_dim  =  Mass fractal dimension
+        cutoff_length  =  Cut-off length
+        background   =  background
+        Ref.:Mildner, Hall,J Phys D Appl Phys(1986), 9, 1535-1545
+        Note I: This model is valid for 1<mass_dim<6.
+        Note II: This model is not in absolute scale.
+        """
+category = "shape-independent"
+
+#             ["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"],
+              ]
+
+
+source = ["lib/J1c.c", "lib/lanczos_gamma.c", "mass_fractal.c"]
+
+demo = dict(scale=1, background=0,
+            radius=10.0,
+            mass_dim=1.9,
+            cutoff_length=2.3)
+
+oldname = 'MassFractalModel'
+oldpars = dict(radius='radius',
+               mass_dim='mass_dim',
+               cutoff_length='co_length')
+
+tests = [[{'radius':        2.0,
+           'mass_dim':      3.3,
+           'cutoff_length': 1.0,
+           }, 0.001, 2.68342996226],
+
+         [{'radius':        1.0,
+           'mass_dim':      1.3,
+           'cutoff_length': 1.0,
+           'background':    0.8,
+           }, 0.001, 1.69747015932],
+
+         [{'radius':        1.0,
+           'mass_dim':      2.3,
+           'cutoff_length': 1.0,
+           'scale':        10.0,
+           }, 0.051, 11.6258202095],
+         ]
Index: sasmodels/models/mass_surface_fractal.c
===================================================================
--- sasmodels/models/mass_surface_fractal.c	(revision 7ed702f2d0d4ecbcb0115970cdabf0d05ac08de0)
+++ sasmodels/models/mass_surface_fractal.c	(revision 7ed702f2d0d4ecbcb0115970cdabf0d05ac08de0)
@@ -0,0 +1,76 @@
+double form_volume(double radius);
+
+double Iq(double q,
+          double mass_dim,
+          double surface_dim,
+          double cluster_rg,
+          double primary_rg);
+
+double Iqxy(double qx, double qy,
+          double mass_dim,
+          double surface_dim,
+          double cluster_rg,
+          double primary_rg);
+
+
+static double _mass_surface_fractal_kernel(double q,
+          double mass_dim,
+          double surface_dim,
+          double cluster_rg,
+          double primary_rg)
+{
+     //computation
+    double tot_dim = 6.0 - surface_dim - mass_dim;
+    mass_dim /= 2.0;
+    tot_dim /= 2.0;
+
+    double rc_norm = cluster_rg * cluster_rg / (3.0 * mass_dim);
+    double rp_norm = primary_rg * primary_rg / (3.0 * tot_dim);
+
+    //x for P
+    double x_val1 = 1.0 +  q * q * rc_norm;
+    double x_val2 = 1.0 +  q * q * rp_norm;
+
+    double inv_form = pow(x_val1, mass_dim) * pow(x_val2, tot_dim);
+
+    //another singular
+    if (inv_form == 0.0) return 0.0;
+
+    double form_factor = 1.0;
+    form_factor /= inv_form;
+
+    return (form_factor);
+}
+double form_volume(double radius){
+
+    return 1.333333333333333*M_PI*radius*radius*radius;
+}
+
+double Iq(double q,
+          double mass_dim,
+          double surface_dim,
+          double cluster_rg,
+          double primary_rg)
+{
+    return _mass_surface_fractal_kernel(q,
+            mass_dim,
+            surface_dim,
+            cluster_rg,
+            primary_rg);
+}
+
+// Iqxy is never called since no orientation or magnetic parameters.
+double Iqxy(double qx, double qy,
+          double mass_dim,
+          double surface_dim,
+          double cluster_rg,
+          double primary_rg)
+{
+    double q = sqrt(qx*qx + qy*qy);
+    return _mass_surface_fractal_kernel(q,
+           mass_dim,
+           surface_dim,
+           cluster_rg,
+           primary_rg);
+}
+
Index: sasmodels/models/mass_surface_fractal.py
===================================================================
--- sasmodels/models/mass_surface_fractal.py	(revision 7ed702f2d0d4ecbcb0115970cdabf0d05ac08de0)
+++ sasmodels/models/mass_surface_fractal.py	(revision 7ed702f2d0d4ecbcb0115970cdabf0d05ac08de0)
@@ -0,0 +1,130 @@
+r"""
+
+A number of natural and commercial processes form high-surface area materials
+as a result of the vapour-phase aggregation of primary particles.
+Examples of such materials include soots, aerosols, and fume or pyrogenic
+silicas. These are all characterised by cluster mass distributions (sometimes
+also cluster size distributions) and internal surfaces that are fractal in
+nature. The scattering from such materials displays two distinct breaks in
+log-log representation, corresponding to the radius-of-gyration of the primary
+particles, $rg$, and the radius-of-gyration of the clusters (aggregates),
+$Rg$. Between these boundaries the scattering follows a power law related to
+the mass fractal dimension, $Dm$, whilst above the high-Q boundary the
+scattering follows a power law related to the surface fractal dimension of
+the primary particles, $Ds$.
+
+Definition
+----------
+
+The scattered intensity I(q) is calculated using a modified
+Ornstein-Zernicke equation
+
+.. math::
+
+    I(q) = scale \times P(q) + background
+
+    P(q) = \left\{ \left[ 1+(q^2a)\right]^{D_m/2} \times
+                   \left[ 1+(q^2b)\right]^{(6-D_s-D_m)/2}
+           \right\}^{-1}
+
+    a = R_{g}^2/(3D_m/2)
+
+    b = r_{g}^2/[-3(D_s+D_m-6)/2]
+
+    scale = scale\_factor \times NV^2 (\rho_{particle} - \rho_{solvent})^2
+
+where $R_g$ is the size of the cluster, $r_g$ is the size of the primary
+particle, $D_s$ is the surface fractal dimension, $D_m$ is the mass fractal
+dimension, $\rho_{solvent}$ is the scattering length density of the solvent,
+and $\rho_{particle}$ is the scattering length density of particles.
+
+.. note::
+
+    The surface ( $D_s$ ) and mass ( $D_m$ ) fractal dimensions are only
+    valid if $0 < surface_dim < 6$ , $0 < mass_dim < 6$ , and
+    $(surface_dim + mass_dim ) < 6$ .
+
+.. figure:: img/mass_surface_fractal_1d.jpg
+
+    1D plot using the default values.
+
+Reference
+---------
+
+P Schmidt, *J Appl. Cryst.*, 24 (1991) 414-435 Equation(19)
+
+A J Hurd, D W Schaefer, J E Martin, *Phys. Rev. A*, 35 (1987) 2361-2364 Equation(2)
+
+"""
+
+from numpy import inf
+
+name = "mass_surface_fractal"
+title = "Mass Surface Fractal model"
+description = """
+        The scattering intensity  I(x) = scale*P(x)*S(x) + background, where
+        p(x)= {[1+(x^2*a)]^(Dm/2) * [1+(x^2*b)]^(6-Ds-Dm)/2}^(-1)
+        a = Rg^2/(3*Dm/2)
+        b = rg^2/(3*(6-Ds-Dm)/2)
+        scale        =  scale factor * N*Volume^2*contrast^2
+        mass_dim       =  Dm (mass fractal dimension)
+        surface_dim  =  Ds
+        cluster_rg  =  Rg
+        primary_rg    =  rg
+        background   =  background
+        Ref: Schmidt, J Appl Cryst, eq(19), (1991), 24, 414-435
+        Hurd, Schaefer, Martin, Phys Rev A, eq(2),(1987),35, 2361-2364
+        Note that 0 < Ds< 6 and 0 < Dm < 6.
+        """
+category = "shape-independent"
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["mass_dim",      "",    1.8, [1e-16, 6.0], "",
+               "Mass fractal dimension"],
+              ["surface_dim",   "",    2.3, [1e-16, 6.0], "",
+               "Surface fractal dimension"],
+              ["cluster_rg", "Ang",   86.7, [0.0, inf], "",
+               "Cluster radius of gyration"],
+              ["primary_rg", "Ang", 4000.,  [0.0, inf], "",
+               "Primary particle radius of gyration"],
+              ]
+
+
+source = ["lib/J1c.c", "mass_surface_fractal.c"]
+
+demo = dict(scale=1, background=0,
+            mass_dim=1.8,
+            surface_dim=2.3,
+            cluster_rg=86.7,
+            primary_rg=4000.0)
+
+oldname = 'MassSurfaceFractal'
+oldpars = dict(radius='radius',
+               mass_dim='mass_dim',
+               surface_dim='surface_dim',
+               cluster_rg='cluster_rg',
+               primary_rg='primary_rg')
+
+tests = [[{'radius':        2.0,
+           'mass_dim':      3.3,
+           'surface_dim':   1.0,
+           'cluster_rg':   90.0,
+           'primary_rg': 4000.0,
+           }, 0.001, 0.18462699016],
+
+         [{'radius':        1.0,
+           'mass_dim':      1.3,
+           'surface_dim':   1.0,
+           'cluster_rg':   90.0,
+           'primary_rg': 2000.0,
+           'background':    0.8,
+           }, 0.001, 1.16539753641],
+
+         [{'radius':        1.0,
+           'mass_dim':      2.3,
+           'surface_dim':   1.0,
+           'cluster_rg':   90.0,
+           'primary_rg': 1000.0,
+           'scale':        10.0,
+           }, 0.051, 0.000169548800377],
+         ]
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 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
@@ -1,30 +1,15 @@
 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);
 
-static double _gamln(double q)
-{
-    // Lanczos approximation to the Gamma function.
-    // Should be refactored out to lib/, if used elsewhere.
-    double x,y,tmp,ser;
-    double coeff[6]=
-        {76.18009172947146,     -86.50532032941677,
-         24.01409824083091,     -1.231739572450155,
-          0.1208650973866179e-2,-0.5395239384953e-5};
-    int j;
+double Iqxy(double qx, double qy,
+            double radius,
+            double surface_dim,
+            double cutoff_length);
 
-    y=x=q;
-    tmp  = x+5.5;
-    tmp -= (x+0.5)*log(tmp);
-    ser  = 1.000000000190015;
-    for (j=0; j<=5; j++) {
-        y+=1.0;
-        ser += coeff[j]/y;
-    }
-    return -tmp+log(2.5066282746310005*ser/x);
-}
-
-static double surface_fractal_kernel(double q,
+static double _surface_fractal_kernel(double q,
     double radius,
     double surface_dim,
@@ -33,10 +18,13 @@
     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)
     mmo = 5.0 - surface_dim;
-    sq  = exp(_gamln(mmo))*sin(-(mmo)*atan(q*cutoff_length));
+    sq  = exp(lanczos_gamma(mmo))*sin(-(mmo)*atan(q*cutoff_length));
     sq *= pow(cutoff_length, mmo);
     sq /= pow((1.0 + (q*cutoff_length)*(q*cutoff_length)),(mmo/2.0));
@@ -59,5 +47,5 @@
     )
 {
-    return surface_fractal_kernel(q, radius, surface_dim, cutoff_length);
+    return _surface_fractal_kernel(q, radius, surface_dim, cutoff_length);
 }
 
@@ -69,5 +57,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 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
@@ -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", "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 d53d3cd773e86c4090a2a101b72ec7579ccbf132)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision 7ed702f2d0d4ecbcb0115970cdabf0d05ac08de0)
@@ -20,5 +20,5 @@
 
     double sn, cn;
-    double st, ct;
+    //double st, ct;
     //const double lower = 0.0;
     //const double upper = 1.0;
@@ -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],
          ]
