Index: explore/angles.py
===================================================================
--- explore/angles.py	(revision ea60e084e3e40eeb3abfc5f728359135acaed020)
+++ explore/angles.py	(revision 20fe0cdaae1d6a1be701c7377a8baaaac269d234)
@@ -229,5 +229,10 @@
 
 if 0:
-    comment("==== asymmetric (old) ====")
+    comment("==== asymmetric (3.x) ====")
     view_inv = Rz(-psi)*Rx(theta)*Ry(-(pi/2 - phi))
     vprint(q_abc, view_inv*q_xy, "reverse view")
+    print("""  existing code
+    cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat;
+    cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat;
+    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat;
+    """)
Index: explore/asymint.py
===================================================================
--- explore/asymint.py	(revision 31eea1f72ee899d439cb11c51d68886b44f27ab2)
+++ explore/asymint.py	(revision 20fe0cdaae1d6a1be701c7377a8baaaac269d234)
@@ -1,2 +1,3 @@
+#!/usr/bin/env python
 """
 Asymmetric shape integration
@@ -6,5 +7,5 @@
 
 import os, sys
-sys.path.insert(0, os.path.dirname(os.path.dirname(__file__)))
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__))))
 
 import numpy as np
@@ -15,16 +16,39 @@
 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
-from sasmodels.special import sas_2J1x_x, sas_sinx_x, sas_3j1x_x
-
-def mp_3j1x_x(x):
-    return 3*(mp.sin(x)/x - mp.cos(x))/(x*x)
-def mp_2J1x_x(x):
-    return 2*mp.j1(x)/x
-def mp_sinx_x(x):
-    return mp.sin(x)/x
+import sasmodels.special as sp
+
+# Need to parse shape early since it determines the kernel function
+# that will be used for the various integrators
+shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1]
+Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2]
+
+class MPenv:
+    sqrt = mp.sqrt
+    exp = mp.exp
+    expm1 = mp.expm1
+    cos = mp.cos
+    sin = mp.sin
+    tan = mp.tan
+    def sas_3j1x_x(x):
+        return 3*(mp.sin(x)/x - mp.cos(x))/(x*x)
+    def sas_2J1x_x(x):
+        return 2*mp.j1(x)/x
+    def sas_sinx_x(x):
+        return mp.sin(x)/x
+    pi = mp.pi
+    mpf = mp.mpf
+
+class NPenv:
+    sqrt = np.sqrt
+    exp = np.exp
+    expm1 = np.expm1
+    cos = np.cos
+    sin = np.sin
+    tan = np.tan
+    sas_3j1x_x = sp.sas_3j1x_x
+    sas_2J1x_x = sp.sas_2J1x_x
+    sas_sinx_x = sp.sas_sinx_x
+    pi = np.pi
+    mpf = float
 
 SLD = 3
@@ -32,102 +56,125 @@
 CONTRAST = SLD - SLD_SOLVENT
 
-def make_parallelepiped(a, b, c):
-    def Fq(qa, qb, qc):
-        siA = sas_sinx_x(0.5*a*qa)
-        siB = sas_sinx_x(0.5*b*qb)
-        siC = sas_sinx_x(0.5*c*qc)
+# Carefully code models so that mpmath will use full precision.  That means:
+#    * wrap inputs in env.mpf
+#    * don't use floating point constants, only integers
+#    * for division, make sure the numerator or denominator is env.mpf
+#    * use env.pi, env.sas_sinx_x, etc. for functions
+def make_parallelepiped(a, b, c, env=NPenv):
+    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
+    def Fq(qa, qb, qc):
+        siA = env.sas_sinx_x(0.5*a*qa/2)
+        siB = env.sas_sinx_x(0.5*b*qb/2)
+        siC = env.sas_sinx_x(0.5*c*qc/2)
         return siA * siB * siC
+    Fq.__doc__ = "parallelepiped a=%g, b=%g c=%g"%(a, b, c)
     volume = a*b*c
-    norm = volume*CONTRAST**2/10**4
-    return norm, Fq
-
-def make_parallelepiped_mp(a, b, c):
-    a, b, c = mp.mpf(a), mp.mpf(b), mp.mpf(c)
-    def Fq(qa, qb, qc):
-        siA = mp_sinx_x(a*qa/2)
-        siB = mp_sinx_x(b*qb/2)
-        siC = mp_sinx_x(c*qc/2)
-        return siA * siB * siC
-    volume = a*b*c
-    norm = (volume*CONTRAST**2)/10000 # mpf since volume=a*b*c is mpf
-    return norm, Fq
-
-def make_triellip(a, b, c):
-    def Fq(qa, qb, qc):
-        qr = sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2)
-        return sas_3j1x_x(qr)
-    volume = 4*pi*a*b*c/3
-    norm = volume*CONTRAST**2/10**4
-    return norm, Fq
-
-def make_triellip_mp(a, b, c):
-    a, b, c = mp.mpf(a), mp.mpf(b), mp.mpf(c)
-    def Fq(qa, qb, qc):
-        qr = mp.sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2)
-        return mp_3j1x_x(qr)
-    volume = (4*mp.pi*a*b*c)/3
-    norm = (volume*CONTRAST**2)/10000  # mpf since mp.pi is mpf
-    return norm, Fq
-
-def make_cylinder(radius, length):
-    def Fq(qa, qb, qc):
-        qab = sqrt(qa**2 + qb**2)
-        return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length)
-    volume = pi*radius**2*length
-    norm = volume*CONTRAST**2/10**4
-    return norm, Fq
-
-def make_cylinder_mp(radius, length):
-    radius, length = mp.mpf(radius), mp.mpf(length)
-    def Fq(qa, qb, qc):
-        qab = mp.sqrt(qa**2 + qb**2)
-        return mp_2J1x_x(qab*radius) * mp_sinx_x((qc*length)/2)
-    volume = mp.pi*radius**2*length
-    norm = (volume*CONTRAST**2)/10000  # mpf since mp.pi is mpf
-    return norm, Fq
-
-def make_sphere(radius):
-    def Fq(qa, qb, qc):
-        q = sqrt(qa**2 + qb**2 + qc**2)
-        return sas_3j1x_x(q*radius)
-    volume = 4*pi*radius**3/3
-    norm = volume*CONTRAST**2/10**4
-    return norm, Fq
-
-def make_sphere_mp(radius):
-    radius = mp.mpf(radius)
-    def Fq(qa, qb, qc):
-        q = mp.sqrt(qa**2 + qb**2 + qc**2)
-        return mp_3j1x_x(q*radius)
-    volume = (4*mp.pi*radius**3)/3
-    norm = (volume*CONTRAST**2)/10000  # mpf since mp.pi is mpf
-    return norm, Fq
-
-shape = 'parallelepiped'
-#shape = 'triellip'
-#shape = 'sphere'
-#shape = 'cylinder'
-if shape == 'cylinder':
+    norm = CONTRAST**2*volume/10000
+    return norm, Fq
+
+def make_triaxial_ellipsoid(a, b, c, env=NPenv):
+    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
+    def Fq(qa, qb, qc):
+        qr = env.sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2)
+        return env.sas_3j1x_x(qr)
+    Fq.__doc__ = "triaxial ellipse minor=%g, major=%g polar=%g"%(a, b, c)
+    volume = 4*env.pi*a*b*c/3
+    norm = CONTRAST**2*volume/10000
+    return norm, Fq
+
+def make_cylinder(radius, length, env=NPenv):
+    radius, length = env.mpf(radius), env.mpf(length)
+    def Fq(qa, qb, qc):
+        qab = env.sqrt(qa**2 + qb**2)
+        return env.sas_2J1x_x(qab*radius) * env.sas_sinx_x((qc*length)/2)
+    Fq.__doc__ = "cylinder radius=%g, length=%g"%(radius, length)
+    volume = env.pi*radius**2*length
+    norm = CONTRAST**2*volume/10000
+    return norm, Fq
+
+def make_sphere(radius, env=NPenv):
+    radius = env.mpf(radius)
+    def Fq(qa, qb, qc):
+        q = env.sqrt(qa**2 + qb**2 + qc**2)
+        return env.sas_3j1x_x(q*radius)
+    Fq.__doc__ = "sphere radius=%g"%(radius, )
+    volume = 4*pi*radius**3
+    norm = CONTRAST**2*volume/10000
+    return norm, Fq
+
+def make_paracrystal(radius, dnn, d_factor, lattice='bcc', env=NPenv):
+    radius, dnn, d_factor = env.mpf(radius), env.mpf(dnn), env.mpf(d_factor)
+    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 fcc(qa, qb, qc):
+        a1 = ( 0. + qb + qc)/2
+        a2 = (-qa + 0. + qc)/2
+        a3 = (-qa + qb + 0.)/2
+        return a1, a2, a3
+    lattice_fn = {'sc': sc, 'bcc': bcc, 'fcc': fcc}[lattice]
+    radius, dnn, d_factor = env.mpf(radius), env.mpf(dnn), env.mpf(d_factor)
+    def Fq(qa, qb, qc):
+        a1, a2, a3 = lattice_fn(qa, qb, qc)
+        # Note: paper says that different directions can have different
+        # distoration factors.  Easy enough to add to the code.
+        # This would definitely break 8-fold symmetry.
+        arg = -(dnn*d_factor)**2*(a1**2 + a2**2 + a3**2)/2
+        exp_arg = env.exp(arg)
+        den = [((exp_arg - 2*env.cos(dnn*a))*exp_arg + 1) for a in (a1, a2, a3)]
+        Sq = -env.expm1(2*arg)**3/(den[0]*den[1]*den[2])
+
+        q = env.sqrt(qa**2 + qb**2 + qc**2)
+        Fq = env.sas_3j1x_x(q*radius)
+        # the kernel computes F(q)**2, but we need S(q)*F(q)**2
+        return env.sqrt(Sq)*Fq
+    Fq.__doc__ = "%s paracrystal a=%g da=%g r=%g"%(lattice, dnn, d_factor, radius)
+    def sphere_volume(r): return 4*env.pi*r**3/3
+    Vf = {
+        'sc': sphere_volume(radius/dnn),
+        'bcc': 2*sphere_volume(env.sqrt(3)/2*radius/dnn),
+        'fcc': 4*sphere_volume(radius/dnn/env.sqrt(2)),
+    }[lattice]
+    volume = sphere_volume(radius)
+    norm = CONTRAST**2*volume*Vf/10000
+    return norm, Fq
+
+if shape == 'sphere':
+    RADIUS = 50  # integer for the sake of mpf
+    NORM, KERNEL = make_sphere(radius=RADIUS)
+    NORM_MP, KERNEL_MP = make_sphere(radius=RADIUS, env=MPenv)
+elif shape == 'cylinder':
     #RADIUS, LENGTH = 10, 100000
     RADIUS, LENGTH = 10, 300  # integer for the sake of mpf
     NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH)
-    NORM_MP, KERNEL_MP = make_cylinder_mp(radius=RADIUS, length=LENGTH)
-elif shape == 'triellip':
+    NORM_MP, KERNEL_MP = make_cylinder(radius=RADIUS, length=LENGTH, env=MPenv)
+elif shape == 'triaxial_ellipsoid':
     #A, B, C = 4450, 14000, 47
     A, B, C = 445, 140, 47  # integer for the sake of mpf
     NORM, KERNEL = make_triellip(A, B, C)
-    NORM_MP, KERNEL_MP = make_triellip_mp(A, B, C)
+    NORM_MP, KERNEL_MP = make_triellip(A, B, C, env=MPenv)
 elif shape == 'parallelepiped':
     #A, B, C = 4450, 14000, 47
     A, B, C = 445, 140, 47  # integer for the sake of mpf
     NORM, KERNEL = make_parallelepiped(A, B, C)
-    NORM_MP, KERNEL_MP = make_parallelepiped_mp(A, B, C)
-elif shape == 'sphere':
-    RADIUS = 50  # integer for the sake of mpf
-    NORM, KERNEL = make_sphere(radius=RADIUS)
-    NORM_MP, KERNEL_MP = make_sphere_mp(radius=RADIUS)
+    NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv)
+elif shape == 'paracrystal':
+    LATTICE = 'bcc'
+    #LATTICE = 'fcc'
+    #LATTICE = 'sc'
+    DNN, D_FACTOR = 220, '0.06'  # mpmath needs to initialize floats from string
+    RADIUS = 40  # integer for the sake of mpf
+    NORM, KERNEL = make_paracrystal(
+        radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE)
+    NORM_MP, KERNEL_MP = make_paracrystal(
+        radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE, env=MPenv)
 else:
     raise ValueError("Unknown shape %r"%shape)
 
+# Note: hardcoded in mp_quad
 THETA_LOW, THETA_HIGH = 0, pi
 PHI_LOW, PHI_HIGH = 0, 2*pi
@@ -145,6 +192,6 @@
 """
 
