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 ]
+    ]
Index: sasmodels/details.py
===================================================================
--- sasmodels/details.py	(revision ccd5f018ad394bf78c2726b92d428f6903233b21)
+++ sasmodels/details.py	(revision f39759c773168cc38c66528d0c9c86a7957986fc)
@@ -15,5 +15,5 @@
 
 import numpy as np  # type: ignore
-from numpy import pi, cos, sin
+from numpy import pi, cos, sin, radians
 
 try:
@@ -29,9 +29,10 @@
 
 try:
-    from typing import List
+    from typing import List, Tuple, Sequence
 except ImportError:
     pass
 else:
     from .modelinfo import ModelInfo
+    from .modelinfo import ParameterTable
 
 
@@ -53,5 +54,5 @@
     coordinates, the total circumference decreases as latitude varies from
     pi r^2 at the equator to 0 at the pole, and the weight associated
-    with a range of phi values needs to be scaled by this circumference.
+    with a range of latitude values needs to be scaled by this circumference.
     This scale factor needs to be updated each time the theta value
     changes.  *theta_par* indicates which of the values in the parameter
@@ -231,5 +232,7 @@
     nvalues = kernel.info.parameters.nvalues
     scalars = [(v[0] if len(v) else np.NaN) for v, w in pairs]
-    values, weights = zip(*pairs[2:npars+2]) if npars else ((),())
+    # skipping scale and background when building values and weights
+    values, weights = zip(*pairs[2:npars+2]) if npars else ((), ())
+    weights = correct_theta_weights(kernel.info.parameters, values, weights)
     length = np.array([len(w) for w in weights])
     offset = np.cumsum(np.hstack((0, length)))
@@ -244,6 +247,32 @@
     return call_details, data, is_magnetic
 
+def correct_theta_weights(parameters, values, weights):
+    # type: (ParameterTable, Sequence[np.ndarray], Sequence[np.ndarray]) -> Sequence[np.ndarray]
+    """
+    If there is a theta parameter, update the weights of that parameter so that
+    the cosine weighting required for polar integration is preserved.  Avoid
+    evaluation strictly at the pole, which would otherwise send the weight to
+    zero.
+
+    Note: values and weights do not include scale and background
+    """
+    # TODO: document code, explaining why scale and background are skipped
+    # given that we don't have scale and background in the list, we
+    # should be calling the variables something other than values and weights
+    # Apparently the parameters.theta_offset similarly skips scale and
+    # and background, so the indexing works out.
+    if parameters.theta_offset >= 0:
+        index = parameters.theta_offset
+        theta = values[index]
+        theta_weight = np.minimum(abs(cos(radians(theta))), 1e-6)
+        # copy the weights list so we can update it
+        weights = list(weights)
+        weights[index] = theta_weight*np.asarray(weights[index])
+        weights = tuple(weights)
+    return weights
+
 
 def convert_magnetism(parameters, values):
+    # type: (ParameterTable, Sequence[np.ndarray]) -> bool
     """
     Convert magnetism values from polar to rectangular coordinates.
@@ -255,5 +284,5 @@
     scale = mag[:,0]
     if np.any(scale):
-        theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180.
+        theta, phi = radians(mag[:, 1]), radians(mag[:, 2])
         cos_theta = cos(theta)
         mag[:, 0] = scale*cos_theta*cos(phi)  # mx
@@ -269,4 +298,12 @@
     """
     Create a mesh grid of dispersion parameters and weights.
+
+    pars is a list of pairs (values, weights), where the values are the
+    individual parameter values at which to evaluate the polydispersity
+    distribution and weights are the weights associated with each value.
+
+    Only the volume parameters should be included in this list.  Orientation
+    parameters do not affect the calculation of effective radius or volume
+    ratio.
 
     Returns [p1,p2,...],w where pj is a vector of values for parameter j
Index: sasmodels/kernel_header.c
===================================================================
--- sasmodels/kernel_header.c	(revision bb4b50990533a15ceda4443988db938e13a28a2b)
+++ sasmodels/kernel_header.c	(revision bb4b50990533a15ceda4443988db938e13a28a2b)
@@ -164,5 +164,5 @@
     SINCOS(phi*M_PI_180, sn, cn); \
     q = sqrt(qx*qx + qy*qy); \
-    cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \
+    cn = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \
     sn = sqrt(1 - cn*cn); \
     } while (0)
Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision bde38b559fcdc618ebafdbb69cf7af6120ea8466)
+++ sasmodels/kernel_iq.c	(revision d4c33d66aa792e691ad5eb91b36d3e854e88f383)
@@ -25,5 +25,5 @@
     int32_t num_weights;        // total length of the weights vector
     int32_t num_active;         // number of non-trivial pd loops
-    int32_t theta_par;          // id of spherical correction variable
+    int32_t theta_par;          // id of spherical correction variable (not used)
 } ProblemDetails;
 
@@ -173,15 +173,4 @@
 
 
-#if MAX_PD>0
-  const int theta_par = details->theta_par;
-  const int fast_theta = (theta_par == p0);
-  const int slow_theta = (theta_par >= 0 && !fast_theta);
-  double spherical_correction = 1.0;
-#else
-  // Note: if not polydisperse the weights cancel and we don't need the
-  // spherical correction.
-  const double spherical_correction = 1.0;
-#endif
-
   int step = pd_start;
 
@@ -220,14 +209,8 @@
 #endif
 #if MAX_PD>0
-  if (slow_theta) { // Theta is not in inner loop
-    spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6);
-  }
   while(i0 < n0) {
     local_values.vector[p0] = v0[i0];
     double weight0 = w0[i0] * weight1;
 //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0);
-    if (fast_theta) { // Theta is in inner loop
-      spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6);
-    }
 #else
     const double weight0 = 1.0;
@@ -244,8 +227,5 @@
       // Note: weight==0 must always be excluded
       if (weight0 > cutoff) {
-        // spherical correction is set at a minimum of 1e-6, otherwise there
-        // would be problems looking at models with theta=90.
-        const double weight = weight0 * spherical_correction;
-        pd_norm += weight * CALL_VOLUME(local_values.table);
+        pd_norm += weight0 * CALL_VOLUME(local_values.table);
 
         #ifdef USE_OPENMP
@@ -304,5 +284,5 @@
 #endif // !MAGNETIC
 //printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0);
-          result[q_index] += weight * scattering;
+          result[q_index] += weight0 * scattering;
         }
       }
Index: sasmodels/kernel_iq.cl
===================================================================
--- sasmodels/kernel_iq.cl	(revision bde38b559fcdc618ebafdbb69cf7af6120ea8466)
+++ sasmodels/kernel_iq.cl	(revision d4c33d66aa792e691ad5eb91b36d3e854e88f383)
@@ -25,5 +25,5 @@
     int32_t num_weights;        // total length of the weights vector
     int32_t num_active;         // number of non-trivial pd loops
-    int32_t theta_par;          // id of spherical correction variable
+    int32_t theta_par;          // id of spherical correction variable (not used)
 } ProblemDetails;
 
@@ -169,15 +169,4 @@
 
 
-#if MAX_PD>0
-  const int theta_par = details->theta_par;
-  const bool fast_theta = (theta_par == p0);
-  const bool slow_theta = (theta_par >= 0 && !fast_theta);
-  double spherical_correction = 1.0;
-#else
-  // Note: if not polydisperse the weights cancel and we don't need the
-  // spherical correction.
-  const double spherical_correction = 1.0;
-#endif
-
   int step = pd_start;
 
@@ -217,14 +206,8 @@
 #endif
 #if MAX_PD>0
-  if (slow_theta) { // Theta is not in inner loop
-    spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6);
-  }
   while(i0 < n0) {
     local_values.vector[p0] = v0[i0];
     double weight0 = w0[i0] * weight1;
 //if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0);
-    if (fast_theta) { // Theta is in inner loop
-      spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6);
-    }
 #else
     const double weight0 = 1.0;
@@ -232,5 +215,4 @@
 
 //if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, local_values.vector[i]); printf("\n"); }
-//if (q_index == 0) printf("sphcor: %g\n", spherical_correction);
 
     #ifdef INVALID
@@ -241,8 +223,5 @@
       // Note: weight==0 must always be excluded
       if (weight0 > cutoff) {
-        // spherical correction is set at a minimum of 1e-6, otherwise there
-        // would be problems looking at models with theta=90.
-        const double weight = weight0 * spherical_correction;
-        pd_norm += weight * CALL_VOLUME(local_values.table);
+        pd_norm += weight0 * CALL_VOLUME(local_values.table);
 
 #if defined(MAGNETIC) && NUM_MAGNETIC > 0
@@ -296,5 +275,5 @@
         const double scattering = CALL_IQ(q, q_index, local_values.table);
 #endif // !MAGNETIC
-        this_result += weight * scattering;
+        this_result += weight0 * scattering;
       }
     }
Index: sasmodels/kernelpy.py
===================================================================
--- sasmodels/kernelpy.py	(revision dbfd47189324c900d6e5cf98a210576a86204ced)
+++ sasmodels/kernelpy.py	(revision d4c33d66aa792e691ad5eb91b36d3e854e88f383)
@@ -197,10 +197,8 @@
 
     pd_norm = 0.0
-    spherical_correction = 1.0
     partial_weight = np.NaN
     weight = np.NaN
 
     p0_par = call_details.pd_par[0]
-    p0_is_theta = (p0_par == call_details.theta_par)
     p0_length = call_details.pd_length[0]
     p0_index = p0_length
@@ -219,14 +217,8 @@
             parameters[pd_par] = pd_value[pd_offset+pd_index]
             partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:])
-            if call_details.theta_par >= 0:
-                cor = sin(pi / 180 * parameters[call_details.theta_par])
-                spherical_correction = max(abs(cor), 1e-6)
             p0_index = loop_index%p0_length
 
         weight = partial_weight * pd_weight[p0_offset + p0_index]
         parameters[p0_par] = pd_value[p0_offset + p0_index]
-        if p0_is_theta:
-            cor = cos(pi/180 * parameters[p0_par])
-            spherical_correction = max(abs(cor), 1e-6)
         p0_index += 1
         if weight > cutoff:
@@ -239,5 +231,4 @@
 
             # update value and norm
-            weight *= spherical_correction
             total += weight * Iq
             pd_norm += weight * form_volume()
Index: sasmodels/model_test.py
===================================================================
--- sasmodels/model_test.py	(revision bb4b50990533a15ceda4443988db938e13a28a2b)
+++ sasmodels/model_test.py	(revision bb4b50990533a15ceda4443988db938e13a28a2b)
@@ -79,7 +79,7 @@
     suite = unittest.TestSuite()
 
-    if models[0] == 'all':
+    if models[0] in core.KINDS:
         skip = models[1:]
-        models = list_models()
+        models = list_models(models[0])
     else:
         skip = []
