Index: sasmodels/models/barbell.c
===================================================================
--- sasmodels/models/barbell.c	(revision 3c97ff07f375d382f26608fa2873fc2ce6733fe0)
+++ sasmodels/models/barbell.c	(revision 50e1e40cf14d4ded28bf96e61a3cfaf35003b899)
@@ -7,23 +7,36 @@
 
 //barbell kernel - same as dumbell
-double _bell_kernel(double q, double h, double bell_radius,
-        double length, double sin_alpha, double cos_alpha);
-double _bell_kernel(double q, double h, double bell_radius,
-        double length, double sin_alpha, double cos_alpha)
+static double
+_bell_kernel(double q, double h, double bell_radius,
+             double half_length, double sin_alpha, double cos_alpha)
 {
+    // translate a point in [-1,1] to a point in [lower,upper]
     const double upper = 1.0;
-    const double lower = -1.0*h/bell_radius;
+    const double lower = h/bell_radius;
+    const double zm = 0.5*(upper-lower);
+    const double zb = 0.5*(upper+lower);
 
+    // cos term in integral is:
+    //    cos (q (R t - h + L/2) cos(alpha))
+    // so turn it into:
+    //    cos (m t + b)
+    // where:
+    //    m = q R cos(alpha)
+    //    b = q(L/2-h) cos(alpha)
+    const double m = q*bell_radius*cos_alpha; // cos argument slope
+    const double b = q*(half_length-h)*cos_alpha; // cos argument intercept
+    const double qrst = q*bell_radius*sin_alpha; // Q*R*sin(theta)
     double total = 0.0;
     for (int i = 0; i < 76; i++){
-        const double t = 0.5*(Gauss76Z[i]*(upper-lower)+upper+lower);
-        const double arg1 = q*cos_alpha*(bell_radius*t+h+length*0.5);
-        const double arg2 = q*bell_radius*sin_alpha*sqrt(1.0-t*t);
-        const double be = (arg2 == 0.0 ? 0.5 :J1(arg2)/arg2);
-        const double Fq = cos(arg1)*(1.0-t*t)*be;
+        const double t = Gauss76Z[i]*zm + zb;
+        const double radical = 1.0 - t*t;
+        const double bj = J1c(qrst*sqrt(radical));
+        const double Fq = cos(m*t + b) * radical * bj;
         total += Gauss76Wt[i] * Fq;
     }
-    const double integral = 0.5*(upper-lower)*total;
-    return 4.0*M_PI*bell_radius*bell_radius*bell_radius*integral;
+    // translate dx in [-1,1] to dx in [lower,upper]
+    const double integral = total*zm;
+    const double bell_Fq = 2*M_PI*cube(bell_radius)*integral;
+    return bell_Fq;
 }
 
@@ -34,90 +47,72 @@
 
     // bell radius should never be less than radius when this is called
-    const double hdist = sqrt(bell_radius*bell_radius - radius*radius);
-    const double p1 = 2.0*bell_radius*bell_radius*bell_radius/3.0;
-    const double p2 = bell_radius*bell_radius*hdist;
-    const double p3 = hdist*hdist*hdist/3.0;
+    const double hdist = sqrt(square(bell_radius) - square(radius));
+    const double p1 = 2.0/3.0*cube(bell_radius);
+    const double p2 = square(bell_radius)*hdist;
+    const double p3 = cube(hdist)/3.0;
 
-    return M_PI*radius*radius*length + 2.0*M_PI*(p1+p2-p3);
+    return M_PI*square(radius)*length + 2.0*M_PI*(p1+p2-p3);
 }
 