-
-def mp_quad(q, shape):
+# 2D integration functions
+def mp_quad_2d(q, shape):
     evals = [0]
     def integrand(theta, phi):
@@ -160,5 +207,5 @@
     return evals[0], Iq
 
-def kernel(q, theta, phi):
+def kernel_2d(q, theta, phi):
     """
     S(q) kernel for paracrystal forms.
@@ -170,5 +217,5 @@
     return NORM*KERNEL(qa, qb, qc)**2
 
-def scipy_dblquad(q):
+def scipy_dblquad_2d(q):
     """
     Compute the integral using scipy dblquad.  This gets the correct answer
@@ -178,5 +225,5 @@
     def integrand(phi, theta):
         evals[0] += 1
-        Zq = kernel(q, theta=theta, phi=phi)
+        Zq = kernel_2d(q, theta=theta, phi=phi)
         return Zq*sin(theta)
     ans = dblquad(integrand, THETA_LOW, THETA_HIGH, lambda x: PHI_LOW, lambda x: PHI_HIGH)[0]
@@ -192,5 +239,5 @@
     def inner(phi, theta):
         evals[0] += 1
-        return kernel(q, theta=theta, phi=phi)
+        return kernel_2d(q, theta=theta, phi=phi)
     def outer(theta):
         Zq = romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(theta,))
