Index: explore/angles.py
===================================================================
--- explore/angles.py	(revision e755d8a01d852e939d7f8af001c9f3a264e8086f)
+++ explore/angles.py	(revision e755d8a01d852e939d7f8af001c9f3a264e8086f)
@@ -0,0 +1,44 @@
+import sympy as sp
+dphi, dpsi, dtheta = sp.var("dphi dpsi dtheta")
+phi, psi, theta = sp.var("phi psi theta")
+x, y, z = sp.var("x y z")
+qx, qy, qz = sp.var("qx qy qz")
+
+def Rx(a):
+    return sp.Matrix([[1, 0, 0], [0, sp.cos(a), sp.sin(a)], [0, -sp.sin(a), sp.cos(a)]])
+def Ry(a):
+    return sp.Matrix([[sp.cos(a), 0, sp.sin(a)], [0, 1, 0], [-sp.sin(a), 0, sp.cos(a)]])
+def Rz(a):
+    return sp.Matrix([[sp.cos(a), sp.sin(a), 0], [-sp.sin(a), sp.cos(a), 0], [0, 0, 1]])
+
+print("==== asymmetric ====")
+q = sp.Matrix([[qx],[qy],[0]])
+p = sp.Matrix([[x],[y],[z]])
+jitter = Rx(dphi)*Ry(dtheta)*Rz(dpsi)
+view = Rz(phi)*Ry(theta)*Rz(psi)
+jitter_inv = Rz(-dpsi)*Ry(-dtheta)*Rx(-dphi)
+view_inv = Rz(-psi)*Ry(-theta)*Rz(-phi)
+print(">> jitter rotation")
+for row in jitter*p: print(row)
+print("\n>> jitter plus view")
+for row in view*jitter*p: print(row)
+print("\n>> view reverse")
+for row in view_inv*q: print(row)
+print("\n>> jitter reverse")
+for row in jitter_inv*view_inv*q: print(row)
+
+print("\n\n==== symmetric ====")
+q = sp.Matrix([[qx],[0],[0]])
+p = sp.Matrix([[x],[y],[z]])
+jitter = Rx(dphi)*Ry(dtheta)
+view = Rz(phi)*Ry(theta)
+jitter_inv = Ry(-dtheta)*Rx(-dphi)
+view_inv = Ry(-theta)*Rz(-phi)
+print(">> jitter rotation")
+for row in jitter*p: print(row)
+print("\n>> jitter plus view")
+for row in view*jitter*p: print(row)
+print("\n>> view reverse")
+for row in view_inv*q: print(row)
+print("\n>> jitter reverse")
+for row in jitter_inv*view_inv*q: print(row)
Index: explore/bccpy.py
===================================================================
--- explore/bccpy.py	(revision 7e0b281e5f175a91044364982fcfbfe802508948)
+++ explore/bccpy.py	(revision 7e0b281e5f175a91044364982fcfbfe802508948)
@@ -0,0 +1,217 @@
+"""
+The current 1D calculations for BCC paracrystal are very wrong at low q, orders
+of magnitude wrong.  The integration fails to capture a very narrow,
+very steep ridge.
+
+Uncomment the plot() line at the bottom of the code to show an image of the
+set of S(q, theta, phi) values that must be summed to form the 1D S(q=0.001)
+for the BCC paracrystal form.  Note particularly that this is a log scale
+image spanning 10 orders of magnitude.  This pattern repeats itself 8 times
+over the entire 4 pi surface integral.
+
+You can explore various integration options by uncommenting more lines.
+Adaptive integration using scpy.integrate.dbsquad is very slow.  Romberg
+didn't even complete in the time I gave it. Accurate brute force calculation
+requires a 4000x4000 grid to get enough precision.
+
+We may need a specialized integrator for low q which can identify and integrate
+the ridges properly.
+
+This program need sasmodels on the path so it is inserted automatically,
+assuming that the explore/bccpy.py is beside sasmodels/special.py in the
+source tree.  Run from the sasmodels directory using:
+
+    python explore/bccpy.py
+"""
+
+from __future__ import print_function, division
+
+import os, sys
+sys.path.insert(0, os.path.dirname(os.path.dirname(__file__)))
+
+import numpy as np
+from numpy import pi, sin, cos, exp, expm1, degrees, log10
+from scipy.integrate import dblquad, simps, romb, romberg
+import pylab
+
+from sasmodels.special import square
+from sasmodels.special import Gauss20Wt, Gauss20Z
+from sasmodels.special import Gauss76Wt, Gauss76Z
+from sasmodels.special import Gauss150Wt, Gauss150Z
+
+Q = 0.001
+DNN = 220.
+D_FACTOR = 0.06
+RADIUS = 40.0
+SLD = 3.0
+SLD_SOLVENT = 6.3
+
+# Note: using Matsuoka 1990; this is different from what
+# is in the sasmodels/models code  (see bcc vs bcc_old).
+# The difference is that the sign of phi and theta seem to be
+# negative in the old vs. the new, yielding a pattern that is
+# swapped left to right and top to bottom.
+def sc(qa, qb, qc):
+    return qa, qb, qc
+
+def bcc(qa, qb, qc):
+    a1 = (+qa + qb + qc)/2
+    a2 = (-qa - qb + qc)/2
+    a3 = (-qa + qb - qc)/2
+    return a1, a2, a3
+
+def bcc_old(qa, qb, qc):
+    a1 = (+qa + qb - qc)/2.0
+    a2 = (+qa - qb + qc)/2.0
+    a3 = (-qa + qb + qc)/2.0
+    return a1, a2, a3
+
+def fcc(qa, qb, qc):
+    a1 = ( 0. + qb + qc)/2
+    a2 = (-qa + 0. + qc)/2
+    a3 = (-qa + qb + 0.)/2
+    return a1, a2, a3
+
+def fcc_old(qa, qb, qc):
+    a1 = ( qa + qb + 0.)/2
+    a2 = ( qa + 0. + qc)/2
+    a3 = ( 0. + qb + qc)/2
+    return a1, a2, a3
+
+KERNEL = bcc
+
+def kernel(q, dnn, d_factor, theta, phi):
+    """
+    S(q) kernel for paracrystal forms.
+    """
+    qab = q*sin(theta)
+    qa = qab*cos(phi)
+    qb = qab*sin(phi)
+    qc = q*cos(theta)
+
+    a1, a2, a3 = KERNEL(qa, qb, qc)
+
+    # Note: paper says that different directions can have different distortion
+    # factors.  Easy enough to add to the code.  This would definitely break
+    # 8-fold symmetry.
+    arg = -0.5*square(dnn*d_factor)*(a1**2 + a2**2 + a3**2)
+    exp_arg = exp(arg)
+    den = [((exp_arg - 2*cos(dnn*a))*exp_arg + 1.0) for a in (a1, a2, a3)]
+    Zq = -expm1(2*arg)**3/np.prod(den, axis=0)
+    return Zq
+
+
+def scipy_dblquad(q=Q, dnn=DNN, d_factor=D_FACTOR):
+    """
+    Compute the integral using scipy dblquad.  This gets the correct answer
+    eventually, but it is slow.
+    """
+    evals = [0]
+    def integrand(theta, phi):
+        evals[0] += 1
+        Zq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=theta, phi=phi)
+        return Zq*sin(theta)
+    ans = dblquad(integrand, 0, pi/2, lambda x: 0, lambda x: pi/2)[0]*8/(4*pi)
+    print("dblquad evals =", evals[0])
+    return ans
+
+
+def scipy_romberg_2d(q=Q, dnn=DNN, d_factor=D_FACTOR):
+    """
+    Compute the integral using romberg integration.  This function does not
+    complete in a reasonable time.  No idea if it is accurate.
+    """
+    def inner(phi, theta):
+        return kernel(q=q, dnn=dnn, d_factor=d_factor, theta=theta, phi=phi)
+    def outer(theta):
+        return romberg(inner, 0, pi/2, divmax=100, args=(theta,))*sin(theta)
+    return romberg(outer, 0, pi/2, divmax=100)*8/(4*pi)
+
+
+def semi_romberg(q=Q, dnn=DNN, d_factor=D_FACTOR, n=100):
+    """
+    Use 1D romberg integration in phi and regular simpsons rule in theta.
+    """
+    evals = [0]
+    def inner(phi, theta):
+        evals[0] += 1
+        return kernel(q=q, dnn=dnn, d_factor=d_factor, theta=theta, phi=phi)
+    theta = np.linspace(0, pi/2, n)
+    f_phi = [romberg(inner, 0, pi/2, divmax=100, args=(t,))
+             for t in theta]
+    ans = simps(sin(theta)*np.array(f_phi), dx=theta[1]-theta[0])
+    print("semi romberg evals =", evals[0])
+    return ans*8/(4*pi)
+
+def gauss_quad(q=Q, dnn=DNN, d_factor=D_FACTOR, n=150):
+    """
+    Compute the integral using gaussian quadrature for n = 20, 76 or 150.
+    """
+    if n == 20:
+        z, w = Gauss20Z, Gauss20Wt
+    elif n == 76:
+        z, w = Gauss76Z, Gauss76Wt
+    else:
+        z, w = Gauss150Z, Gauss150Wt
+    theta = pi/4*(z + 1)
+    phi = pi/4*(z + 1)
+    Atheta, Aphi = np.meshgrid(theta, phi)
+    Aw = w[None, :] * w[:, None]
+    sin_theta = np.fmax(abs(sin(Atheta)), 1e-6)
+    Zq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi)
+    print("gauss %d evals ="%n, n**2)
+    return np.sum(Zq*Aw*sin_theta)*8/(4*pi)
+
+
+def gridded_integrals(q=0.001, dnn=DNN, d_factor=D_FACTOR, n=300):
+    """
+    Compute the integral on a regular grid using rectangular, trapezoidal,
+    simpsons, and romberg integration.  Romberg integration requires that
+    the grid be of size n = 2**k + 1.
+    """
+    theta = np.linspace(0, pi/2, n)
+    phi = np.linspace(0, pi/2, n)
+    Atheta, Aphi = np.meshgrid(theta, phi)
+    Zq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi)
+    Zq *= abs(sin(Atheta))
+    dx, dy = theta[1]-theta[0], phi[1]-phi[0]
+    print("rect", n, np.sum(Zq)*dx*dy*8/(4*pi))
+    print("trapz", n, np.trapz(np.trapz(Zq, dx=dx), dx=dy)*8/(4*pi))
+    print("simpson", n, simps(simps(Zq, dx=dx), dx=dy)*8/(4*pi))
+    print("romb", n, romb(romb(Zq, dx=dx), dx=dy)*8/(4*pi))
+    print("gridded %d evals ="%n, n**2)
+
+def plot(q=0.001, dnn=DNN, d_factor=D_FACTOR, n=300):
+    """
+    Plot the 2D surface that needs to be integrated in order to compute
+    the BCC S(q) at a particular q, dnn and d_factor.  *n* is the number
+    of points in the grid.
+    """
+    theta = np.linspace(0, pi, n)
+    phi = np.linspace(0, 2*pi, n)
+    #theta = np.linspace(0, pi/2, n)
+    #phi = np.linspace(0, pi/2, n)
+    Atheta, Aphi = np.meshgrid(theta, phi)
+    Zq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi)
+    Zq *= abs(sin(Atheta))
+    pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Zq, 1.e-6)))
+    pylab.axis('tight')
+    pylab.title("%s Z(q) for q=%g, dnn=%g d_factor=%g"
+                % (KERNEL.__name__, q, dnn, d_factor))
+    pylab.xlabel("theta (degrees)")
+    pylab.ylabel("phi (degrees)")
+    cbar = pylab.colorbar()
+    cbar.set_label('log10 S(q)')
+    pylab.show()
+
+if __name__ == "__main__":
+    #print("gauss", 20, gauss_quad(n=20))
+    #print("gauss", 76, gauss_quad(n=76))
+    #print("gauss", 150, gauss_quad(n=150))
+    #print("dblquad", scipy_dblquad())
+    #print("semi romberg", semi_romberg())
+    #gridded_integrals(n=2**8+1)
+    #gridded_integrals(n=2**10+1)
+    #gridded_integrals(n=2**13+1)
+    #print("romberg", scipy_romberg())
+    plot(n=400)
Index: explore/jitter.py
===================================================================
--- explore/jitter.py	(revision 85190c2070ee71c354d018a96403e7d1e6045492)
+++ explore/jitter.py	(revision d4c33d66aa792e691ad5eb91b36d3e854e88f383)
@@ -3,13 +3,14 @@
 dispersity and possible replacement algorithms.
 """
+import sys
+
 import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot
 import matplotlib.pyplot as plt
 from matplotlib.widgets import Slider, CheckButtons
 from matplotlib import cm
-
 import numpy as np
 from numpy import pi, cos, sin, sqrt, exp, degrees, radians
 
-def draw_beam(ax):
+def draw_beam(ax, view=(0, 0)):
     #ax.plot([0,0],[0,0],[1,-1])
     #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8)
@@ -22,14 +23,18 @@
     x = r*np.outer(np.cos(u), np.ones_like(v))
     y = r*np.outer(np.sin(u), np.ones_like(v))
-    z = np.outer(np.ones_like(u), v)
+    z = 1.3*np.outer(np.ones_like(u), v)
+
+    theta, phi = view
+    shape = x.shape
+    points = np.matrix([x.flatten(), y.flatten(), z.flatten()])
+    points = Rz(phi)*Ry(theta)*points
+    x, y, z = [v.reshape(shape) for v in points]
 
     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5)
 
-def draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi):
-    size=[0.1, 0.4, 1.0]
-    view=[theta, phi, psi]
-    shimmy=[0,0,0]
-    #draw_shape = draw_parallelepiped
-    draw_shape = draw_ellipsoid
+def draw_jitter(ax, view, jitter):
+    size = [0.1, 0.4, 1.0]
+    draw_shape = draw_parallelepiped
+    #draw_shape = draw_ellipsoid
 
     #np.random.seed(10)
@@ -64,4 +69,5 @@
         [ 1,  1,  1],
     ]
+    dtheta, dphi, dpsi = jitter
     if dtheta == 0:
         cloud = [v for v in cloud if v[0] == 0]
@@ -70,8 +76,8 @@
     if dpsi == 0:
         cloud = [v for v in cloud if v[2] == 0]
-    draw_shape(ax, size, view, shimmy, steps=100, alpha=0.8)
+    draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8)
     for point in cloud:
-        shimmy=[dtheta*point[0], dphi*point[1], dpsi*point[2]]
-        draw_shape(ax, size, view, shimmy, alpha=0.8)
+        delta = [dtheta*point[0], dphi*point[1], dpsi*point[2]]
+        draw_shape(ax, size, view, delta, alpha=0.8)
     for v in 'xyz':
         a, b, c = size
@@ -80,9 +86,6 @@
         getattr(ax, v+'axis').label.set_text(v)
 
-def draw_ellipsoid(ax, size, view, shimmy, steps=25, alpha=1):
+def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1):
     a,b,c = size
-    theta, phi, psi = view
-    dtheta, dphi, dpsi = shimmy
-
     u = np.linspace(0, 2 * np.pi, steps)
     v = np.linspace(0, np.pi, steps)
@@ -90,18 +93,19 @@
     y = b*np.outer(np.sin(u), np.sin(v))
     z = c*np.outer(np.ones_like(u), np.cos(v))
-
-    shape = x.shape
-    points = np.matrix([x.flatten(),y.flatten(),z.flatten()])
-    points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points
-    points = Rz(phi)*Ry(theta)*Rz(psi)*points
-    x,y,z = [v.reshape(shape) for v in points]
+    x, y, z = transform_xyz(view, jitter, x, y, z)
 
     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha)
 
-def draw_parallelepiped(ax, size, view, shimmy, alpha=1):
+    draw_labels(ax, view, jitter, [
+         ('c+', [ 0, 0, c], [ 1, 0, 0]),
+         ('c-', [ 0, 0,-c], [ 0, 0,-1]),
+         ('a+', [ a, 0, 0], [ 0, 0, 1]),
+         ('a-', [-a, 0, 0], [ 0, 0,-1]),
+         ('b+', [ 0, b, 0], [-1, 0, 0]),
+         ('b-', [ 0,-b, 0], [-1, 0, 0]),
+    ])
+
+def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1):
     a,b,c = size
-    theta, phi, psi = view
-    dtheta, dphi, dpsi = shimmy
-
     x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1])
     y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1])
@@ -118,53 +122,40 @@
     ])
 
-    points = np.matrix([x,y,z])
-    points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points
-    points = Rz(phi)*Ry(theta)*Rz(psi)*points
-
-    x,y,z = [np.array(v).flatten() for v in points]
+    x, y, z = transform_xyz(view, jitter, x, y, z)
     ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha)
 
-def draw_sphere(ax, radius=10., steps=100):
-    u = np.linspace(0, 2 * np.pi, steps)
-    v = np.linspace(0, np.pi, steps)
-
-    x = radius * np.outer(np.cos(u), np.sin(v))
-    y = radius * np.outer(np.sin(u), np.sin(v))
-    z = radius * np.outer(np.ones(np.size(u)), np.cos(v))
-    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w')
-
-def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'):
-    theta_center = radians(theta)
-    phi_center = radians(phi)
-    flow_center = radians(flow)
-    dtheta = radians(dtheta)
-    dphi = radians(dphi)
-
-    # 10 point 3-sigma gaussian weights
-    t = np.linspace(-3., 3., 11)
+    draw_labels(ax, view, jitter, [
+         ('c+', [ 0, 0, c], [ 1, 0, 0]),
+         ('c-', [ 0, 0,-c], [ 0, 0,-1]),
+         ('a+', [ a, 0, 0], [ 0, 0, 1]),
+         ('a-', [-a, 0, 0], [ 0, 0,-1]),
+         ('b+', [ 0, b, 0], [-1, 0, 0]),
+         ('b-', [ 0,-b, 0], [-1, 0, 0]),
+    ])
+
+def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gauss'):
+    theta, phi, psi = view
+    dtheta, dphi, dpsi = jitter
     if dist == 'gauss':
+        t = np.linspace(-3, 3, n)
         weights = exp(-0.5*t**2)
     elif dist == 'rect':
+        t = np.linspace(0, 1, n)
         weights = np.ones_like(t)
     else:
         raise ValueError("expected dist to be 'gauss' or 'rect'")
-    theta = dtheta*t
-    phi = dphi*t
-
-    x = radius * np.outer(cos(phi), cos(theta))
-    y = radius * np.outer(sin(phi), cos(theta))
-    z = radius * np.outer(np.ones_like(phi), sin(theta))
-    #w = np.outer(weights, weights*abs(cos(dtheta*t)))
-    w = np.outer(weights, weights*abs(cos(theta)))
-
-    x, y, z, w = [v.flatten() for v in (x,y,z,w)]
-    x, y, z = rotate(x, y, z, phi_center, theta_center, flow_center)
-
-    ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.)
-
-def rotate(x, y, z, phi, theta, psi):
-    R = Rz(phi)*Ry(theta)*Rz(psi)
-    p = np.vstack([x,y,z])
-    return R*p
+
+    # mesh in theta, phi formed by rotating z
+    z = np.matrix([[0], [0], [radius]])
+    points = np.hstack([Rx(phi_i)*Ry(theta_i)*z
+                        for theta_i in dtheta*t
+                        for phi_i in dphi*t])
+    # rotate relative to beam
+    points = orient_relative_to_beam(view, points)
+
+    w = np.outer(weights, weights)
+
+    x, y, z = [np.array(v).flatten() for v in points]
+    ax.scatter(x, y, z, c=w.flatten(), marker='o', vmin=0., vmax=1.)
 
 def Rx(angle):
@@ -188,4 +179,44 @@
          [0., 0., 1.]]
     return np.matrix(R)
+
+def transform_xyz(view, jitter, x, y, z):
+    x, y, z = [np.asarray(v) for v in (x, y, z)]
+    shape = x.shape
+    points = np.matrix([x.flatten(),y.flatten(),z.flatten()])
+    points = apply_jitter(jitter, points)
+    points = orient_relative_to_beam(view, points)
+    x, y, z = [np.array(v).reshape(shape) for v in points]
+    return x, y, z
+
+def apply_jitter(jitter, points):
+    dtheta, dphi, dpsi = jitter
+    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points
+    return points
+
+def orient_relative_to_beam(view, points):
+    theta, phi, psi = view
+    points = Rz(phi)*Ry(theta)*Rz(psi)*points
+    return points
+
+def draw_labels(ax, view, jitter, text):
+    labels, locations, orientations = zip(*text)
+    px, py, pz = zip(*locations)
+    dx, dy, dz = zip(*orientations)
+
+    px, py, pz = transform_xyz(view, jitter, px, py, pz)
+    dx, dy, dz = transform_xyz(view, jitter, dx, dy, dz)
+
+    for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)):
+        zdir = np.asarray(zdir).flatten()
+        ax.text(p[0], p[1], p[2], label, zdir=zdir)
+
+def draw_sphere(ax, radius=10., steps=100):
+    u = np.linspace(0, 2 * np.pi, steps)
+    v = np.linspace(0, np.pi, steps)
+
+    x = radius * np.outer(np.cos(u), np.sin(v))
+    y = radius * np.outer(np.sin(u), np.sin(v))
+    z = radius * np.outer(np.ones(np.size(u)), np.cos(v))
+    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w')
 
 def main():
@@ -206,4 +237,5 @@
 
     axcolor = 'lightgoldenrodyellow'
+
     axtheta  = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor)
     axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor)
@@ -212,4 +244,5 @@
     sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi)
     spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi)
+
     axdtheta  = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor)
     axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor)
@@ -217,14 +250,14 @@
     sdtheta = Slider(axdtheta, 'dTheta', 0, 30, valinit=dtheta)
     sdphi = Slider(axdphi, 'dPhi', 0, 30, valinit=dphi)
-    sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dphi)
+    sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dpsi)
 
     def update(val, axis=None):
-        theta, phi, psi = stheta.val, sphi.val, spsi.val
-        dtheta, dphi, dpsi = sdtheta.val, sdphi.val, sdpsi.val
+        view = stheta.val, sphi.val, spsi.val
+        jitter = sdtheta.val, sdphi.val, sdpsi.val
         ax.cla()
-        draw_beam(ax)
-        draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi)
-        #if not axis.startswith('d'):
-        #    ax.view_init(elev=theta, azim=phi)
+        draw_beam(ax, (0, 0))
+        draw_jitter(ax, view, jitter)
+        #draw_jitter(ax, view, (0,0,0))
+        draw_mesh(ax, view, jitter)
         plt.gcf().canvas.draw()
 
Index: explore/sc.c
===================================================================
--- explore/sc.c	(revision 64ca16350cce3cfe6ae7ba9b085022c0f44cc6e0)
+++ explore/sc.c	(revision 64ca16350cce3cfe6ae7ba9b085022c0f44cc6e0)
@@ -0,0 +1,148 @@
+static double
+sc_Zq(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.
+    #if 1  // SC
+    const double a1 = qa;
+    const double a2 = qb;
+    const double a3 = qc;
+    #elif 1 // BCC
+    const double a1 = (+qa + qb + qc)/2.;
+    const double a2 = (-qa - qb + qc)/2.;
+    const double a3 = (-qa + qb - qc)/2.;
+    #elif 1 // FCC
+    const double a1 = ( qa + qb)/2.0;
+    const double a2 = (-qa + qc)/2.0;
+    const double a3 = (-qa + qb)/2.0;
+    #endif
+
+    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 Zq = -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 Zq;
+}
+
+// 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 form = sc_Zq(qa, qb, qc, dnn, d_factor);
+            inner_sum += form;
+        }
+        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 Zq;
+if (sym > 0.) {
+    Zq = outer_sum/M_PI_2;
+} else {
+    Zq = outer_sum/(4.0*M_PI);
+}
+
+    //return Zq;
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+    return sc_volume_fraction(radius, dnn) * Pq * Zq;
+}
+
+
+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 Zq = sc_Zq(qa, qb, qc, dnn, d_factor);
+    return Zq;
+    return sc_volume_fraction(radius, dnn) * Pq * Zq;
+}
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 ]
+    ]
