Index: sasmodels/models/barbell.c
===================================================================
--- sasmodels/models/barbell.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
+++ sasmodels/models/barbell.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -10,6 +10,6 @@
 //barbell kernel - same as dumbell
 static double
-_bell_kernel(double q, double h, double radius_bell,
-             double half_length, double sin_alpha, double cos_alpha)
+_bell_kernel(double qab, double qc, double h, double radius_bell,
+             double half_length)
 {
     // translate a point in [-1,1] to a point in [lower,upper]
@@ -26,12 +26,12 @@
     //    m = q R cos(alpha)
     //    b = q(L/2-h) cos(alpha)
-    const double m = q*radius_bell*cos_alpha; // cos argument slope
-    const double b = q*(half_length-h)*cos_alpha; // cos argument intercept
-    const double qrst = q*radius_bell*sin_alpha; // Q*R*sin(theta)
+    const double m = radius_bell*qc; // cos argument slope
+    const double b = (half_length-h)*qc; // cos argument intercept
+    const double qab_r = radius_bell*qab; // Q*R*sin(theta)
     double total = 0.0;
     for (int i = 0; i < 76; i++){
         const double t = Gauss76Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
-        const double bj = sas_2J1x_x(qrst*sqrt(radical));
+        const double bj = sas_2J1x_x(qab_r*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
         total += Gauss76Wt[i] * Fq;
@@ -44,11 +44,10 @@
 
 static double
-_fq(double q, double h,
-    double radius_bell, double radius, double half_length,
-    double sin_alpha, double cos_alpha)
+_fq(double qab, double qc, double h,
+    double radius_bell, double radius, double half_length)
 {
-    const double bell_fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, 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 bell_fq = _bell_kernel(qab, qc, h, radius_bell, half_length);
+    const double bj = sas_2J1x_x(radius*qab);
+    const double si = sas_sinx_x(half_length*qc);
     const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si;
     const double Aq = bell_fq + cyl_fq;
@@ -84,5 +83,5 @@
         double sin_alpha, cos_alpha; // slots to hold sincos function output
         SINCOS(alpha, sin_alpha, cos_alpha);
-        const double Aq = _fq(q, h, radius_bell, radius, half_length, sin_alpha, cos_alpha);
+        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length);
         total += Gauss76Wt[i] * Aq * Aq * sin_alpha;
     }
@@ -103,7 +102,9 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
     const double h = -sqrt(square(radius_bell) - square(radius));
-    const double Aq = _fq(q, h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha);
+    const double Aq = _fq(qab, qc, h, radius_bell, radius, 0.5*length);
 
     // Multiply by contrast^2 and convert to cm-1
Index: sasmodels/models/bcc_paracrystal.c
===================================================================
--- sasmodels/models/bcc_paracrystal.c	(revision 50beefeac18f2f952bf0863526ec9bbd3009e29c)
+++ sasmodels/models/bcc_paracrystal.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,89 +1,93 @@
-double form_volume(double radius);
-double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld);
-double Iqxy(double qx, double qy, double dnn,
-    double d_factor, double radius,double sld, double solvent_sld,
-    double theta, double phi, double psi);
+static double
+_sq_bcc(double qa, double qb, double qc, double dnn, double d_factor)
+{
+    // Rewriting equations for efficiency, accuracy and readability, and so
+    // code is reusable between 1D and 2D models.
+    const double a1 = +qa - qc + qb;
+    const double a2 = +qa + qc - qb;
+    const double a3 = -qa + qc + qb;
 
-double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
-double _BCCeval(double Theta, double Phi, double temp1, double temp3);
-double _sphereform(double q, double radius, double sld, double solvent_sld);
+    const double half_dnn = 0.5*dnn;
+    const double arg = 0.5*square(half_dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+
+#if 1
+    // Numerator: (1 - exp(a)^2)^3
+    //         => (-(exp(2a) - 1))^3
+    //         => -expm1(2a)^3
+    // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2)
+    //         => exp(a)^2 - 2 cos(xk) exp(a) + 1
+    //         => (exp(a) - 2 cos(xk)) * exp(a) + 1
+    const double exp_arg = exp(-arg);
+    const double Sq = -cube(expm1(-2.0*arg))
+        / ( ((exp_arg - 2.0*cos(half_dnn*a1))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(half_dnn*a2))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(half_dnn*a3))*exp_arg + 1.0));
+#else
+    // Alternate form, which perhaps is more approachable
+    const double sinh_qd = sinh(arg);
+    const double cosh_qd = cosh(arg);
+    const double Sq = sinh_qd/(cosh_qd - cos(half_dnn*a1))
+                    * sinh_qd/(cosh_qd - cos(half_dnn*a2))
+                    * sinh_qd/(cosh_qd - cos(half_dnn*a3));
+#endif
+
+    return Sq;
+}
 
 
-double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) {
-
-	const double Da = d_factor*dnn;
-	const double temp1 = q*q*Da*Da;
-	const double temp3 = q*dnn;
-
-	double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI);
-	return(retVal);
+// occupied volume fraction calculated from lattice symmetry and sphere radius
+static double
+_bcc_volume_fraction(double radius, double dnn)
+{
+    return 2.0*sphere_volume(sqrt(0.75)*radius/dnn);
 }
 
-double _BCCeval(double Theta, double Phi, double temp1, double temp3) {
-
-	double result;
-	double sin_theta,cos_theta,sin_phi,cos_phi;
-	SINCOS(Theta, sin_theta, cos_theta);
-	SINCOS(Phi, sin_phi, cos_phi);
-
-	const double temp6 =  sin_theta;
-	const double temp7 =  sin_theta*cos_phi + sin_theta*sin_phi + cos_theta;
-	const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta;
-	const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta;
-
-	const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9));
-	result = cube(1.0 - (temp10*temp10))*temp6
-	    / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10)
-	      * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10)
-	      * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10));
-
-	return (result);
-}
-
-double form_volume(double radius){
+static double
+form_volume(double radius)
+{
     return sphere_volume(radius);
 }
 
 
-double Iq(double q, double dnn,
+static double Iq(double q, double dnn,
   double d_factor, double radius,
-  double sld, double solvent_sld){
+  double sld, double solvent_sld)
+{
+    // translate a point in [-1,1] to a point in [0, 2 pi]
+    const double phi_m = M_PI;
+    const double phi_b = M_PI;
+    // translate a point in [-1,1] to a point in [0, pi]
+    const double theta_m = M_PI_2;
+    const double theta_b = M_PI_2;
 
-	//Volume fraction calculated from lattice symmetry and sphere radius
-	const double s1 = dnn/sqrt(0.75);
-	const double latticescale = 2.0*sphere_volume(radius/s1);
+    double outer_sum = 0.0;
+    for(int i=0; i<150; i++) {
+        double inner_sum = 0.0;
+        const double theta = Gauss150Z[i]*theta_m + theta_b;
+        double sin_theta, cos_theta;
+        SINCOS(theta, sin_theta, cos_theta);
+        const double qc = q*cos_theta;
+        const double qab = q*sin_theta;
+        for(int j=0;j<150;j++) {
+            const double phi = Gauss150Z[j]*phi_m + phi_b;
+            double sin_phi, cos_phi;
+            SINCOS(phi, sin_phi, cos_phi);
+            const double qa = qab*cos_phi;
+            const double qb = qab*sin_phi;
+            const double fq = _sq_bcc(qa, qb, qc, dnn, d_factor);
+            inner_sum += Gauss150Wt[j] * fq;
+        }
+        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
+        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+    }
+    outer_sum *= theta_m;
+    const double Sq = outer_sum/(4.0*M_PI);
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
 
-    const double va = 0.0;
-    const double vb = 2.0*M_PI;
-    const double vaj = 0.0;
-    const double vbj = M_PI;
-
-    double summ = 0.0;
-    double answer = 0.0;
-	for(int i=0; i<150; i++) {
-		//setup inner integral over the ellipsoidal cross-section
-		double summj=0.0;
-		const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;		//the outer dummy is phi
-		for(int j=0;j<150;j++) {
-			//20 gauss points for the inner integral
-			double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;		//the inner dummy is theta
-			double yyy = Gauss150Wt[j] * _BCC_Integrand(q,dnn,d_factor,ztheta,zphi);
-			summj += yyy;
-		}
-		//now calculate the value of the inner integral
-		double answer = (vbj-vaj)/2.0*summj;
-
-		//now calculate outer integral
-		summ = summ+(Gauss150Wt[i] * answer);
-	}		//final scaling is done at the end of the function, after the NT_FP64 case
-
-	answer = (vb-va)/2.0*summ;
-	answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale;
-
-    return answer;
+    return _bcc_volume_fraction(radius, dnn) * Pq * Sq;
 }
 
 
-double Iqxy(double qx, double qy,
+static double Iqxy(double qx, double qy,
     double dnn, double d_factor, double radius,
     double sld, double solvent_sld,
@@ -92,20 +96,11 @@
     double q, zhat, yhat, xhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
 
-    const double a1 = +xhat - zhat + yhat;
-    const double a2 = +xhat + zhat - yhat;
-    const double a3 = -xhat + zhat + yhat;
-
-    const double qd = 0.5*q*dnn;
-    const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3);
-    const double tanh_qd = tanh(arg);
-    const double cosh_qd = cosh(arg);
-    const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*a2)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*a3)/cosh_qd);
-
-    const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq;
-    //the occupied volume of the lattice
-    const double lattice_scale = 2.0*sphere_volume(sqrt(0.75)*radius/dnn);
-    return lattice_scale * Fq;
+    q = sqrt(qa*qa + qb*qb + qc*qc);
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+    const double Sq = _sq_bcc(qa, qb, qc, dnn, d_factor);
+    return _bcc_volume_fraction(radius, dnn) * Pq * Sq;
 }