@@ -200,5 +247,5 @@
 
 
-def semi_romberg(q, n=100):
+def semi_romberg_2d(q, n=100):
     """
     Use 1D romberg integration in phi and regular simpsons rule in theta.
@@ -207,5 +254,5 @@
     def inner(phi, theta):
         evals[0] += 1
-        return kernel(q, theta=theta, phi=phi)
+        return kernel_2d(q, theta=theta, phi=phi)
     theta = np.linspace(THETA_LOW, THETA_HIGH, n)
     Zq = [romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(t,)) for t in theta]
@@ -213,16 +260,9 @@
     return evals[0], ans*SCALE/(4*pi)
 
-def gauss_quad(q, n=150):
+def gauss_quad_2d(q, 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
-    elif n == 150:
-        z, w = Gauss150Z, Gauss150Wt
-    else:
-        z, w = leggauss(n)
+    z, w = leggauss(n)
     theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW
     phi = (PHI_HIGH-PHI_LOW)*(z + 1)/2 + PHI_LOW
@@ -230,5 +270,5 @@
     Aw = w[None, :] * w[:, None]
     sin_theta = abs(sin(Atheta))
-    Zq = kernel(q=q, theta=Atheta, phi=Aphi)
+    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
     # change from [-1,1] x [-1,1] range to [0, pi] x [0, 2 pi] range
     dxdy_stretch = (THETA_HIGH-THETA_LOW)/2 * (PHI_HIGH-PHI_LOW)/2
@@ -236,6 +276,5 @@
     return n**2, Iq
 
-
-def gridded_integrals(q, n=300):
+def gridded_2d(q, n=300):
     """
     Compute the integral on a regular grid using rectangular, trapezoidal,