-double Iq(double q, double sld,
-        double solvent_sld,
-        double bell_radius,
-        double radius,
-        double length)
+double Iq(double q, double sld, double solvent_sld,
+          double bell_radius, double radius, double length)
 {
-    double sn, cn; // slots to hold sincos function output
+    // Exclude invalid inputs.
+    if (bell_radius < radius) return NAN;
+    const double h = -sqrt(bell_radius*bell_radius - radius*radius);
+    const double half_length = 0.5*length;
 
-    if (bell_radius < radius) return NAN;
-
-    const double lower = 0.0;
-    const double upper = M_PI_2;
-    const double h = sqrt(bell_radius*bell_radius-radius*radius);
+    // translate a point in [-1,1] to a point in [0, pi/2]
+    const double zm = M_PI_4;
+    const double zb = M_PI_4;
     double total = 0.0;
     for (int i = 0; i < 76; i++){
-        const double alpha= 0.5*(Gauss76Z[i]*(upper-lower) + upper + lower);
-        SINCOS(alpha, sn, cn);
+        const double alpha= Gauss76Z[i]*zm + zb;
+        double sin_alpha, cos_alpha; // slots to hold sincos function output
+        SINCOS(alpha, sin_alpha, cos_alpha);
 
-        const double bell_Fq = _bell_kernel(q, h, bell_radius, length, sn, cn);
+        const double bell_Fq = _bell_kernel(q, h, bell_radius, half_length, sin_alpha, cos_alpha);
+        const double bj = J1c(q*radius*sin_alpha);
+        const double si = sinc(q*half_length*cos_alpha);
+        const double cyl_Fq = M_PI*radius*radius*length*bj*si;
+        const double Aq = bell_Fq + cyl_Fq;
+        total += Gauss76Wt[i] * Aq * Aq * sin_alpha;
+    }
+    // translate dx in [-1,1] to dx in [lower,upper]
+    const double form = total*zm;
 
-        const double arg1 = q*length*0.5*cn;
-        const double arg2 = q*radius*sn;
-        // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1
-        const double be = (arg2 == 0.0 ? 0.5 :J1(arg2)/arg2);
-        const double si = (arg1 == 0.0 ? 1.0 :sin(arg1)/arg1);
-        const double cyl_Fq = M_PI*radius*radius*length*si*2.0*be;
-
-        const double Aq = cyl_Fq + bell_Fq;
-        total += Gauss76Wt[i] * Aq * Aq * sn;
-    }
-
-    const double form = total*(upper-lower)*0.5;
-
-    //Contrast and volume normalization
+    //Contrast
     const double s = (sld - solvent_sld);
-    return form*1.0e-4*s*s; //form_volume(bell_radius,radius,length);
+    return 1.0e-4 * s * s * form;
 }
 
 
-
 double Iqxy(double qx, double qy,
-        double sld,
-        double solvent_sld,
-        double bell_radius,
-        double radius,
-        double length,
-        double theta,
-        double phi)
+        double sld, double solvent_sld,
+        double bell_radius, double radius, double length,
+        double theta, double phi)
 {
-     double sn, cn; // slots to hold sincos function output
+    // Compute angle alpha between q and the cylinder axis
+    double sn, cn; // slots to hold sincos function output
+    SINCOS(theta*M_PI_180, sn, cn);
+    const double q = sqrt(qx*qx+qy*qy);
+    const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);
+    const double alpha = acos(cos_val); // rod angle relative to q
 
     // Exclude invalid inputs.
     if (bell_radius < radius) return NAN;
+    const double h = -sqrt(square(bell_radius) - square(radius));
+    const double half_length = 0.5*length;
 