Index: sasmodels/models/bcc_paracrystal.py
===================================================================
--- sasmodels/models/bcc_paracrystal.py	(revision 69e1afc57270f67b2a0d7157a6d263fc2db98cc0)
+++ sasmodels/models/bcc_paracrystal.py	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -81,5 +81,5 @@
 .. figure:: img/parallelepiped_angle_definition.png
 
-    Orientation of the crystal with respect to the scattering plane, when 
+    Orientation of the crystal with respect to the scattering plane, when
     $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis).
 
Index: sasmodels/models/capped_cylinder.c
===================================================================
--- sasmodels/models/capped_cylinder.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
+++ sasmodels/models/capped_cylinder.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -14,8 +14,8 @@
 //   radius_cap is the radius of the lens
 //   length is the cylinder length, or the separation between the lens halves
-//   alpha is the angle of the cylinder wrt q.
+//   theta is the angle of the cylinder wrt q.
 static double
-_cap_kernel(double q, double h, double radius_cap,
-                      double half_length, double sin_alpha, double cos_alpha)
+_cap_kernel(double qab, double qc, double h, double radius_cap,
+            double half_length)
 {
     // translate a point in [-1,1] to a point in [lower,upper]
@@ -26,18 +26,18 @@
 
     // cos term in integral is:
-    //    cos (q (R t - h + L/2) cos(alpha))
+    //    cos (q (R t - h + L/2) cos(theta))
     // 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*radius_cap*cos_alpha; // cos argument slope
-    const double b = q*(half_length-h)*cos_alpha; // cos argument intercept
-    const double qrst = q*radius_cap*sin_alpha; // Q*R*sin(theta)
+    //    m = q R cos(theta)
+    //    b = q(L/2-h) cos(theta)
+    const double m = radius_cap*qc; // cos argument slope
+    const double b = (half_length-h)*qc; // cos argument intercept
+    const double qab_r = radius_cap*qab; // Q*R*sin(theta)
     double total = 0.0;
     for (int i=0; i<76 ;i++) {
         const double t = Gauss76Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
-        const double bj = sas_2J1x_x(qrst*sqrt(radical));
+        const double bj = sas_2J1x_x(qab_r*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
         total += Gauss76Wt[i] * Fq;
@@ -50,10 +50,9 @@
 
 static double
-_fq(double q, double h, double radius_cap, double radius, double half_length,
-    double sin_alpha, double cos_alpha)
+_fq(double qab, double qc, double h, double radius_cap, double radius, double half_length)
 {
-    const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, 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 cap_Fq = _cap_kernel(qab, qc, h, radius_cap, half_length);
+    const double bj = sas_2J1x_x(radius*qab);
+    const double si = sas_sinx_x(half_length*qc);
     const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si;
     const double Aq = cap_Fq + cyl_Fq;
@@ -101,11 +100,12 @@
     double total = 0.0;
     for (int i=0; i<76 ;i++) {
-        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 Aq = _fq(q, h, radius_cap, radius, half_length, sin_alpha, cos_alpha);
-        // sin_alpha for spherical coord integration
-        total += Gauss76Wt[i] * Aq * Aq * sin_alpha;
+        const double theta = Gauss76Z[i]*zm + zb;
+        double sin_theta, cos_theta; // slots to hold sincos function output
+        SINCOS(theta, sin_theta, cos_theta);
+        const double qab = q*sin_theta;
+        const double qc = q*cos_theta;
+        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length);
+        // scale by sin_theta for spherical coord integration
+        total += Gauss76Wt[i] * Aq * Aq * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -125,7 +125,9 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
     const double h = sqrt(radius_cap*radius_cap - radius*radius);
-    const double Aq = _fq(q, h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha);
+    const double Aq = _fq(qab, qc, h, radius_cap, radius, 0.5*length);
 
     // Multiply by contrast^2 and convert to cm-1
Index: sasmodels/models/core_shell_bicelle.c
===================================================================
--- sasmodels/models/core_shell_bicelle.c	(revision b260926c42c2ae859a32c20b98935378d51fab4d)
+++ sasmodels/models/core_shell_bicelle.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,56 +1,31 @@
-double form_volume(double radius, double thick_rim, double thick_face, double length);
-double Iq(double q,
-          double radius,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld);
-
-
-double Iqxy(double qx, double qy,
-          double radius,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld,
-          double theta,
-          double phi);
-
-
-double form_volume(double radius, double thick_rim, double thick_face, double length)
+static double
+form_volume(double radius, double thick_rim, double thick_face, double length)
 {
-    return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face);
+    return M_PI*square(radius+thick_rim)*(length+2.0*thick_face);
 }
 
 static double
-bicelle_kernel(double q,
-              double rad,
-              double radthick,
-              double facthick,
-              double halflength,
-              double rhoc,
-              double rhoh,
-              double rhor,
-              double rhosolv,
-              double sin_alpha,
-              double cos_alpha)
+bicelle_kernel(double qab,
+    double qc,
+    double radius,
+    double thick_radius,
+    double thick_face,
+    double halflength,
+    double sld_core,
+    double sld_face,
+    double sld_rim,
+    double sld_solvent)
 {
-    const double dr1 = rhoc-rhoh;
-    const double dr2 = rhor-rhosolv;
-    const double dr3 = rhoh-rhor;
-    const double vol1 = M_PI*square(rad)*2.0*(halflength);
-    const double vol2 = M_PI*square(rad+radthick)*2.0*(halflength+facthick);
-    const double vol3 = M_PI*square(rad)*2.0*(halflength+facthick);
+    const double dr1 = sld_core-sld_face;
+    const double dr2 = sld_rim-sld_solvent;
+    const double dr3 = sld_face-sld_rim;
+    const double vol1 = M_PI*square(radius)*2.0*(halflength);
+    const double vol2 = M_PI*square(radius+thick_radius)*2.0*(halflength+thick_face);
+    const double vol3 = M_PI*square(radius)*2.0*(halflength+thick_face);
 
-    const double be1 = sas_2J1x_x(q*(rad)*sin_alpha);
-    const double be2 = sas_2J1x_x(q*(rad+radthick)*sin_alpha);
-    const double si1 = sas_sinx_x(q*(halflength)*cos_alpha);
-    const double si2 = sas_sinx_x(q*(halflength+facthick)*cos_alpha);
+    const double be1 = sas_2J1x_x((radius)*qab);
+    const double be2 = sas_2J1x_x((radius+thick_radius)*qab);
+    const double si1 = sas_sinx_x((halflength)*qc);
+    const double si2 = sas_sinx_x((halflength+thick_face)*qc);
 
     const double t = vol1*dr1*si1*be1 +
@@ -58,20 +33,17 @@
                      vol3*dr3*si2*be1;
 
-    const double retval = t*t;
-
-    return retval;
-
+    return t;
 }
 
 static double
-bicelle_integration(double q,
-                   double rad,
-                   double radthick,
-                   double facthick,
-                   double length,
-                   double rhoc,
-                   double rhoh,
-                   double rhor,
-                   double rhosolv)
+Iq(double q,
+    double radius,
+    double thick_radius,
+    double thick_face,
+    double length,
+    double sld_core,
+    double sld_face,
+    double sld_rim,
+    double sld_solvent)
 {
     // set up the integration end points
@@ -79,82 +51,40 @@
     const double halflength = 0.5*length;
 
-    double summ = 0.0;
+    double total = 0.0;
     for(int i=0;i<N_POINTS_76;i++) {
-        double alpha = (Gauss76Z[i] + 1.0)*uplim;
-        double sin_alpha, cos_alpha; // slots to hold sincos function output
-        SINCOS(alpha, sin_alpha, cos_alpha);
-        double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick,
-                             halflength, rhoc, rhoh, rhor, rhosolv,
-                             sin_alpha, cos_alpha);
-        summ += yyy*sin_alpha;
+        double theta = (Gauss76Z[i] + 1.0)*uplim;
+        double sin_theta, cos_theta; // slots to hold sincos function output
+        SINCOS(theta, sin_theta, cos_theta);
+        double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face,
+                                   halflength, sld_core, sld_face, sld_rim, sld_solvent);
+        total += Gauss76Wt[i]*fq*fq*sin_theta;
     }
 
     // calculate value of integral to return
-    double answer = uplim*summ;
-    return answer;
+    double answer = total*uplim;
+    return 1.0e-4*answer;
 }
 
 static double
-bicelle_kernel_2d(double qx, double qy,
-          double radius,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld,
-          double theta,
-          double phi)
+Iqxy(double qx, double qy,
+    double radius,
+    double thick_rim,
+    double thick_face,
+    double length,
+    double core_sld,
+    double face_sld,
+    double rim_sld,
+    double solvent_sld,
+    double theta,
+    double phi)
 {
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
-    double answer = bicelle_kernel(q, radius, thick_rim, thick_face,
+    double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face,
                            0.5*length, core_sld, face_sld, rim_sld,
-                           solvent_sld, sin_alpha, cos_alpha);
-    return 1.0e-4*answer;
+                           solvent_sld);
+    return 1.0e-4*fq*fq;
 }