Index: sasmodels/models/_py_parallelepiped.py
===================================================================
--- sasmodels/models/_py_parallelepiped.py	(revision 597878a0f8b1b1b38cfce28323dfa620bdbaea50)
+++ sasmodels/models/_py_parallelepiped.py	(revision 597878a0f8b1b1b38cfce28323dfa620bdbaea50)
@@ -0,0 +1,316 @@
+# parallelepiped model
+# Note: model title and parameter table are inserted automatically
+r"""
+The form factor is normalized by the particle volume.
+For information about polarised and magnetic scattering, see
+the :ref:`magnetism` documentation.
+
+Definition
+----------
+
+ This model calculates the scattering from a rectangular parallelepiped
+ (\:numref:`parallelepiped-image`\).
+ If you need to apply polydispersity, see also :ref:`rectangular-prism`.
+
+.. _parallelepiped-image:
+
+
+.. figure:: img/parallelepiped_geometry.jpg
+
+   Parallelepiped with the corresponding definition of sides.
+
+.. note::
+
+   The edge of the solid used to have to satisfy the condition that $A < B < C$.
+   After some improvements to the effective radius calculation, used with
+   an S(Q), it is beleived that this is no longer the case.
+
+The 1D scattering intensity $I(q)$ is calculated as:
+
+.. Comment by Miguel Gonzalez:
+   I am modifying the original text because I find the notation a little bit
+   confusing. I think that in most textbooks/papers, the notation P(Q) is
+   used for the form factor (adim, P(Q=0)=1), although F(q) seems also to
+   be used. But here (as for many other models), P(q) is used to represent
+   the scattering intensity (in cm-1 normally). It would be good to agree on
+   a common notation.
+
+.. math::
+
+    I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2
+           \left< P(q, \alpha) \right> + \text{background}
+
+where the volume $V = A B C$, the contrast is defined as
+$\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$,
+$P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented
+at an angle $\alpha$ (angle between the long axis C and $\vec q$),
+and the averaging $\left<\ldots\right>$ is applied over all orientations.
+
+Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the
+form factor is given by (Mittelbach and Porod, 1961)
+
+.. math::
+
+    P(q, \alpha) = \int_0^1 \phi_Q\left(\mu \sqrt{1-\sigma^2},a\right)
+        \left[S(\mu c \sigma/2)\right]^2 d\sigma
+
+with
+
+.. math::
+
+    \phi_Q(\mu,a) &= \int_0^1
+        \left\{S\left[\frac{\mu}{2}\cos\left(\frac{\pi}{2}u\right)\right]
+               S\left[\frac{\mu a}{2}\sin\left(\frac{\pi}{2}u\right)\right]
+               \right\}^2 du
+
+    S(x) &= \frac{\sin x}{x}
+
+    \mu &= qB
+
+The scattering intensity per unit volume is returned in units of |cm^-1|.
+
+NB: The 2nd virial coefficient of the parallelepiped is calculated based on
+the averaged effective radius, after appropriately sorting the three
+dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and
+length $(= C)$ values, and used as the effective radius for
+$S(q)$ when $P(q) \cdot S(q)$ is applied.
+
+To provide easy access to the orientation of the parallelepiped, we define
+three angles $\theta$, $\phi$ and $\Psi$. The definition of $\theta$ and
+$\phi$ is the same as for the cylinder model (see also figures below).
+
+.. Comment by Miguel Gonzalez:
+   The following text has been commented because I think there are two
+   mistakes. Psi is the rotational angle around C (but I cannot understand
+   what it means against the q plane) and psi=0 corresponds to a||x and b||y.
+
+   The angle $\Psi$ is the rotational angle around the $C$ axis against
+   the $q$ plane. For example, $\Psi = 0$ when the $B$ axis is parallel
+   to the $x$-axis of the detector.
+
+The angle $\Psi$ is the rotational angle around the $C$ axis.
+For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the $B$ axis
+oriented parallel to the y-axis of the detector with $A$ along the z-axis.
+For other $\theta$, $\phi$ values, the parallelepiped has to be first rotated
+$\theta$ degrees around $z$ and $\phi$ degrees around $y$,
+before doing a final rotation of $\Psi$ degrees around the resulting $C$ to
+obtain the final orientation of the parallelepiped.
+For example, for $\theta = 0$ and $\phi = 90$, we have that $\Psi = 0$
+corresponds to $A$ along $x$ and $B$ along $y$,
+while for $\theta = 90$ and $\phi = 0$, $\Psi = 0$ corresponds to
+$A$ along $z$ and $B$ along $x$.
+
+.. _parallelepiped-orientation:
+
+.. figure:: img/parallelepiped_angle_definition.png
+
+    Definition of the angles for oriented parallelepiped, shown with $A<B<C$.
+
+.. figure:: img/parallelepiped_angle_projection.png
+
+    Examples of the angles for an oriented parallelepiped against the
+    detector plane.
+
+On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will
+appear. These are actually rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped, perpendicular to the $a$ x $c$ and $b$ x $c$ faces. 
+(When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.) The third orientation distribution, in $\psi$, is 
+about the $c$ axis of the particle, perpendicular to the $a$ x $b$ face. Some experimentation may be required to 
+understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation 
+distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.) 
+
+    
+For a given orientation of the parallelepiped, the 2D form factor is
+calculated as
+
+.. math::
+
+    P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}{2}qA\cos\alpha)}\right]^2
+                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}{2}qB\cos\beta)}\right]^2
+                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}{2}qC\cos\gamma)}\right]^2
+
+with
+
+.. math::
+
+    \cos\alpha &= \hat A \cdot \hat q,
+
+    \cos\beta  &= \hat B \cdot \hat q,
+
+    \cos\gamma &= \hat C \cdot \hat q
+
+and the scattering intensity as:
+
+.. math::
+
+    I(q_x, q_y) = \frac{\text{scale}}{V} V^2 \Delta\rho^2 P(q_x, q_y)
+            + \text{background}
+
+.. Comment by Miguel Gonzalez:
+   This reflects the logic of the code, as in parallelepiped.c the call
+   to _pkernel returns $P(q_x, q_y)$ and then this is multiplied by
+   $V^2 * (\Delta \rho)^2$. And finally outside parallelepiped.c it will be
+   multiplied by scale, normalized by $V$ and the background added. But
+   mathematically it makes more sense to write
+   $I(q_x, q_y) = \text{scale} V \Delta\rho^2 P(q_x, q_y) + \text{background}$,
+   with scale being the volume fraction.
+
+
+Validation
+----------
+
+Validation of the code was done by comparing the output of the 1D calculation
+to the angular average of the output of a 2D calculation over all possible
+angles.
+
+
+References
+----------
+
+P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
+
+R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
+
+Authorship and Verification
+----------------------------
+
+* **Author:** This model is based on form factor calculations implemented
+    in a c-library provided by the NIST Center for Neutron Research (Kline, 2006).
+* **Last Modified by:**  Paul Kienzle **Date:** April 05, 2017
+* **Last Reviewed by:**  Richard Heenan **Date:** April 06, 2017
+
+"""
+
+import numpy as np
+from numpy import pi, inf, sqrt, sin, cos, radians
+from sasmodels.special import sas_sinx_x, Gauss76Wt, Gauss76Z, square, ORIENT_ASYMMETRIC
+
+name = "parallelepiped"
+title = "Rectangular parallelepiped with uniform scattering length density."
+description = """
+    I(q)= scale*V*(sld - sld_solvent)^2*P(q,alpha)+background
+        P(q,alpha) = integral from 0 to 1 of ...
+           phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma
+        with
+            phi(mu,a) = integral from 0 to 1 of ..
+            (S((/mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du
+            S(x) = sin(x)/x
+            mu = q*B
+        V: Volume of the rectangular parallelepiped
+        alpha: angle between the long axis of the
+            parallelepiped and the q-vector for 1D
+"""
+category = "shape:parallelepiped"
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
+               "Parallelepiped scattering length density"],
+              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
+               "Solvent scattering length density"],
+              ["length_a", "Ang", 35, [0, inf], "volume",
+               "Shorter side of the parallelepiped"],
+              ["length_b", "Ang", 75, [0, inf], "volume",
+               "Second side of the parallelepiped"],
+              ["length_c", "Ang", 400, [0, inf], "volume",
+               "Larger side of the parallelepiped"],
+              ["theta", "degrees", 60, [-360, 360], "orientation",
+               "c axis to beam angle"],
+              ["phi", "degrees", 60, [-360, 360], "orientation",
+               "rotation about beam"],
+              ["psi", "degrees", 60, [-360, 360], "orientation",
+               "rotation about c axis"],
+             ]
+
+def form_volume(length_a, length_b, length_c):
+    return length_a * length_b * length_c
+             
+def Iq(q,
+        sld, solvent_sld,
+        length_a, length_b, length_c):
+
+    mu = 0.5 * q * length_b
+    
+    # Scale sides by B
+    a_scaled = length_a / length_b
+    c_scaled = length_c / length_b
+        
+    # outer integral (with gauss points), integration limits = 0, 1
+    outer_total = 0 #initialize integral
+
+    for i in range(76):
+        sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
+        mu_proj = mu * sqrt(1.0-sigma*sigma)
+
+        # inner integral (with gauss points), integration limits = 0, 1
+        # corresponding to angles from 0 to pi/2.
+        inner_total = 0.0
+        for j in range(76):
+            uu = 0.5 * ( Gauss76Z[j] + 1.0 );
+            sin_uu, cos_uu = sin(0.5*pi*uu), cos(0.5*pi*uu)
+            si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled)
+            si2 = sas_sinx_x(mu_proj * cos_uu)
+            inner_total += Gauss76Wt[j] * square(si1 * si2)
+        inner_total *= 0.5
+
+        si = sas_sinx_x(mu * c_scaled * sigma)
+        outer_total += Gauss76Wt[i] * inner_total * si * si
+    outer_total *= 0.5
+
+    # Multiply by contrast^2 and convert from [1e-12 A-1] to [cm-1]
+    V = form_volume(length_a, length_b, length_c)
+    drho = (sld-solvent_sld)
+    return 1.0e-4 * square(drho * V) * outer_total
+Iq.vectorized = True
+
+def Iqxy(qx, qy,
+        sld, solvent_sld,
+        length_a, length_b, length_c,
+        theta, phi, psi):
+    q, xhat, yhat, zhat = ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi)
+    siA = sas_sinx_x(0.5*length_a*q*xhat)
+    siB = sas_sinx_x(0.5*length_b*q*yhat)
+    siC = sas_sinx_x(0.5*length_c*q*zhat)
+    V = form_volume(length_a, length_b, length_c)
+    drho = (sld - solvent_sld)
+    form = V * drho * siA * siB * siC
+    # Square and convert from [1e-12 A-1] to [cm-1]
+    return 1.0e-4 * form * form
+Iqxy.vectorized = True
+
+
+def ER(length_a, length_b, length_c):
+    """
+    Return effective radius (ER) for P(q)*S(q)
+    """
+    # now that axes can be in any size order, need to sort a,b,c where a~b and c is either much smaller
+    # or much larger
+    abc = np.vstack((length_a, length_b, length_c))
+    abc = np.sort(abc, axis=0)
+    selector = (abc[1] - abc[0]) > (abc[2] - abc[1])
+    length = np.where(selector, abc[0], abc[2])
+    # surface average radius (rough approximation)
+    radius = np.sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi)
+
+    ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius))
+    return 0.5 * (ddd) ** (1. / 3.)
+
+# VR defaults to 1.0
+
+# parameters for demo
+demo = dict(scale=1, background=0,
+            sld=6.3, sld_solvent=1.0,
+            length_a=35, length_b=75, length_c=400,
+            theta=45, phi=30, psi=15,
+            length_a_pd=0.1, length_a_pd_n=10,
+            length_b_pd=0.1, length_b_pd_n=1,
+            length_c_pd=0.1, length_c_pd_n=1,
+            theta_pd=10, theta_pd_n=1,
+            phi_pd=10, phi_pd_n=1,
+            psi_pd=10, psi_pd_n=10)
+# rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models
+qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
+tests = [[{}, 0.2, 0.17758004974],
+         [{}, [0.2], [0.17758004974]],
+         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0089517140475],
+         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0089517140475]],
+        ]
+del qx, qy  # not necessary to delete, but cleaner
Index: sasmodels/models/barbell.c
===================================================================
--- sasmodels/models/barbell.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
+++ sasmodels/models/barbell.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -10,6 +10,6 @@
 //barbell kernel - same as dumbell
 static double
-_bell_kernel(double q, double h, double radius_bell,
-             double half_length, double sin_alpha, double cos_alpha)
+_bell_kernel(double qab, double qc, double h, double radius_bell,
+             double half_length)
 {
     // translate a point in [-1,1] to a point in [lower,upper]
@@ -26,12 +26,12 @@
     //    m = q R cos(alpha)
     //    b = q(L/2-h) cos(alpha)
-    const double m = q*radius_bell*cos_alpha; // cos argument slope
-    const double b = q*(half_length-h)*cos_alpha; // cos argument intercept
-    const double qrst = q*radius_bell*sin_alpha; // Q*R*sin(theta)
+    const double m = radius_bell*qc; // cos argument slope
+    const double b = (half_length-h)*qc; // cos argument intercept
+    const double qab_r = radius_bell*qab; // Q*R*sin(theta)
     double total = 0.0;
     for (int i = 0; i < 76; i++){
         const double t = Gauss76Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
-        const double bj = sas_2J1x_x(qrst*sqrt(radical));
+        const double bj = sas_2J1x_x(qab_r*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
         total += Gauss76Wt[i] * Fq;
@@ -44,11 +44,10 @@
 
 static double
-_fq(double q, double h,
-    double radius_bell, double radius, double half_length,
-    double sin_alpha, double cos_alpha)
+_fq(double qab, double qc, double h,
+    double radius_bell, double radius, double half_length)
 {
-    const double bell_fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha);
-    const double bj = sas_2J1x_x(q*radius*sin_alpha);
-    const double si = sas_sinx_x(q*half_length*cos_alpha);
+    const double bell_fq = _bell_kernel(qab, qc, h, radius_bell, half_length);
+    const double bj = sas_2J1x_x(radius*qab);
+    const double si = sas_sinx_x(half_length*qc);
     const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si;
     const double Aq = bell_fq + cyl_fq;
@@ -84,5 +83,5 @@
         double sin_alpha, cos_alpha; // slots to hold sincos function output
         SINCOS(alpha, sin_alpha, cos_alpha);
-        const double Aq = _fq(q, h, radius_bell, radius, half_length, sin_alpha, cos_alpha);
+        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length);
         total += Gauss76Wt[i] * Aq * Aq * sin_alpha;
     }
@@ -103,7 +102,9 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
     const double h = -sqrt(square(radius_bell) - square(radius));
-    const double Aq = _fq(q, h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha);
+    const double Aq = _fq(qab, qc, h, radius_bell, radius, 0.5*length);
 
     // Multiply by contrast^2 and convert to cm-1
Index: sasmodels/models/bcc_paracrystal.c
===================================================================
--- sasmodels/models/bcc_paracrystal.c	(revision 50beefeac18f2f952bf0863526ec9bbd3009e29c)
+++ sasmodels/models/bcc_paracrystal.c	(revision 7e0b281e5f175a91044364982fcfbfe802508948)
@@ -1,89 +1,95 @@
-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);
+static double
+bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
+{
+#if 0  // Equations as written in Matsuoka
+    const double a1 = (+qa + qb + qc)/2.0;
+    const double a2 = (-qa - qb + qc)/2.0;
+    const double a3 = (-qa + qb - qc)/2.0;
+#else
+    const double a1 = (+qa + qb - qc)/2.0;
+    const double a2 = (+qa - qb + qc)/2.0;
+    const double a3 = (-qa + qb + qc)/2.0;
+#endif
 
-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);
+#if 1
+    // Numerator: (1 - exp(a)^2)^3
+    //         => (-(exp(2a) - 1))^3
+    //         => -expm1(2a)^3
+    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a))
+    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1)
+    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1)
+    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    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));
+#else
+    // Alternate form, which perhaps is more approachable
+    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    const double sinh_qd = sinh(arg);
+    const double cosh_qd = cosh(arg);
+    const double Zq = sinh_qd/(cosh_qd - cos(dnn*a1))
+                    * sinh_qd/(cosh_qd - cos(dnn*a2))
+                    * sinh_qd/(cosh_qd - cos(dnn*a3));
+#endif
+
+    return Zq;
+}
 
 
-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);
+// occupied volume fraction calculated from lattice symmetry and sphere radius
+static double
+bcc_volume_fraction(double radius, double dnn)
+{
+    return 2.0*sphere_volume(sqrt(0.75)*radius/dnn);
 }
 
