Index: sasmodels/models/bcc.c
===================================================================
--- sasmodels/models/bcc.c	(revision 7ed702f2d0d4ecbcb0115970cdabf0d05ac08de0)
+++ sasmodels/models/bcc.c	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
@@ -40,12 +40,4 @@
 }
 
-double _sphereform(double q, double radius, double sld, double solvent_sld){
-    const double qr = q*radius;
-    double sn, cn;
-    SINCOS(qr, sn, cn);
-    const double bes = (qr == 0.0 ? 1.0 : 3.0*(sn-qr*cn)/(qr*qr*qr));
-    const double fq = bes * (sld - solvent_sld)*form_volume(radius);
-    return 1.0e-4*fq*fq;
-}
 
 double form_volume(double radius){
@@ -87,5 +79,5 @@
 
 	answer = (vb-va)/2.0*summ;
-	answer = answer*_sphereform(q,radius,sld,solvent_sld)*latticescale;
+	answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale;
 
     return answer;
@@ -174,5 +166,5 @@
 
   // Use SphereForm directly from libigor
-  answer = _sphereform(q,radius,sld,solvent_sld)*Zq*latticescale;
+  answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale;
 
   return answer;
Index: sasmodels/models/bcc.py
===================================================================
--- sasmodels/models/bcc.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
+++ sasmodels/models/bcc.py	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
@@ -132,5 +132,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/J1.c", "lib/gauss150.c", "bcc.c"]
+source = ["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 eeb8bacf06af69341f80d5e76d7b10a7239c4e67)
+++ sasmodels/models/fcc.c	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
@@ -7,5 +7,4 @@
 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
 double _FCCeval(double Theta, double Phi, double temp1, double temp3);
-double _sphereform(double q, double radius, double sld, double solvent_sld);
 
 
@@ -39,13 +38,4 @@
 
 	return (result);
-}
-
-double _sphereform(double q, double radius, double sld, double solvent_sld){
-    const double qr = q*radius;
-    double sn, cn;
-    SINCOS(qr, sn, cn);
-    const double bes = (qr == 0.0 ? 1.0 : 3.0*(sn-qr*cn)/(qr*qr*qr));
-    const double fq = bes * (sld - solvent_sld)*form_volume(radius);
-    return 1.0e-4*fq*fq;
 }
 
@@ -88,5 +78,5 @@
 
 	answer = (vb-va)/2.0*summ;
-	answer = answer*_sphereform(q,radius,sld,solvent_sld)*latticescale;
+	answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale;
 
     return answer;
@@ -175,5 +165,5 @@
 
   // Use SphereForm directly from libigor
-  answer = _sphereform(q,radius,sld,solvent_sld)*Zq*latticescale;
+  answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale;
 
   return answer;
Index: sasmodels/models/fcc.py
===================================================================
--- sasmodels/models/fcc.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
+++ sasmodels/models/fcc.py	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
@@ -125,5 +125,5 @@
              ]
 