-
-double Iq(double q,
-          double radius,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld)
-{
-    double intensity = bicelle_integration(q, radius, thick_rim, thick_face,
-                       length, core_sld, face_sld, rim_sld, solvent_sld);
-    return intensity*1.0e-4;
-}
-
-
-double Iqxy(double qx, double qy,
-          double radius,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld,
-          double theta,
-          double phi)
-{
-    double intensity = bicelle_kernel_2d(qx, qy,
-                      radius,
-                      thick_rim,
-                      thick_face,
-                      length,
-                      core_sld,
-                      face_sld,
-                      rim_sld,
-                      solvent_sld,
-                      theta,
-                      phi);
-
-    return intensity;
-}
Index: sasmodels/models/core_shell_bicelle_elliptical.c
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical.c	(revision dedcf34f547ba57a5e8ca3c5e9892ede11d1fcaf)
+++ sasmodels/models/core_shell_bicelle_elliptical.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -17,29 +17,21 @@
         double thick_face,
         double length,
-        double rhoc,
-        double rhoh,
-        double rhor,
-        double rhosolv)
+        double sld_core,
+        double sld_face,
+        double sld_rim,
+        double sld_solvent)
 {
-    double si1,si2,be1,be2;
      // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle
      // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle
-     //    const double uplim = M_PI_4;
     const double halfheight = 0.5*length;
-    //const double va = 0.0;
-    //const double vb = 1.0;
-    // inner integral limits
-    //const double vaj=0.0;
-    //const double vbj=M_PI;
-
     const double r_major = r_minor * x_core;
     const double r2A = 0.5*(square(r_major) + square(r_minor));
     const double r2B = 0.5*(square(r_major) - square(r_minor));
-    const double dr1 = (rhoc-rhoh)   *M_PI*r_minor*r_major*(2.0*halfheight);;
-    const double dr2 = (rhor-rhosolv)*M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);
-    const double dr3 = (rhoh-rhor)   *M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
-    //const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight);
-    //const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);
-    //const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
+    const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight);
+    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);
+    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
+    const double dr1 = vol1*(sld_core-sld_face);
+    const double dr2 = vol2*(sld_rim-sld_solvent);
+    const double dr3 = vol3*(sld_face-sld_rim);
 
     //initialize integral
@@ -47,25 +39,27 @@
     for(int i=0;i<76;i++) {
         //setup inner integral over the ellipsoidal cross-section
-        // since we generate these lots of times, why not store them somewhere?
-        //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
-        const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0;
-        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
-        double inner_sum=0;
-        double sinarg1 = q*halfheight*cos_alpha;
-        double sinarg2 = q*(halfheight+thick_face)*cos_alpha;
-        si1 = sas_sinx_x(sinarg1);
-        si2 = sas_sinx_x(sinarg2);
+        //const double va = 0.0;
+        //const double vb = 1.0;
+        //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
+        const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0;
+        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
+        const double qab = q*sin_theta;
+        const double qc = q*cos_theta;
+        const double si1 = sas_sinx_x(halfheight*qc);
+        const double si2 = sas_sinx_x((halfheight+thick_face)*qc);
+        double inner_sum=0.0;
         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)
-            //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
-            const double beta = ( Gauss76Z[j] +1.0)*M_PI_2;
-            const double rr = sqrt(r2A - r2B*cos(beta));
-            double besarg1 = q*rr*sin_alpha;
-            double besarg2 = q*(rr+thick_rim)*sin_alpha;
-            be1 = sas_2J1x_x(besarg1);
-            be2 = sas_2J1x_x(besarg2);
-            inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 +
-                                              dr2*si2*be2 +
-                                              dr3*si2*be1);
+            // inner integral limits
+            //const double vaj=0.0;
+            //const double vbj=M_PI;
+            //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
+            const double phi = ( Gauss76Z[j] +1.0)*M_PI_2;
+            const double rr = sqrt(r2A - r2B*cos(phi));
+            const double be1 = sas_2J1x_x(rr*qab);
+            const double be2 = sas_2J1x_x((rr+thick_rim)*qab);
+            const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1;
+
+            inner_sum += Gauss76Wt[j] * fq * fq;
         }
         //now calculate outer integral
@@ -83,18 +77,21 @@
           double thick_face,
           double length,
-          double rhoc,
-          double rhoh,
-          double rhor,
-          double rhosolv,
+          double sld_core,
+          double sld_face,
+          double sld_rim,
+          double sld_solvent,
           double theta,
           double phi,
           double psi)
 {
-       // THIS NEEDS TESTING
     double q, xhat, yhat, zhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
-    const double dr1 = rhoc-rhoh;
-    const double dr2 = rhor-rhosolv;
-    const double dr3 = rhoh-rhor;
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
+
+    const double dr1 = sld_core-sld_face;
+    const double dr2 = sld_rim-sld_solvent;
+    const double dr3 = sld_face-sld_rim;
     const double r_major = r_minor*x_core;
     const double halfheight = 0.5*length;
@@ -104,13 +101,12 @@
 
     // Compute effective radius in rotated coordinates
-    const double r_hat = sqrt(square(r_major*xhat) + square(r_minor*yhat));
-    const double rshell_hat = sqrt(square((r_major+thick_rim)*xhat)
-                                   + square((r_minor+thick_rim)*yhat));
-    const double be1 = sas_2J1x_x( q*r_hat );
-    const double be2 = sas_2J1x_x( q*rshell_hat );
-    const double si1 = sas_sinx_x( q*halfheight*zhat );
-    const double si2 = sas_sinx_x( q*(halfheight + thick_face)*zhat );
-    const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1);
-    return 1.0e-4 * Aq;
+    const double qr_hat = sqrt(square(r_major*qa) + square(r_minor*qb));
+    const double qrshell_hat = sqrt(square((r_major+thick_rim)*qa)
+                                   + square((r_minor+thick_rim)*qb));
+    const double be1 = sas_2J1x_x( qr_hat );
+    const double be2 = sas_2J1x_x( qrshell_hat );
+    const double si1 = sas_sinx_x( halfheight*qc );
+    const double si2 = sas_sinx_x( (halfheight + thick_face)*qc );
+    const double fq = vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1;
+    return 1.0e-4 * fq*fq;
 }
-
Index: sasmodels/models/core_shell_cylinder.c
===================================================================
--- sasmodels/models/core_shell_cylinder.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
+++ sasmodels/models/core_shell_cylinder.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,23 +1,18 @@
-double form_volume(double radius, double thickness, double length);
-double Iq(double q, double core_sld, double shell_sld, double solvent_sld,
-    double radius, double thickness, double length);
-double Iqxy(double qx, double qy, double core_sld, double shell_sld, double solvent_sld,
-    double radius, double thickness, double length, double theta, double phi);
-
 // vd = volume * delta_rho