-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 = cube(1.0 - (temp10*temp10))*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){
+static double
+form_volume(double radius)
+{
     return sphere_volume(radius);
 }
 
 
-double Iq(double q, double dnn,
+static double Iq(double q, double dnn,
   double d_factor, double radius,
-  double sld, double solvent_sld){
+  double sld, double solvent_sld)
+{
+    // translate a point in [-1,1] to a point in [0, 2 pi]
+    const double phi_m = M_PI;
+    const double phi_b = M_PI;
+    // translate a point in [-1,1] to a point in [0, pi]
+    const double theta_m = M_PI_2;
+    const double theta_b = M_PI_2;
 
-	//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 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 form = bcc_Zq(qa, qb, qc, dnn, d_factor);
+            inner_sum += Gauss150Wt[j] * form;
+        }
+        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
+        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+    }
+    outer_sum *= theta_m;
+    const double Zq = outer_sum/(4.0*M_PI);
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+    return bcc_volume_fraction(radius, dnn) * Pq * Zq;
 }
 
 
-double Iqxy(double qx, double qy,
+static double Iqxy(double qx, double qy,
     double dnn, double d_factor, double radius,
     double sld, double solvent_sld,
@@ -92,20 +98,11 @@
     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;
 
-    const double a1 = +xhat - zhat + yhat;
-    const double a2 = +xhat + zhat - yhat;
-    const double a3 = -xhat + zhat + yhat;
-
-    const double qd = 0.5*q*dnn;
-    const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3);
-    const double tanh_qd = tanh(arg);
-    const double cosh_qd = cosh(arg);
-    const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*a2)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*a3)/cosh_qd);
-
-    const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq;
-    //the occupied volume of the lattice
-    const double lattice_scale = 2.0*sphere_volume(sqrt(0.75)*radius/dnn);
-    return lattice_scale * Fq;
+    q = sqrt(qa*qa + qb*qb + qc*qc);
+    const double Zq = bcc_Zq(qa, qb, qc, dnn, d_factor);
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+    return bcc_volume_fraction(radius, dnn) * Pq * Zq;
 }
Index: sasmodels/models/bcc_paracrystal.py
===================================================================
--- sasmodels/models/bcc_paracrystal.py	(revision 69e1afc57270f67b2a0d7157a6d263fc2db98cc0)
+++ sasmodels/models/bcc_paracrystal.py	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -81,5 +81,5 @@
 .. figure:: img/parallelepiped_angle_definition.png
 
-    Orientation of the crystal with respect to the scattering plane, when 
+    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).
 
Index: sasmodels/models/capped_cylinder.c
===================================================================
--- sasmodels/models/capped_cylinder.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
+++ sasmodels/models/capped_cylinder.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -14,8 +14,8 @@
 //   radius_cap is the radius of the lens
 //   length is the cylinder length, or the separation between the lens halves
-//   alpha is the angle of the cylinder wrt q.
+//   theta is the angle of the cylinder wrt q.
 static double
-_cap_kernel(double q, double h, double radius_cap,
-                      double half_length, double sin_alpha, double cos_alpha)
+_cap_kernel(double qab, double qc, double h, double radius_cap,
+            double half_length)
 {
     // translate a point in [-1,1] to a point in [lower,upper]
@@ -26,18 +26,18 @@
 
     // cos term in integral is:
-    //    cos (q (R t - h + L/2) cos(alpha))
+    //    cos (q (R t - h + L/2) cos(theta))
     // so turn it into:
     //    cos (m t + b)
     // where:
-    //    m = q R cos(alpha)
-    //    b = q(L/2-h) cos(alpha)
-    const double m = q*radius_cap*cos_alpha; // cos argument slope
-    const double b = q*(half_length-h)*cos_alpha; // cos argument intercept
-    const double qrst = q*radius_cap*sin_alpha; // Q*R*sin(theta)
+    //    m = q R cos(theta)
+    //    b = q(L/2-h) cos(theta)
+    const double m = radius_cap*qc; // cos argument slope
+    const double b = (half_length-h)*qc; // cos argument intercept
+    const double qab_r = radius_cap*qab; // Q*R*sin(theta)
     double total = 0.0;
     for (int i=0; i<76 ;i++) {
         const double t = Gauss76Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
-        const double bj = sas_2J1x_x(qrst*sqrt(radical));
+        const double bj = sas_2J1x_x(qab_r*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
         total += Gauss76Wt[i] * Fq;
@@ -50,10 +50,9 @@
 
 static double
-_fq(double q, double h, double radius_cap, double radius, double half_length,
-    double sin_alpha, double cos_alpha)
+_fq(double qab, double qc, double h, double radius_cap, double radius, double half_length)
 {
-    const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha);
-    const double bj = sas_2J1x_x(q*radius*sin_alpha);
-    const double si = sas_sinx_x(q*half_length*cos_alpha);
+    const double cap_Fq = _cap_kernel(qab, qc, h, radius_cap, half_length);
+    const double bj = sas_2J1x_x(radius*qab);
+    const double si = sas_sinx_x(half_length*qc);
     const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si;
     const double Aq = cap_Fq + cyl_Fq;
@@ -101,11 +100,12 @@
     double total = 0.0;
     for (int i=0; i<76 ;i++) {
-        const double alpha= Gauss76Z[i]*zm + zb;
-        double sin_alpha, cos_alpha; // slots to hold sincos function output
-        SINCOS(alpha, sin_alpha, cos_alpha);
-
-        const double Aq = _fq(q, h, radius_cap, radius, half_length, sin_alpha, cos_alpha);
-        // sin_alpha for spherical coord integration
-        total += Gauss76Wt[i] * Aq * Aq * sin_alpha;
+        const double theta = Gauss76Z[i]*zm + zb;
+        double sin_theta, cos_theta; // slots to hold sincos function output
+        SINCOS(theta, sin_theta, cos_theta);
+        const double qab = q*sin_theta;
+        const double qc = q*cos_theta;
+        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length);
+        // scale by sin_theta for spherical coord integration
+        total += Gauss76Wt[i] * Aq * Aq * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -125,7 +125,9 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
     const double h = sqrt(radius_cap*radius_cap - radius*radius);
-    const double Aq = _fq(q, h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha);
+    const double Aq = _fq(qab, qc, h, radius_cap, radius, 0.5*length);
 
     // Multiply by contrast^2 and convert to cm-1
Index: sasmodels/models/core_shell_bicelle.c
===================================================================
--- sasmodels/models/core_shell_bicelle.c	(revision b260926c42c2ae859a32c20b98935378d51fab4d)
+++ sasmodels/models/core_shell_bicelle.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,56 +1,31 @@
-double form_volume(double radius, double thick_rim, double thick_face, double length);
-double Iq(double q,
-          double radius,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld);
-
-
-double Iqxy(double qx, double qy,
-          double radius,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld,
-          double theta,
-          double phi);
-
-
-double form_volume(double radius, double thick_rim, double thick_face, double length)
+static double
+form_volume(double radius, double thick_rim, double thick_face, double length)
 {
-    return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face);
+    return M_PI*square(radius+thick_rim)*(length+2.0*thick_face);
 }
 
 static double
-bicelle_kernel(double q,
-              double rad,
-              double radthick,
-              double facthick,
-              double halflength,
-              double rhoc,
-              double rhoh,
-              double rhor,
-              double rhosolv,
-              double sin_alpha,
-              double cos_alpha)
+bicelle_kernel(double qab,
+    double qc,
+    double radius,
+    double thick_radius,
+    double thick_face,
+    double halflength,
+    double sld_core,
+    double sld_face,
+    double sld_rim,
+    double sld_solvent)
 {
-    const double dr1 = rhoc-rhoh;
-    const double dr2 = rhor-rhosolv;
-    const double dr3 = rhoh-rhor;
-    const double vol1 = M_PI*square(rad)*2.0*(halflength);
-    const double vol2 = M_PI*square(rad+radthick)*2.0*(halflength+facthick);
-    const double vol3 = M_PI*square(rad)*2.0*(halflength+facthick);
+    const double dr1 = sld_core-sld_face;
+    const double dr2 = sld_rim-sld_solvent;
+    const double dr3 = sld_face-sld_rim;
+    const double vol1 = M_PI*square(radius)*2.0*(halflength);
+    const double vol2 = M_PI*square(radius+thick_radius)*2.0*(halflength+thick_face);
+    const double vol3 = M_PI*square(radius)*2.0*(halflength+thick_face);
 
-    const double be1 = sas_2J1x_x(q*(rad)*sin_alpha);
-    const double be2 = sas_2J1x_x(q*(rad+radthick)*sin_alpha);
-    const double si1 = sas_sinx_x(q*(halflength)*cos_alpha);
-    const double si2 = sas_sinx_x(q*(halflength+facthick)*cos_alpha);
+    const double be1 = sas_2J1x_x((radius)*qab);
+    const double be2 = sas_2J1x_x((radius+thick_radius)*qab);
+    const double si1 = sas_sinx_x((halflength)*qc);
+    const double si2 = sas_sinx_x((halflength+thick_face)*qc);
 
     const double t = vol1*dr1*si1*be1 +
@@ -58,20 +33,17 @@
                      vol3*dr3*si2*be1;
 
-    const double retval = t*t;
-
-    return retval;
-
+    return t;
 }
 
 static double
-bicelle_integration(double q,
-                   double rad,
-                   double radthick,
-                   double facthick,
-                   double length,
-                   double rhoc,
-                   double rhoh,
-                   double rhor,
-                   double rhosolv)
+Iq(double q,
+    double radius,
+    double thick_radius,
+    double thick_face,
+    double length,
+    double sld_core,
+    double sld_face,
+    double sld_rim,
+    double sld_solvent)
 {
     // set up the integration end points
@@ -79,82 +51,40 @@
     const double halflength = 0.5*length;
 
-    double summ = 0.0;
+    double total = 0.0;
     for(int i=0;i<N_POINTS_76;i++) {
-        double alpha = (Gauss76Z[i] + 1.0)*uplim;
-        double sin_alpha, cos_alpha; // slots to hold sincos function output
-        SINCOS(alpha, sin_alpha, cos_alpha);
-        double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick,
-                             halflength, rhoc, rhoh, rhor, rhosolv,
-                             sin_alpha, cos_alpha);
-        summ += yyy*sin_alpha;
+        double theta = (Gauss76Z[i] + 1.0)*uplim;
+        double sin_theta, cos_theta; // slots to hold sincos function output
+        SINCOS(theta, sin_theta, cos_theta);
+        double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face,
+                                   halflength, sld_core, sld_face, sld_rim, sld_solvent);
+        total += Gauss76Wt[i]*fq*fq*sin_theta;
     }
 
     // calculate value of integral to return
-    double answer = uplim*summ;
-    return answer;
+    double answer = total*uplim;
+    return 1.0e-4*answer;
 }
 
 static double
-bicelle_kernel_2d(double qx, double qy,
-          double radius,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld,
-          double theta,
-          double phi)
+Iqxy(double qx, double qy,
+    double radius,
+    double thick_rim,
+    double thick_face,
+    double length,
+    double core_sld,
+    double face_sld,
+    double rim_sld,
+    double solvent_sld,
+    double theta,
+    double phi)
 {
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
-    double answer = bicelle_kernel(q, radius, thick_rim, thick_face,
+    double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face,
                            0.5*length, core_sld, face_sld, rim_sld,
-                           solvent_sld, sin_alpha, cos_alpha);
-    return 1.0e-4*answer;
+                           solvent_sld);
+    return 1.0e-4*fq*fq;
 }
-
-double Iq(double q,
-          double radius,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld)
-{
-    double intensity = bicelle_integration(q, radius, thick_rim, thick_face,
-                       length, core_sld, face_sld, rim_sld, solvent_sld);
-    return intensity*1.0e-4;
-}
-
-
-double Iqxy(double qx, double qy,
-          double radius,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld,
-          double theta,
-          double phi)
-{
-    double intensity = bicelle_kernel_2d(qx, qy,
-                      radius,
-                      thick_rim,
-                      thick_face,
-                      length,
-                      core_sld,
-                      face_sld,
-                      rim_sld,
-                      solvent_sld,
-                      theta,
-                      phi);
-
-    return intensity;
-}
Index: sasmodels/models/core_shell_bicelle_elliptical.c
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical.c	(revision dedcf34f547ba57a5e8ca3c5e9892ede11d1fcaf)
+++ sasmodels/models/core_shell_bicelle_elliptical.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -17,29 +17,21 @@
         double thick_face,
         double length,
-        double rhoc,
-        double rhoh,
-        double rhor,
-        double rhosolv)
+        double sld_core,
+        double sld_face,
+        double sld_rim,
+        double sld_solvent)
 {
-    double si1,si2,be1,be2;
      // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle
      // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle
-     //    const double uplim = M_PI_4;
     const double halfheight = 0.5*length;
-    //const double va = 0.0;
-    //const double vb = 1.0;
-    // inner integral limits
-    //const double vaj=0.0;
-    //const double vbj=M_PI;
-
     const double r_major = r_minor * x_core;
     const double r2A = 0.5*(square(r_major) + square(r_minor));
     const double r2B = 0.5*(square(r_major) - square(r_minor));
-    const double dr1 = (rhoc-rhoh)   *M_PI*r_minor*r_major*(2.0*halfheight);;
-    const double dr2 = (rhor-rhosolv)*M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);
-    const double dr3 = (rhoh-rhor)   *M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
-    //const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight);
-    //const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);
-    //const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
+    const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight);
+    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);
+    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
+    const double dr1 = vol1*(sld_core-sld_face);
+    const double dr2 = vol2*(sld_rim-sld_solvent);
+    const double dr3 = vol3*(sld_face-sld_rim);
 
     //initialize integral