@@ -246,5 +285,5 @@
     phi = np.linspace(PHI_LOW, PHI_HIGH, n)
     Atheta, Aphi = np.meshgrid(theta, phi)
-    Zq = kernel(q=q, theta=Atheta, phi=Aphi)
+    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
     Zq *= abs(sin(Atheta))
     dx, dy = theta[1]-theta[0], phi[1]-phi[0]
@@ -254,5 +293,5 @@
     print("romb-%d"%n, n**2, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi))
 
-def plot(q, n=300):
+def plot_2d(q, n=300):
     """
     Plot the 2D surface that needs to be integrated in order to compute
@@ -263,9 +302,9 @@
     phi = np.linspace(PHI_LOW, PHI_HIGH, n)
     Atheta, Aphi = np.meshgrid(theta, phi)
-    Zq = kernel(q=q, theta=Atheta, phi=Aphi)
+    Zq = kernel_2d(q=q, 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" % (KERNEL.__name__, q))
+    pylab.title("%s I(q,t) sin(t) for q=%g" % (KERNEL.__doc__, q))
     pylab.xlabel("theta (degrees)")
     pylab.ylabel("phi (degrees)")
@@ -274,23 +313,24 @@
     pylab.show()
 
-if __name__ == "__main__":
-    Qstr = '0.005'
-    #Qstr = '0.8'
-    #Qstr = '0.0003'
+def main(Qstr):
     Q = float(Qstr)
     if shape == 'sphere':
         print("exact", NORM*sas_3j1x_x(Q*RADIUS)**2)
-    print("gauss-20", *gauss_quad(Q, n=20))
-    print("gauss-76", *gauss_quad(Q, n=76))
-    print("gauss-150", *gauss_quad(Q, n=150))
-    print("gauss-500", *gauss_quad(Q, n=500))
-    print("dblquad", *scipy_dblquad(Q))
-    print("semi-romberg-100", *semi_romberg(Q, n=100))
-    print("romberg", *scipy_romberg_2d(Q))
-    #gridded_integrals(Q, n=2**8+1)
-    gridded_integrals(Q, n=2**10+1)
-    #gridded_integrals(Q, n=2**13+1)
-    #gridded_integrals(Q, n=2**15+1)
-    with mp.workprec(100):
-        print("mpmath", *mp_quad(mp.mpf(Qstr), shape))
-    #plot(Q, n=200)
+    print("gauss-20", *gauss_quad_2d(Q, n=20))
+    print("gauss-76", *gauss_quad_2d(Q, n=76))
+    print("gauss-150", *gauss_quad_2d(Q, n=150))
+    print("gauss-500", *gauss_quad_2d(Q, n=500))
+    #gridded_2d(Q, n=2**8+1)
+    gridded_2d(Q, n=2**10+1)
+    #gridded_2d(Q, n=2**13+1)
+    #gridded_2d(Q, n=2**15+1)
+    if shape != 'paracrystal':  # adaptive forms are too slow!
+        print("dblquad", *scipy_dblquad_2d(Q))
+        print("semi-romberg-100", *semi_romberg_2d(Q, n=100))
+        print("romberg", *scipy_romberg_2d(Q))
+        with mp.workprec(100):
+            print("mpmath", *mp_quad_2d(mp.mpf(Qstr), shape))
+    plot_2d(Q, n=200)
+
+if __name__ == "__main__":
+    main(Qstr)
Index: plore/bccpy.py
===================================================================
--- explore/bccpy.py	(revision 7e0b281e5f175a91044364982fcfbfe802508948)
+++ 	(revision )
@@ -1,217 +1,0 @@
-"""
-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/symint.py
===================================================================
--- explore/symint.py	(revision e47b06b53a04280d71f38038efd95d90cb05a2f8)
+++ explore/symint.py	(revision 20fe0cdaae1d6a1be701c7377a8baaaac269d234)
@@ -10,6 +10,6 @@
 import numpy as np
 from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10