-// besarg = q * R * sin(alpha)
-// siarg = q * L/2 * cos(alpha)
-double _cyl(double vd, double besarg, double siarg);
-double _cyl(double vd, double besarg, double siarg)
+// besarg = q * R * sin(theta)
+// siarg = q * L/2 * cos(theta)
+static double _cyl(double vd, double besarg, double siarg)
 {
     return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg);
 }
 
-double form_volume(double radius, double thickness, double length)
+static double
+form_volume(double radius, double thickness, double length)
 {
-    return M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness);
+    return M_PI*square(radius+thickness)*(length+2.0*thickness);
 }
 
-double Iq(double q,
+static double
+Iq(double q,
     double core_sld,
     double shell_sld,
@@ -28,21 +23,22 @@
 {
     // precalculate constants
-    const double core_qr = q*radius;
-    const double core_qh = q*0.5*length;
+    const double core_r = radius;
+    const double core_h = 0.5*length;
     const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld);
-    const double shell_qr = q*(radius + thickness);
-    const double shell_qh = q*(0.5*length + thickness);
+    const double shell_r = (radius + thickness);
+    const double shell_h = (0.5*length + thickness);
     const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld);
     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;
-        double sn, cn;
-        const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2);
-        SINCOS(alpha, sn, cn);
-        const double fq = _cyl(core_vd, core_qr*sn, core_qh*cn)
-            + _cyl(shell_vd, shell_qr*sn, shell_qh*cn);
-        total += Gauss76Wt[i] * fq * fq * sn;
+        // translate a point in [-1,1] to a point in [0, pi/2]
+        //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
+        double sin_theta, cos_theta;
+        const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4;
+        SINCOS(theta, sin_theta,  cos_theta);
+        const double qab = q*sin_theta;
+        const double qc = q*cos_theta;
+        const double fq = _cyl(core_vd, core_r*qab, core_h*qc)
+            + _cyl(shell_vd, shell_r*qab, shell_h*qc);
+        total += Gauss76Wt[i] * fq * fq * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -64,14 +60,16 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
-    const double core_qr = q*radius;
-    const double core_qh = q*0.5*length;
+    const double core_r = radius;
+    const double core_h = 0.5*length;
     const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld);
-    const double shell_qr = q*(radius + thickness);
-    const double shell_qh = q*(0.5*length + thickness);
+    const double shell_r = (radius + thickness);
+    const double shell_h = (0.5*length + thickness);
     const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld);
 
-    const double fq = _cyl(core_vd, core_qr*sin_alpha, core_qh*cos_alpha)
-        + _cyl(shell_vd, shell_qr*sin_alpha, shell_qh*cos_alpha);
+    const double fq = _cyl(core_vd, core_r*qab, core_h*qc)
+        + _cyl(shell_vd, shell_r*qab, shell_h*qc);
     return 1.0e-4 * fq * fq;
 }
Index: sasmodels/models/core_shell_ellipsoid.c
===================================================================
--- sasmodels/models/core_shell_ellipsoid.c	(revision 0a3d9b2c366d93b98a02e5c9077f26098e630bd0)
+++ sasmodels/models/core_shell_ellipsoid.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,32 +1,34 @@
-double form_volume(double radius_equat_core,
-                   double polar_core,
-                   double equat_shell,
-                   double polar_shell);
-double Iq(double q,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld);
 
+// Converted from Igor function gfn4, using the same pattern as ellipsoid
+// for evaluating the parts of the integral.
+//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN
+//                       BY (53) & (58-59) IN CHEN AND
+//                       KOTLARCHYK REFERENCE
+//
+//       <OBLATE ELLIPSOID>
+static double
+_cs_ellipsoid_kernel(double qab, double qc,
+    double equat_core, double polar_core,
+    double equat_shell, double polar_shell,
+    double sld_core_shell, double sld_shell_solvent)
+{
+    const double qr_core = sqrt(square(equat_core*qab) + square(polar_core*qc));
+    const double si_core = sas_3j1x_x(qr_core);
+    const double volume_core = M_4PI_3*equat_core*equat_core*polar_core;
+    const double fq_core = si_core*volume_core*sld_core_shell;
 
-double Iqxy(double qx, double qy,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double theta,
-          double phi);
+    const double qr_shell = sqrt(square(equat_shell*qab) + square(polar_shell*qc));
+    const double si_shell = sas_3j1x_x(qr_shell);
+    const double volume_shell = M_4PI_3*equat_shell*equat_shell*polar_shell;
+    const double fq_shell = si_shell*volume_shell*sld_shell_solvent;
 
+    return fq_core + fq_shell;
+}
 
-double form_volume(double radius_equat_core,
-                   double x_core,
-                   double thick_shell,
-                   double x_polar_shell)
+static double
+form_volume(double radius_equat_core,
+    double x_core,
+    double thick_shell,
+    double x_polar_shell)
 {
     const double equat_shell = radius_equat_core + thick_shell;
@@ -37,20 +39,15 @@
 
 static double
-core_shell_ellipsoid_xt_kernel(double q,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld)
+Iq(double q,
+    double radius_equat_core,
+    double x_core,
+    double thick_shell,
+    double x_polar_shell,
+    double core_sld,
+    double shell_sld,
+    double solvent_sld)
 {
-    const double lolim = 0.0;
-    const double uplim = 1.0;
-
-
-    const double delpc = core_sld - shell_sld; //core - shell
-    const double delps = shell_sld - solvent_sld; //shell - solvent
-
+    const double sld_core_shell = core_sld - shell_sld;
+    const double sld_shell_solvent = shell_sld - solvent_sld;
 
     const double polar_core = radius_equat_core*x_core;
@@ -58,34 +55,42 @@
     const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell;
 
-    double summ = 0.0;	 //initialize intergral
+    // translate from [-1, 1] => [0, 1]
+    const double m = 0.5;
+    const double b = 0.5;
+    double total = 0.0;     //initialize intergral
     for(int i=0;i<76;i++) {
-        double zi = 0.5*( Gauss76Z[i]*(uplim-lolim) + uplim + lolim );
-        double yyy = gfn4(zi, radius_equat_core, polar_core, equat_shell,
-                          polar_shell, delpc, delps, q);
-        summ += Gauss76Wt[i] * yyy;
+        const double cos_theta = Gauss76Z[i]*m + b;
+        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
+        double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta,
+            radius_equat_core, polar_core,
+            equat_shell, polar_shell,
+            sld_core_shell, sld_shell_solvent);
+        total += Gauss76Wt[i] * fq * fq;
     }
-    summ *= 0.5*(uplim-lolim);
+    total *= m;
 
     // convert to [cm-1]
-    return 1.0e-4 * summ;
+    return 1.0e-4 * total;
 }
 
 static double
-core_shell_ellipsoid_xt_kernel_2d(double qx, double qy,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double theta,
-          double phi)
+Iqxy(double qx, double qy,
+    double radius_equat_core,
+    double x_core,
+    double thick_shell,
+    double x_polar_shell,
+    double core_sld,
+    double shell_sld,
+    double solvent_sld,
+    double theta,
+    double phi)
 {
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
-    const double sldcs = core_sld - shell_sld;
-    const double sldss = shell_sld- solvent_sld;
+    const double sld_core_shell = core_sld - shell_sld;
+    const double sld_shell_solvent = shell_sld - solvent_sld;
 
     const double polar_core = radius_equat_core*x_core;
@@ -93,65 +98,10 @@
     const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell;
 
-    // Call the IGOR library function to get the kernel:
-    // MUST use gfn4 not gf2 because of the def of params.
-    double answer = gfn4(cos_alpha,
-                  radius_equat_core,
-                  polar_core,
-                  equat_shell,
-                  polar_shell,
-                  sldcs,
-                  sldss,
-                  q);
+    double fq = _cs_ellipsoid_kernel(qab, qc,
+                  radius_equat_core, polar_core,
+                  equat_shell, polar_shell,
+                  sld_core_shell, sld_shell_solvent);
 
     //convert to [cm-1]
-    answer *= 1.0e-4;
-
-    return answer;
+    return 1.0e-4 * fq * fq;
 }