@@ -47,25 +39,27 @@
     for(int i=0;i<76;i++) {
         //setup inner integral over the ellipsoidal cross-section
-        // since we generate these lots of times, why not store them somewhere?
-        //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
-        const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0;
-        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
-        double inner_sum=0;
-        double sinarg1 = q*halfheight*cos_alpha;
-        double sinarg2 = q*(halfheight+thick_face)*cos_alpha;
-        si1 = sas_sinx_x(sinarg1);
-        si2 = sas_sinx_x(sinarg2);
+        //const double va = 0.0;
+        //const double vb = 1.0;
+        //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
+        const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0;
+        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
+        const double qab = q*sin_theta;
+        const double qc = q*cos_theta;
+        const double si1 = sas_sinx_x(halfheight*qc);
+        const double si2 = sas_sinx_x((halfheight+thick_face)*qc);
+        double inner_sum=0.0;
         for(int j=0;j<76;j++) {
             //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe)
-            //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
-            const double beta = ( Gauss76Z[j] +1.0)*M_PI_2;
-            const double rr = sqrt(r2A - r2B*cos(beta));
-            double besarg1 = q*rr*sin_alpha;
-            double besarg2 = q*(rr+thick_rim)*sin_alpha;
-            be1 = sas_2J1x_x(besarg1);
-            be2 = sas_2J1x_x(besarg2);
-            inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 +
-                                              dr2*si2*be2 +
-                                              dr3*si2*be1);
+            // inner integral limits
+            //const double vaj=0.0;
+            //const double vbj=M_PI;
+            //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
+            const double phi = ( Gauss76Z[j] +1.0)*M_PI_2;
+            const double rr = sqrt(r2A - r2B*cos(phi));
+            const double be1 = sas_2J1x_x(rr*qab);
+            const double be2 = sas_2J1x_x((rr+thick_rim)*qab);
+            const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1;
+
+            inner_sum += Gauss76Wt[j] * fq * fq;
         }
         //now calculate outer integral
@@ -83,18 +77,21 @@
           double thick_face,
           double length,
-          double rhoc,
-          double rhoh,
-          double rhor,
-          double rhosolv,
+          double sld_core,
+          double sld_face,
+          double sld_rim,
+          double sld_solvent,
           double theta,
           double phi,
           double psi)
 {
-       // THIS NEEDS TESTING
     double q, xhat, yhat, zhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
-    const double dr1 = rhoc-rhoh;
-    const double dr2 = rhor-rhosolv;
-    const double dr3 = rhoh-rhor;
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
+
+    const double dr1 = sld_core-sld_face;
+    const double dr2 = sld_rim-sld_solvent;
+    const double dr3 = sld_face-sld_rim;
     const double r_major = r_minor*x_core;
     const double halfheight = 0.5*length;
@@ -104,13 +101,12 @@
 
     // Compute effective radius in rotated coordinates
-    const double r_hat = sqrt(square(r_major*xhat) + square(r_minor*yhat));
-    const double rshell_hat = sqrt(square((r_major+thick_rim)*xhat)
-                                   + square((r_minor+thick_rim)*yhat));
-    const double be1 = sas_2J1x_x( q*r_hat );
-    const double be2 = sas_2J1x_x( q*rshell_hat );
-    const double si1 = sas_sinx_x( q*halfheight*zhat );
-    const double si2 = sas_sinx_x( q*(halfheight + thick_face)*zhat );
-    const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1);
-    return 1.0e-4 * Aq;
+    const double qr_hat = sqrt(square(r_major*qa) + square(r_minor*qb));
+    const double qrshell_hat = sqrt(square((r_major+thick_rim)*qa)
+                                   + square((r_minor+thick_rim)*qb));
+    const double be1 = sas_2J1x_x( qr_hat );
+    const double be2 = sas_2J1x_x( qrshell_hat );
+    const double si1 = sas_sinx_x( halfheight*qc );
+    const double si2 = sas_sinx_x( (halfheight + thick_face)*qc );
+    const double fq = vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1;
+    return 1.0e-4 * fq*fq;
 }
-
Index: sasmodels/models/core_shell_cylinder.c
===================================================================
--- sasmodels/models/core_shell_cylinder.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
+++ sasmodels/models/core_shell_cylinder.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,23 +1,18 @@
-double form_volume(double radius, double thickness, double length);
-double Iq(double q, double core_sld, double shell_sld, double solvent_sld,
-    double radius, double thickness, double length);
-double Iqxy(double qx, double qy, double core_sld, double shell_sld, double solvent_sld,
-    double radius, double thickness, double length, double theta, double phi);
-
 // vd = volume * delta_rho
-// besarg = q * R * sin(alpha)
-// siarg = q * L/2 * cos(alpha)
-double _cyl(double vd, double besarg, double siarg);
-double _cyl(double vd, double besarg, double siarg)
+// besarg = q * R * sin(theta)
+// siarg = q * L/2 * cos(theta)
+static double _cyl(double vd, double besarg, double siarg)
 {
     return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg);
 }
 
-double form_volume(double radius, double thickness, double length)
+static double
+form_volume(double radius, double thickness, double length)
 {
-    return M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness);
+    return M_PI*square(radius+thickness)*(length+2.0*thickness);
 }
 
-double Iq(double q,
+static double
+Iq(double q,
     double core_sld,
     double shell_sld,
@@ -28,21 +23,22 @@
 {
     // precalculate constants
-    const double core_qr = q*radius;
-    const double core_qh = q*0.5*length;
+    const double core_r = radius;
+    const double core_h = 0.5*length;
     const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld);
-    const double shell_qr = q*(radius + thickness);
-    const double shell_qh = q*(0.5*length + thickness);
+    const double shell_r = (radius + thickness);
+    const double shell_h = (0.5*length + thickness);
     const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld);
     double total = 0.0;
-    // double lower=0, upper=M_PI_2;
     for (int i=0; i<76 ;i++) {
-        // translate a point in [-1,1] to a point in [lower,upper]
-        //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
-        double sn, cn;
-        const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2);
-        SINCOS(alpha, sn, cn);
-        const double fq = _cyl(core_vd, core_qr*sn, core_qh*cn)
-            + _cyl(shell_vd, shell_qr*sn, shell_qh*cn);
-        total += Gauss76Wt[i] * fq * fq * sn;
+        // translate a point in [-1,1] to a point in [0, pi/2]
+        //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
+        double sin_theta, cos_theta;
+        const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4;
+        SINCOS(theta, sin_theta,  cos_theta);
+        const double qab = q*sin_theta;
+        const double qc = q*cos_theta;
+        const double fq = _cyl(core_vd, core_r*qab, core_h*qc)
+            + _cyl(shell_vd, shell_r*qab, shell_h*qc);
+        total += Gauss76Wt[i] * fq * fq * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -64,14 +60,16 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
-    const double core_qr = q*radius;
-    const double core_qh = q*0.5*length;
+    const double core_r = radius;
+    const double core_h = 0.5*length;
     const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld);
-    const double shell_qr = q*(radius + thickness);
-    const double shell_qh = q*(0.5*length + thickness);
+    const double shell_r = (radius + thickness);
+    const double shell_h = (0.5*length + thickness);
     const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld);
 
-    const double fq = _cyl(core_vd, core_qr*sin_alpha, core_qh*cos_alpha)
-        + _cyl(shell_vd, shell_qr*sin_alpha, shell_qh*cos_alpha);
+    const double fq = _cyl(core_vd, core_r*qab, core_h*qc)
+        + _cyl(shell_vd, shell_r*qab, shell_h*qc);
     return 1.0e-4 * fq * fq;
 }
Index: sasmodels/models/core_shell_ellipsoid.c
===================================================================
--- sasmodels/models/core_shell_ellipsoid.c	(revision 0a3d9b2c366d93b98a02e5c9077f26098e630bd0)
+++ sasmodels/models/core_shell_ellipsoid.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,32 +1,34 @@
-double form_volume(double radius_equat_core,
-                   double polar_core,
-                   double equat_shell,
-                   double polar_shell);
-double Iq(double q,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld);
 
+// Converted from Igor function gfn4, using the same pattern as ellipsoid
+// for evaluating the parts of the integral.
+//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN
+//                       BY (53) & (58-59) IN CHEN AND
+//                       KOTLARCHYK REFERENCE
+//
+//       <OBLATE ELLIPSOID>
+static double
+_cs_ellipsoid_kernel(double qab, double qc,
+    double equat_core, double polar_core,
+    double equat_shell, double polar_shell,
+    double sld_core_shell, double sld_shell_solvent)
+{
+    const double qr_core = sqrt(square(equat_core*qab) + square(polar_core*qc));
+    const double si_core = sas_3j1x_x(qr_core);
+    const double volume_core = M_4PI_3*equat_core*equat_core*polar_core;
+    const double fq_core = si_core*volume_core*sld_core_shell;
 
-double Iqxy(double qx, double qy,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double theta,
-          double phi);
+    const double qr_shell = sqrt(square(equat_shell*qab) + square(polar_shell*qc));
+    const double si_shell = sas_3j1x_x(qr_shell);
+    const double volume_shell = M_4PI_3*equat_shell*equat_shell*polar_shell;
+    const double fq_shell = si_shell*volume_shell*sld_shell_solvent;
 
+    return fq_core + fq_shell;
+}
 
-double form_volume(double radius_equat_core,
-                   double x_core,
-                   double thick_shell,
-                   double x_polar_shell)
+static double
+form_volume(double radius_equat_core,
+    double x_core,
+    double thick_shell,
+    double x_polar_shell)
 {
     const double equat_shell = radius_equat_core + thick_shell;
@@ -37,20 +39,15 @@
 
 static double
-core_shell_ellipsoid_xt_kernel(double q,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld)
+Iq(double q,
+    double radius_equat_core,
+    double x_core,
+    double thick_shell,
+    double x_polar_shell,
+    double core_sld,
+    double shell_sld,
+    double solvent_sld)
 {
-    const double lolim = 0.0;
-    const double uplim = 1.0;
-
-
-    const double delpc = core_sld - shell_sld; //core - shell
-    const double delps = shell_sld - solvent_sld; //shell - solvent
-
+    const double sld_core_shell = core_sld - shell_sld;
+    const double sld_shell_solvent = shell_sld - solvent_sld;
 
     const double polar_core = radius_equat_core*x_core;
@@ -58,34 +55,42 @@
     const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell;
 
-    double summ = 0.0;	 //initialize intergral
+    // translate from [-1, 1] => [0, 1]
+    const double m = 0.5;
+    const double b = 0.5;
+    double total = 0.0;     //initialize intergral
     for(int i=0;i<76;i++) {
-        double zi = 0.5*( Gauss76Z[i]*(uplim-lolim) + uplim + lolim );
-        double yyy = gfn4(zi, radius_equat_core, polar_core, equat_shell,
-                          polar_shell, delpc, delps, q);
-        summ += Gauss76Wt[i] * yyy;
+        const double cos_theta = Gauss76Z[i]*m + b;
+        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
+        double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta,
+            radius_equat_core, polar_core,
+            equat_shell, polar_shell,
+            sld_core_shell, sld_shell_solvent);
+        total += Gauss76Wt[i] * fq * fq;
     }
-    summ *= 0.5*(uplim-lolim);
+    total *= m;
 
     // convert to [cm-1]
-    return 1.0e-4 * summ;
+    return 1.0e-4 * total;
 }
 
 static double
-core_shell_ellipsoid_xt_kernel_2d(double qx, double qy,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double theta,
-          double phi)
+Iqxy(double qx, double qy,
+    double radius_equat_core,
+    double x_core,
+    double thick_shell,
+    double x_polar_shell,
+    double core_sld,
+    double shell_sld,
+    double solvent_sld,
+    double theta,
+    double phi)
 {
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
-    const double sldcs = core_sld - shell_sld;
-    const double sldss = shell_sld- solvent_sld;
+    const double sld_core_shell = core_sld - shell_sld;
+    const double sld_shell_solvent = shell_sld - solvent_sld;
 
     const double polar_core = radius_equat_core*x_core;
@@ -93,65 +98,10 @@
     const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell;
 
-    // Call the IGOR library function to get the kernel:
-    // MUST use gfn4 not gf2 because of the def of params.
-    double answer = gfn4(cos_alpha,
-                  radius_equat_core,
-                  polar_core,
-                  equat_shell,
-                  polar_shell,
-                  sldcs,
-                  sldss,
-                  q);
+    double fq = _cs_ellipsoid_kernel(qab, qc,
+                  radius_equat_core, polar_core,
+                  equat_shell, polar_shell,
+                  sld_core_shell, sld_shell_solvent);
 
     //convert to [cm-1]
-    answer *= 1.0e-4;
-
-    return answer;
+    return 1.0e-4 * fq * fq;
 }
-
-double Iq(double q,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld)
-{
-    double intensity = core_shell_ellipsoid_xt_kernel(q,
-           radius_equat_core,
-           x_core,
-           thick_shell,
-           x_polar_shell,
-           core_sld,
-           shell_sld,
-           solvent_sld);
-
-    return intensity;
-}
-
-
-double Iqxy(double qx, double qy,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double theta,
-          double phi)
-{
-    double intensity = core_shell_ellipsoid_xt_kernel_2d(qx, qy,
-                       radius_equat_core,
-                       x_core,
-                       thick_shell,
-                       x_polar_shell,
-                       core_sld,
-                       shell_sld,
-                       solvent_sld,
-                       theta,
-                       phi);
-
-    return intensity;
-}
Index: sasmodels/models/core_shell_ellipsoid.py
===================================================================
--- sasmodels/models/core_shell_ellipsoid.py	(revision 9802ab3945a154b871264715dcb463d4ec644afa)
+++ sasmodels/models/core_shell_ellipsoid.py	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -25,6 +25,6 @@
 ellipsoid. This may have some undesirable effects if the aspect ratio of the
 ellipsoid is large (ie, if $X << 1$ or $X >> 1$ ), when the $S(q)$