-    // Compute angle alpha between q and the cylinder axis
-    SINCOS(theta*M_PI_180, sn, cn);
-    // # The following correction factor exists in sasview, but it can't be
-    // # right, so we are leaving it out for now.
-    const double q = sqrt(qx*qx+qy*qy);
-    const double cos_val = cn*cos(phi*M_PI_180)*qx + sn*qy;
-    const double alpha = acos(cos_val); // rod angle relative to q
-    SINCOS(alpha, sn, cn);
-
-    const double h = sqrt(bell_radius*bell_radius - radius*radius); // negative h
-    const double bell_Fq = _bell_kernel(q, h, bell_radius, length, sn, cn)/sn;
-
-    const double besarg = q*radius*sn;
-    const double siarg = q*0.5*length*cn;
-    // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1
-    const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg);
-    const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg);
-    const double cyl_Fq = M_PI*radius*radius*length*2.0*bj*si;
-
-    // Volume weighted average F(q)
+    double sin_alpha, cos_alpha; // slots to hold sincos function output
+    SINCOS(alpha, sin_alpha, cos_alpha);
+    const double bell_Fq = _bell_kernel(q, h, bell_radius, half_length, sin_alpha, cos_alpha);
+    const double bj = J1c(q*radius*sin_alpha);
+    const double si = sinc(q*half_length*cos_alpha);
+    const double cyl_Fq = M_PI*radius*radius*length*bj*si;
     const double Aq = cyl_Fq + bell_Fq;
 
-    // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1
+    // Multiply by contrast^2 and convert to cm-1
     const double s = (sld - solvent_sld);
-    return 1.0e-4 * Aq * Aq * s * s; // form_volume(radius, cap_radius, length);
+    return 1.0e-4 * square(s * Aq);
 }
Index: sasmodels/models/barbell.py
===================================================================
--- sasmodels/models/barbell.py	(revision dcdf29d606cfaa2f3bb9d578ac394c6f167fa360)
+++ sasmodels/models/barbell.py	(revision 50e1e40cf14d4ded28bf96e61a3cfaf35003b899)
@@ -121,5 +121,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/J1.c", "lib/gauss76.c", "barbell.c"]
+source = ["lib/J1c.c", "lib/gauss76.c", "barbell.c"]
 
 # parameters for demo
Index: sasmodels/models/capped_cylinder.c
===================================================================
--- sasmodels/models/capped_cylinder.c	(revision 139c528292ae3514b98f9222f04775881a46a7e1)
+++ sasmodels/models/capped_cylinder.c	(revision 50e1e40cf14d4ded28bf96e61a3cfaf35003b899)
@@ -13,15 +13,14 @@
 //   length is the cylinder length, or the separation between the lens halves
 //   alpha is the angle of the cylinder wrt q.
-double _cap_kernel(double q, double h, double cap_radius, double length,
-                 double sin_alpha, double cos_alpha);
-double _cap_kernel(double q, double h, double cap_radius, double length,
-                 double sin_alpha, double cos_alpha)
+static double
+_cap_kernel(double q, double h, double cap_radius,
+                      double half_length, double sin_alpha, double cos_alpha)
 {
-    // For speed, we are pre-calculating terms which are constant over the
-    // kernel.
+    // translate a point in [-1,1] to a point in [lower,upper]
     const double upper = 1.0;
     const double lower = h/cap_radius; // integral lower bound
     const double zm = 0.5*(upper-lower);
     const double zb = 0.5*(upper+lower);
+
     // cos term in integral is:
     //    cos (q (R t - h + L/2) cos(alpha))
@@ -32,21 +31,18 @@
     //    b = q(L/2-h) cos(alpha)
     const double m = q*cap_radius*cos_alpha; // cos argument slope
-    const double b = q*(0.5*length-h)*cos_alpha; // cos argument intercept
+    const double b = q*(half_length-h)*cos_alpha; // cos argument intercept
     const double qrst = q*cap_radius*sin_alpha; // Q*R*sin(theta)
     double total = 0.0;
     for (int i=0; i<76 ;i++) {
-        // translate a point in [-1,1] to a point in [lower,upper]
-        //const double t = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
         const double t = Gauss76Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
-        const double arg = qrst*sqrt(radical); // cap bessel function arg
-        const double be = (arg == 0.0 ? 0.5 : J1(arg)/arg);
-        const double Fq = cos(m*t + b) * radical * be;
+        const double bj = J1c(qrst*sqrt(radical));
+        const double Fq = cos(m*t + b) * radical * bj;
         total += Gauss76Wt[i] * Fq;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
-    //const double form = (upper-lower)/2.0*total;
-    const double integral = 0.5*(upper-lower)*total;
-    return 4.0*M_PI*cap_radius*cap_radius*cap_radius*integral;
+    const double integral = total*zm;
+    const double cap_Fq = 2*M_PI*cube(cap_radius)*integral;
+    return cap_Fq;
 }
 
@@ -66,5 +62,5 @@
     // The first part is clearly V_cyl.  The second part requires some work:
     //    (R-h)^2 => h_c^2
-    //    (2R+h) => 2R+ h_c-h_c + h => 2R + (R-h)-hc + h => 3R-h_c
+    //    (2R+h) => 2R+ h_c-h_c + h => 2R + (R-h)-h_c + h => 3R-h_c
     // And so:
     //    2/3 pi (R-h)^2 (2R + h) => 2/3 pi h_c^2 (3 r_cap - h_c)
@@ -77,93 +73,67 @@
     //        = pi (r^2 (L+hc) + hc^3/3)
     const double hc = cap_radius - sqrt(cap_radius*cap_radius - radius*radius);
-    return M_PI*(radius*radius*(length+hc) + 0.333333333333333*hc*hc*hc);
+    return M_PI*(radius*radius*(length+hc) + hc*hc*hc/3.0);
 }
 