-
-double Iq(double q,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld)
-{
-    double intensity = core_shell_ellipsoid_xt_kernel(q,
-           radius_equat_core,
-           x_core,
-           thick_shell,
-           x_polar_shell,
-           core_sld,
-           shell_sld,
-           solvent_sld);
-
-    return intensity;
-}
-
-
-double Iqxy(double qx, double qy,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double theta,
-          double phi)
-{
-    double intensity = core_shell_ellipsoid_xt_kernel_2d(qx, qy,
-                       radius_equat_core,
-                       x_core,
-                       thick_shell,
-                       x_polar_shell,
-                       core_sld,
-                       shell_sld,
-                       solvent_sld,
-                       theta,
-                       phi);
-
-    return intensity;
-}
Index: sasmodels/models/core_shell_ellipsoid.py
===================================================================
--- sasmodels/models/core_shell_ellipsoid.py	(revision 9802ab3945a154b871264715dcb463d4ec644afa)
+++ sasmodels/models/core_shell_ellipsoid.py	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -25,6 +25,6 @@
 ellipsoid. This may have some undesirable effects if the aspect ratio of the
 ellipsoid is large (ie, if $X << 1$ or $X >> 1$ ), when the $S(q)$
-- which assumes spheres - will not in any case be valid.  Generating a 
-custom product model will enable separate effective volume fraction and effective 
+- which assumes spheres - will not in any case be valid.  Generating a
+custom product model will enable separate effective volume fraction and effective
 radius in the $S(q)$.
 
@@ -44,11 +44,11 @@
 
 .. math::
-    \begin{align}    
+    \begin{align}
     F(q,\alpha) = &f(q,radius\_equat\_core,radius\_equat\_core.x\_core,\alpha) \\
     &+ f(q,radius\_equat\_core + thick\_shell,radius\_equat\_core.x\_core + thick\_shell.x\_polar\_shell,\alpha)
-    \end{align} 
+    \end{align}
 
 where
- 
+
 .. math::
 
@@ -77,5 +77,5 @@
    F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}
 
-For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data, 
+For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,
 see the :ref:`elliptical-cylinder` model for further information.
 
@@ -139,6 +139,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sas_3j1x_x.c", "lib/gfn.c", "lib/gauss76.c",
-          "core_shell_ellipsoid.c"]
+source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"]
 
 def ER(radius_equat_core, x_core, thick_shell, x_polar_shell):
Index: sasmodels/models/core_shell_parallelepiped.c
===================================================================
--- sasmodels/models/core_shell_parallelepiped.c	(revision 92dfe0c540274c8b800a9929c3f1e1a7d4cba9dc)
+++ sasmodels/models/core_shell_parallelepiped.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,3 +1,3 @@
-double form_volume(double length_a, double length_b, double length_c, 
+double form_volume(double length_a, double length_b, double length_c,
                    double thick_rim_a, double thick_rim_b, double thick_rim_c);
 double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld,
@@ -9,10 +9,10 @@
             double thick_rim_c, double theta, double phi, double psi);
 
-double form_volume(double length_a, double length_b, double length_c, 
+double form_volume(double length_a, double length_b, double length_c,
                    double thick_rim_a, double thick_rim_b, double thick_rim_c)
 {
     //return length_a * length_b * length_c;
-    return length_a * length_b * length_c + 
-           2.0 * thick_rim_a * length_b * length_c + 
+    return length_a * length_b * length_c +
+           2.0 * thick_rim_a * length_b * length_c +
            2.0 * thick_rim_b * length_a * length_c +
            2.0 * thick_rim_c * length_a * length_b;
@@ -34,14 +34,14 @@
     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled
     // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
-    
+
     const double mu = 0.5 * q * length_b;
-    
+
     //calculate volume before rescaling (in original code, but not used)
-    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);        
-    //double vol = length_a * length_b * length_c + 
-    //       2.0 * thick_rim_a * length_b * length_c + 
+    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
+    //double vol = length_a * length_b * length_c +
+    //       2.0 * thick_rim_a * length_b * length_c +
     //       2.0 * thick_rim_b * length_a * length_c +
     //       2.0 * thick_rim_c * length_a * length_b;
-    
+
     // Scale sides by B
     const double a_scaled = length_a / length_b;
@@ -101,6 +101,6 @@
             //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors
             // This is the function called by csparallelepiped_analytical_2D_scaled,
-            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c        
-            
+            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c
+
             //  correct FF : sum of square of phase factors
             inner_total += Gauss76Wt[j] * form * form;
@@ -136,4 +136,7 @@
     double q, zhat, yhat, xhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
 
     // cspkernel in csparallelepiped recoded here
@@ -160,11 +163,11 @@
     double tc = length_a + 2.0*thick_rim_c;
     //handle arg=0 separately, as sin(t)/t -> 1 as t->0
-    double siA = sas_sinx_x(0.5*q*length_a*xhat);
-    double siB = sas_sinx_x(0.5*q*length_b*yhat);
-    double siC = sas_sinx_x(0.5*q*length_c*zhat);
-    double siAt = sas_sinx_x(0.5*q*ta*xhat);
-    double siBt = sas_sinx_x(0.5*q*tb*yhat);
-    double siCt = sas_sinx_x(0.5*q*tc*zhat);
-    
+    double siA = sas_sinx_x(0.5*length_a*qa);
+    double siB = sas_sinx_x(0.5*length_b*qb);
+    double siC = sas_sinx_x(0.5*length_c*qc);
+    double siAt = sas_sinx_x(0.5*ta*qa);
+    double siBt = sas_sinx_x(0.5*tb*qb);
+    double siCt = sas_sinx_x(0.5*tc*qc);
+
 
     // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed
@@ -174,5 +177,5 @@
                + drB*siA*(siBt-siB)*siC*V2
                + drC*siA*siB*(siCt*siCt-siC)*V3);
-   
+
     return 1.0e-4 * f * f;
 }
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
+++ sasmodels/models/cylinder.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,37 +1,31 @@
-double form_volume(double radius, double length);
-double fq(double q, double sn, double cn,double radius, double length);
-double orient_avg_1D(double q, double radius, double length);
-double Iq(double q, double sld, double solvent_sld, double radius, double length);
-double Iqxy(double qx, double qy, double sld, double solvent_sld,
-    double radius, double length, double theta, double phi);
-
 #define INVALID(v) (v.radius<0 || v.length<0)
 
-double form_volume(double radius, double length)
+static double
+form_volume(double radius, double length)
 {
     return M_PI*radius*radius*length;
 }
 
-double fq(double q, double sn, double cn, double radius, double length)
+static double
+fq(double qab, double qc, double radius, double length)
 {
-    // precompute qr and qh to save time in the loop
-    const double qr = q*radius;
-    const double qh = q*0.5*length; 
-    return sas_2J1x_x(qr*sn) * sas_sinx_x(qh*cn);
+    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length);
 }
 
