Index: explore/sc.c
===================================================================
--- explore/sc.c	(revision fdd56a1decdae24f11cdccb11ac7ffade9812587)
+++ explore/sc.c	(revision fdd56a1decdae24f11cdccb11ac7ffade9812587)
@@ -0,0 +1,137 @@
+static double
+_sq_sc(double qa, double qb, double qc, double dnn, double d_factor)
+{
+    // Rewriting equations for efficiency, accuracy and readability, and so
+    // code is reusable between 1D and 2D models.
+    const double a1 = qa;
+    const double a2 = qb;
+    const double a3 = qc;
+
+    const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+
+    // Numerator: (1 - exp(a)^2)^3
+    //         => (-(exp(2a) - 1))^3
+    //         => -expm1(2a)^3
+    // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2)
+    //         => exp(a)^2 - 2 cos(xk) exp(a) + 1
+    //         => (exp(a) - 2 cos(xk)) * exp(a) + 1
+    const double exp_arg = exp(-arg);
+    const double Sq = -cube(expm1(-2.0*arg))
+        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0));
+
+    return Sq;
+}
+
+// occupied volume fraction calculated from lattice symmetry and sphere radius
+static double
+_sc_volume_fraction(double radius, double dnn)
+{
+    return sphere_volume(radius/dnn);
+}
+
+static double
+form_volume(double radius)
+{
+    return sphere_volume(radius);
+}
+
+
+static double Iq(double q, double dnn,
+  double d_factor, double radius,
+  double sld, double solvent_sld,
+  double n, double sym)
+{
+double phi_m, phi_b, theta_m, theta_b;
+if (sym>0.) {
+    // translate a point in [-1,1] to a point in [0, 2 pi]
+    phi_m = M_PI_4;
+    phi_b = M_PI_4;
+    // translate a point in [-1,1] to a point in [0, pi]
+    theta_m = M_PI_4;
+    theta_b = M_PI_4;
+} else {
+    // translate a point in [-1,1] to a point in [0, 2 pi]
+    phi_m = M_PI;
+    phi_b = M_PI;
+    // translate a point in [-1,1] to a point in [0, pi]
+    theta_m = M_PI_2;
+    theta_b = M_PI_2;
+}
+
+#if 0
+    double outer_sum = 0.0;
+    for(int i=0; i<150; i++) {
+        double inner_sum = 0.0;
+        const double theta = Gauss150Z[i]*theta_m + theta_b;
+        double sin_theta, cos_theta;
+        SINCOS(theta, sin_theta, cos_theta);
+        const double qc = q*cos_theta;
+        const double qab = q*sin_theta;
+        for(int j=0;j<150;j++) {
+            const double phi = Gauss150Z[j]*phi_m + phi_b;
+            double sin_phi, cos_phi;
+            SINCOS(phi, sin_phi, cos_phi);
+            const double qa = qab*cos_phi;
+            const double qb = qab*sin_phi;
+            const double fq = _sq_sc(qa, qb, qc, dnn, d_factor);
+            inner_sum += Gauss150Wt[j] * fq;
+        }
+        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
+        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+    }
+    outer_sum *= theta_m;
+#else
+    double outer_sum = 0.0;
+    for(int i=0; i<(int)n; i++) {
+        double inner_sum = 0.0;
+        const double theta = (i*2./n-1.)*theta_m + theta_b;
+        double sin_theta, cos_theta;
+        SINCOS(theta, sin_theta, cos_theta);
+        const double qc = q*cos_theta;
+        const double qab = q*sin_theta;
+        for(int j=0;j<(int)n;j++) {
+            const double phi = (j*2./n-1.)*phi_m + phi_b;
+            double sin_phi, cos_phi;
+            SINCOS(phi, sin_phi, cos_phi);
+            const double qa = qab*cos_phi;
+            const double qb = qab*sin_phi;
+            const double fq = _sq_sc(qa, qb, qc, dnn, d_factor);
+            inner_sum += fq;
+        }
+        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
+        outer_sum += inner_sum * sin_theta;
+    }
+    outer_sum *= theta_m/(n*n);
+#endif
+double Sq;
+if (sym > 0.) {
+    Sq = outer_sum/M_PI_2;
+} else {
+    Sq = outer_sum/(4.0*M_PI);
+}
+
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+
+    return _sc_volume_fraction(radius, dnn) * Pq * Sq;
+}
+
+
+static double Iqxy(double qx, double qy,
+    double dnn, double d_factor, double radius,
+    double sld, double solvent_sld,
+    double n, double sym,
+    double theta, double phi, double psi)
+{
+    double q, zhat, yhat, xhat;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
+
+    q = sqrt(qa*qa + qb*qb + qc*qc);
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+    const double Sq = _sq_sc(qa, qb, qc, dnn, d_factor);
+    return _sc_volume_fraction(radius, dnn) * Pq * Sq;
+}
Index: explore/sc.py
===================================================================
--- explore/sc.py	(revision fdd56a1decdae24f11cdccb11ac7ffade9812587)
+++ explore/sc.py	(revision fdd56a1decdae24f11cdccb11ac7ffade9812587)
@@ -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) = \text{scale}\frac{V_\text{lattice}P(q)Z(q)}{V_p} + \text{background}
+
+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 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_\text{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/parallelepiped_angle_definition.png
+
+    Orientation of the crystal with respect to the scattering plane, when
+    $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis).
+
+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],  "sld",         "Sphere scattering length density"],
+              ["sld_solvent", "1e-6/Ang^2",  6.3, [0.0, inf],  "sld",         "Solvent scattering length density"],
+              ["n", "", 150, [0,inf], "", "integration order"],
+              ["sym", "", 1, [0,1], "", "use symmetry"],
+              ["theta",       "degrees",    0,    [-360, 360], "orientation", "c axis to beam angle"],
+              ["phi",         "degrees",    0,    [-360, 360], "orientation", "rotation about beam"],
+              ["psi",         "degrees",    0,    [-360, 360], "orientation", "rotation about c axis"]
+             ]
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c", "lib/gauss150.c", "sc.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)
+
+tests = [
+    # Accuracy tests based on content in test/utest_extra_models.py, 2d tests added April 10, 2017
+    [{}, 0.001, 10.3048],
+    [{}, 0.215268, 0.00814889],
+    [{}, 0.414467, 0.001313289],
+    [{'theta':10.0,'phi':20,'psi':30.0},(0.045,-0.035),18.0397138402 ],
+    [{'theta':10.0,'phi':20,'psi':30.0},(0.023,0.045),0.0177333171285 ]
+    ]