+from numpy.polynomial.legendre import leggauss
 from scipy.integrate import dblquad, simps, romb, romberg
-from scipy.special.orthogonal import p_roots
 import pylab
 
@@ -21,5 +21,5 @@
 
 SLD = 3.0
-SLD_SOLVENT = 6.3
+SLD_SOLVENT = 6
 CONTRAST = SLD - SLD_SOLVENT
 
@@ -29,5 +29,5 @@
     cylinder.__doc__ = "cylinder radius=%g, length=%g"%(radius, length)
     volume = pi*radius**2*length
-    norm = 1e-4*volume*CONTRAST**2
+    norm = CONTRAST**2*volume/10000
     return norm, cylinder
 
@@ -37,5 +37,5 @@
     long_cylinder.__doc__ = "long cylinder radius=%g, length=%g"%(radius, length)
     volume = pi*radius**2*length
-    norm = 1e-4*volume*CONTRAST**2*pi/length
+    norm = CONTRAST**2*volume/10000*pi/length
     return long_cylinder
 
@@ -46,5 +46,5 @@
     sphere.__doc__ = "sphere radius=%g"%(radius,)
     volume = 4*pi*radius**3/3
-    norm = 1e-4*volume*CONTRAST**2
+    norm = CONTRAST**2*volume/10000
     return norm, sphere
 