-double orient_avg_1D(double q, double radius, double length)
+static double
+orient_avg_1D(double q, double radius, double 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; 
+    const double zb = M_PI_4;
 
     double total = 0.0;
     for (int i=0; i<76 ;i++) {
-        const double alpha = Gauss76Z[i]*zm + zb;
-        double sn, cn; // slots to hold sincos function output
-        // alpha(theta,phi) the projection of the cylinder on the detector plane
-        SINCOS(alpha, sn, cn);
-        total += Gauss76Wt[i] * square( fq(q, sn, cn, radius, length) ) * sn;
+        const double theta = Gauss76Z[i]*zm + zb;
+        double sin_theta, cos_theta; // slots to hold sincos function output
+        // theta (theta,phi) the projection of the cylinder on the detector plane
+        SINCOS(theta , sin_theta, cos_theta);
+        const double form = fq(q*sin_theta, q*cos_theta, radius, length);
+        total += Gauss76Wt[i] * form * form * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -39,5 +33,6 @@
 }
 
-double Iq(double q,
+static double
+Iq(double q,
     double sld,
     double solvent_sld,
@@ -49,6 +44,6 @@
 }
 
-
-double Iqxy(double qx, double qy,
+static double
+Iqxy(double qx, double qy,
     double sld,
     double solvent_sld,
@@ -60,7 +55,8 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
-    //printf("sn: %g cn: %g\n", sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
     const double s = (sld-solvent_sld) * form_volume(radius, length);
-    const double form = fq(q, sin_alpha, cos_alpha, radius, length);
+    const double form = fq(qab, qc, radius, length);
     return 1.0e-4 * square(s * form);
 }
Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision 3b571ae4950c5a99059d2ce8c74e80d0eb908a8f)
+++ sasmodels/models/ellipsoid.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,13 +1,10 @@
-double form_volume(double radius_polar, double radius_equatorial);
-double Iq(double q, double sld, double sld_solvent, double radius_polar, double radius_equatorial);
-double Iqxy(double qx, double qy, double sld, double sld_solvent,
-    double radius_polar, double radius_equatorial, double theta, double phi);
-
-double form_volume(double radius_polar, double radius_equatorial)
+static double
+form_volume(double radius_polar, double radius_equatorial)
 {
     return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial;
 }
 
-double Iq(double q,
+static  double
+Iq(double q,
     double sld,
     double sld_solvent,
@@ -41,5 +38,6 @@
 }
 
-double Iqxy(double qx, double qy,
+static double
+Iqxy(double qx, double qy,
     double sld,
     double sld_solvent,
@@ -51,10 +49,11 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
-    const double r = sqrt(square(radius_equatorial*sin_alpha)
-                          + square(radius_polar*cos_alpha));
-    const double f = sas_3j1x_x(q*r);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
+
+    const double qr = sqrt(square(radius_equatorial*qab) + square(radius_polar*qc));
+    const double f = sas_3j1x_x(qr);
     const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial);
 
     return 1.0e-4 * square(f * s);
 }
-
Index: sasmodels/models/elliptical_cylinder.c
===================================================================
--- sasmodels/models/elliptical_cylinder.c	(revision 61104c829b100869b42c1afd75f30ed30a934cd1)
+++ sasmodels/models/elliptical_cylinder.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,10 +1,3 @@
-double form_volume(double radius_minor, double r_ratio, double length);
-double Iq(double q, double radius_minor, double r_ratio, double length,
-          double sld, double solvent_sld);
-double Iqxy(double qx, double qy, double radius_minor, double r_ratio, double length,
-            double sld, double solvent_sld, double theta, double phi, double psi);
-
-
-double
+static double
 form_volume(double radius_minor, double r_ratio, double length)
 {
@@ -12,5 +5,5 @@
 }
 
-double
+static double
 Iq(double q, double radius_minor, double r_ratio, double length,
    double sld, double solvent_sld)
@@ -61,5 +54,5 @@
 
 
-double
+static double
 Iqxy(double qx, double qy,
      double radius_minor, double r_ratio, double length,
@@ -69,10 +62,13 @@
     double q, xhat, yhat, zhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
 
     // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2)
     // Given:    radius_major = r_ratio * radius_minor
-    const double r = radius_minor*sqrt(square(r_ratio*xhat) + square(yhat));
-    const double be = sas_2J1x_x(q*r);
-    const double si = sas_sinx_x(q*zhat*0.5*length);
+    const double qr = radius_minor*sqrt(square(r_ratio*qa) + square(qb));
+    const double be = sas_2J1x_x(qr);
+    const double si = sas_sinx_x(qc*0.5*length);
     const double Aq = be * si;
     const double delrho = sld - solvent_sld;
Index: sasmodels/models/fcc_paracrystal.c
===================================================================
--- sasmodels/models/fcc_paracrystal.c	(revision 50beefeac18f2f952bf0863526ec9bbd3009e29c)
+++ sasmodels/models/fcc_paracrystal.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,89 +1,84 @@
-double form_volume(double radius);
-double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld);
-double Iqxy(double qx, double qy, double dnn,
-    double d_factor, double radius,double sld, double solvent_sld,
-    double theta, double phi, double psi);
+static double
+_sq_fcc(double qa, double qb, double qc, double dnn, double d_factor)
+{
+    // Rewriting equations for efficiency, accuracy and readability, and so
+    // code is reusable between 1D and 2D models.
+    const double a1 = qb + qa;
+    const double a2 = qa + qc;
+    const double a3 = qb + qc;
 
-double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
-double _FCCeval(double Theta, double Phi, double temp1, double temp3);
+    const double half_dnn = 0.5*dnn;
+    const double arg = 0.5*square(half_dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+
+    // Numerator: (1 - exp(a)^2)^3
+    //         => (-(exp(2a) - 1))^3
+    //         => -expm1(2a)^3
+    // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2)
+    //         => exp(a)^2 - 2 cos(xk) exp(a) + 1
+    //         => (exp(a) - 2 cos(xk)) * exp(a) + 1
+    const double exp_arg = exp(-arg);
+    const double Sq = -cube(expm1(-2.0*arg))
+        / ( ((exp_arg - 2.0*cos(half_dnn*a1))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(half_dnn*a2))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(half_dnn*a3))*exp_arg + 1.0));
+
+    return Sq;
+}
 
 
-double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) {
-
-	const double Da = d_factor*dnn;
-	const double temp1 = q*q*Da*Da;
-	const double temp3 = q*dnn;
-
-	double retVal = _FCCeval(theta,phi,temp1,temp3)/(4.0*M_PI);
-	return(retVal);
+// occupied volume fraction calculated from lattice symmetry and sphere radius
+static double
+_fcc_volume_fraction(double radius, double dnn)
+{
+    return 4.0*sphere_volume(M_SQRT1_2*radius/dnn);
 }
 
-double _FCCeval(double Theta, double Phi, double temp1, double temp3) {
-
-	double result;
-        double sin_theta,cos_theta,sin_phi,cos_phi;
-	SINCOS(Theta, sin_theta, cos_theta);
-	SINCOS(Phi, sin_phi, cos_phi);
-
-	const double temp6 =  sin_theta;
-	const double temp7 =  sin_theta*sin_phi + cos_theta;
-	const double temp8 = -sin_theta*cos_phi + cos_theta;
-	const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi;
-
-	const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
-	result = cube(1.0-(temp10*temp10))*temp6
-	    / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10)
-	      * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10)
-	      * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10));
-
-	return (result);
-}
-
-double form_volume(double radius){
+static double
+form_volume(double radius)
+{
     return sphere_volume(radius);
 }
 
 
-double Iq(double q, double dnn,
+static double Iq(double q, double dnn,
   double d_factor, double radius,
-  double sld, double solvent_sld){
+  double sld, double solvent_sld)
+{
+    // translate a point in [-1,1] to a point in [0, 2 pi]
+    const double phi_m = M_PI;
+    const double phi_b = M_PI;
+    // translate a point in [-1,1] to a point in [0, pi]
+    const double theta_m = M_PI_2;
+    const double theta_b = M_PI_2;
 
-	//Volume fraction calculated from lattice symmetry and sphere radius
-	const double s1 = dnn*sqrt(2.0);
-	const double latticescale = 4.0*sphere_volume(radius/s1);
+    double outer_sum = 0.0;
+    for(int i=0; i<150; i++) {
+        double inner_sum = 0.0;
+        const double theta = Gauss150Z[i]*theta_m + theta_b;
+        double sin_theta, cos_theta;
+        SINCOS(theta, sin_theta, cos_theta);
+        const double qc = q*cos_theta;
+        const double qab = q*sin_theta;
+        for(int j=0;j<150;j++) {
+            const double phi = Gauss150Z[j]*phi_m + phi_b;
+            double sin_phi, cos_phi;
+            SINCOS(phi, sin_phi, cos_phi);
+            const double qa = qab*cos_phi;
+            const double qb = qab*sin_phi;
+            const double fq = _sq_fcc(qa, qb, qc, dnn, d_factor);
+            inner_sum += Gauss150Wt[j] * fq;
+        }
+        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
+        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+    }
+    outer_sum *= theta_m;
+    const double Sq = outer_sum/(4.0*M_PI);
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
 
-    const double va = 0.0;
-    const double vb = 2.0*M_PI;
-    const double vaj = 0.0;
-    const double vbj = M_PI;
-
-    double summ = 0.0;
-    double answer = 0.0;
-	for(int i=0; i<150; i++) {
-		//setup inner integral over the ellipsoidal cross-section
-		double summj=0.0;
-		const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;		//the outer dummy is phi
-		for(int j=0;j<150;j++) {
-			//20 gauss points for the inner integral
-			double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;		//the inner dummy is theta
-			double yyy = Gauss150Wt[j] * _FCC_Integrand(q,dnn,d_factor,ztheta,zphi);
-			summj += yyy;
-		}
-		//now calculate the value of the inner integral
-		double answer = (vbj-vaj)/2.0*summj;
-
-		//now calculate outer integral
-		summ = summ+(Gauss150Wt[i] * answer);
-	}		//final scaling is done at the end of the function, after the NT_FP64 case
-
-	answer = (vb-va)/2.0*summ;
-	answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale;
-
-    return answer;
+    return _fcc_volume_fraction(radius, dnn) * Pq * Sq;
+}
 
 
-}
-
-double Iqxy(double qx, double qy,
+static double Iqxy(double qx, double qy,
     double dnn, double d_factor, double radius,
     double sld, double solvent_sld,
@@ -92,21 +87,11 @@
     double q, zhat, yhat, xhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
 
-    const double a1 = yhat + xhat;
-    const double a2 = xhat + zhat;
-    const double a3 = yhat + zhat;
-    const double qd = 0.5*q*dnn;
-    const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3);
-    const double tanh_qd = tanh(arg);
-    const double cosh_qd = cosh(arg);
-    const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*a2)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*a3)/cosh_qd);
-
-    //if (isnan(Zq)) printf("q:(%g,%g) qd: %g a1: %g a2: %g a3: %g arg: %g\n", qx, qy, qd, a1, a2, a3, arg);
-
-    const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq;
-    //the occupied volume of the lattice
-    const double lattice_scale = 4.0*sphere_volume(M_SQRT1_2*radius/dnn);
-    return lattice_scale * Fq;
+    q = sqrt(qa*qa + qb*qb + qc*qc);
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+    const double Sq = _sq_fcc(qa, qb, qc, dnn, d_factor);
+    return _fcc_volume_fraction(radius, dnn) * Pq * Sq;
 }
