Index: explore/bccpy.py
===================================================================
--- explore/bccpy.py	(revision cb038a27e078c34759b4e3be747478bb31d93f9b)
+++ explore/bccpy.py	(revision 7e0b281e5f175a91044364982fcfbfe802508948)
@@ -47,4 +47,38 @@
 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):
     """
@@ -56,23 +90,14 @@
     qc = q*cos(theta)
 
-    if 0: # sc
-        a1, a2, a3 = qa, qb, qc
-        dcos = dnn
-    if 1: # bcc
-        a1 = +qa - qc + qb
-        a2 = +qa + qc - qb
-        a3 = -qa + qc + qb
-        dcos = dnn/2
-    if 0: # fcc
-        a1 = qb + qa
-        a2 = qa + qc
-        a3 = qb + qc
-        dcos = dnn/2
-
-    arg = 0.5*square(dnn*d_factor)*(a1**2 + a2**2 + a3**2)
-    exp_arg = exp(-arg)
-    den = [((exp_arg - 2*cos(dcos*a))*exp_arg + 1.0) for a in (a1, a2, a3)]
-    Sq = -expm1(-2*arg)**3/np.prod(den, axis=0)
-    return Sq
+    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
 
 
@@ -85,6 +110,6 @@
     def integrand(theta, phi):
         evals[0] += 1
-        Sq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=theta, phi=phi)
-        return Sq*sin(theta)
+        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])
@@ -134,7 +159,7 @@
     Aw = w[None, :] * w[:, None]
     sin_theta = np.fmax(abs(sin(Atheta)), 1e-6)
-    Sq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi)
+    Zq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi)
     print("gauss %d evals ="%n, n**2)
-    return np.sum(Sq*Aw*sin_theta)*8/(4*pi)
+    return np.sum(Zq*Aw*sin_theta)*8/(4*pi)
 
 
@@ -148,11 +173,11 @@
     phi = np.linspace(0, pi/2, n)
     Atheta, Aphi = np.meshgrid(theta, phi)
-    Sq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi)
-    Sq *= abs(sin(Atheta))
+    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(Sq)*dx*dy*8/(4*pi))
-    print("trapz", n, np.trapz(np.trapz(Sq, dx=dx), dx=dy)*8/(4*pi))
-    print("simpson", n, simps(simps(Sq, dx=dx), dx=dy)*8/(4*pi))
-    print("romb", n, romb(romb(Sq, dx=dx), dx=dy)*8/(4*pi))
+    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)
 
@@ -168,9 +193,10 @@
     #phi = np.linspace(0, pi/2, n)
     Atheta, Aphi = np.meshgrid(theta, phi)
-    Sq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi)
-    Sq *= abs(sin(Atheta))
-    pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Sq, 1.e-6)))
+    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("BCC S(q) for q=%g, dnn=%g d_factor=%g" % (q, dnn, d_factor))
+    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)")
Index: explore/sc.c
===================================================================
--- explore/sc.c	(revision fdd56a1decdae24f11cdccb11ac7ffade9812587)
+++ explore/sc.c	(revision 7e0b281e5f175a91044364982fcfbfe802508948)
@@ -1,4 +1,4 @@
 static double
-_sq_sc(double qa, double qb, double qc, double dnn, double d_factor)
+sc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
 {
     // Rewriting equations for efficiency, accuracy and readability, and so
@@ -8,5 +8,5 @@
     const double a3 = qc;
 
-    const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
 
     // Numerator: (1 - exp(a)^2)^3
@@ -16,16 +16,16 @@
     //         => exp(a)^2 - 2 cos(xk) exp(a) + 1
     //         => (exp(a) - 2 cos(xk)) * exp(a) + 1
-    const double exp_arg = exp(-arg);
-    const double Sq = -cube(expm1(-2.0*arg))
+    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 Sq;
+    return Zq;
 }
 
 // occupied volume fraction calculated from lattice symmetry and sphere radius
 static double
-_sc_volume_fraction(double radius, double dnn)
+sc_volume_fraction(double radius, double dnn)
 {
     return sphere_volume(radius/dnn);
@@ -98,6 +98,6 @@
             const double qa = qab*cos_phi;
             const double qb = qab*sin_phi;
-            const double fq = _sq_sc(qa, qb, qc, dnn, d_factor);
-            inner_sum += fq;
+            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
@@ -106,14 +106,14 @@
     outer_sum *= theta_m/(n*n);
 #endif
-double Sq;
+double Zq;
 if (sym > 0.) {
-    Sq = outer_sum/M_PI_2;
+    Zq = outer_sum/M_PI_2;
 } else {
-    Sq = outer_sum/(4.0*M_PI);
+    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 * Sq;
+    return sc_volume_fraction(radius, dnn) * Pq * Zq;
 }
 
@@ -133,5 +133,5 @@
     q = sqrt(qa*qa + qb*qb + qc*qc);
     const double Pq = sphere_form(q, radius, sld, solvent_sld);
-    const double Sq = _sq_sc(qa, qb, qc, dnn, d_factor);
-    return _sc_volume_fraction(radius, dnn) * Pq * Sq;
+    const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor);
+    return sc_volume_fraction(radius, dnn) * Pq * Zq;
 }
Index: sasmodels/models/bcc_paracrystal.c
===================================================================
--- sasmodels/models/bcc_paracrystal.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
+++ sasmodels/models/bcc_paracrystal.c	(revision 7e0b281e5f175a91044364982fcfbfe802508948)
@@ -1,13 +1,14 @@
 static double
-_sq_bcc(double qa, double qb, double qc, double dnn, double d_factor)
+bcc_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.
-    const double a1 = +qa - qc + qb;
-    const double a2 = +qa + qc - qb;
-    const double a3 = -qa + qc + qb;
-
-    const double half_dnn = 0.5*dnn;
-    const double arg = 0.5*square(half_dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+#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
 
 #if 1
@@ -15,22 +16,24 @@
     //         => (-(exp(2a) - 1))^3
     //         => -expm1(2a)^3
-    // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2)
-    //         => exp(a)^2 - 2 cos(xk) exp(a) + 1
-    //         => (exp(a) - 2 cos(xk)) * exp(a) + 1
-    const double exp_arg = exp(-arg);
-    const double Sq = -cube(expm1(-2.0*arg))
-        / ( ((exp_arg - 2.0*cos(half_dnn*a1))*exp_arg + 1.0)
-          * ((exp_arg - 2.0*cos(half_dnn*a2))*exp_arg + 1.0)
-          * ((exp_arg - 2.0*cos(half_dnn*a3))*exp_arg + 1.0));
+    // 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 Sq = sinh_qd/(cosh_qd - cos(half_dnn*a1))
-                    * sinh_qd/(cosh_qd - cos(half_dnn*a2))
-                    * sinh_qd/(cosh_qd - cos(half_dnn*a3));
+    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 Sq;
+    return Zq;
 }
 
@@ -38,5 +41,5 @@
 // occupied volume fraction calculated from lattice symmetry and sphere radius
 static double
-_bcc_volume_fraction(double radius, double dnn)
+bcc_volume_fraction(double radius, double dnn)
 {
     return 2.0*sphere_volume(sqrt(0.75)*radius/dnn);
@@ -75,6 +78,6 @@
             const double qa = qab*cos_phi;
             const double qb = qab*sin_phi;
-            const double fq = _sq_bcc(qa, qb, qc, dnn, d_factor);
-            inner_sum += Gauss150Wt[j] * fq;
+            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
@@ -82,8 +85,7 @@
     }
     outer_sum *= theta_m;
-    const double Sq = outer_sum/(4.0*M_PI);
+    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 * Sq;
+    return bcc_volume_fraction(radius, dnn) * Pq * Zq;
 }
 
@@ -101,6 +103,6 @@
 
     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);
-    const double Sq = _sq_bcc(qa, qb, qc, dnn, d_factor);
-    return _bcc_volume_fraction(radius, dnn) * Pq * Sq;
+    return bcc_volume_fraction(radius, dnn) * Pq * Zq;
 }
Index: sasmodels/models/fcc_paracrystal.c
===================================================================
--- sasmodels/models/fcc_paracrystal.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
+++ sasmodels/models/fcc_paracrystal.c	(revision 7e0b281e5f175a91044364982fcfbfe802508948)
@@ -1,27 +1,29 @@
 static double
-_sq_fcc(double qa, double qb, double qc, double dnn, double d_factor)
+fcc_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.
-    const double a1 = qb + qa;
-    const double a2 = qa + qc;
-    const double a3 = qb + qc;
-
-    const double half_dnn = 0.5*dnn;
-    const double arg = 0.5*square(half_dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+#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
 
     // Numerator: (1 - exp(a)^2)^3
     //         => (-(exp(2a) - 1))^3
     //         => -expm1(2a)^3
-    // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2)
-    //         => exp(a)^2 - 2 cos(xk) exp(a) + 1
-    //         => (exp(a) - 2 cos(xk)) * exp(a) + 1
-    const double exp_arg = exp(-arg);
-    const double Sq = -cube(expm1(-2.0*arg))
-        / ( ((exp_arg - 2.0*cos(half_dnn*a1))*exp_arg + 1.0)
-          * ((exp_arg - 2.0*cos(half_dnn*a2))*exp_arg + 1.0)
-          * ((exp_arg - 2.0*cos(half_dnn*a3))*exp_arg + 1.0));
+    // 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 Sq;
+    return Zq;
 }
 
@@ -29,5 +31,5 @@
 // occupied volume fraction calculated from lattice symmetry and sphere radius
 static double
-_fcc_volume_fraction(double radius, double dnn)
+fcc_volume_fraction(double radius, double dnn)
 {
     return 4.0*sphere_volume(M_SQRT1_2*radius/dnn);
@@ -66,6 +68,6 @@
             const double qa = qab*cos_phi;
             const double qb = qab*sin_phi;
-            const double fq = _sq_fcc(qa, qb, qc, dnn, d_factor);
-            inner_sum += Gauss150Wt[j] * fq;
+            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
@@ -73,8 +75,8 @@
     }
     outer_sum *= theta_m;
-    const double Sq = outer_sum/(4.0*M_PI);
+    const double Zq = outer_sum/(4.0*M_PI);
     const double Pq = sphere_form(q, radius, sld, solvent_sld);
 
-    return _fcc_volume_fraction(radius, dnn) * Pq * Sq;
+    return fcc_volume_fraction(radius, dnn) * Pq * Zq;
 }
 
@@ -93,5 +95,5 @@
     q = sqrt(qa*qa + qb*qb + qc*qc);
     const double Pq = sphere_form(q, radius, sld, solvent_sld);
-    const double Sq = _sq_fcc(qa, qb, qc, dnn, d_factor);
-    return _fcc_volume_fraction(radius, dnn) * Pq * Sq;
+    const double Zq = fcc_Zq(qa, qb, qc, dnn, d_factor);
+    return fcc_volume_fraction(radius, dnn) * Pq * Zq;
 }
Index: sasmodels/models/sc_paracrystal.c
===================================================================
--- sasmodels/models/sc_paracrystal.c	(revision 2a0b2b133086c070f586ead5c4a343a7e57ecad5)
+++ sasmodels/models/sc_paracrystal.c	(revision 7e0b281e5f175a91044364982fcfbfe802508948)
@@ -1,31 +1,28 @@
 static double
-_sq_sc(double qa, double qb, double qc, double dnn, double d_factor)
+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.
     const double a1 = qa;
     const double a2 = qb;
     const double a3 = qc;
 
-    const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
-
     // Numerator: (1 - exp(a)^2)^3
     //         => (-(exp(2a) - 1))^3
     //         => -expm1(2a)^3
-    // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2)
-    //         => exp(a)^2 - 2 cos(xk) exp(a) + 1
-    //         => (exp(a) - 2 cos(xk)) * exp(a) + 1
-    const double exp_arg = exp(-arg);
-    const double Sq = -cube(expm1(-2.0*arg))
+    // 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 Sq;
+    return Zq;
 }
 
 // occupied volume fraction calculated from lattice symmetry and sphere radius
 static double
-_sc_volume_fraction(double radius, double dnn)
+sc_volume_fraction(double radius, double dnn)
 {
     return sphere_volume(radius/dnn);
@@ -65,6 +62,6 @@
             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;
+            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
@@ -72,8 +69,8 @@
     }
     outer_sum *= theta_m;
-    const double Sq = outer_sum/M_PI_2;
+    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 * Sq;
+    return sc_volume_fraction(radius, dnn) * Pq * Zq;
 }
 
@@ -92,5 +89,5 @@
     q = sqrt(qa*qa + qb*qb + qc*qc);
     const double Pq = sphere_form(q, radius, sld, solvent_sld);
-    const double Sq = _sq_sc(qa, qb, qc, dnn, d_factor);
-    return _sc_volume_fraction(radius, dnn) * Pq * Sq;
+    const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor);
+    return sc_volume_fraction(radius, dnn) * Pq * Zq;
 }