-source = ["lib/J1.c", "lib/gauss150.c", "fcc.c"]
+source = ["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 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
@@ -0,0 +1,26 @@
+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;
+}
+
+// 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 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
@@ -0,0 +1,228 @@
+double form_volume(double radius);
+
+double Iq(double q,
+          double dnn,
+          double d_factor,
+          double radius,
+          double sphere_sld,
+          double solvent_sld);
+
+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);
+
+double form_volume(double radius)
+{
+    return 1.333333333333333*M_PI*radius*radius*radius;
+}
+
+static double
+sc_eval(double theta, double phi, double temp3, double temp4, double temp5)
+{
+    double cnt, snt;
+    SINCOS(theta, cnt, snt);
+
+    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);
+}
+
+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 = pow( 1.0-exp(-1.0*temp1) ,3);
+	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);
+	integrand *= 2.0/M_PI;
+
+	return(integrand);
+}
+
+static
+double sc_crystal_kernel(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.0; //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
+	const double latticeScale = (4.0/3.0)*M_PI*(radius*radius*radius)/pow(dnn,3);
+
+	//answer *= sphere_form_paracrystal(q, radius,contrast)*latticeScale;
+	answer *= sphere_form(q, radius, sphere_sld, solvent_sld)*latticeScale;
+
+	return answer;
+}
+
+static
+double sc_crystal_kernel_2d(double q, double q_x, double q_y,
+          double dnn,
+          double d_factor,
+          double radius,
+          double sphere_sld,
+          double solvent_sld,
+          double theta,
+          double phi,
+          double psi)
+{
+    //convert angle degree to radian
+    theta = theta * M_PI_180;
+    phi = phi * M_PI_180;
+    psi = psi * M_PI_180;
+
+    const double qda_2 = pow(q*d_factor*dnn,2.0);
+
+    double snt, cnt;
+    SINCOS(theta, snt, cnt);
+
+    double snp, cnp;
+    SINCOS(phi, snp, cnp);
+
+    double sns, cns;
+    SINCOS(psi, sns, cns);
+
+    /// Angles here are respect to detector coordinate instead of against
+    //  q coordinate(PRB 36, 3, 1754)
+    // a3 axis orientation
+
+    const double a3_x = cnt * cnp;
+    const double a3_y = snt;
+
+    // Compute the angle btw vector q and the a3 axis
+    double cos_val_a3 = a3_x*q_x + a3_y*q_y;
+
+    // a1 axis orientation
+    const double a1_x = -cnp*sns * snt+snp*cns;
+    const double a1_y = sns*cnt;
+
+    double cos_val_a1 = a1_x*q_x + a1_y*q_y;
+
+    // a2 axis orientation
+    const double a2_x = -snt*cns*cnp-sns*snp;
+    const double a2_y = cnt*cns;
+
+    // a2 axis
+    const double cos_val_a2 =  a2_x*q_x + a2_y*q_y;
+
+    // The following test should always pass
+    if (fabs(cos_val_a3)>1.0) {
+        //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+        cos_val_a3 = 1.0;
+    }
+    if (fabs(cos_val_a1)>1.0) {
+        //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+        cos_val_a1 = 1.0;
+    }
+    if (fabs(cos_val_a2)>1.0) {
+        //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+        cos_val_a3 = 1.0;
+    }
+
+    const double a3_dot_q = dnn*q*cos_val_a3;
+    const double a1_dot_q = dnn*q*cos_val_a1;
+    const double a2_dot_q = dnn*q*cos_val_a2;
+
+    // Call Zq=Z1*Z2*Z3
+    double Zq = (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a1_dot_q)+exp(-qda_2));
+    Zq *= (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a2_dot_q)+exp(-qda_2));
+    Zq *= (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a3_dot_q)+exp(-qda_2));
+
+    // Use SphereForm directly from libigor
+    double answer = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq;
+
+    //consider scales
+    const double latticeScale = (4.0/3.0)*M_PI*(radius*radius*radius)/pow(dnn,3.0);
+    answer *= latticeScale;
+
+    return answer;
+}
+
+double Iq(double q,
+          double dnn,
+          double d_factor,
+          double radius,
+          double sphere_sld,
+          double solvent_sld)
+{
+    return sc_crystal_kernel(q,
+              dnn,
+              d_factor,
+              radius,
+              sphere_sld,
+              solvent_sld);
+}
+
+// Iqxy is never called since no orientation or magnetic parameters.
+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)
+{
+    double q = sqrt(qx*qx + qy*qy);
+
+
+    return sc_crystal_kernel_2d(q, qx/q, qy/q,
+                  dnn,
+                  d_factor,
+                  radius,
+                  sphere_sld,
+                  solvent_sld,
+                  theta,
+                  phi,
+                  psi);
+
+}
+
Index: sasmodels/models/sc_crystal.py
===================================================================
--- sasmodels/models/sc_crystal.py	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
+++ sasmodels/models/sc_crystal.py	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
@@ -0,0 +1,168 @@
+r"""
+Calculates the scattering from a **simple cubic lattice** with
+paracrystalline distortion. Thermal vibrations are considered to be
+negligible, and the size of the paracrystal is infinitely large.
+Paracrystalline distortion is assumed to be isotropic and characterized
+by a Gaussian distribution.
+
+Definition
+----------
+
+The scattering intensity $I(q)$ is calculated as
+
+.. math::
+
+    I(q) = \frac{scale}{V_p}V_{lattice}P(q)Z(q)
+
+where scale is the volume fraction of spheres, $V_p$ is the volume of
+the primary particle, $V_{lattice}$ is a volume correction for the crystal
+structure, $P(q)$ is the form factor of the sphere (normalized), and
+$Z(q)$ is the paracrystalline structure factor for a simple cubic structure.
+
+Equation (16) of the 1987 reference is used to calculate $Z(q)$, using
+equations (13)-(15) from the 1987 paper for Z1, Z2, and Z3.
+
+The lattice correction (the occupied volume of the lattice) for a simple
+cubic structure of particles of radius *R* and nearest neighbor separation *D* is
+
+.. math::
+
+    V_{lattice}=\frac{4\pi}{3}\frac{R^3}{D^3}
+
+The distortion factor (one standard deviation) of the paracrystal is included
+in the calculation of $Z(q)$
+
+.. math::
+
+    \Delta a = gD
+
+where *g* is a fractional distortion based on the nearest neighbor distance.
+
+The simple cubic lattice is
+
+.. figure:: img/sc_crystal_fig1.jpg
+
+For a crystal, diffraction peaks appear at reduced q-values given by
+
+.. math::
+
+    \frac{qD}{2\pi} = \sqrt{h^2+k^2+l^2}
+
+where for a simple cubic lattice any h, k, l are allowed and none are
+forbidden. Thus the peak positions correspond to (just the first 5)
+
+.. math::
+
+    \begin{align*}
+    q/q_0 \quad & \quad 1
+                & \sqrt{2} \quad
+                & \quad  \sqrt{3} \quad
+                & \sqrt{4} \quad
+                & \quad \sqrt{5}\quad \\
+    Indices \quad & (100)
+                  & \quad (110) \quad
+                  & \quad (111)
+                  & (200) \quad
+                  & \quad (210)
+    \end{align*}
+
+.. note::
+
+    The calculation of *Z(q)* is a double numerical integral that must be
+    carried out with a high density of points to properly capture the sharp
+    peaks of the paracrystalline scattering.
+    So be warned that the calculation is SLOW. Go get some coffee.
+    Fitting of any experimental data must be resolution smeared for any
+    meaningful fit. This makes a triple integral. Very, very slow.
+    Go get lunch!
+
+This example dataset is produced using 200 data points,
+$q_{min} = 0.01A^{-1}, q_{max} = 0.1A^{-1}$ and the above default values.
+
+.. figure:: img/sc_crystal_1d.jpg
+
+    1D plot in the linear scale using the default values (w/200 data point).
+
+The 2D (Anisotropic model) is based on the reference below where *I(q)* is
+approximated for 1d scattering. Thus the scattering pattern for 2D may not
+be accurate. Note that we are not responsible for any incorrectness of the 2D
+model computation.
+
+.. figure:: img/sc_crystal_fig2.jpg
+.. figure:: img/sc_crystal_fig3.jpg
+
+    2D plot using the default values (w/200X200 pixels).
+
+Reference
+---------
+Hideki Matsuoka et. al. *Physical Review B,* 36 (1987) 1754-1765
+(Original Paper)
+
+Hideki Matsuoka et. al. *Physical Review B,* 41 (1990) 3854 -3856
+(Corrections to FCC and BCC lattice structure calculation)
+
+"""
+
+from numpy import inf
+
+name = "sc_crystal"
+title = "Simple cubic lattice with paracrystalline distortion"
+description = """
+        P(q)=(scale/Vp)*V_lattice*P(q)*Z(q)+bkg where scale is the volume
+        fraction of sphere,
+        Vp = volume of the primary particle,
+        V_lattice = volume correction for
+        for the crystal structure,
+        P(q)= form factor of the sphere (normalized),
+        Z(q)= paracrystalline structure factor
+        for a simple cubic structure.
+        [Simple Cubic ParaCrystal Model]
+        Parameters;
+        scale: volume fraction of spheres
+        bkg:background, R: radius of sphere
+        dnn: Nearest neighbor distance
+        d_factor: Paracrystal distortion factor
+        radius: radius of the spheres
+        sldSph: SLD of the sphere
+        sldSolv: SLD of the solvent
+        """
+category = "shape:paracrystal"
+single = False
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["dnn",         "Ang",       220.0,  [0.0, inf],  "",            "Nearest neighbor distance"],
+              ["d_factor",    "",            0.06, [-inf, inf], "",            "Paracrystal distortion factor"],
+              ["radius",      "Ang",        40.0,  [0.0, inf],  "volume",      "Radius of sphere"],
+              ["sphere_sld",  "1e-6/Ang^2",  3.0,  [0.0, inf],  "",            "Sphere scattering length density"],
+              ["solvent_sld", "1e-6/Ang^2",  6.3,  [0.0, inf],  "",            "Solvent scattering length density"],
+              ["theta",       "degrees",     0.0,  [-inf, inf], "orientation", "Orientation of the a1 axis w/respect incoming beam"],
+              ["phi",         "degrees",     0.0,  [-inf, inf], "orientation", "Orientation of the a2 in the plane of the detector"],
+              ["psi",         "degrees",     0.0,  [-inf, inf], "orientation", "Orientation of the a3 in the plane of the detector"],
+             ]
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["lib/sphere_form.c", "lib/gauss150.c", "sc_crystal.c"]
+
+demo = dict(scale=1, background=0,
+            dnn=220.0,
+            d_factor=0.06,
+            radius=40.0,
+            sphere_sld=3.0,
+            solvent_sld=6.3,
+            theta=0.0,
+            phi=0.0,
+            psi=0.0)
+
+oldname = 'SCCrystalModel'
+
+oldpars = dict(sphere_sld='sldSph',
+               solvent_sld='sldSolv')
+
+tests = [
+    # Accuracy tests based on content in test/utest_extra_models.py
+    [{}, 0.001, 10.3038],
+    [{}, 0.215268, 0.00714889],
+    [{}, (0.414467), 0.000313289]
+    ]
+
+