-double Iq(double q,
-    double sld,
-    double solvent_sld,
-    double radius,
-    double cap_radius,
-    double length)
+double Iq(double q, double sld, double solvent_sld,
+          double radius, double cap_radius, double length)
 {
-    double sn, cn; // slots to hold sincos function output
-
     // Exclude invalid inputs.
     if (cap_radius < radius) return NAN;
+    const double h = sqrt(cap_radius*cap_radius - radius*radius);
+    const double half_length = 0.5*length;
 
-    const double lower = 0.0;
-    const double upper = M_PI_2;
-    const double h = sqrt(cap_radius*cap_radius - radius*radius); // negative h
+    // translate a point in [-1,1] to a point in [0, pi/2]
+    const double zm = M_PI_4;
+    const double zb = M_PI_4;
     double total = 0.0;
     for (int i=0; i<76 ;i++) {
-        // translate a point in [-1,1] to a point in [lower,upper]
-        const double alpha= 0.5*(Gauss76Z[i]*(upper-lower) + upper + lower);
-        SINCOS(alpha, sn, cn);
+        const double alpha= Gauss76Z[i]*zm + zb;
+        double sin_alpha, cos_alpha; // slots to hold sincos function output
+        SINCOS(alpha, sin_alpha, cos_alpha);
 
-        const double cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn);
-
-        // The following is CylKernel() / sin(alpha), but we are doing it in place
-        // to avoid sin(alpha)/sin(alpha) for alpha = 0.  It is also a teensy bit
-        // faster since we don't multiply and divide sin(alpha).
-        const double besarg = q*radius*sn;
-        const double siarg = q*0.5*length*cn;
-        // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1
-        const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg);
-        const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg);
-        const double cyl_Fq = M_PI*radius*radius*length*2.0*bj*si;
-
-        // Volume weighted average F(q)
-        const double Aq = cyl_Fq + cap_Fq;
-        total += Gauss76Wt[i] * Aq * Aq * sn; // sn for spherical coord integration
+        const double cap_Fq = _cap_kernel(q, h, cap_radius, half_length, sin_alpha, cos_alpha);
+        const double bj = J1c(q*radius*sin_alpha);
+        const double si = sinc(q*half_length*cos_alpha);
+        const double cyl_Fq = M_PI*radius*radius*length*bj*si;
+        const double Aq = cap_Fq + cyl_Fq;
+        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; // sin_alpha for spherical coord integration
     }
     // translate dx in [-1,1] to dx in [lower,upper]