-- which assumes spheres - will not in any case be valid.  Generating a 
-custom product model will enable separate effective volume fraction and effective 
+- which assumes spheres - will not in any case be valid.  Generating a
+custom product model will enable separate effective volume fraction and effective
 radius in the $S(q)$.
 
@@ -44,11 +44,11 @@
 
 .. math::
-    \begin{align}    
+    \begin{align}
     F(q,\alpha) = &f(q,radius\_equat\_core,radius\_equat\_core.x\_core,\alpha) \\
     &+ f(q,radius\_equat\_core + thick\_shell,radius\_equat\_core.x\_core + thick\_shell.x\_polar\_shell,\alpha)
-    \end{align} 
+    \end{align}
 
 where
- 
+
 .. math::
 
@@ -77,5 +77,5 @@
    F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}
 
-For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data, 
+For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,
 see the :ref:`elliptical-cylinder` model for further information.
 
@@ -139,6 +139,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sas_3j1x_x.c", "lib/gfn.c", "lib/gauss76.c",
-          "core_shell_ellipsoid.c"]
+source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"]
 
 def ER(radius_equat_core, x_core, thick_shell, x_polar_shell):
Index: sasmodels/models/core_shell_parallelepiped.c
===================================================================
--- sasmodels/models/core_shell_parallelepiped.c	(revision 92dfe0c540274c8b800a9929c3f1e1a7d4cba9dc)
+++ sasmodels/models/core_shell_parallelepiped.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,3 +1,3 @@
-double form_volume(double length_a, double length_b, double length_c, 
+double form_volume(double length_a, double length_b, double length_c,
                    double thick_rim_a, double thick_rim_b, double thick_rim_c);
 double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld,
@@ -9,10 +9,10 @@
             double thick_rim_c, double theta, double phi, double psi);
 
-double form_volume(double length_a, double length_b, double length_c, 
+double form_volume(double length_a, double length_b, double length_c,
                    double thick_rim_a, double thick_rim_b, double thick_rim_c)
 {
     //return length_a * length_b * length_c;
-    return length_a * length_b * length_c + 
-           2.0 * thick_rim_a * length_b * length_c + 
+    return length_a * length_b * length_c +
+           2.0 * thick_rim_a * length_b * length_c +
            2.0 * thick_rim_b * length_a * length_c +
            2.0 * thick_rim_c * length_a * length_b;
@@ -34,14 +34,14 @@
     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled
     // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
-    
+
     const double mu = 0.5 * q * length_b;
-    
+
     //calculate volume before rescaling (in original code, but not used)
-    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);        
-    //double vol = length_a * length_b * length_c + 
-    //       2.0 * thick_rim_a * length_b * length_c + 
+    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
+    //double vol = length_a * length_b * length_c +
+    //       2.0 * thick_rim_a * length_b * length_c +
     //       2.0 * thick_rim_b * length_a * length_c +
     //       2.0 * thick_rim_c * length_a * length_b;
-    
+
     // Scale sides by B
     const double a_scaled = length_a / length_b;
@@ -101,6 +101,6 @@
             //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors
             // This is the function called by csparallelepiped_analytical_2D_scaled,
-            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c        
-            
+            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c
+
             //  correct FF : sum of square of phase factors
             inner_total += Gauss76Wt[j] * form * form;
@@ -136,4 +136,7 @@
     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;
 
     // cspkernel in csparallelepiped recoded here
@@ -160,11 +163,11 @@
     double tc = length_a + 2.0*thick_rim_c;
     //handle arg=0 separately, as sin(t)/t -> 1 as t->0
-    double siA = sas_sinx_x(0.5*q*length_a*xhat);
-    double siB = sas_sinx_x(0.5*q*length_b*yhat);
-    double siC = sas_sinx_x(0.5*q*length_c*zhat);
-    double siAt = sas_sinx_x(0.5*q*ta*xhat);
-    double siBt = sas_sinx_x(0.5*q*tb*yhat);
-    double siCt = sas_sinx_x(0.5*q*tc*zhat);
-    
+    double siA = sas_sinx_x(0.5*length_a*qa);
+    double siB = sas_sinx_x(0.5*length_b*qb);
+    double siC = sas_sinx_x(0.5*length_c*qc);
+    double siAt = sas_sinx_x(0.5*ta*qa);
+    double siBt = sas_sinx_x(0.5*tb*qb);
+    double siCt = sas_sinx_x(0.5*tc*qc);
+
 
     // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed
@@ -174,5 +177,5 @@
                + drB*siA*(siBt-siB)*siC*V2
                + drC*siA*siB*(siCt*siCt-siC)*V3);
-   
+
     return 1.0e-4 * f * f;
 }
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
+++ sasmodels/models/cylinder.c	(revision b34fc77253f928aabc59bd17cc4a0132796579d2)
@@ -1,37 +1,31 @@
-double form_volume(double radius, double length);
-double fq(double q, double sn, double cn,double radius, double length);
-double orient_avg_1D(double q, double radius, double length);
-double Iq(double q, double sld, double solvent_sld, double radius, double length);
-double Iqxy(double qx, double qy, double sld, double solvent_sld,
-    double radius, double length, double theta, double phi);
-
 #define INVALID(v) (v.radius<0 || v.length<0)
 
-double form_volume(double radius, double length)
+static double
+form_volume(double radius, double length)
 {
     return M_PI*radius*radius*length;
 }
 
-double fq(double q, double sn, double cn, double radius, double length)
+static double
+fq(double qab, double qc, double radius, double length)
 {
-    // precompute qr and qh to save time in the loop
-    const double qr = q*radius;
-    const double qh = q*0.5*length; 
-    return sas_2J1x_x(qr*sn) * sas_sinx_x(qh*cn);
+    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length);
 }
 
-double orient_avg_1D(double q, double radius, double length)
+static double
+orient_avg_1D(double q, double radius, double length)
 {
     // translate a point in [-1,1] to a point in [0, pi/2]
     const double zm = M_PI_4;
-    const double zb = M_PI_4; 
+    const double zb = M_PI_4;
 
     double total = 0.0;
     for (int i=0; i<76 ;i++) {
-        const double alpha = Gauss76Z[i]*zm + zb;
-        double sn, cn; // slots to hold sincos function output
-        // alpha(theta,phi) the projection of the cylinder on the detector plane
-        SINCOS(alpha, sn, cn);
-        total += Gauss76Wt[i] * square( fq(q, sn, cn, radius, length) ) * sn;
+        const double theta = Gauss76Z[i]*zm + zb;
+        double sin_theta, cos_theta; // slots to hold sincos function output
+        // theta (theta,phi) the projection of the cylinder on the detector plane
+        SINCOS(theta , sin_theta, cos_theta);
+        const double form = fq(q*sin_theta, q*cos_theta, radius, length);
+        total += Gauss76Wt[i] * form * form * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -39,5 +33,6 @@
 }
 
-double Iq(double q,
+static double
+Iq(double q,
     double sld,
     double solvent_sld,
@@ -49,6 +44,6 @@
 }
 
-
-double Iqxy(double qx, double qy,
+static double
+Iqxy(double qx, double qy,
     double sld,
     double solvent_sld,
@@ -60,7 +55,9 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
-    //printf("sn: %g cn: %g\n", sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
+
     const double s = (sld-solvent_sld) * form_volume(radius, length);
-    const double form = fq(q, sin_alpha, cos_alpha, radius, length);
+    const double form = fq(qab, qc, radius, length);
     return 1.0e-4 * square(s * form);
 }
Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision 3b571ae4950c5a99059d2ce8c74e80d0eb908a8f)
+++ sasmodels/models/ellipsoid.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,13 +1,10 @@
-double form_volume(double radius_polar, double radius_equatorial);
-double Iq(double q, double sld, double sld_solvent, double radius_polar, double radius_equatorial);
-double Iqxy(double qx, double qy, double sld, double sld_solvent,
-    double radius_polar, double radius_equatorial, double theta, double phi);
-
-double form_volume(double radius_polar, double radius_equatorial)
+static double
+form_volume(double radius_polar, double radius_equatorial)
 {
     return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial;
 }
 
-double Iq(double q,
+static  double
+Iq(double q,
     double sld,
     double sld_solvent,
@@ -41,5 +38,6 @@
 }
 
-double Iqxy(double qx, double qy,
+static double
+Iqxy(double qx, double qy,
     double sld,
     double sld_solvent,
@@ -51,10 +49,11 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
-    const double r = sqrt(square(radius_equatorial*sin_alpha)
-                          + square(radius_polar*cos_alpha));
-    const double f = sas_3j1x_x(q*r);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
+
+    const double qr = sqrt(square(radius_equatorial*qab) + square(radius_polar*qc));
+    const double f = sas_3j1x_x(qr);
     const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial);
 
     return 1.0e-4 * square(f * s);
 }
-
Index: sasmodels/models/elliptical_cylinder.c
===================================================================
--- sasmodels/models/elliptical_cylinder.c	(revision 61104c829b100869b42c1afd75f30ed30a934cd1)
+++ sasmodels/models/elliptical_cylinder.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,10 +1,3 @@
-double form_volume(double radius_minor, double r_ratio, double length);
-double Iq(double q, double radius_minor, double r_ratio, double length,
-          double sld, double solvent_sld);
-double Iqxy(double qx, double qy, double radius_minor, double r_ratio, double length,
-            double sld, double solvent_sld, double theta, double phi, double psi);
-
-
-double
+static double
 form_volume(double radius_minor, double r_ratio, double length)
 {
@@ -12,5 +5,5 @@
 }
 
-double
+static double
 Iq(double q, double radius_minor, double r_ratio, double length,
    double sld, double solvent_sld)
@@ -61,5 +54,5 @@
 
 
-double
+static double
 Iqxy(double qx, double qy,
      double radius_minor, double r_ratio, double length,
@@ -69,10 +62,13 @@
     double q, xhat, yhat, zhat;
     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;
 
     // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2)
     // Given:    radius_major = r_ratio * radius_minor
-    const double r = radius_minor*sqrt(square(r_ratio*xhat) + square(yhat));
-    const double be = sas_2J1x_x(q*r);
-    const double si = sas_sinx_x(q*zhat*0.5*length);
+    const double qr = radius_minor*sqrt(square(r_ratio*qa) + square(qb));
+    const double be = sas_2J1x_x(qr);
+    const double si = sas_sinx_x(qc*0.5*length);
     const double Aq = be * si;
     const double delrho = sld - solvent_sld;
Index: sasmodels/models/fcc_paracrystal.c
===================================================================
--- sasmodels/models/fcc_paracrystal.c	(revision 50beefeac18f2f952bf0863526ec9bbd3009e29c)
+++ sasmodels/models/fcc_paracrystal.c	(revision 7e0b281e5f175a91044364982fcfbfe802508948)
@@ -1,89 +1,86 @@
-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);
+static double
+fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
+{
+#if 0  // Equations as written in Matsuoka
+    const double a1 = ( qa + qb)/2.0;
+    const double a2 = (-qa + qc)/2.0;
+    const double a3 = (-qa + qb)/2.0;
+#else
+    const double a1 = ( qa + qb)/2.0;
+    const double a2 = ( qa + qc)/2.0;
+    const double a3 = ( qb + qc)/2.0;
+#endif
 
-double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
-double _FCCeval(double Theta, double Phi, double temp1, double temp3);
+    // Numerator: (1 - exp(a)^2)^3
+    //         => (-(exp(2a) - 1))^3
+    //         => -expm1(2a)^3
+    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a))
+    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1)
+    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1)
+    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    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;
+}
 
 
-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);
+// occupied volume fraction calculated from lattice symmetry and sphere radius
+static double
+fcc_volume_fraction(double radius, double dnn)
+{
+    return 4.0*sphere_volume(M_SQRT1_2*radius/dnn);
 }
 
-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 = cube(1.0-(temp10*temp10))*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){
+static double
+form_volume(double radius)
+{
     return sphere_volume(radius);
 }
 
 
-double Iq(double q, double dnn,
+static double Iq(double q, double dnn,
   double d_factor, double radius,
-  double sld, double solvent_sld){
+  double sld, double solvent_sld)
+{
+    // translate a point in [-1,1] to a point in [0, 2 pi]
+    const double phi_m = M_PI;
+    const double phi_b = M_PI;
+    // translate a point in [-1,1] to a point in [0, pi]
+    const double theta_m = M_PI_2;
+    const double theta_b = M_PI_2;
 
-	//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);
+    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 form = fcc_Zq(qa, qb, qc, dnn, d_factor);
+            inner_sum += Gauss150Wt[j] * form;
+        }
+        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
+        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+    }
+    outer_sum *= theta_m;
+    const double Zq = outer_sum/(4.0*M_PI);
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
 
-    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;
+    return fcc_volume_fraction(radius, dnn) * Pq * Zq;
+}
 
 
-}
-
-double Iqxy(double qx, double qy,
+static double Iqxy(double qx, double qy,
     double dnn, double d_factor, double radius,
     double sld, double solvent_sld,
@@ -92,21 +89,11 @@
     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;
 
-    const double a1 = yhat + xhat;
-    const double a2 = xhat + zhat;
-    const double a3 = yhat + zhat;
-    const double qd = 0.5*q*dnn;
-    const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3);
-    const double tanh_qd = tanh(arg);
-    const double cosh_qd = cosh(arg);
-    const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*a2)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*a3)/cosh_qd);
-
-    //if (isnan(Zq)) printf("q:(%g,%g) qd: %g a1: %g a2: %g a3: %g arg: %g\n", qx, qy, qd, a1, a2, a3, arg);
-
-    const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq;
-    //the occupied volume of the lattice
-    const double lattice_scale = 4.0*sphere_volume(M_SQRT1_2*radius/dnn);
-    return lattice_scale * Fq;
+    q = sqrt(qa*qa + qb*qb + qc*qc);
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+    const double Zq = fcc_Zq(qa, qb, qc, dnn, d_factor);
+    return fcc_volume_fraction(radius, dnn) * Pq * Zq;
 }