Index: sasmodels/models/hollow_cylinder.c
===================================================================
--- sasmodels/models/hollow_cylinder.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
+++ sasmodels/models/hollow_cylinder.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,7 +1,7 @@
 double form_volume(double radius, double thickness, double length);
 double Iq(double q, double radius, double thickness, double length, double sld,
-	double solvent_sld);
+    double solvent_sld);
 double Iqxy(double qx, double qy, double radius, double thickness, double length, double sld,
-	double solvent_sld, double theta, double phi);
+    double solvent_sld, double theta, double phi);
 
 //#define INVALID(v) (v.radius_core >= v.radius)
@@ -14,17 +14,15 @@
 }
 
-
 static double
-_hollow_cylinder_kernel(double q,
-    double radius, double thickness, double length, double sin_val, double cos_val)
+_fq(double qab, double qc,
+    double radius, double thickness, double length)
 {
-    const double qs = q*sin_val;
-    const double lam1 = sas_2J1x_x((radius+thickness)*qs);
-    const double lam2 = sas_2J1x_x(radius*qs);
+    const double lam1 = sas_2J1x_x((radius+thickness)*qab);
+    const double lam2 = sas_2J1x_x(radius*qab);
     const double gamma_sq = square(radius/(radius+thickness));
-    //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 = sas_sinx_x(0.5*q*length*cos_val);
+    //Note: lim_{thickness -> 0} psi = sas_J0(radius*qab)
+    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab)
+    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq);    //SRK 10/19/00
+    const double t2 = sas_sinx_x(0.5*length*qc);
     return psi*t2;
 }
@@ -43,13 +41,13 @@
 {
     const double lower = 0.0;
-    const double upper = 1.0;		//limits of numerical integral
+    const double upper = 1.0;        //limits of numerical integral
 
-    double summ = 0.0;			//initialize intergral
+    double summ = 0.0;            //initialize intergral
     for (int i=0;i<76;i++) {
-        const double cos_val = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper );
-        const double sin_val = sqrt(1.0 - cos_val*cos_val);
-        const double inter = _hollow_cylinder_kernel(q, radius, thickness, length,
-                                                     sin_val, cos_val);
-        summ += Gauss76Wt[i] * inter * inter;
+        const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper );
+        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
+        const double form = _fq(q*sin_theta, q*cos_theta,
+                                radius, thickness, length);
+        summ += Gauss76Wt[i] * form * form;
     }
 
@@ -66,9 +64,10 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
-    const double Aq = _hollow_cylinder_kernel(q, radius, thickness, length,
-        sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
+
+    const double form = _fq(qab, qc, radius, thickness, length);
 
     const double vol = form_volume(radius, thickness, length);
-    return _hollow_cylinder_scaling(Aq*Aq, solvent_sld-sld, vol);
+    return _hollow_cylinder_scaling(form*form, solvent_sld-sld, vol);
 }
-
Index: sasmodels/models/parallelepiped.c
===================================================================
--- sasmodels/models/parallelepiped.c	(revision d6050806338b4bcae5a997acb08f4b1e7e14dd36)
+++ sasmodels/models/parallelepiped.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -20,9 +20,9 @@
 {
     const double mu = 0.5 * q * length_b;
-    
+
     // Scale sides by B
     const double a_scaled = length_a / length_b;
     const double c_scaled = length_c / length_b;
-        
+
     // outer integral (with gauss points), integration limits = 0, 1
     double outer_total = 0; //initialize integral
@@ -69,8 +69,11 @@
     double q, xhat, yhat, zhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
 
-    const double siA = sas_sinx_x(0.5*length_a*q*xhat);
-    const double siB = sas_sinx_x(0.5*length_b*q*yhat);
-    const double siC = sas_sinx_x(0.5*length_c*q*zhat);
+    const double siA = sas_sinx_x(0.5*length_a*qa);
+    const double siB = sas_sinx_x(0.5*length_b*qb);
+    const double siC = sas_sinx_x(0.5*length_c*qc);
     const double V = form_volume(length_a, length_b, length_c);
     const double drho = (sld - solvent_sld);
Index: sasmodels/models/sc_paracrystal.c
===================================================================
--- sasmodels/models/sc_paracrystal.c	(revision 50beefeac18f2f952bf0863526ec9bbd3009e29c)
+++ sasmodels/models/sc_paracrystal.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,128 +1,96 @@
-double form_volume(double radius);
+static double
+_sq_sc(double qa, double qb, double qc, double dnn, double d_factor)
+{
+    // Rewriting equations for efficiency, accuracy and readability, and so
+    // code is reusable between 1D and 2D models.
+    const double a1 = qa;
+    const double a2 = qb;
+    const double a3 = qc;
 
-double Iq(double q,
-          double dnn,
-          double d_factor,
-          double radius,
-          double sphere_sld,
-          double solvent_sld);
+    const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
 
-double Iqxy(double qx, double qy,
-            double dnn,
-            double d_factor,
-            double radius,
-            double sphere_sld,
-            double solvent_sld,
-            double theta,
-            double phi,
-            double psi);
+    // Numerator: (1 - exp(a)^2)^3
+    //         => (-(exp(2a) - 1))^3
+    //         => -expm1(2a)^3
+    // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2)
+    //         => exp(a)^2 - 2 cos(xk) exp(a) + 1
+    //         => (exp(a) - 2 cos(xk)) * exp(a) + 1
+    const double exp_arg = exp(-arg);
+    const double Sq = -cube(expm1(-2.0*arg))
+        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0));
 
-double form_volume(double radius)
+    return Sq;
+}
+
+// occupied volume fraction calculated from lattice symmetry and sphere radius
+static double
+_sc_volume_fraction(double radius, double dnn)
+{
+    return sphere_volume(radius/dnn);
+}
+
+static double
+form_volume(double radius)
 {
     return sphere_volume(radius);
 }
 
-static double
-sc_eval(double theta, double phi, double temp3, double temp4, double temp5)
+
+static double Iq(double q, double dnn,
+  double d_factor, double radius,
+  double sld, double solvent_sld)
 {
-    double cnt, snt;
-    SINCOS(theta, cnt, snt);
+    // translate a point in [-1,1] to a point in [0, 2 pi]
+    const double phi_m = M_PI_4;
+    const double phi_b = M_PI_4;
+    // translate a point in [-1,1] to a point in [0, pi]
+    const double theta_m = M_PI_4;
+    const double theta_b = M_PI_4;
 
-    double cnp, snp;
-    SINCOS(phi, cnp, snp);
 
-	double temp6 = snt;
-	double temp7 = -1.0*temp3*snt*cnp;
-	double temp8 = temp3*snt*snp;
-	double temp9 = temp3*cnt;
-	double result = temp6/((1.0-temp4*cos((temp7))+temp5)*
-	                       (1.0-temp4*cos((temp8))+temp5)*
-	                       (1.0-temp4*cos((temp9))+temp5));
-	return (result);
+    double outer_sum = 0.0;
+    for(int i=0; i<150; i++) {
+        double inner_sum = 0.0;
+        const double theta = Gauss150Z[i]*theta_m + theta_b;
+        double sin_theta, cos_theta;
+        SINCOS(theta, sin_theta, cos_theta);
+        const double qc = q*cos_theta;
+        const double qab = q*sin_theta;
+        for(int j=0;j<150;j++) {
+            const double phi = Gauss150Z[j]*phi_m + phi_b;
+            double sin_phi, cos_phi;
+            SINCOS(phi, sin_phi, cos_phi);
+            const double qa = qab*cos_phi;
+            const double qb = qab*sin_phi;
+            const double fq = _sq_sc(qa, qb, qc, dnn, d_factor);
+            inner_sum += Gauss150Wt[j] * fq;
+        }
+        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
+        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+    }
+    outer_sum *= theta_m;
+    const double Sq = outer_sum/M_PI_2;
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+
+    return _sc_volume_fraction(radius, dnn) * Pq * Sq;
 }
 