-    const double form = total * 0.5*(upper-lower);
+    const double form = total * zm;
 
-    // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1
-    // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
-    // The additional volume factor is for polydisperse volume normalization.
+    // Contrast
     const double s = (sld - solvent_sld);
-    return 1.0e-4 * form * s * s; // form_volume(radius, cap_radius, length);
+    return 1.0e-4 * s * s * form;
 }
 
 
 double Iqxy(double qx, double qy,
-    double sld,
-    double solvent_sld,
-    double radius,
-    double cap_radius,
-    double length,
-    double theta,
-    double phi)
+    double sld, double solvent_sld, double radius,
+    double cap_radius, double length,
+    double theta, double phi)
 {
-    double sn, cn; // slots to hold sincos function output
+    // Compute angle alpha between q and the cylinder axis
+    double sn, cn;
+    SINCOS(theta*M_PI_180, sn, cn);
+    const double q = sqrt(qx*qx+qy*qy);
+    const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);
+    const double alpha = acos(cos_val); // rod angle relative to q
 
     // Exclude invalid inputs.
     if (cap_radius < radius) return NAN;
+    const double h = sqrt(cap_radius*cap_radius - radius*radius);
+    const double half_length = 0.5*length;
 
-    // Compute angle alpha between q and the cylinder axis
-    SINCOS(theta*M_PI_180, sn, cn);
-    // # The following correction factor exists in sasview, but it can't be
-    // # right, so we are leaving it out for now.
-    const double q = sqrt(qx*qx+qy*qy);
-    const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);
-    const double alpha = acos(cos_val); // rod angle relative to q
-    SINCOS(alpha, sn, cn);
+    double sin_alpha, cos_alpha; // slots to hold sincos function output
+    SINCOS(alpha, sin_alpha, cos_alpha);
+    const double cap_Fq = _cap_kernel(q, h, cap_radius, half_length, sin_alpha, cos_alpha);
+    const double bj = J1c(q*radius*sin_alpha);
+    const double si = sinc(q*half_length*cos_alpha);
+    const double cyl_Fq = M_PI*radius*radius*length*bj*si;
+    const double Aq = cap_Fq + cyl_Fq;
 
-    const double h = sqrt(cap_radius*cap_radius - radius*radius); // negative h
-    const double cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn);
-
-    const double besarg = q*radius*sn;
-    const double siarg = q*0.5*length*cn;
-    // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1
-    const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg);
-    const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg);
-    const double cyl_Fq = M_PI*radius*radius*length*2.0*bj*si;
-
-    // Volume weighted average F(q)
-    const double Aq = cyl_Fq + cap_Fq;
-
-    // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1
+    // Multiply by contrast^2 and convert to cm-1
     const double s = (sld - solvent_sld);
-    return 1.0e-4 * Aq * Aq * s * s; // form_volume(radius, cap_radius, length);
+    return 1.0e-4 * square(s * Aq);
 }
Index: sasmodels/models/capped_cylinder.py
===================================================================
--- sasmodels/models/capped_cylinder.py	(revision dcdf29d606cfaa2f3bb9d578ac394c6f167fa360)
+++ sasmodels/models/capped_cylinder.py	(revision 50e1e40cf14d4ded28bf96e61a3cfaf35003b899)
@@ -147,5 +147,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/J1.c", "lib/gauss76.c", "capped_cylinder.c"]
+source = ["lib/J1c.c", "lib/gauss76.c", "capped_cylinder.c"]
 
 demo = dict(scale=1, background=0,
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision 3c97ff07f375d382f26608fa2873fc2ce6733fe0)
+++ sasmodels/models/cylinder.c	(revision 50e1e40cf14d4ded28bf96e61a3cfaf35003b899)
@@ -3,15 +3,4 @@
 double Iqxy(double qx, double qy, double sld, double solvent_sld,
     double radius, double length, double theta, double phi);
