Index: sasmodels/models/bcc.c
===================================================================
--- sasmodels/models/bcc.c	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
+++ sasmodels/models/bcc.c	(revision ad90df9f3ed42019113506f1366fe4b851b5f389)
@@ -40,7 +40,6 @@
 }
 
-
 double form_volume(double radius){
-    return 1.333333333333333*M_PI*radius*radius*radius;
+    return sphere_volume(radius);
 }
 
@@ -52,5 +51,5 @@
 	//Volume fraction calculated from lattice symmetry and sphere radius
 	const double s1 = dnn/sqrt(0.75);
-	const double latticescale = 2.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1);
+	const double latticescale = 2.0*sphere_volume(radius/s1);
 
     const double va = 0.0;
@@ -113,5 +112,5 @@
 
   //the occupied volume of the lattice
-  const double latticescale = 2.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1);
+  const double latticescale = 2.0*sphere_volume(radius/s1);
   // q vector
   //q_z = 0.0; // for SANS; assuming qz is negligible
Index: sasmodels/models/bcc.py
===================================================================
--- sasmodels/models/bcc.py	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
+++ sasmodels/models/bcc.py	(revision ad90df9f3ed42019113506f1366fe4b851b5f389)
@@ -132,5 +132,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/J1.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc.c"]
+source = ["lib/sph_j1c.c", "lib/J1.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc.c"]
 
 # parameters for demo
Index: sasmodels/models/fcc.c
===================================================================
--- sasmodels/models/fcc.c	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
+++ sasmodels/models/fcc.c	(revision ad90df9f3ed42019113506f1366fe4b851b5f389)
@@ -41,5 +41,5 @@
 
 double form_volume(double radius){
-    return 1.333333333333333*M_PI*radius*radius*radius;
+    return sphere_volume(radius);
 }
 
@@ -51,5 +51,5 @@
 	//Volume fraction calculated from lattice symmetry and sphere radius
 	const double s1 = dnn*sqrt(2.0);
-	const double latticescale = 4.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1);
+	const double latticescale = 4.0*sphere_volume(radius/s1);
 
     const double va = 0.0;
@@ -112,5 +112,5 @@
 
   //the occupied volume of the lattice
-  const double latticescale = 2.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1);
+  const double latticescale = 2.0*sphere_volume(radius/s1);
   // q vector
   // q_z = 0.0; // for SANS; assuming qz is negligible
Index: sasmodels/models/fcc.py
===================================================================
--- sasmodels/models/fcc.py	(revision c0ccea8f892e537b9ac3a573d9883875986f4b15)
+++ sasmodels/models/fcc.py	(revision ad90df9f3ed42019113506f1366fe4b851b5f389)
@@ -127,5 +127,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/J1.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc.c"]
+source = ["lib/sph_j1c.c", "lib/J1.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc.c"]
 
 # parameters for demo
Index: sasmodels/models/lib/sphere_form.c
===================================================================
--- sasmodels/models/lib/sphere_form.c	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
+++ sasmodels/models/lib/sphere_form.c	(revision ad90df9f3ed42019113506f1366fe4b851b5f389)
@@ -1,26 +1,13 @@
-double
+inline double
+sphere_volume(double radius)
+{
+    return M_4PI_3*cube(radius);
+}
+
+inline double
 sphere_form(double q, double radius, double sld, double solvent_sld)
 {
-    const double qr = q*radius;
-    double sn, cn;
-    SINCOS(qr, sn, cn);
-    const double vol = 1.333333333333333*M_PI*radius*radius*radius;
-    const double bes = (qr == 0.0 ? 1.0 : 3.0*(sn-qr*cn)/(qr*qr*qr));
-    const double fq = bes * (sld - solvent_sld)*vol;
-    return 1.0e-4*fq*fq;
+    const double fq = sphere_volume(radius) * sph_j1c(q*radius);
+    const double contrast = (sld - solvent_sld);
+    return 1.0e-4*square(contrast * fq);
 }
-
-// Do we even need the paracrystal form?
-// The only difference is the return value at r=0
-double
-sphere_form_paracrystal(double q, double radius, double delrho)
-{
-    const double qr = q*radius;
-    const double vol = 1.333333333333333*M_PI*radius*radius*radius;
-	const double f = vol*delrho*delrho*1.0e-4;
-    double sn, cn;
-    SINCOS(qr, sn, cn);
-	const double bes = (qr == 0.0 ? f : 3.0*(sn-qr*cn)/(qr*qr*qr));
-	const double fq = bes*delrho*vol;
-	return fq*fq * 1.0e-4;
-}
Index: sasmodels/models/sc_crystal.c
===================================================================
--- sasmodels/models/sc_crystal.c	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
+++ sasmodels/models/sc_crystal.c	(revision ad90df9f3ed42019113506f1366fe4b851b5f389)
@@ -20,5 +20,5 @@
 double form_volume(double radius)
 {
-    return 1.333333333333333*M_PI*radius*radius*radius;
+    return sphere_volume(radius);
 }
 
@@ -95,7 +95,7 @@
 
 	//Volume fraction calculated from lattice symmetry and sphere radius
-	const double latticeScale = (4.0/3.0)*M_PI*(radius*radius*radius)/pow(dnn,3);
-
-	//answer *= sphere_form_paracrystal(q, radius,contrast)*latticeScale;
+	// 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;
 
@@ -180,5 +180,5 @@
 
     //consider scales
-    const double latticeScale = (4.0/3.0)*M_PI*(radius*radius*radius)/pow(dnn,3.0);
+    const double latticeScale = sphere_volume(radius/dnn);
     answer *= latticeScale;
 
Index: sasmodels/models/sc_crystal.py
===================================================================
--- sasmodels/models/sc_crystal.py	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
+++ sasmodels/models/sc_crystal.py	(revision ad90df9f3ed42019113506f1366fe4b851b5f389)
@@ -142,5 +142,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sphere_form.c", "lib/gauss150.c", "sc_crystal.c"]
+source = ["lib/sph_j1c.c", "lib/sphere_form.c", "lib/gauss150.c", "sc_crystal.c"]
 
 demo = dict(scale=1, background=0,
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision 50e1e40cf14d4ded28bf96e61a3cfaf35003b899)
+++ sasmodels/models/sphere.py	(revision ad90df9f3ed42019113506f1366fe4b851b5f389)
@@ -78,17 +78,14 @@
              ]
 
-source = ["lib/sph_j1c.c"]
+source = ["lib/sph_j1c.c", "lib/sphere_form.c"]
 
 # No volume normalization despite having a volume parameter
 # This should perhaps be volume normalized?
 form_volume = """
-    return M_4PI_3*cube(radius);
+    return sphere_volume(radius);
     """
 
 Iq = """
-    const double qr = q*radius;
-    const double bes = sph_j1c(qr);
-    const double fq = bes * (sld - solvent_sld) * form_volume(radius);
-    return 1.0e-4*fq*fq;
+    return sphere_form(q, radius, sld, solvent_sld);
     """
 