Index: sasmodels/models/hollow_cylinder.c
===================================================================
--- sasmodels/models/hollow_cylinder.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
+++ sasmodels/models/hollow_cylinder.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,7 +1,7 @@
 double form_volume(double radius, double thickness, double length);
 double Iq(double q, double radius, double thickness, double length, double sld,
-	double solvent_sld);
+    double solvent_sld);
 double Iqxy(double qx, double qy, double radius, double thickness, double length, double sld,
-	double solvent_sld, double theta, double phi);
+    double solvent_sld, double theta, double phi);
 
 //#define INVALID(v) (v.radius_core >= v.radius)
@@ -14,17 +14,15 @@
 }
 
-
 static double
-_hollow_cylinder_kernel(double q,
-    double radius, double thickness, double length, double sin_val, double cos_val)
+_fq(double qab, double qc,
+    double radius, double thickness, double length)
 {
-    const double qs = q*sin_val;
-    const double lam1 = sas_2J1x_x((radius+thickness)*qs);
-    const double lam2 = sas_2J1x_x(radius*qs);
+    const double lam1 = sas_2J1x_x((radius+thickness)*qab);
+    const double lam2 = sas_2J1x_x(radius*qab);
     const double gamma_sq = square(radius/(radius+thickness));
-    //Note: lim_{thickness -> 0} psi = sas_J0(radius*qs)
-    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qs)
-    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq);	//SRK 10/19/00
-    const double t2 = sas_sinx_x(0.5*q*length*cos_val);
+    //Note: lim_{thickness -> 0} psi = sas_J0(radius*qab)
+    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab)
+    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq);    //SRK 10/19/00
+    const double t2 = sas_sinx_x(0.5*length*qc);
     return psi*t2;
 }
@@ -43,13 +41,13 @@
 {
     const double lower = 0.0;
-    const double upper = 1.0;		//limits of numerical integral
+    const double upper = 1.0;        //limits of numerical integral
 
-    double summ = 0.0;			//initialize intergral
+    double summ = 0.0;            //initialize intergral
     for (int i=0;i<76;i++) {
-        const double cos_val = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper );
-        const double sin_val = sqrt(1.0 - cos_val*cos_val);
-        const double inter = _hollow_cylinder_kernel(q, radius, thickness, length,
-                                                     sin_val, cos_val);
-        summ += Gauss76Wt[i] * inter * inter;
+        const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper );
+        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
+        const double form = _fq(q*sin_theta, q*cos_theta,
+                                radius, thickness, length);
+        summ += Gauss76Wt[i] * form * form;
     }
 
@@ -66,9 +64,10 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
-    const double Aq = _hollow_cylinder_kernel(q, radius, thickness, length,
-        sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
+
+    const double form = _fq(qab, qc, radius, thickness, length);
 
     const double vol = form_volume(radius, thickness, length);
-    return _hollow_cylinder_scaling(Aq*Aq, solvent_sld-sld, vol);
+    return _hollow_cylinder_scaling(form*form, solvent_sld-sld, vol);
 }
-
Index: sasmodels/models/parallelepiped.c
===================================================================
--- sasmodels/models/parallelepiped.c	(revision d6050806338b4bcae5a997acb08f4b1e7e14dd36)
+++ sasmodels/models/parallelepiped.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -20,9 +20,9 @@
 {
     const double mu = 0.5 * q * length_b;
-    
+
     // Scale sides by B
     const double a_scaled = length_a / length_b;
     const double c_scaled = length_c / length_b;
-        
+
     // outer integral (with gauss points), integration limits = 0, 1
     double outer_total = 0; //initialize integral
@@ -69,8 +69,11 @@
     double q, xhat, yhat, zhat;
     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;
 
-    const double siA = sas_sinx_x(0.5*length_a*q*xhat);
-    const double siB = sas_sinx_x(0.5*length_b*q*yhat);
-    const double siC = sas_sinx_x(0.5*length_c*q*zhat);
+    const double siA = sas_sinx_x(0.5*length_a*qa);
+    const double siB = sas_sinx_x(0.5*length_b*qb);
+    const double siC = sas_sinx_x(0.5*length_c*qc);
     const double V = form_volume(length_a, length_b, length_c);
     const double drho = (sld - solvent_sld);
Index: sasmodels/models/sc_paracrystal.c
===================================================================
--- sasmodels/models/sc_paracrystal.c	(revision 50beefeac18f2f952bf0863526ec9bbd3009e29c)
+++ sasmodels/models/sc_paracrystal.c	(revision 7e0b281e5f175a91044364982fcfbfe802508948)
@@ -1,128 +1,93 @@
-double form_volume(double radius);
+static double
+sc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
+{
+    const double a1 = qa;
+    const double a2 = qb;
+    const double a3 = qc;
 
-double Iq(double q,
-          double dnn,
-          double d_factor,
-          double radius,
-          double sphere_sld,
-          double solvent_sld);
+    // Numerator: (1 - exp(a)^2)^3
+    //         => (-(exp(2a) - 1))^3
+    //         => -expm1(2a)^3
+    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a))
+    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1)
+    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1)
+    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    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));
 
-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);
+    return Zq;
+}
 
-double form_volume(double radius)
+// 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
-sc_eval(double theta, double phi, double temp3, double temp4, double temp5)
+
+static double Iq(double q, double dnn,
+  double d_factor, double radius,
+  double sld, double solvent_sld)
 {
-    double cnt, snt;
-    SINCOS(theta, cnt, snt);
+    // translate a point in [-1,1] to a point in [0, 2 pi]
+    const double phi_m = M_PI_4;
+    const double phi_b = M_PI_4;
+    // translate a point in [-1,1] to a point in [0, pi]
+    const double theta_m = M_PI_4;
+    const double theta_b = M_PI_4;
 
-    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);
+    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 form = sc_Zq(qa, qb, qc, dnn, d_factor);
+            inner_sum += Gauss150Wt[j] * form;
+        }
+        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
+        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+    }
+    outer_sum *= theta_m;
+    const double Zq = outer_sum/M_PI_2;
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+
+    return sc_volume_fraction(radius, dnn) * Pq * Zq;
 }
 
-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 = cube(-expm1(-temp1));
-	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)/M_PI_2;
-
-	return(integrand);
-}
-
-double Iq(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; //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;
-}
-
-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)
+static 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 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;
 
-    const double qd = q*dnn;
-    const double arg = 0.5*square(qd*d_factor);
-    const double tanh_qd = tanh(arg);
-    const double cosh_qd = cosh(arg);
-    const double Zq = tanh_qd/(1. - cos(qd*zhat)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*yhat)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*xhat)/cosh_qd);
-
-    const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq;
-    //the occupied volume of the lattice
-    const double lattice_scale = sphere_volume(radius/dnn);
-    return lattice_scale * Fq;
+    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 sc_volume_fraction(radius, dnn) * Pq * Zq;
 }
Index: sasmodels/models/sc_paracrystal.py
===================================================================
--- sasmodels/models/sc_paracrystal.py	(revision 69e1afc57270f67b2a0d7157a6d263fc2db98cc0)
+++ sasmodels/models/sc_paracrystal.py	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -152,8 +152,6 @@
     [{}, 0.001, 10.3048],
     [{}, 0.215268, 0.00814889],
-    [{}, (0.414467), 0.001313289],
+    [{}, 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 ]
     ]
-
-
Index: sasmodels/models/stacked_disks.c
===================================================================
--- sasmodels/models/stacked_disks.c	(revision 19f996bbf3615f0161b71f00777e5f16155bef70)
+++ sasmodels/models/stacked_disks.c	(revision b34fc77253f928aabc59bd17cc4a0132796579d2)
@@ -1,4 +1,5 @@
 static double stacked_disks_kernel(
-    double q,
+    double qab,
+    double qc,
     double halfheight,
     double thick_layer,
@@ -9,6 +10,4 @@
     double layer_sld,
     double solvent_sld,
-    double sin_alpha,
-    double cos_alpha,
     double d)
 
@@ -20,9 +19,9 @@
     // zi is the dummy variable for the integration (x in Feigin's notation)
 
-    const double besarg1 = q*radius*sin_alpha;
-    //const double besarg2 = q*radius*sin_alpha;
+    const double besarg1 = radius*qab;
+    //const double besarg2 = radius*qab;
 
-    const double sinarg1 = q*halfheight*cos_alpha;
-    const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha;
+    const double sinarg1 = halfheight*qc;
+    const double sinarg2 = (halfheight+thick_layer)*qc;
 
     const double be1 = sas_2J1x_x(besarg1);
@@ -43,5 +42,5 @@
 
     // loop for the structure factor S(q)
-    double qd_cos_alpha = q*d*cos_alpha;
+    double qd_cos_alpha = d*qc;
     //d*cos_alpha is the projection of d onto q (in other words the component
     //of d that is parallel to q.
@@ -84,5 +83,5 @@
         double sin_alpha, cos_alpha; // slots to hold sincos function output
         SINCOS(zi, sin_alpha, cos_alpha);
-        double yyy = stacked_disks_kernel(q,
+        double yyy = stacked_disks_kernel(q*sin_alpha, q*cos_alpha,
                            halfheight,
                            thick_layer,
@@ -93,6 +92,4 @@
                            layer_sld,
                            solvent_sld,
-                           sin_alpha,
-                           cos_alpha,
                            d);
         summ += Gauss76Wt[i] * yyy * sin_alpha;
@@ -152,11 +149,13 @@
     double phi)
 {
-    int n_stacking = (int)(fp_n_stacking + 0.5);
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
+    int n_stacking = (int)(fp_n_stacking + 0.5);
     double d = 2.0 * thick_layer + thick_core;
     double halfheight = 0.5*thick_core;
-    double answer = stacked_disks_kernel(q,
+    double answer = stacked_disks_kernel(qab, qc,
                      halfheight,
                      thick_layer,
@@ -167,6 +166,4 @@
                      layer_sld,
                      solvent_sld,
-                     sin_alpha,
-                     cos_alpha,
                      d);
 
Index: sasmodels/models/triaxial_ellipsoid.c
===================================================================
--- sasmodels/models/triaxial_ellipsoid.c	(revision 68dd6a95741269738dcbfb051d93fe9f26b89c43)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
@@ -1,17 +1,12 @@
-double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar);
-double Iq(double q, double sld, double sld_solvent,
-    double radius_equat_minor, double radius_equat_major, double radius_polar);
-double Iqxy(double qx, double qy, double sld, double sld_solvent,
-    double radius_equat_minor, double radius_equat_major, double radius_polar, double theta, double phi, double psi);
-
 //#define INVALID(v) (v.radius_equat_minor > v.radius_equat_major || v.radius_equat_major > v.radius_polar)
 
-
-double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar)
+static double
+form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar)
 {
     return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar;
 }
 
-double Iq(double q,
+static double
+Iq(double q,
     double sld,
     double sld_solvent,
@@ -45,9 +40,11 @@
     // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi
     const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI);
-    const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar);
-    return 1.0e-4 * s * s * fqsq;
+    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar);
+    const double drho = (sld - sld_solvent);
+    return 1.0e-4 * square(vol*drho) * fqsq;
 }
 
-double Iqxy(double qx, double qy,
+static double
+Iqxy(double qx, double qy,
     double sld,
     double sld_solvent,
@@ -61,12 +58,15 @@
     double q, xhat, yhat, zhat;
     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;
 
-    const double r = sqrt(square(radius_equat_minor*xhat)
-                          + square(radius_equat_major*yhat)
-                          + square(radius_polar*zhat));
-    const double fq = sas_3j1x_x(q*r);
-    const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar);
+    const double qr = sqrt(square(radius_equat_minor*qa)
+                           + square(radius_equat_major*qb)
+                           + square(radius_polar*qc));
+    const double fq = sas_3j1x_x(qr);
+    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar);
+    const double drho = (sld - sld_solvent);
 
-    return 1.0e-4 * square(s * fq);
+    return 1.0e-4 * square(vol * drho * fq);
 }