@@ -53,5 +53,5 @@
 
 
-def kernel(q, theta):
+def kernel_1d(q, theta):
     """
     S(q) kernel for paracrystal forms.
@@ -61,17 +61,15 @@
     return NORM*KERNEL(qab, qc)**2
 
-
-def gauss_quad(q, n=150):
+def gauss_quad_1d(q, n=150):
     """
     Compute the integral using gaussian quadrature for n = 20, 76 or 150.
     """
-    z, w = p_roots(n)
+    z, w = leggauss(n)
     theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW
     sin_theta = abs(sin(theta))
-    Zq = kernel(q=q, theta=theta)
+    Zq = kernel_1d(q=q, theta=theta)
     return np.sum(Zq*w*sin_theta)*(THETA_HIGH-THETA_LOW)/2
 
-
-def gridded_integrals(q, n=300):
+def gridded_1d(q, n=300):
     """
     Compute the integral on a regular grid using rectangular, trapezoidal,
@@ -80,13 +78,13 @@
     """
     theta = np.linspace(THETA_LOW, THETA_HIGH, n)
-    Zq = kernel(q=q, theta=theta)
+    Zq = kernel_1d(q=q, theta=theta)
     Zq *= abs(sin(theta))
     dx = theta[1]-theta[0]
-    print("rect", n, np.sum(Zq)*dx*SCALE)
-    print("trapz", n, np.trapz(Zq, dx=dx)*SCALE)
-    print("simpson", n, simps(Zq, dx=dx)*SCALE)
-    print("romb", n, romb(Zq, dx=dx)*SCALE)
+    print("rect-%d"%n, np.sum(Zq)*dx*SCALE)
+    print("trapz-%d"%n, np.trapz(Zq, dx=dx)*SCALE)
+    print("simpson-%d"%n, simps(Zq, dx=dx)*SCALE)
+    print("romb-%d"%n, romb(Zq, dx=dx)*SCALE)
 
-def scipy_romberg(q):
+def scipy_romberg_1d(q):
     """
     Compute the integral using romberg integration.  This function does not
@@ -96,16 +94,15 @@
     def outer(theta):
         evals[0] += 1
-        return kernel(q, theta=theta)*abs(sin(theta))
+        return kernel_1d(q, theta=theta)*abs(sin(theta))
     result = romberg(outer, THETA_LOW, THETA_HIGH, divmax=100)*SCALE
     print("scipy romberg", evals[0], result)
 
