Index: sasmodels/models/barbell.py
===================================================================
--- sasmodels/models/barbell.py	(revision 26141cba42ea53a3ca3b54e60b579f5525dd583d)
+++ sasmodels/models/barbell.py	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
@@ -101,6 +101,6 @@
 # pylint: disable=bad-whitespace, line-too-long
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["sld",         "4e-6/Ang^2",   4, [-inf, inf], "",            "Barbell scattering length density"],
-              ["solvent_sld", "1e-6/Ang^2",   1, [-inf, inf], "",            "Solvent scattering length density"],
+parameters = [["sld",         "1e-6/Ang^2",   4, [-inf, inf], "",            "Barbell scattering length density"],
+              ["sld_solvent", "1e-6/Ang^2",   1, [-inf, inf], "",            "Solvent scattering length density"],
               ["bell_radius", "Ang",         40, [0, inf],    "volume",      "Spherical bell radius"],
               ["radius",      "Ang",         20, [0, inf],    "volume",      "Cylindrical bar radius"],
@@ -115,5 +115,5 @@
 # parameters for demo
 demo = dict(scale=1, background=0,
-            sld=6, solvent_sld=1,
+            sld=6, sld_solvent=1,
             bell_radius=40, radius=20, length=400,
             theta=60, phi=60,
@@ -128,4 +128,4 @@
 oldname = 'BarBellModel'
 oldpars = dict(sld='sld_barbell',
-               solvent_sld='sld_solv', bell_radius='rad_bell',
+               sld_solvent='sld_solv', bell_radius='rad_bell',
                radius='rad_bar', length='len_bar')
Index: smodels/models/bcc.c
===================================================================
--- sasmodels/models/bcc.c	(revision ad90df9f3ed42019113506f1366fe4b851b5f389)
+++ 	(revision )
@@ -1,170 +1,0 @@
-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);
-
-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);
-
-
-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);
-}
-
-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 = pow(1.0-(temp10*temp10),3)*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){
-    return sphere_volume(radius);
-}
-
-
-double Iq(double q, double dnn,
-  double d_factor, double radius,
-  double sld, double solvent_sld){
-
-	//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);
-
-    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;
-
-
-}
-
-
-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 b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
-  //double q_z;
-  double cos_val_b3, cos_val_b2, cos_val_b1;
-  double a1_dot_q, a2_dot_q,a3_dot_q;
-  double answer;
-  double Zq, Fkq, Fkq_2;
-
-  //convert to q and make scaled values
-  double q = sqrt(qx*qx+qy*qy);
-  double q_x = qx/q;
-  double q_y = qy/q;
-
-  //convert angle degree to radian
-  theta = theta * M_PI_180;
-  phi = phi * M_PI_180;
-  psi = psi * M_PI_180;
-
-  const double Da = d_factor*dnn;
-  const double s1 = dnn/sqrt(0.75);
-
-
-  //the occupied volume of the lattice
-  const double latticescale = 2.0*sphere_volume(radius/s1);
-  // q vector
-  //q_z = 0.0; // for SANS; assuming qz is negligible
-  /// Angles here are respect to detector coordinate
-  ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
-    // b3 axis orientation
-    b3_x = cos(theta) * cos(phi);
-    b3_y = sin(theta);
-    //b3_z = -cos(theta) * sin(phi);
-    cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z;
-
-    //alpha = acos(cos_val_b3);
-    // b1 axis orientation
-    b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
-    b1_y = sin(psi)*cos(theta);
-    cos_val_b1 = b1_x*q_x + b1_y*q_y;
-    // b2 axis orientation
-    b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
-  	b2_y = cos(theta)*cos(psi);
-    cos_val_b2 = b2_x*q_x + b2_y*q_y;
-
-    // The following test should always pass
-    if (fabs(cos_val_b3)>1.0) {
-      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
-      cos_val_b3 = 1.0;
-    }
-    if (fabs(cos_val_b2)>1.0) {
-      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
-      cos_val_b2 = 1.0;
-    }
-    if (fabs(cos_val_b1)>1.0) {
-      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
-      cos_val_b1 = 1.0;
-    }
-    // Compute the angle btw vector q and the a3 axis
-    a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3);
-
-    // a1 axis
-    a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1);
-
-    // a2 axis
-    a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2);
-
-
-    // Get Fkq and Fkq_2
-    Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
-    Fkq_2 = Fkq*Fkq;
-    // Call Zq=Z1*Z2*Z3
-    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
-    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
-    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2);
-
-  // Use SphereForm directly from libigor
-  answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale;
-
-  return answer;
- }
Index: smodels/models/bcc.py
===================================================================
--- sasmodels/models/bcc.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
+++ 	(revision )
@@ -1,143 +1,0 @@
-#bcc paracrystal model
-#note model title and parameter table are automatically inserted
-#note - calculation requires double precision
-r"""
-Calculates the scattering from a **body-centered 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{\text{scale}}{V_p} V_\text{lattice} P(q) Z(q)
-
-
-where *scale* is the volume fraction of spheres, $V_p$ is the volume of the
-primary particle, $V_\text{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 body-centered cubic structure.
-
-Equation (1) of the 1990 reference is used to calculate $Z(q)$, using
-equations (29)-(31) from the 1987 paper for $Z1$, $Z2$, and $Z3$.
-
-The lattice correction (the occupied volume of the lattice) for a
-body-centered cubic structure of particles of radius $R$ and nearest neighbor
-separation $D$ is
-
-.. math::
-
-    V_\text{lattice} = \frac{16\pi}{3} \frac{R^3}{\left(D\sqrt{2}\right)^3}
-
-
-The distortion factor (one standard deviation) of the paracrystal is included
-in the calculation of $Z(q)$
-
-.. math::
-
-    \Delta a = g D
-
-where $g$ is a fractional distortion based on the nearest neighbor distance.
-
-
-.. figure:: img/bcc_geometry.jpg
-
-    Body-centered cubic lattice.
-
-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 body-centered cubic lattice, only reflections where
-$(h + k + l) = \text{even}$ are allowed and reflections where
-$(h + k + l) = \text{odd}$ are forbidden. Thus the peak positions
-correspond to (just the first 5)
-
-.. math::
-
-    \begin{array}{lccccc}
-    q/q_o          &   1   & \sqrt{2} & \sqrt{3} & \sqrt{4} & \sqrt{5} \\
-    \text{Indices} & (110) &    (200) & (211)    & (220)    & (310)    \\
-    \end{array}
-
-**NB**: 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,
-*qmin* = 0.001 |Ang^-1|, *qmax* = 0.1 |Ang^-1| and the above default values.
-
-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/bcc_angle_definition.png
-
-    Orientation of the crystal with respect to the scattering plane.
-
-References
-----------
-
-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 = "bcc_paracrystal"
-title = "Body-centred cubic lattic with paracrystalline distortion"
-description = """
-    Calculates the scattering from a **body-centered 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.
-    """
-category = "shape:paracrystal"
-
-single = False
-
-# pylint: disable=bad-whitespace, line-too-long
-#             ["name", "units", default, [lower, upper], "type","description" ],
-parameters = [["dnn",         "Ang",       220,    [-inf, inf], "",            "Nearest neighbour distance"],
-              ["d_factor",    "",            0.06, [-inf, inf], "",            "Paracrystal distortion factor"],
-              ["radius",      "Ang",        40,    [0, inf],    "volume",      "Particle radius"],
-              ["sld",         "1e-6/Ang^2",  4,    [-inf, inf], "",            "Particle scattering length density"],
-              ["solvent_sld", "1e-6/Ang^2",  1,    [-inf, inf], "",            "Solvent scattering length density"],
-              ["theta",       "degrees",    60,    [-inf, inf], "orientation", "In plane angle"],
-              ["phi",         "degrees",    60,    [-inf, inf], "orientation", "Out of plane angle"],
-              ["psi",         "degrees",    60,    [-inf, inf], "orientation", "Out of plane angle"]
-             ]
-# pylint: enable=bad-whitespace, line-too-long
-
-source = ["lib/sph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc.c"]
-
-# parameters for demo
-demo = dict(
-    scale=1, background=0,
-    dnn=220, d_factor=0.06, sld=4, solvent_sld=1,
-    radius=40,
-    theta=60, phi=60, psi=60,
-    radius_pd=.2, radius_pd_n=2,
-    theta_pd=15, theta_pd_n=0,
-    phi_pd=15, phi_pd_n=0,
-    psi_pd=15, psi_pd_n=0,
-    )
-
-# For testing against the old sasview models, include the converted parameter
-# names and the target sasview model name.
-oldname = 'BCCrystalModel'
-oldpars = dict(sld='sldSph', solvent_sld='sldSolv')
Index: sasmodels/models/bcc_paracrystal.py
===================================================================
--- sasmodels/models/bcc_paracrystal.py	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
+++ sasmodels/models/bcc_paracrystal.py	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
@@ -0,0 +1,143 @@
+#bcc paracrystal model
+#note model title and parameter table are automatically inserted
+#note - calculation requires double precision
+r"""
+Calculates the scattering from a **body-centered 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{\text{scale}}{V_p} V_\text{lattice} P(q) Z(q)
+
+
+where *scale* is the volume fraction of spheres, $V_p$ is the volume of the
+primary particle, $V_\text{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 body-centered cubic structure.
+
+Equation (1) of the 1990 reference is used to calculate $Z(q)$, using
+equations (29)-(31) from the 1987 paper for $Z1$, $Z2$, and $Z3$.
+
+The lattice correction (the occupied volume of the lattice) for a
+body-centered cubic structure of particles of radius $R$ and nearest neighbor
+separation $D$ is
+
+.. math::
+
+    V_\text{lattice} = \frac{16\pi}{3} \frac{R^3}{\left(D\sqrt{2}\right)^3}
+
+
+The distortion factor (one standard deviation) of the paracrystal is included
+in the calculation of $Z(q)$
+
+.. math::
+
+    \Delta a = g D
+
+where $g$ is a fractional distortion based on the nearest neighbor distance.
+
+
+.. figure:: img/bcc_geometry.jpg
+
+    Body-centered cubic lattice.
+
+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 body-centered cubic lattice, only reflections where
+$(h + k + l) = \text{even}$ are allowed and reflections where
+$(h + k + l) = \text{odd}$ are forbidden. Thus the peak positions
+correspond to (just the first 5)
+
+.. math::
+
+    \begin{array}{lccccc}
+    q/q_o          &   1   & \sqrt{2} & \sqrt{3} & \sqrt{4} & \sqrt{5} \\
+    \text{Indices} & (110) &    (200) & (211)    & (220)    & (310)    \\
+    \end{array}
+
+**NB**: 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,
+*qmin* = 0.001 |Ang^-1|, *qmax* = 0.1 |Ang^-1| and the above default values.
+
+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/bcc_angle_definition.png
+
+    Orientation of the crystal with respect to the scattering plane.
+
+References
+----------
+
+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 = "bcc_paracrystal"
+title = "Body-centred cubic lattic with paracrystalline distortion"
+description = """
+    Calculates the scattering from a **body-centered 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.
+    """
+category = "shape:paracrystal"
+
+single = False
+
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type","description" ],
+parameters = [["dnn",         "Ang",       220,    [-inf, inf], "",            "Nearest neighbour distance"],
+              ["d_factor",    "",            0.06, [-inf, inf], "",            "Paracrystal distortion factor"],
+              ["radius",      "Ang",        40,    [0, inf],    "volume",      "Particle radius"],
+              ["sld",         "1e-6/Ang^2",  4,    [-inf, inf], "",            "Particle scattering length density"],
+              ["sld_solvent", "1e-6/Ang^2",  1,    [-inf, inf], "",            "Solvent scattering length density"],
+              ["theta",       "degrees",    60,    [-inf, inf], "orientation", "In plane angle"],
+              ["phi",         "degrees",    60,    [-inf, inf], "orientation", "Out of plane angle"],
+              ["psi",         "degrees",    60,    [-inf, inf], "orientation", "Out of plane angle"]
+             ]
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["lib/sph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc_paracrystal_kernel.c"]
+
+# parameters for demo
+demo = dict(
+    scale=1, background=0,
+    dnn=220, d_factor=0.06, sld=4, sld_solvent=1,
+    radius=40,
+    theta=60, phi=60, psi=60,
+    radius_pd=.2, radius_pd_n=2,
+    theta_pd=15, theta_pd_n=0,
+    phi_pd=15, phi_pd_n=0,
+    psi_pd=15, psi_pd_n=0,
+    )
+
+# For testing against the old sasview models, include the converted parameter
+# names and the target sasview model name.
+oldname = 'BCCrystalModel'
+oldpars = dict(sld='sldSph', sld_solvent='sldSolv')
Index: sasmodels/models/bcc_paracrystal_kernel.c
===================================================================
--- sasmodels/models/bcc_paracrystal_kernel.c	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
+++ sasmodels/models/bcc_paracrystal_kernel.c	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
@@ -0,0 +1,170 @@
+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);
+
+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);
+
+
+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);
+}
+
+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 = pow(1.0-(temp10*temp10),3)*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){
+    return sphere_volume(radius);
+}
+
+
+double Iq(double q, double dnn,
+  double d_factor, double radius,
+  double sld, double solvent_sld){
+
+	//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);
+
+    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;
+
+
+}
+
+
+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 b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
+  //double q_z;
+  double cos_val_b3, cos_val_b2, cos_val_b1;
+  double a1_dot_q, a2_dot_q,a3_dot_q;
+  double answer;
+  double Zq, Fkq, Fkq_2;
+
+  //convert to q and make scaled values
+  double q = sqrt(qx*qx+qy*qy);
+  double q_x = qx/q;
+  double q_y = qy/q;
+
+  //convert angle degree to radian
+  theta = theta * M_PI_180;
+  phi = phi * M_PI_180;
+  psi = psi * M_PI_180;
+
+  const double Da = d_factor*dnn;
+  const double s1 = dnn/sqrt(0.75);
+
+
+  //the occupied volume of the lattice
+  const double latticescale = 2.0*sphere_volume(radius/s1);
+  // q vector
+  //q_z = 0.0; // for SANS; assuming qz is negligible
+  /// Angles here are respect to detector coordinate
+  ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
+    // b3 axis orientation
+    b3_x = cos(theta) * cos(phi);
+    b3_y = sin(theta);
+    //b3_z = -cos(theta) * sin(phi);
+    cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z;
+
+    //alpha = acos(cos_val_b3);
+    // b1 axis orientation
+    b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
+    b1_y = sin(psi)*cos(theta);
+    cos_val_b1 = b1_x*q_x + b1_y*q_y;
+    // b2 axis orientation
+    b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
+  	b2_y = cos(theta)*cos(psi);
+    cos_val_b2 = b2_x*q_x + b2_y*q_y;
+
+    // The following test should always pass
+    if (fabs(cos_val_b3)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b3 = 1.0;
+    }
+    if (fabs(cos_val_b2)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b2 = 1.0;
+    }
+    if (fabs(cos_val_b1)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b1 = 1.0;
+    }
+    // Compute the angle btw vector q and the a3 axis
+    a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3);
+
+    // a1 axis
+    a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1);
+
+    // a2 axis
+    a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2);
+
+
+    // Get Fkq and Fkq_2
+    Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
+    Fkq_2 = Fkq*Fkq;
+    // Call Zq=Z1*Z2*Z3
+    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2);
+
+  // Use SphereForm directly from libigor
+  answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale;
+
+  return answer;
+ }
Index: smodels/models/fcc.c
===================================================================
--- sasmodels/models/fcc.c	(revision ad90df9f3ed42019113506f1366fe4b851b5f389)
+++ 	(revision )
@@ -1,170 +1,0 @@
-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);
-
-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 _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);
-}
-
-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 = pow((1.0-(temp10*temp10)),3)*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){
-    return sphere_volume(radius);
-}
-
-
-double Iq(double q, double dnn,
-  double d_factor, double radius,
-  double sld, double solvent_sld){
-
-	//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);
-
-    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;
-
-
-}
-
-
-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 b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
-  // double q_z;
-  double cos_val_b3, cos_val_b2, cos_val_b1;
-  double a1_dot_q, a2_dot_q,a3_dot_q;
-  double answer;
-  double Zq, Fkq, Fkq_2;
-
-  //convert to q and make scaled values
-  double q = sqrt(qx*qx+qy*qy);
-  double q_x = qx/q;
-  double q_y = qy/q;
-
-  //convert angle degree to radian
-  theta = theta * M_PI_180;
-  phi = phi * M_PI_180;
-  psi = psi * M_PI_180;
-
-  const double Da = d_factor*dnn;
-  const double s1 = dnn/sqrt(0.75);
-
-
-  //the occupied volume of the lattice
-  const double latticescale = 2.0*sphere_volume(radius/s1);
-  // q vector
-  // q_z = 0.0; // for SANS; assuming qz is negligible
-  /// Angles here are respect to detector coordinate
-  ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
-    // b3 axis orientation
-    b3_x = cos(theta) * cos(phi);
-    b3_y = sin(theta);
-    //b3_z = -cos(theta) * sin(phi);
-    cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z;
-
-    //alpha = acos(cos_val_b3);
-    // b1 axis orientation
-    b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
-    b1_y = sin(psi)*cos(theta);
-    cos_val_b1 = b1_x*q_x + b1_y*q_y;
-    // b2 axis orientation
-    b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
-  	b2_y = cos(theta)*cos(psi);
-    cos_val_b2 = b2_x*q_x + b2_y*q_y;
-
-    // The following test should always pass
-    if (fabs(cos_val_b3)>1.0) {
-      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
-      cos_val_b3 = 1.0;
-    }
-    if (fabs(cos_val_b2)>1.0) {
-      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
-      cos_val_b2 = 1.0;
-    }
-    if (fabs(cos_val_b1)>1.0) {
-      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
-      cos_val_b1 = 1.0;
-    }
-    // Compute the angle btw vector q and the a3 axis
-    a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3);
-
-    // a1 axis
-    a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1);
-
-    // a2 axis
-    a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2);
-
-
-    // Get Fkq and Fkq_2
-    Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
-    Fkq_2 = Fkq*Fkq;
-    // Call Zq=Z1*Z2*Z3
-    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
-    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
-    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2);
-
-  // Use SphereForm directly from libigor
-  answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale;
-
-  return answer;
- }
Index: smodels/models/fcc.py
===================================================================
--- sasmodels/models/fcc.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
+++ 	(revision )
@@ -1,134 +1,0 @@
-#fcc paracrystal model
-#note model title and parameter table are automatically inserted
-#note - calculation requires double precision
-r"""
-Calculates the scattering from a **face-centered 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{\text{scale}}{V_p} V_\text{lattice} P(q) Z(q)
-
-where *scale* is the volume fraction of spheres, $V_p$ is the volume of
-the primary particle, $V_\text{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 face-centered cubic structure.
-
-Equation (1) of the 1990 reference is used to calculate $Z(q)$, using
-equations (23)-(25) from the 1987 paper for $Z1$, $Z2$, and $Z3$.
-
-The lattice correction (the occupied volume of the lattice) for a
-face-centered cubic structure of particles of radius $R$ and nearest
-neighbor separation $D$ is
-
-.. math::
-
-   V_\text{lattice} = \frac{16\pi}{3}\frac{R^3}{\left(D\sqrt{2}\right)^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.
-
-.. figure:: img/fcc_geometry.jpg
-
-    Face-centered cubic lattice.
-
-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 face-centered cubic lattice $h, k , l$ all odd or all
-even are allowed and reflections where $h, k, l$ are mixed odd/even
-are forbidden. Thus the peak positions correspond to (just the first 5)
-
-.. math::
-
-    \begin{array}{cccccc}
-    q/q_0 & 1 & \sqrt{4/3} & \sqrt{8/3} & \sqrt{11/3} & \sqrt{4} \\
-    \text{Indices} & (111)  & (200) & (220) & (311) & (222)
-    \end{array}
-
-**NB**: 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!
-
-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/bcc_angle_definition.png
-
-    Orientation of the crystal with respect to the scattering plane.
-
-References
-----------
-
-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 = "fcc_paracrystal"
-title = "Face-centred cubic lattic with paracrystalline distortion"
-description = """
-    Calculates the scattering from a **face-centered 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.
-    """
-category = "shape:paracrystal"
-
-single = False
-
-# pylint: disable=bad-whitespace, line-too-long
-#             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["dnn", "Ang", 220, [-inf, inf], "", "Nearest neighbour distance"],
-              ["d_factor", "", 0.06, [-inf, inf], "", "Paracrystal distortion factor"],
-              ["radius", "Ang", 40, [0, inf], "volume", "Particle radius"],
-              ["sld", "1e-6/Ang^2", 4, [-inf, inf], "", "Particle scattering length density"],
-              ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "", "Solvent scattering length density"],
-              ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"],
-              ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"],
-              ["psi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"]
-             ]
-# pylint: enable=bad-whitespace, line-too-long
-
-source = ["lib/sph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc.c"]
-
-# parameters for demo
-demo = dict(scale=1, background=0,
-            dnn=220, d_factor=0.06, sld=4, solvent_sld=1,
-            radius=40,
-            theta=60, phi=60, psi=60,
-            radius_pd=.2, radius_pd_n=0.2,
-            theta_pd=15, theta_pd_n=0,
-            phi_pd=15, phi_pd_n=0,
-            psi_pd=15, psi_pd_n=0,
-           )
-
-# For testing against the old sasview models, include the converted parameter
-# names and the target sasview model name.
-oldname = 'FCCrystalModel'
-oldpars = dict(sld='sldSph', solvent_sld='sldSolv')
Index: sasmodels/models/fcc_paracrystal.py
===================================================================
--- sasmodels/models/fcc_paracrystal.py	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
+++ sasmodels/models/fcc_paracrystal.py	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
@@ -0,0 +1,134 @@
+#fcc paracrystal model
+#note model title and parameter table are automatically inserted
+#note - calculation requires double precision
+r"""
+Calculates the scattering from a **face-centered 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{\text{scale}}{V_p} V_\text{lattice} P(q) Z(q)
+
+where *scale* is the volume fraction of spheres, $V_p$ is the volume of
+the primary particle, $V_\text{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 face-centered cubic structure.
+
+Equation (1) of the 1990 reference is used to calculate $Z(q)$, using
+equations (23)-(25) from the 1987 paper for $Z1$, $Z2$, and $Z3$.
+
+The lattice correction (the occupied volume of the lattice) for a
+face-centered cubic structure of particles of radius $R$ and nearest
+neighbor separation $D$ is
+
+.. math::
+
+   V_\text{lattice} = \frac{16\pi}{3}\frac{R^3}{\left(D\sqrt{2}\right)^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.
+
+.. figure:: img/fcc_geometry.jpg
+
+    Face-centered cubic lattice.
+
+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 face-centered cubic lattice $h, k , l$ all odd or all
+even are allowed and reflections where $h, k, l$ are mixed odd/even
+are forbidden. Thus the peak positions correspond to (just the first 5)
+
+.. math::
+
+    \begin{array}{cccccc}
+    q/q_0 & 1 & \sqrt{4/3} & \sqrt{8/3} & \sqrt{11/3} & \sqrt{4} \\
+    \text{Indices} & (111)  & (200) & (220) & (311) & (222)
+    \end{array}
+
+**NB**: 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!
+
+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/bcc_angle_definition.png
+
+    Orientation of the crystal with respect to the scattering plane.
+
+References
+----------
+
+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 = "fcc_paracrystal"
+title = "Face-centred cubic lattic with paracrystalline distortion"
+description = """
+    Calculates the scattering from a **face-centered 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.
+    """
+category = "shape:paracrystal"
+
+single = False
+
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["dnn", "Ang", 220, [-inf, inf], "", "Nearest neighbour distance"],
+              ["d_factor", "", 0.06, [-inf, inf], "", "Paracrystal distortion factor"],
+              ["radius", "Ang", 40, [0, inf], "volume", "Particle radius"],
+              ["sld", "1e-6/Ang^2", 4, [-inf, inf], "", "Particle scattering length density"],
+              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "", "Solvent scattering length density"],
+              ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"],
+              ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"],
+              ["psi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"]
+             ]
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["lib/sph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc_paracrystal_kernel.c"]
+
+# parameters for demo
+demo = dict(scale=1, background=0,
+            dnn=220, d_factor=0.06, sld=4, sld_solvent=1,
+            radius=40,
+            theta=60, phi=60, psi=60,
+            radius_pd=.2, radius_pd_n=0.2,
+            theta_pd=15, theta_pd_n=0,
+            phi_pd=15, phi_pd_n=0,
+            psi_pd=15, psi_pd_n=0,
+           )
+
+# For testing against the old sasview models, include the converted parameter
+# names and the target sasview model name.
+oldname = 'FCCrystalModel'
+oldpars = dict(sld='sldSph', sld_solvent='sldSolv')
Index: sasmodels/models/fcc_paracrystal_kernel.c
===================================================================
--- sasmodels/models/fcc_paracrystal_kernel.c	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
+++ sasmodels/models/fcc_paracrystal_kernel.c	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
@@ -0,0 +1,170 @@
+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);
+
+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 _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);
+}
+
+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 = pow((1.0-(temp10*temp10)),3)*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){
+    return sphere_volume(radius);
+}
+
+
+double Iq(double q, double dnn,
+  double d_factor, double radius,
+  double sld, double solvent_sld){
+
+	//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);
+
+    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;
+
+
+}
+
+
+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 b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
+  // double q_z;
+  double cos_val_b3, cos_val_b2, cos_val_b1;
+  double a1_dot_q, a2_dot_q,a3_dot_q;
+  double answer;
+  double Zq, Fkq, Fkq_2;
+
+  //convert to q and make scaled values
+  double q = sqrt(qx*qx+qy*qy);
+  double q_x = qx/q;
+  double q_y = qy/q;
+
+  //convert angle degree to radian
+  theta = theta * M_PI_180;
+  phi = phi * M_PI_180;
+  psi = psi * M_PI_180;
+
+  const double Da = d_factor*dnn;
+  const double s1 = dnn/sqrt(0.75);
+
+
+  //the occupied volume of the lattice
+  const double latticescale = 2.0*sphere_volume(radius/s1);
+  // q vector
+  // q_z = 0.0; // for SANS; assuming qz is negligible
+  /// Angles here are respect to detector coordinate
+  ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
+    // b3 axis orientation
+    b3_x = cos(theta) * cos(phi);
+    b3_y = sin(theta);
+    //b3_z = -cos(theta) * sin(phi);
+    cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z;
+
+    //alpha = acos(cos_val_b3);
+    // b1 axis orientation
+    b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
+    b1_y = sin(psi)*cos(theta);
+    cos_val_b1 = b1_x*q_x + b1_y*q_y;
+    // b2 axis orientation
+    b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
+  	b2_y = cos(theta)*cos(psi);
+    cos_val_b2 = b2_x*q_x + b2_y*q_y;
+
+    // The following test should always pass
+    if (fabs(cos_val_b3)>1.0) {
+      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b3 = 1.0;
+    }
+    if (fabs(cos_val_b2)>1.0) {
+      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b2 = 1.0;
+    }
+    if (fabs(cos_val_b1)>1.0) {
+      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b1 = 1.0;
+    }
+    // Compute the angle btw vector q and the a3 axis
+    a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3);
+
+    // a1 axis
+    a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1);
+
+    // a2 axis
+    a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2);
+
+
+    // Get Fkq and Fkq_2
+    Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
+    Fkq_2 = Fkq*Fkq;
+    // Call Zq=Z1*Z2*Z3
+    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2);
+
+  // Use SphereForm directly from libigor
+  answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale;
+
+  return answer;
+ }
Index: smodels/models/sc_crystal.c
===================================================================
--- sasmodels/models/sc_crystal.c	(revision ad90df9f3ed42019113506f1366fe4b851b5f389)
+++ 	(revision )
@@ -1,228 +1,0 @@
-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 sphere_volume(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
-	// 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;
-}
-
-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 = sphere_volume(radius/dnn);
-    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: smodels/models/sc_crystal.py
===================================================================
--- sasmodels/models/sc_crystal.py	(revision a84a0cac83cecece77ff644df4f2552ef1452b62)
+++ 	(revision )
@@ -1,159 +1,0 @@
-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_geometry.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::
-    :nowrap:
-
-    \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!
-
-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_angle_definition.jpg
-
-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/sph_j1c.c", "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.3048],
-    [{}, 0.215268, 0.00814889],
-    [{}, (0.414467), 0.001313289]
-    ]
-
-
Index: sasmodels/models/sc_paracrystal.py
===================================================================
--- sasmodels/models/sc_paracrystal.py	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
+++ sasmodels/models/sc_paracrystal.py	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
@@ -0,0 +1,159 @@
+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_geometry.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::
+    :nowrap:
+
+    \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!
+
+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_angle_definition.jpg
+
+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_paracrystal"
+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"],
+              ["sld",  "1e-6/Ang^2",  3.0,  [0.0, inf],  "",            "Sphere scattering length density"],
+              ["sld_solvent", "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/sph_j1c.c", "lib/sphere_form.c", "lib/gauss150.c", "sc_paracrystal_kernel.c"]
+
+demo = dict(scale=1, background=0,
+            dnn=220.0,
+            d_factor=0.06,
+            radius=40.0,
+            sld=3.0,
+            sld_solvent=6.3,
+            theta=0.0,
+            phi=0.0,
+            psi=0.0)
+
+oldname = 'SCCrystalModel'
+
+oldpars = dict(sld='sldSph',
+               sld_solvent='sldSolv')
+
+tests = [
+    # Accuracy tests based on content in test/utest_extra_models.py
+    [{}, 0.001, 10.3048],
+    [{}, 0.215268, 0.00814889],
+    [{}, (0.414467), 0.001313289]
+    ]
+
+
Index: sasmodels/models/sc_paracrystal_kernel.c
===================================================================
--- sasmodels/models/sc_paracrystal_kernel.c	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
+++ sasmodels/models/sc_paracrystal_kernel.c	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
@@ -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 sphere_volume(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
+	// 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;
+}
+
+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 = sphere_volume(radius/dnn);
+    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);
+
+}
+