-static double
-sc_integrand(double dnn, double d_factor, double qq, double xx, double yy)
-{
-    //Function to calculate integrand values for simple cubic structure
 
-	double da = d_factor*dnn;
-	double temp1 = qq*qq*da*da;
-	double temp2 = cube(-expm1(-temp1));
-	double temp3 = qq*dnn;
-	double temp4 = 2.0*exp(-0.5*temp1);
-	double temp5 = exp(-1.0*temp1);
-
-	double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5)/M_PI_2;
-
-	return(integrand);
-}
-
-double Iq(double q,
-          double dnn,
-          double d_factor,
-          double radius,
-          double sphere_sld,
-          double solvent_sld)
-{
-	const double va = 0.0;
-	const double vb = M_PI_2; //orientation average, outer integral
-
-    double summ=0.0;
-    double answer=0.0;
-
-	for(int i=0;i<150;i++) {
-		//setup inner integral over the ellipsoidal cross-section
-		double summj=0.0;
-		double zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;
-		for(int j=0;j<150;j++) {
-			//150 gauss points for the inner integral
-			double zij = ( Gauss150Z[j]*(vb-va) + va + vb )/2.0;
-			double tmp = Gauss150Wt[j] * sc_integrand(dnn,d_factor,q,zi,zij);
-			summj += tmp;
-		}
-		//now calculate the value of the inner integral
-		answer = (vb-va)/2.0*summj;
-
-		//now calculate outer integral
-		double tmp = Gauss150Wt[i] * answer;
-		summ += tmp;
-	}		//final scaling is done at the end of the function, after the NT_FP64 case
-
-	answer = (vb-va)/2.0*summ;
-
-	//Volume fraction calculated from lattice symmetry and sphere radius
-	// NB: 4/3 pi r^3 / dnn^3 = 4/3 pi(r/dnn)^3
-	const double latticeScale = sphere_volume(radius/dnn);
-
-	answer *= sphere_form(q, radius, sphere_sld, solvent_sld)*latticeScale;
-
-	return answer;
-}
-
-double Iqxy(double qx, double qy,
-          double dnn,
-          double d_factor,
-          double radius,
-          double sphere_sld,
-          double solvent_sld,
-          double theta,
-          double phi,
-          double psi)
+static double Iqxy(double qx, double qy,
+    double dnn, double d_factor, double radius,
+    double sld, double solvent_sld,
+    double theta, double phi, double psi)
 {
     double q, zhat, yhat, xhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
 
-    const double qd = q*dnn;
-    const double arg = 0.5*square(qd*d_factor);
-    const double tanh_qd = tanh(arg);
-    const double cosh_qd = cosh(arg);
-    const double Zq = tanh_qd/(1. - cos(qd*zhat)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*yhat)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*xhat)/cosh_qd);
-
-    const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq;
-    //the occupied volume of the lattice
-    const double lattice_scale = sphere_volume(radius/dnn);
-    return lattice_scale * Fq;
+    q = sqrt(qa*qa + qb*qb + qc*qc);
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+    const double Sq = _sq_sc(qa, qb, qc, dnn, d_factor);
+    return _sc_volume_fraction(radius, dnn) * Pq * Sq;
 }
Index: sasmodels/models/sc_paracrystal.py
===================================================================
--- sasmodels/models/sc_paracrystal.py	(revision 69e1afc57270f67b2a0d7157a6d263fc2db98cc0)
+++ sasmodels/models/sc_paracrystal.py	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -152,8 +152,6 @@
     [{}, 0.001, 10.3048],
     [{}, 0.215268, 0.00814889],
-    [{}, (0.414467), 0.001313289],
+    [{}, 0.414467, 0.001313289],
     [{'theta':10.0,'phi':20,'psi':30.0},(0.045,-0.035),18.0397138402 ],
     [{'theta':10.0,'phi':20,'psi':30.0},(0.023,0.045),0.0177333171285 ]
     ]
-
-
Index: sasmodels/models/stacked_disks.c
===================================================================
--- sasmodels/models/stacked_disks.c	(revision 19f996bbf3615f0161b71f00777e5f16155bef70)
+++ sasmodels/models/stacked_disks.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,4 +1,5 @@
 static double stacked_disks_kernel(
-    double q,
+    double qab,
+    double qc,
     double halfheight,
     double thick_layer,
@@ -9,6 +10,4 @@
     double layer_sld,
     double solvent_sld,
-    double sin_alpha,
-    double cos_alpha,
     double d)
 
@@ -20,9 +19,9 @@
     // zi is the dummy variable for the integration (x in Feigin's notation)
 
-    const double besarg1 = q*radius*sin_alpha;
-    //const double besarg2 = q*radius*sin_alpha;
+    const double besarg1 = radius*qab;
+    //const double besarg2 = radius*qab;
 
-    const double sinarg1 = q*halfheight*cos_alpha;
-    const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha;
+    const double sinarg1 = halfheight*qc;
+    const double sinarg2 = (halfheight+thick_layer)*qc;
 
     const double be1 = sas_2J1x_x(besarg1);
@@ -43,5 +42,5 @@
 
     // loop for the structure factor S(q)
-    double qd_cos_alpha = q*d*cos_alpha;
+    double qd_cos_alpha = d*qc;
     //d*cos_alpha is the projection of d onto q (in other words the component
     //of d that is parallel to q.
@@ -84,5 +83,5 @@
         double sin_alpha, cos_alpha; // slots to hold sincos function output
         SINCOS(zi, sin_alpha, cos_alpha);
-        double yyy = stacked_disks_kernel(q,
+        double yyy = stacked_disks_kernel(q*sin_alpha, q*cos_alpha,
                            halfheight,
                            thick_layer,
@@ -93,6 +92,4 @@
                            layer_sld,
                            solvent_sld,
-                           sin_alpha,
-                           cos_alpha,
                            d);
         summ += Gauss76Wt[i] * yyy * sin_alpha;
@@ -155,8 +152,10 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
     double d = 2.0 * thick_layer + thick_core;
     double halfheight = 0.5*thick_core;
-    double answer = stacked_disks_kernel(q,
+    double answer = stacked_disks_kernel(qab, qc,
                      halfheight,
                      thick_layer,
@@ -167,6 +166,4 @@
                      layer_sld,
                      solvent_sld,
-                     sin_alpha,
-                     cos_alpha,
                      d);
 
Index: sasmodels/models/triaxial_ellipsoid.c
===================================================================
--- sasmodels/models/triaxial_ellipsoid.c	(revision 68dd6a95741269738dcbfb051d93fe9f26b89c43)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,17 +1,12 @@
-double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar);
-double Iq(double q, double sld, double sld_solvent,
-    double radius_equat_minor, double radius_equat_major, double radius_polar);
-double Iqxy(double qx, double qy, double sld, double sld_solvent,
-    double radius_equat_minor, double radius_equat_major, double radius_polar, double theta, double phi, double psi);
-
 //#define INVALID(v) (v.radius_equat_minor > v.radius_equat_major || v.radius_equat_major > v.radius_polar)
 
-
-double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar)
+static double
+form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar)
 {
     return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar;
 }
 
-double Iq(double q,
+static double
+Iq(double q,
     double sld,
     double sld_solvent,
@@ -45,9 +40,11 @@
     // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi
     const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI);
-    const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar);
-    return 1.0e-4 * s * s * fqsq;
+    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar);
+    const double drho = (sld - sld_solvent);
+    return 1.0e-4 * square(vol*drho) * fqsq;
 }
 
-double Iqxy(double qx, double qy,
+static double
+Iqxy(double qx, double qy,
     double sld,
     double sld_solvent,
@@ -61,12 +58,15 @@
     double q, xhat, yhat, zhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
 
-    const double r = sqrt(square(radius_equat_minor*xhat)
-                          + square(radius_equat_major*yhat)
-                          + square(radius_polar*zhat));
-    const double fq = sas_3j1x_x(q*r);
-    const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar);
+    const double qr = sqrt(square(radius_equat_minor*qa)
+                           + square(radius_equat_major*qb)
+                           + square(radius_polar*qc));
+    const double fq = sas_3j1x_x(qr);
+    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar);
+    const double drho = (sld - sld_solvent);
 
-    return 1.0e-4 * square(s * fq);
+    return 1.0e-4 * square(vol * drho * fq);
 }
-