-def plot(q, n=300):
+def plot_1d(q, 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.
+    Plot the function that needs to be integrated in order to compute
+    the I(q) at a particular q.  *n* is the number of points in the grid.
     """
     theta = np.linspace(THETA_LOW, THETA_HIGH, n)
-    Zq = kernel(q=q, theta=theta)
+    Zq = kernel_1d(q=q, theta=theta)
     Zq *= abs(sin(theta))
     pylab.semilogy(degrees(theta), np.fmax(Zq, 1.e-6), label="Q=%g"%q)
@@ -116,5 +113,5 @@
 def Iq_trapz(q, n):
     theta = np.linspace(THETA_LOW, THETA_HIGH, n)
-    Zq = kernel(q=q, theta=theta)
+    Zq = kernel_1d(q=q, theta=theta)
     Zq *= abs(sin(theta))
     dx = theta[1]-theta[0]
@@ -142,13 +139,13 @@
     Q = 0.386
     for n in (20, 76, 150, 300, 1000): #, 10000, 30000):
-        print("gauss", n, gauss_quad(Q, n=n))
+        print("gauss-%d"%n, gauss_quad(Q, n=n))
     for k in (8, 10, 13, 16, 19):
-        gridded_integrals(Q, n=2**k+1)
+        gridded_1d(Q, n=2**k+1)
     #print("inf cyl", 0, long_cyl(Q))
     #scipy_romberg(Q)
 
-    plot(0.386, n=2000)
-    plot(0.5, n=2000)
-    plot(0.8, n=2000)
+    plot_1d(0.386, n=2000)
+    plot_1d(0.5, n=2000)
+    plot_1d(0.8, n=2000)
     pylab.legend()
     pylab.figure()
Index: explore/transform_angles.py
===================================================================
--- explore/transform_angles.py	(revision 20fe0cdaae1d6a1be701c7377a8baaaac269d234)
+++ explore/transform_angles.py	(revision 20fe0cdaae1d6a1be701c7377a8baaaac269d234)
@@ -0,0 +1,57 @@
+#!/usr/bin/env python
+from __future__ import print_function, division
+
+import sys
+
+import numpy as np
+from numpy import pi, cos, sin, sqrt, exp, degrees, radians
+from scipy.optimize import fmin
+
+# Definition of rotation matrices comes from wikipedia:
+#    https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations
+def Rx(angle):
+    """Construct a matrix to rotate points about *x* by *angle* degrees."""
+    a = radians(angle)
+    R = [[1, 0, 0],
+         [0, +cos(a), -sin(a)],
+         [0, +sin(a), +cos(a)]]
+    return np.matrix(R)
+
+def Ry(angle):
+    """Construct a matrix to rotate points about *y* by *angle* degrees."""
+    a = radians(angle)
+    R = [[+cos(a), 0, +sin(a)],
+         [0, 1, 0],
+         [-sin(a), 0, +cos(a)]]
+    return np.matrix(R)
+
+def Rz(angle):
+    """Construct a matrix to rotate points about *z* by *angle* degrees."""
+    a = radians(angle)
+    R = [[+cos(a), -sin(a), 0],
+         [+sin(a), +cos(a), 0],
+         [0, 0, 1]]
+    return np.matrix(R)
+
+
+def transform_angles(theta, phi, psi, qx=0.1, qy=0.1):
+    Rold = Rz(-psi)*Rx(theta)*Ry(-(90 - phi))
+    cost = lambda p: np.linalg.norm(Rz(-p[2])*Ry(-p[0])*Rz(-p[1]) - Rold)
+    result = fmin(cost, (theta, phi, psi))
+    theta_p, phi_p, psi_p = result
+    Rnew = Rz(-psi_p)*Ry(-theta_p)*Rz(-phi_p)
+
+    print("old: theta, phi, psi =", ", ".join(str(v) for v in (theta, phi, psi)))
+    print("new: theta, phi, psi =", ", ".join(str(v) for v in result))
+    try:
+        point = np.matrix([qx, qy, [0]*len(qx)])
+    except TypeError:
+        point = np.matrix([[qx],[qy],[0]])
+    for p in point.T:
+        print("q abc old for", p, (Rold*p.T).T)
+        print("q abc new for", p, (Rnew*p.T).T)
+
+if __name__ == "__main__":
+    theta, phi, psi = (float(v) for v in sys.argv[1:])
+    #transform_angles(theta, phi, psi)
+    transform_angles(theta, phi, psi, qx=-0.017, qy=0.035)
Index: sasmodels/model_test.py
===================================================================
--- sasmodels/model_test.py	(revision 65314f75d3d74de18ada1d98440acf06ab8c4ce3)
+++ sasmodels/model_test.py	(revision 20fe0cdaae1d6a1be701c7377a8baaaac269d234)
@@ -202,4 +202,5 @@
                 ]
             tests = smoke_tests
+            #tests = []
             if self.info.tests is not None:
                 tests += self.info.tests