-
-// twovd = 2 * volume * delta_rho
-// besarg = q * R * sin(alpha)
-// siarg = q * L/2 * cos(alpha)
-double _cyl(double besarg, double siarg);
-double _cyl(double besarg, double siarg)
-{
-    const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg);
-    const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg);
-    return si*bj;
-}
 
 double form_volume(double radius, double length)
@@ -26,23 +15,25 @@
     double length)
 {
+    // TODO: return NaN if radius<0 or length<0?
+    // precompute qr and qh to save time in the loop
     const double qr = q*radius;
     const double qh = q*0.5*length;
+
+    // translate a point in [-1,1] to a point in [0, pi/2]
+    const double zm = M_PI_4;
+    const double zb = M_PI_4;
+
     double total = 0.0;
-    // double lower=0, upper=M_PI_2;
     for (int i=0; i<76 ;i++) {
-        // translate a point in [-1,1] to a point in [lower,upper]
-        //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
-        const double alpha = M_PI_4*(Gauss76Z[i] + 1.0);
+        const double alpha = Gauss76Z[i]*zm + zb;
         double sn, cn;
         SINCOS(alpha, sn, cn);
-        // For a bit of efficiency, we are moving the 2 V delta rho constant
-        // factor, 2Vdrho, out of the loop, so this is fq/2Vdrho rather than fq.
-        const double fq = _cyl(qr*sn, qh*cn);
-        total += Gauss76Wt[i] * fq * fq * sn;
+        const double fq = sinc(qh*cn) * J1c(qr*sn);
+        total += Gauss76Wt[i] * fq*fq * sn;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
-    //const double form = (upper-lower)/2.0*total;
-    const double twoVdrho = 2.0*(sld-solvent_sld)*form_volume(radius, length);
-    return 1.0e-4 * twoVdrho * twoVdrho * total * M_PI_4;
+    const double form = total*zm;
+    const double s = (sld - solvent_sld) * form_volume(radius, length);
+    return 1.0e-4 * s * s * form;
 }
 
@@ -56,21 +47,16 @@
     double phi)
 {
-    // TODO: check that radius<0 and length<0 give zero scattering.
-    // This should be the case since the polydispersity weight vector should
-    // be zero length, and this function never called.
+    // TODO: return NaN if radius<0 or length<0?
     double sn, cn; // slots to hold sincos function output
 
     // Compute angle alpha between q and the cylinder axis
     SINCOS(theta*M_PI_180, sn, cn);
-    const double q = sqrt(qx*qx+qy*qy);
+    const double q = sqrt(qx*qx + qy*qy);
     const double cos_val = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q);
     const double alpha = acos(cos_val);
+
     SINCOS(alpha, sn, cn);
-    //sn = sqrt(1.0 - cos_val*cos_val);
-    //sn = 1.0 - 0.5*cos_val*cos_val;  // if cos_val is very small
-    //cn = cos_val;
-
-    const double twovd = 2.0*(sld-solvent_sld)*form_volume(radius, length);
-    const double fq = twovd * _cyl(q*radius*sn, q*0.5*length*cn);
-    return 1.0e-4 * fq * fq;
+    const double fq = sinc(q*0.5*length*cn) * J1c(q*radius*sn);
+    const double s = (sld-solvent_sld) * form_volume(radius, length);
+    return 1.0e-4 * square(s * fq);
 }
Index: sasmodels/models/cylinder.py
===================================================================
--- sasmodels/models/cylinder.py	(revision 0c3a22666528577f65591c031dbd899b43980331)
+++ sasmodels/models/cylinder.py	(revision 50e1e40cf14d4ded28bf96e61a3cfaf35003b899)
@@ -141,5 +141,5 @@
              ]
 