-
Index: sasmodels/special.py
===================================================================
--- sasmodels/special.py	(revision 859567e703d85ffc9cd3c2aea784e5ef9dfdc901)
+++ sasmodels/special.py	(revision 859567e703d85ffc9cd3c2aea784e5ef9dfdc901)
@@ -0,0 +1,804 @@
+"""
+Special Functions
+.................
+
+The C code follows the C99 standard, with the usual math functions,
+as defined in
+`OpenCL <https://www.khronos.org/registry/cl/sdk/1.1/docs/man/xhtml/mathFunctions.html>`_.
+This includes the following:
+
+    M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E:
+        $\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler's constant $e$
+    exp, log, pow(x,y), expm1, sqrt:
+        Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\sqrt{x}$.
+        The function expm1(x) is accurate across all $x$, including $x$
+        very close to zero.
+    sin, cos, tan, asin, acos, atan:
+        Trigonometry functions and inverses, operating on radians.
+    sinh, cosh, tanh, asinh, acosh, atanh:
+        Hyperbolic trigonometry functions.
+    atan2(y,x):
+        Angle from the $x$\ -axis to the point $(x,y)$, which is equal to
+        $\tan^{-1}(y/x)$ corrected for quadrant.  That is, if $x$ and $y$ are
+        both negative, then atan2(y,x) returns a value in quadrant III where
+        atan(y/x) would return a value in quadrant I. Similarly for
+        quadrants II and IV when $x$ and $y$ have opposite sign.
+    fmin(x,y), fmax(x,y), trunc, rint:
+        Floating point functions.  rint(x) returns the nearest integer.
+    NAN:
+        NaN, Not a Number, $0/0$.  Use isnan(x) to test for NaN.  Note that
+        you cannot use :code:`x == NAN` to test for NaN values since that
+        will always return false.  NAN does not equal NAN!
+    INFINITY:
+        $\infty, 1/0$.  Use isinf(x) to test for infinity, or isfinite(x)
+        to test for finite and not NaN.
+    erf, erfc, tgamma, lgamma:  **do not use**
+        Special functions that should be part of the standard, but are missing
+        or inaccurate on some platforms. Use sas_erf, sas_erfc and sas_gamma
+        instead (see below). Note: lgamma(x) has not yet been tested.
+
+Some non-standard constants and functions are also provided:
+
+    M_PI_180, M_4PI_3:
+        $\frac{\pi}{180}$, $\frac{4\pi}{3}$
+    SINCOS(x, s, c):
+        Macro which sets s=sin(x) and c=cos(x). The variables *c* and *s*
+        must be declared first.
+    square(x):
+        $x^2$
+    cube(x):
+        $x^3$
+    sas_sinx_x(x):
+        $\sin(x)/x$, with limit $\sin(0)/0 = 1$.
+    powr(x, y):
+        $x^y$ for $x \ge 0$; this is faster than general $x^y$ on some GPUs.
+    pown(x, n):
+        $x^n$ for $n$ integer; this is faster than general $x^n$ on some GPUs.
+    FLOAT_SIZE:
+        The number of bytes in a floating point value.  Even though all
+        variables are declared double, they may be converted to single
+        precision float before running. If your algorithm depends on
+        precision (which is not uncommon for numerical algorithms), use
+        the following::
+
+            #if FLOAT_SIZE>4
+            ... code for double precision ...
+            #else
+            ... code for single precision ...
+            #endif
+    SAS_DOUBLE:
+        A replacement for :code:`double` so that the declared variable will
+        stay double precision; this should generally not be used since some
+        graphics cards do not support double precision.  There is no provision
+        for forcing a constant to stay double precision.
+
+The following special functions and scattering calculations are defined in
+`sasmodels/models/lib <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib>`_.
+These functions have been tuned to be fast and numerically stable down
+to $q=0$ even in single precision.  In some cases they work around bugs
+which appear on some platforms but not others, so use them where needed.
+Add the files listed in :code:`source = ["lib/file.c", ...]` to your *model.py*
+file in the order given, otherwise these functions will not be available.
+
+    polevl(x, c, n):
+        Polynomial evaluation $p(x) = \sum_{i=0}^n c_i x^i$ using Horner's
+        method so it is faster and more accurate.
+
+        $c = \{c_n, c_{n-1}, \ldots, c_0 \}$ is the table of coefficients,
+        sorted from highest to lowest.
+
+        :code:`source = ["lib/polevl.c", ...]` (`link to code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/polevl.c>`_)
+
+    p1evl(x, c, n):
+        Evaluation of normalized polynomial $p(x) = x^n + \sum_{i=0}^{n-1} c_i x^i$
+        using Horner's method so it is faster and more accurate.
+
+        $c = \{c_{n-1}, c_{n-2} \ldots, c_0 \}$ is the table of coefficients,
+        sorted from highest to lowest.
+
+        :code:`source = ["lib/polevl.c", ...]`
+        (`link to code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/polevl.c>`_)
+
+    sas_gamma(x):
+        Gamma function $\text{sas_gamma}(x) = \Gamma(x)$.
+
+        The standard math function, tgamma(x) is unstable for $x < 1$
+        on some platforms.
+
+        :code:`source = ["lib/sasgamma.c", ...]`
+        (`link to code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gamma.c>`_)
+
+    sas_erf(x), sas_erfc(x):
+        Error function
+        $\text{sas_erf}(x) = \frac{2}{\sqrt\pi}\int_0^x e^{-t^2}\,dt$
+        and complementary error function
+        $\text{sas_erfc}(x) = \frac{2}{\sqrt\pi}\int_x^{\infty} e^{-t^2}\,dt$.
+
+        The standard math functions erf(x) and erfc(x) are slower and broken
+        on some platforms.
+
+        :code:`source = ["lib/polevl.c", "lib/sas_erf.c", ...]`
+        (`link to error functions' code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_erf.c>`_)
+
+    sas_J0(x):
+        Bessel function of the first kind $\text{sas_J0}(x)=J_0(x)$ where
+        $J_0(x) = \frac{1}{\pi}\int_0^\pi \cos(x\sin(\tau))\,d\tau$.
+
+        The standard math function j0(x) is not available on all platforms.
+
+        :code:`source = ["lib/polevl.c", "lib/sas_J0.c", ...]`
+        (`link to Bessel function's code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_J0.c>`_)
+
+    sas_J1(x):
+        Bessel function of the first kind  $\text{sas_J1}(x)=J_1(x)$ where
+        $J_1(x) = \frac{1}{\pi}\int_0^\pi \cos(\tau - x\sin(\tau))\,d\tau$.
+
+        The standard math function j1(x) is not available on all platforms.
+
+        :code:`source = ["lib/polevl.c", "lib/sas_J1.c", ...]`
+        (`link to Bessel function's code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_J1.c>`_)
+
+    sas_JN(n, x):
+        Bessel function of the first kind and integer order $n$:
+        $\text{sas_JN}(n, x)=J_n(x)$ where
+        $J_n(x) = \frac{1}{\pi}\int_0^\pi \cos(n\tau - x\sin(\tau))\,d\tau$.
+        If $n$ = 0 or 1, it uses sas_J0(x) or sas_J1(x), respectively.
+
+        The standard math function jn(n, x) is not available on all platforms.
+
+        :code:`source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c", ...]`
+        (`link to Bessel function's code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_JN.c>`_)
+
+    sas_Si(x):
+        Sine integral $\text{Si}(x) = \int_0^x \tfrac{\sin t}{t}\,dt$.
+
+        This function uses Taylor series for small and large arguments:
+
+        For large arguments,
+
+        .. math::
+
+             \text{Si}(x) \sim \frac{\pi}{2}
+             - \frac{\cos(x)}{x}\left(1 - \frac{2!}{x^2} + \frac{4!}{x^4} - \frac{6!}{x^6} \right)
+             - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right)
+
+        For small arguments,
+
+        .. math::
+
+           \text{Si}(x) \sim x
+           - \frac{x^3}{3\times 3!} + \frac{x^5}{5 \times 5!} - \frac{x^7}{7 \times 7!}
+           + \frac{x^9}{9\times 9!} - \frac{x^{11}}{11\times 11!}
+
+        :code:`source = ["lib/Si.c", ...]`
+        (`link to code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/Si.c>`_)
+
+    sas_3j1x_x(x):
+        Spherical Bessel form
+        $\text{sph_j1c}(x) = 3 j_1(x)/x = 3 (\sin(x) - x \cos(x))/x^3$,
+        with a limiting value of 1 at $x=0$, where $j_1(x)$ is the spherical
+        Bessel function of the first kind and first order.
+
+        This function uses a Taylor series for small $x$ for numerical accuracy.
+
+        :code:`source = ["lib/sas_3j1x_x.c", ...]`
+        (`link to code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_3j1x_x.c>`_)
+
+
+    sas_2J1x_x(x):
+        Bessel form $\text{sas_J1c}(x) = 2 J_1(x)/x$, with a limiting value
+        of 1 at $x=0$, where $J_1(x)$ is the Bessel function of first kind
+        and first order.
+
+        :code:`source = ["lib/polevl.c", "lib/sas_J1.c", ...]`
+        (`link to Bessel form's code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_J1.c>`_)
+
+
+    Gauss76Z[i], Gauss76Wt[i]:
+        Points $z_i$ and weights $w_i$ for 76-point Gaussian quadrature, respectively,
+        computing $\int_{-1}^1 f(z)\,dz \approx \sum_{i=1}^{76} w_i\,f(z_i)$.
+
+        Similar arrays are available in :code:`gauss20.c` for 20-point
+        quadrature and in :code:`gauss150.c` for 150-point quadrature.
+
+        :code:`source = ["lib/gauss76.c", ...]`
+        (`link to code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/gauss76.c>`_)
+
+
+"""
+import numpy as np
+import scipy.special
+
+# Functions to add to our standard set
+from numpy import degrees, radians
+
+# C99 standard math library functions
+M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E = np.pi, np.pi/2, np.pi/4, np.sqrt(0.5), np.e
+from numpy import exp, log, power as pow, expm1, sqrt
+from numpy import sin, cos, tan, arcsin as asin, arccos as acos, arctan as atan
+from numpy import sinh, cosh, tanh, arcsinh as asinh, arccosh as acosh, arctanh as atanh
+from numpy import arctan2 as atan2
+from numpy import fmin, fmax, trunc, rint
+from numpy import NAN, inf as INFINITY
+
+# erf, erfc, tgamma, lgamma  **do not use**
+
+# Rotations of q
+def ORIENT_SYMMETRIC(qx, qy, theta, phi):
+    q = sqrt(qx*qx + qy*qy)
+    q = sqrt(qx*qx + qy*qy)
+    sin_phi, cos_phi = sin(radians(phi)), cos(radians(phi))
+    cn = (cn*qx + sn*qy)/q * sin(radians(theta))
+    cn[q==0.] = 1.
+    sn = sqrt(1 - cn**2)
+    return q, cn, sn
+
+def ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi):
+    q = sqrt(qx*qx + qy*qy)
+    qxhat = qx/q
+    qyhat = qy/q
+    sin_theta, cos_theta = sin(radians(theta)), cos(radians(theta))
+    sin_phi, cos_phi = sin(radians(phi)), cos(radians(phi))
+    sin_psi, cos_psi = sin(radians(psi)), cos(radians(psi))
+
+    xhat = (qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi)
+          + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi))
+    yhat = (qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi)
+          + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi))
+    zhat = (qxhat*(-sin_theta*cos_phi)
+          + qyhat*(-sin_theta*sin_phi))
+    return q, xhat, yhat, zhat
+
+
+# non-standard constants and functions
+
+M_PI_180, M_4PI_3 = M_PI/180, 4*M_PI/3
+
+# can't do SINCOS in python; use "s, c = SINCOS(x)" instead
+def SINCOS(x): return sin(x), cos(x)
+
+def square(x): return x*x
+
+def cube(x): return x*x*x
+
+from numpy import sinc as _sinc
+def sas_sinx_x(x): return _sinc(x/M_PI)
+def powr(x, y): return x**y
+def pown(x, n): return x**n
+
+FLOAT_SIZE = 8
+
+def polevl(x, c, n): return np.polyval(c[:n], x)
+def p1evl(x, c, n): return np.polyval(np.hstack(([1.], c))[:n], x)
+
+from scipy.special import gamma as sas_gamma
+from scipy.special import erf as sas_erf
+from scipy.special import erfc as sas_erfc
+from scipy.special import j0 as sas_J0
+from scipy.special import j1 as sas_J1
+from scipy.special import jn as sas_JN
+
+# missing sas_Si
+def sas_Si(x):
+    return scipy.special.sici(x)[0]
+
+def sas_3j1x_x(x):
+    retvalue = 3*sas_j1(x)/x
+    retvalue[x==0] = 1.
+
+def sas_2J1x_x(x):
+    retvalue = 2*sas_J1(x)/x
+    retvalue[x==0] = 1.
+
+
+# Gaussians
+
+Gauss20Wt = np.array([
+	.0176140071391521,
+	.0406014298003869,
+	.0626720483341091,
+	.0832767415767047,
+	.10193011981724,
+	.118194531961518,
+	.131688638449177,
+	.142096109318382,
+	.149172986472604,
+	.152753387130726,
+	.152753387130726,
+	.149172986472604,
+	.142096109318382,
+	.131688638449177,
+	.118194531961518,
+	.10193011981724,
+	.0832767415767047,
+	.0626720483341091,
+	.0406014298003869,
+	.0176140071391521
+])
+
+Gauss20Z = np.array([
+	-.993128599185095,
+	-.963971927277914,
+	-.912234428251326,
+	-.839116971822219,
+	-.746331906460151,
+	-.636053680726515,
+	-.510867001950827,
+	-.37370608871542,
+	-.227785851141645,
+	-.076526521133497,
+	.0765265211334973,
+	.227785851141645,
+	.37370608871542,
+	.510867001950827,
+	.636053680726515,
+	.746331906460151,
+	.839116971822219,
+	.912234428251326,
+	.963971927277914,
+	.993128599185095
+])
+
+Gauss76Wt = np.array([
+	.00126779163408536,		#0
+	.00294910295364247,
+	.00462793522803742,
+	.00629918049732845,
+	.00795984747723973,
+	.00960710541471375,
+	.0112381685696677,
+	.0128502838475101,
+	.0144407317482767,
+	.0160068299122486,
+	.0175459372914742,		#10
+	.0190554584671906,
+	.020532847967908,
+	.0219756145344162,
+	.0233813253070112,
+	.0247476099206597,
+	.026072164497986,
+	.0273527555318275,
+	.028587223650054,
+	.029773487255905,
+	.0309095460374916,		#20
+	.0319934843404216,
+	.0330234743977917,
+	.0339977794120564,
+	.0349147564835508,
+	.0357728593807139,
+	.0365706411473296,
+	.0373067565423816,
+	.0379799643084053,
+	.0385891292645067,
+	.0391332242205184,		#30
+	.0396113317090621,
+	.0400226455325968,
+	.040366472122844,
+	.0406422317102947,
+	.0408494593018285,
+	.040987805464794,
+	.0410570369162294,
+	.0410570369162294,
+	.040987805464794,
+	.0408494593018285,		#40
+	.0406422317102947,
+	.040366472122844,
+	.0400226455325968,
+	.0396113317090621,
+	.0391332242205184,
+	.0385891292645067,
+	.0379799643084053,
+	.0373067565423816,
+	.0365706411473296,
+	.0357728593807139,		#50
+	.0349147564835508,
+	.0339977794120564,
+	.0330234743977917,
+	.0319934843404216,
+	.0309095460374916,
+	.029773487255905,
+	.028587223650054,
+	.0273527555318275,
+	.026072164497986,
+	.0247476099206597,		#60
+	.0233813253070112,
+	.0219756145344162,
+	.020532847967908,
+	.0190554584671906,
+	.0175459372914742,
+	.0160068299122486,
+	.0144407317482767,
+	.0128502838475101,
+	.0112381685696677,
+	.00960710541471375,		#70
+	.00795984747723973,
+	.00629918049732845,
+	.00462793522803742,
+	.00294910295364247,
+	.00126779163408536		#75 (indexed from 0)
+])
+
+Gauss76Z = np.array([
+	-.999505948362153,		#0
+	-.997397786355355,
+	-.993608772723527,
+	-.988144453359837,
+	-.981013938975656,
+	-.972229228520377,
+	-.961805126758768,
+	-.949759207710896,
+	-.936111781934811,
+	-.92088586125215,
+	-.904107119545567,		#10
+	-.885803849292083,
+	-.866006913771982,
+	-.844749694983342,
+	-.822068037328975,
+	-.7980001871612,
+	-.77258672828181,
+	-.74587051350361,
+	-.717896592387704,
+	-.688712135277641,
+	-.658366353758143,		#20
+	-.626910417672267,
+	-.594397368836793,
+	-.560882031601237,
+	-.526420920401243,
+	-.491072144462194,
+	-.454895309813726,
+	-.417951418780327,
+	-.380302767117504,
+	-.342012838966962,
+	-.303146199807908,		#30
+	-.263768387584994,
+	-.223945802196474,
+	-.183745593528914,
+	-.143235548227268,
+	-.102483975391227,
+	-.0615595913906112,
+	-.0205314039939986,
+	.0205314039939986,
+	.0615595913906112,
+	.102483975391227,			#40
+	.143235548227268,
+	.183745593528914,
+	.223945802196474,
+	.263768387584994,
+	.303146199807908,
+	.342012838966962,
+	.380302767117504,
+	.417951418780327,
+	.454895309813726,
+	.491072144462194,		#50
+	.526420920401243,
+	.560882031601237,
+	.594397368836793,
+	.626910417672267,
+	.658366353758143,
+	.688712135277641,
+	.717896592387704,
+	.74587051350361,
+	.77258672828181,
+	.7980001871612,	#60
+	.822068037328975,
+	.844749694983342,
+	.866006913771982,
+	.885803849292083,
+	.904107119545567,
+	.92088586125215,
+	.936111781934811,
+	.949759207710896,
+	.961805126758768,
+	.972229228520377,		#70
+	.981013938975656,
+	.988144453359837,
+	.993608772723527,
+	.997397786355355,
+	.999505948362153		#75
+])
+
+Gauss150Z = np.array([
+  	-0.9998723404457334,
+  	-0.9993274305065947,
+  	-0.9983473449340834,
+  	-0.9969322929775997,
+  	-0.9950828645255290,
+  	-0.9927998590434373,
+  	-0.9900842691660192,
+  	-0.9869372772712794,
+  	-0.9833602541697529,
+  	-0.9793547582425894,
+  	-0.9749225346595943,
+  	-0.9700655145738374,
+  	-0.9647858142586956,
+  	-0.9590857341746905,
+  	-0.9529677579610971,
+  	-0.9464345513503147,
+  	-0.9394889610042837,
+  	-0.9321340132728527,
+  	-0.9243729128743136,
+  	-0.9162090414984952,
+  	-0.9076459563329236,
+  	-0.8986873885126239,
+  	-0.8893372414942055,
+  	-0.8795995893549102,
+  	-0.8694786750173527,
+  	-0.8589789084007133,
+  	-0.8481048644991847,
+  	-0.8368612813885015,
+  	-0.8252530581614230,
+  	-0.8132852527930605,
+  	-0.8009630799369827,
+  	-0.7882919086530552,
+  	-0.7752772600680049,
+  	-0.7619248049697269,
+  	-0.7482403613363824,
+  	-0.7342298918013638,
+  	-0.7198995010552305,
+  	-0.7052554331857488,
+  	-0.6903040689571928,
+  	-0.6750519230300931,
+  	-0.6595056411226444,
+  	-0.6436719971150083,
+  	-0.6275578900977726,
+  	-0.6111703413658551,
+  	-0.5945164913591590,
+  	-0.5776035965513142,
+  	-0.5604390262878617,
+  	-0.5430302595752546,
+  	-0.5253848818220803,
+  	-0.5075105815339176,
+  	-0.4894151469632753,
+  	-0.4711064627160663,
+  	-0.4525925063160997,
+  	-0.4338813447290861,
+  	-0.4149811308476706,
+  	-0.3959000999390257,
+  	-0.3766465660565522,
+  	-0.3572289184172501,
+  	-0.3376556177463400,
+  	-0.3179351925907259,
+  	-0.2980762356029071,
+  	-0.2780873997969574,
+  	-0.2579773947782034,
+  	-0.2377549829482451,
+  	-0.2174289756869712,
+  	-0.1970082295132342,
+  	-0.1765016422258567,
+  	-0.1559181490266516,
+  	-0.1352667186271445,
+  	-0.1145563493406956,
+  	-0.0937960651617229,
+  	-0.0729949118337358,
+  	-0.0521619529078925,
+  	-0.0313062657937972,
+  	-0.0104369378042598,
+  	0.0104369378042598,
+  	0.0313062657937972,
+  	0.0521619529078925,
+  	0.0729949118337358,
+  	0.0937960651617229,
+  	0.1145563493406956,
+  	0.1352667186271445,
+  	0.1559181490266516,
+  	0.1765016422258567,
+  	0.1970082295132342,
+  	0.2174289756869712,
+  	0.2377549829482451,
+  	0.2579773947782034,
+  	0.2780873997969574,
+  	0.2980762356029071,
+  	0.3179351925907259,
+  	0.3376556177463400,
+  	0.3572289184172501,
+  	0.3766465660565522,
+  	0.3959000999390257,
+  	0.4149811308476706,
+  	0.4338813447290861,
+  	0.4525925063160997,
+  	0.4711064627160663,
+  	0.4894151469632753,
+  	0.5075105815339176,
+  	0.5253848818220803,
+  	0.5430302595752546,
+  	0.5604390262878617,
+  	0.5776035965513142,
+  	0.5945164913591590,
+  	0.6111703413658551,
+  	0.6275578900977726,
+  	0.6436719971150083,
+  	0.6595056411226444,
+  	0.6750519230300931,
+  	0.6903040689571928,
+  	0.7052554331857488,
+  	0.7198995010552305,
+  	0.7342298918013638,
+  	0.7482403613363824,
+  	0.7619248049697269,
+  	0.7752772600680049,
+  	0.7882919086530552,
+  	0.8009630799369827,
+  	0.8132852527930605,
+  	0.8252530581614230,
+  	0.8368612813885015,
+  	0.8481048644991847,
+  	0.8589789084007133,
+  	0.8694786750173527,
+  	0.8795995893549102,
+  	0.8893372414942055,
+  	0.8986873885126239,
+  	0.9076459563329236,
+  	0.9162090414984952,
+  	0.9243729128743136,
+  	0.9321340132728527,
+  	0.9394889610042837,
+  	0.9464345513503147,
+  	0.9529677579610971,
+  	0.9590857341746905,
+  	0.9647858142586956,
+  	0.9700655145738374,
+  	0.9749225346595943,
+  	0.9793547582425894,
+  	0.9833602541697529,
+  	0.9869372772712794,
+  	0.9900842691660192,
+  	0.9927998590434373,
+  	0.9950828645255290,
+  	0.9969322929775997,
+  	0.9983473449340834,
+  	0.9993274305065947,
+  	0.9998723404457334
+])
+
+Gauss150Wt = np.array([
+  	0.0003276086705538,
+  	0.0007624720924706,
+  	0.0011976474864367,
+  	0.0016323569986067,
+  	0.0020663664924131,
+  	0.0024994789888943,
+  	0.0029315036836558,
+  	0.0033622516236779,
+  	0.0037915348363451,
+  	0.0042191661429919,
+  	0.0046449591497966,
+  	0.0050687282939456,
+  	0.0054902889094487,
+  	0.0059094573005900,
+  	0.0063260508184704,
+  	0.0067398879387430,
+  	0.0071507883396855,
+  	0.0075585729801782,
+  	0.0079630641773633,
+  	0.0083640856838475,
+  	0.0087614627643580,
+  	0.0091550222717888,
+  	0.0095445927225849,
+  	0.0099300043714212,
+  	0.0103110892851360,
+  	0.0106876814158841,
+  	0.0110596166734735,
+  	0.0114267329968529,
+  	0.0117888704247183,
+  	0.0121458711652067,
+  	0.0124975796646449,
+  	0.0128438426753249,
+  	0.0131845093222756,
+  	0.0135194311690004,
+  	0.0138484622795371,
+  	0.0141714592928592,
+  	0.0144882814685445,
+  	0.0147987907597169,
+  	0.0151028518701744,
+  	0.0154003323133401,
+  	0.0156911024699895,
+  	0.0159750356447283,
+  	0.0162520081211971,
+  	0.0165218992159766,
+  	0.0167845913311726,
+  	0.0170399700056559,
+  	0.0172879239649355,
+  	0.0175283451696437,
+  	0.0177611288626114,
+  	0.0179861736145128,
+  	0.0182033813680609,
+  	0.0184126574807331,
+  	0.0186139107660094,
+  	0.0188070535331042,
+  	0.0189920016251754,
+  	0.0191686744559934,
+  	0.0193369950450545,
+  	0.0194968900511231,
+  	0.0196482898041878,
+  	0.0197911283358190,
+  	0.0199253434079123,
+  	0.0200508765398072,
+  	0.0201676730337687,
+  	0.0202756819988200,
+  	0.0203748563729175,
+  	0.0204651529434560,
+  	0.0205465323660984,
+  	0.0206189591819181,
+  	0.0206824018328499,
+  	0.0207368326754401,
+  	0.0207822279928917,
+  	0.0208185680053983,
+  	0.0208458368787627,
+  	0.0208640227312962,
+  	0.0208731176389954,
+  	0.0208731176389954,
+  	0.0208640227312962,
+  	0.0208458368787627,
+  	0.0208185680053983,
+  	0.0207822279928917,
+  	0.0207368326754401,
+  	0.0206824018328499,
+  	0.0206189591819181,
+  	0.0205465323660984,
+  	0.0204651529434560,
+  	0.0203748563729175,
+  	0.0202756819988200,
+  	0.0201676730337687,
+  	0.0200508765398072,
+  	0.0199253434079123,
+  	0.0197911283358190,
+  	0.0196482898041878,
+  	0.0194968900511231,
+  	0.0193369950450545,
+  	0.0191686744559934,
+  	0.0189920016251754,
+  	0.0188070535331042,
+  	0.0186139107660094,
+  	0.0184126574807331,
+  	0.0182033813680609,
+  	0.0179861736145128,
+  	0.0177611288626114,
+  	0.0175283451696437,
+  	0.0172879239649355,
+  	0.0170399700056559,
+  	0.0167845913311726,
+  	0.0165218992159766,
+  	0.0162520081211971,
+  	0.0159750356447283,
+  	0.0156911024699895,
+  	0.0154003323133401,
+  	0.0151028518701744,
+  	0.0147987907597169,
+  	0.0144882814685445,
+  	0.0141714592928592,
+  	0.0138484622795371,
+  	0.0135194311690004,
+  	0.0131845093222756,
+  	0.0128438426753249,
+  	0.0124975796646449,
+  	0.0121458711652067,
+  	0.0117888704247183,
+  	0.0114267329968529,
+  	0.0110596166734735,
+  	0.0106876814158841,
+  	0.0103110892851360,
+  	0.0099300043714212,
+  	0.0095445927225849,
+  	0.0091550222717888,
+  	0.0087614627643580,
+  	0.0083640856838475,
+  	0.0079630641773633,
+  	0.0075585729801782,
+  	0.0071507883396855,
+  	0.0067398879387430,
+  	0.0063260508184704,
+  	0.0059094573005900,
+  	0.0054902889094487,
+  	0.0050687282939456,
+  	0.0046449591497966,
+  	0.0042191661429919,
+  	0.0037915348363451,
+  	0.0033622516236779,
+  	0.0029315036836558,
+  	0.0024994789888943,
+  	0.0020663664924131,
+  	0.0016323569986067,
+  	0.0011976474864367,
+  	0.0007624720924706,
+  	0.0003276086705538
+])
Index: sasmodels/weights.py
===================================================================
--- sasmodels/weights.py	(revision f1a8811c147abcdc8d3fb84ee7e8a25fa91be625)
+++ sasmodels/weights.py	(revision d4c33d66aa792e691ad5eb91b36d3e854e88f383)
@@ -55,4 +55,8 @@
         """
         sigma = self.width * center if relative else self.width
+        if not relative:
+            # For orientation, the jitter is relative to 0 not the angle
+            #center = 0
+            pass
         if sigma == 0 or self.npts < 2:
             if lb <= center <= ub:
