Index: doc/guide/pd/polydispersity.rst
===================================================================
--- doc/guide/pd/polydispersity.rst	(revision f4ae8c46b94ae66b74ad2c8cfe09975d2caa7227)
+++ doc/guide/pd/polydispersity.rst	(revision 29afc50fdd56895a3c14d9c1e09b002d087c956a)
@@ -28,9 +28,13 @@
 sigmas $N_\sigma$ to include from the tails of the distribution, and the
 number of points used to compute the average. The center of the distribution
-is set by the value of the model parameter.
-
-Volume parameters have polydispersity *PD* (not to be confused with a
-molecular weight distributions in polymer science), but orientation parameters
-use angular distributions of width $\sigma$.
+is set by the value of the model parameter. The meaning of a polydispersity 
+parameter *PD* (not to be confused with a molecular weight distributions 
+in polymer science) in a model depends on the type of parameter it is being 
+applied too.
+
+The distribution width applied to *volume* (ie, shape-describing) parameters 
+is relative to the center value such that $\sigma = \mathrm{PD} \cdot \bar x$. 
+However, the distribution width applied to *orientation* (ie, angle-describing) 
+parameters is just $\sigma = \mathrm{PD}$.
 
 $N_\sigma$ determines how far into the tails to evaluate the distribution,
@@ -69,4 +73,8 @@
 or angular orientations, use the Gaussian or Boltzmann distributions.
 
+If applying polydispersion to parameters describing angles, use the Uniform 
+distribution. Beware of using distributions that are always positive (eg, the 
+Lognormal) because angles can be negative!
+
 The array distribution allows a user-defined distribution to be applied.
 
@@ -215,8 +223,8 @@
 The polydispersity in sasmodels is given by
 
-.. math:: \text{PD} = p = \sigma / x_\text{med}
-
-The mean value of the distribution is given by $\bar x = \exp(\mu+ p^2/2)$
-and the peak value by $\max x = \exp(\mu - p^2)$.
+.. math:: \text{PD} = \sigma = p / x_\text{med}
+
+The mean value of the distribution is given by $\bar x = \exp(\mu+ \sigma^2/2)$
+and the peak value by $\max x = \exp(\mu - \sigma^2)$.
 
 The variance (the square of the standard deviation) of the *lognormal*
@@ -232,5 +240,5 @@
 .. figure:: pd_lognormal.jpg
 
-    Lognormal distribution.
+    Lognormal distribution for PD=0.1.
 
 For further information on the Lognormal distribution see:
@@ -334,2 +342,3 @@
 | 2017-05-08 Paul Kienzle
 | 2018-03-20 Steve King
+| 2018-04-04 Steve King
Index: explore/asymint.py
===================================================================
--- explore/asymint.py	(revision a1c32c237d9d08ea8ba721dbaf09ecf6b862c812)
+++ explore/asymint.py	(revision c4621699875d34774d863620c13cccb42524c8bf)
@@ -30,5 +30,5 @@
 import numpy as np
 import mpmath as mp
-from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10
+from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10, arccos
 from numpy.polynomial.legendre import leggauss
 from scipy.integrate import dblquad, simps, romb, romberg
@@ -41,4 +41,6 @@
 shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1]
 Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2]
+
+DTYPE = 'd'
 
 class MPenv:
@@ -72,5 +74,6 @@
     sas_sinx_x = staticmethod(sp.sas_sinx_x)
     pi = np.pi
-    mpf = staticmethod(float)
+    #mpf = staticmethod(float)
+    mpf = staticmethod(lambda x: np.array(x, DTYPE))
 
 SLD = 3
@@ -220,7 +223,11 @@
     #A, B, C = 4450, 14000, 47
     #A, B, C = 445, 140, 47  # integer for the sake of mpf
-    A, B, C = 6800, 114, 1380
-    DA, DB, DC = 2300, 21, 58
+    A, B, C = 114, 1380, 6800
+    DA, DB, DC = 21, 58, 2300
     SLDA, SLDB, SLDC = "5", "-0.3", "11.5"
+    ## default parameters from sasmodels
+    #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2
+    ## swap A-B-C to C-B-A
+    #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA
     #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3
     #SLD_SOLVENT,CONTRAST = 0, 4
@@ -233,4 +240,5 @@
         DA, DC = DC, DA
         SLDA, SLDC = SLDC, SLDA
+    #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC)
     NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC)
     NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv)
@@ -348,4 +356,25 @@
     return n**2, Iq
 
+def gauss_quad_usub(q, n=150, dtype=DTYPE):
+    """
+    Compute the integral using gaussian quadrature for n = 20, 76 or 150.
+
+    Use *u = sin theta* substitution, and restrict integration over a single
+    quadrant for shapes that are mirror symmetric about AB, AC and BC planes.
+
+    Note that this doesn't work for fcc/bcc paracrystals, which instead step
+    over the entire 4 pi surface uniformly in theta-phi.
+    """
+    z, w = leggauss(n)
+    cos_theta = 0.5 * (z + 1)
+    theta = arccos(cos_theta)
+    phi = pi/2*(0.5 * (z + 1))
+    Atheta, Aphi = np.meshgrid(theta, phi)
+    Aw = w[None, :] * w[:, None]
+    q, Atheta, Aphi, Aw = [np.asarray(v, dtype=dtype) for v in (q, Atheta, Aphi, Aw)]
+    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
+    Iq = np.sum(Zq*Aw)*0.25
+    return n**2, Iq
+
 def gridded_2d(q, n=300):
     """
@@ -395,4 +424,7 @@
     print("gauss-1025", *gauss_quad_2d(Q, n=1025))
     print("gauss-2049", *gauss_quad_2d(Q, n=2049))
+    print("gauss-20 usub", *gauss_quad_usub(Q, n=20))
+    print("gauss-76 usub", *gauss_quad_usub(Q, n=76))
+    print("gauss-150 usub", *gauss_quad_usub(Q, n=150))
     #gridded_2d(Q, n=2**8+1)
     gridded_2d(Q, n=2**10+1)
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
+++ sasmodels/compare.py	(revision 5770493cd59e15962cc1fb80022d0920b561cf63)
@@ -752,8 +752,8 @@
     """
     for k in range(opts['sets']):
-        if k > 1:
+        if k > 0:
             # print a separate seed for each dataset for better reproducibility
             new_seed = np.random.randint(1000000)
-            print("Set %d uses -random=%i"%(k+1, new_seed))
+            print("=== Set %d uses -random=%i ==="%(k+1, new_seed))
             np.random.seed(new_seed)
         opts['pars'] = parse_pars(opts, maxdim=maxdim)
@@ -762,4 +762,7 @@
         result = run_models(opts, verbose=True)
         if opts['plot']:
+            if opts['is2d'] and k > 0:
+                import matplotlib.pyplot as plt
+                plt.figure()
             limits = plot_models(opts, result, limits=limits, setnum=k)
         if opts['show_weights']:
@@ -1329,16 +1332,19 @@
 
     # Evaluate preset parameter expressions
+    # Note: need to replace ':' with '_' in parameter names and expressions
+    # in order to support math on magnetic parameters.
     context = MATH.copy()
     context['np'] = np
-    context.update(pars)
+    context.update((k.replace(':', '_'), v) for k, v in pars.items())
     context.update((k, v) for k, v in presets.items() if isinstance(v, float))
+    #for k,v in sorted(context.items()): print(k, v)
     for k, v in presets.items():
         if not isinstance(v, float) and not k.endswith('_type'):
-            presets[k] = eval(v, context)
+            presets[k] = eval(v.replace(':', '_'), context)
     context.update(presets)
-    context.update((k, v) for k, v in presets2.items() if isinstance(v, float))
+    context.update((k.replace(':', '_'), v) for k, v in presets2.items() if isinstance(v, float))
     for k, v in presets2.items():
         if not isinstance(v, float) and not k.endswith('_type'):
-            presets2[k] = eval(v, context)
+            presets2[k] = eval(v.replace(':', '_'), context)
 
     # update parameters with presets
Index: sasmodels/details.py
===================================================================
--- sasmodels/details.py	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
+++ sasmodels/details.py	(revision 885753a925cb7e3423bdd0db6db395557a6fb1c6)
@@ -243,5 +243,5 @@
     offset = np.cumsum(np.hstack((0, length)))
     call_details = make_details(kernel.info, length, offset[:-1], offset[-1])
-    # Pad value array to a 32 value boundaryd
+    # Pad value array to a 32 value boundary
     data_len = nvalues + 2*sum(len(v) for v in dispersity)
     extra = (32 - data_len%32)%32
@@ -250,4 +250,5 @@
     is_magnetic = convert_magnetism(kernel.info.parameters, data)
     #call_details.show()
+    #print("data", data)
     return call_details, data, is_magnetic
 
@@ -296,11 +297,10 @@
     mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues]
     mag = mag.reshape(-1, 3)
-    scale = mag[:, 0]
-    if np.any(scale):
+    if np.any(mag[:, 0] != 0.0):
+        M0 = mag[:, 0].copy()
         theta, phi = radians(mag[:, 1]), radians(mag[:, 2])
-        cos_theta = cos(theta)
-        mag[:, 0] = scale*cos_theta*cos(phi)  # mx
-        mag[:, 1] = scale*sin(theta)  # my
-        mag[:, 2] = -scale*cos_theta*sin(phi)  # mz
+        mag[:, 0] = +M0*cos(theta)*cos(phi)  # mx
+        mag[:, 1] = +M0*sin(theta) # my
+        mag[:, 2] = -M0*cos(theta)*sin(phi)  # mz
         return True
     else:
Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
+++ sasmodels/kernel_iq.c	(revision dc6f601ca60f7d630a484562e056a94ed48bae82)
@@ -79,19 +79,22 @@
 //     du * (m_sigma_y + 1j*m_sigma_z);
 // weights for spin crosssections: dd du real, ud real, uu, du imag, ud imag
-static void set_spin_weights(double in_spin, double out_spin, double spins[4])
+static void set_spin_weights(double in_spin, double out_spin, double weight[6])
 {
   in_spin = clip(in_spin, 0.0, 1.0);
   out_spin = clip(out_spin, 0.0, 1.0);
-  spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd
-  spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du real
-  spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud real
-  spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu
-  spins[4] = spins[1]; // du imag
-  spins[5] = spins[2]; // ud imag
+  // Note: sasview 3.1 scaled all slds by sqrt(weight) and assumed that
+  //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...)
+  // which is likely to be the case for simple models.
+  weight[0] = sqrt((1.0-in_spin) * (1.0-out_spin)); // dd
+  weight[1] = sqrt((1.0-in_spin) * out_spin);       // du.real
+  weight[2] = sqrt(in_spin * (1.0-out_spin));       // ud.real
+  weight[3] = sqrt(in_spin * out_spin);             // uu
+  weight[4] = weight[1]; // du.imag
+  weight[5] = weight[2]; // ud.imag
 }
 
 // Compute the magnetic sld
 static double mag_sld(
-  const unsigned int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag
+  const unsigned int xs, // 0=dd, 1=du.real, 2=ud.real, 3=uu, 4=du.imag, 5=ud.imag
   const double qx, const double qy,
   const double px, const double py,
@@ -106,16 +109,16 @@
       case 0: // uu => sld - D M_perpx
           return sld - px*perp;
-      case 1: // ud real => -D M_perpy
+      case 1: // ud.real => -D M_perpy
           return py*perp;
-      case 2: // du real => -D M_perpy
+      case 2: // du.real => -D M_perpy
           return py*perp;
-      case 3: // dd real => sld + D M_perpx
+      case 3: // dd => sld + D M_perpx
           return sld + px*perp;
     }
   } else {
     if (xs== 4) {
-      return -mz;  // ud imag => -D M_perpz
+      return -mz;  // du.imag => +D M_perpz
     } else { // index == 5
-      return mz;   // du imag => D M_perpz
+      return +mz;  // ud.imag => -D M_perpz
     }
   }