-source = ["lib/J1.c", "lib/gauss76.c", "cylinder.c"]
+source = ["lib/J1c.c", "lib/gauss76.c", "cylinder.c"]
 
 def ER(radius, length):
Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
+++ sasmodels/models/ellipsoid.c	(revision 50e1e40cf14d4ded28bf96e61a3cfaf35003b899)
@@ -17,5 +17,5 @@
 double form_volume(double rpolar, double requatorial)
 {
-    return 1.333333333333333*M_PI*rpolar*requatorial*requatorial;
+    return M_4PI_3*rpolar*requatorial*requatorial;
 }
 
@@ -26,16 +26,17 @@
     double requatorial)
 {
-    //const double lower = 0.0;
-    //const double upper = 1.0;
+    // translate a point in [-1,1] to a point in [0, 1]
+    const double zm = 0.5;
+    const double zb = 0.5;
     double total = 0.0;
     for (int i=0;i<76;i++) {
         //const double sin_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2;
-        const double sin_alpha = 0.5*(Gauss76Z[i] + 1.0);
+        const double sin_alpha = Gauss76Z[i]*zm + zb;
         total += Gauss76Wt[i] * _ellipsoid_kernel(q, rpolar, requatorial, sin_alpha);
     }
-    //const double form = (upper-lower)/2*total;
-    const double form = 0.5*total;
+    // translate dx in [-1,1] to dx in [lower,upper]
+    const double form = total*zm;
     const double s = (sld - solvent_sld) * form_volume(rpolar, requatorial);
-    return 1.0e-4 * form * s * s;
+    return 1.0e-4 * s * s * form;
 }
 
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision c6915512ab98df810af8905ba552572737ceb9fa)
+++ sasmodels/models/sphere.py	(revision 50e1e40cf14d4ded28bf96e61a3cfaf35003b899)
@@ -83,5 +83,5 @@
 # This should perhaps be volume normalized?
 form_volume = """
-    return 1.333333333333333*M_PI*radius*radius*radius;
+    return M_4PI_3*cube(radius);
     """
 
Index: sasmodels/models/triaxial_ellipsoid.c
===================================================================
--- sasmodels/models/triaxial_ellipsoid.c	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision 50e1e40cf14d4ded28bf96e61a3cfaf35003b899)
@@ -20,7 +20,7 @@
 
     double sn, cn;
-    //double st, ct;
-    //const double lower = 0.0;
-    //const double upper = 1.0;
+    // translate a point in [-1,1] to a point in [0, 1]
+    const double zm = 0.5;
+    const double zb = 0.5;
     double outer = 0.0;
     for (int i=0;i<76;i++) {
@@ -34,6 +34,6 @@
         double inner = 0.0;
         for (int j=0;j<76;j++) {
-            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 ysq = square(Gauss76Z[j]*zm + zb);
+            const double t = q*sqrt(acosx2 + bsinx2*(1.0-ysq) + c2*ysq);
             const double fq = sph_j1c(t);
             inner += Gauss76Wt[j] * fq * fq ;
@@ -41,8 +41,8 @@
         outer += Gauss76Wt[i] * 0.5 * inner;
     }
-    //const double fq2 = (upper-lower)/2*outer;
-    const double fq2 = 0.5*outer;
+    // translate dx in [-1,1] to dx in [lower,upper]
+    const double fqsq = outer*zm;
     const double s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar);
-    return 1.0e-4 * fq2 * s * s;
+    return 1.0e-4 * s * s * fqsq;
 }
 
@@ -62,5 +62,4 @@
     double sphi, cphi;
     double spsi, cpsi;
-    double st, ct;
 
     const double q = sqrt(qx*qx + qy*qy);
@@ -76,9 +75,8 @@
                           + req_major*req_major*cmu*cmu
                           + rpolar*rpolar*calpha*calpha);
-    SINCOS(t, st, ct);
-    const double fq = ( t==0.0 ? 1.0 : 3.0*(st-t*ct)/(t*t*t) );
+    const double fq = sph_j1c(t);
     const double s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar);
 
-    return 1.0e-4 * fq * fq * s * s;
+    return 1.0e-4 * square(s * fq);
 }
 