@@ -312,7 +315,7 @@
   //     up_angle = values[NUM_PARS+4];
   // TODO: could precompute more magnetism parameters before calling the kernel.
-  double spins[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill
+  double xs_weights[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill
   double cos_mspin, sin_mspin;
-  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], spins);
+  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], xs_weights);
   SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin);
 #endif // MAGNETIC
@@ -661,5 +664,5 @@
             // loop over uu, ud real, du real, dd, ud imag, du imag
             for (unsigned int xs=0; xs<6; xs++) {
-              const double xs_weight = spins[xs];
+              const double xs_weight = xs_weights[xs];
               if (xs_weight > 1.e-8) {
                 // Since the cross section weight is significant, set the slds
@@ -674,4 +677,6 @@
                   local_values.vector[sld_index] =
                     mag_sld(xs, qx, qy, px, py, values[sld_index+2], mx, my, mz);
+//if (q_index==0) printf("%d: (qx,qy)=(%g,%g) xs=%d sld%d=%g p=(%g,%g) m=(%g,%g,%g)\n",
+//  q_index, qx, qy, xs, sk, local_values.vector[sld_index], px, py, mx, my, mz);
                 }
                 scattering += xs_weight * CALL_KERNEL();
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision 3221de0ed9b3480163246ceb0503432126bd5134)
+++ sasmodels/sasview_model.py	(revision 7c3fb150f0b78da8ddc4f5fc61815995e60a2486)
@@ -593,8 +593,5 @@
             # Check whether we have a list of ndarrays [qx,qy]
             qx, qy = qdist
-            if not self._model_info.parameters.has_2d:
-                return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2))
-            else:
-                return self.calculate_Iq(qx, qy)
+            return self.calculate_Iq(qx, qy)
 
         elif isinstance(qdist, np.ndarray):
@@ -677,5 +674,6 @@
         call_details, values, is_magnetic = make_kernel_args(calculator, pairs)
         #call_details.show()
-        #print("pairs", pairs)
+        #print("================ parameters ==================")
+        #for p, v in zip(parameters.call_parameters, pairs): print(p.name, v[0])
         #for k, p in enumerate(self._model_info.parameters.call_parameters):
         #    print(k, p.name, *pairs[k])
@@ -721,5 +719,5 @@
 
     def set_dispersion(self, parameter, dispersion):
-        # type: (str, weights.Dispersion) -> Dict[str, Any]
+        # type: (str, weights.Dispersion) -> None
         """
         Set the dispersion object for a model parameter
@@ -744,5 +742,6 @@
             self.dispersion[parameter] = dispersion.get_pars()
         else:
-            raise ValueError("%r is not a dispersity or orientation parameter")
+            raise ValueError("%r is not a dispersity or orientation parameter"
+                             % parameter)
 
     def _dispersion_mesh(self):
@@ -871,6 +870,20 @@
     CylinderModel().evalDistribution([0.1, 0.1])
 
+def magnetic_demo():
+    Model = _make_standard_model('sphere')
+    model = Model()
+    model.setParam('M0:sld', 8)
+    q = np.linspace(-0.35, 0.35, 500)
+    qx, qy = np.meshgrid(q, q)
+    result = model.calculate_Iq(qx.flatten(), qy.flatten())
+    result = result.reshape(qx.shape)
+
+    import pylab
+    pylab.imshow(np.log(result + 0.001))
+    pylab.show()
+
 if __name__ == "__main__":
     print("cylinder(0.1,0.1)=%g"%test_cylinder())
+    #magnetic_demo()
     #test_product()
     #test_structure_factor()
