Index: sasmodels/guyou.py
===================================================================
--- sasmodels/guyou.py	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
+++ sasmodels/guyou.py	(revision b3703f51bafe67afb8acefe986ad7466720f755b)
@@ -195,9 +195,10 @@
         plt.plot(x, y, 'g')
 
-    for long in range(-limit, limit+1, step):
-        x, y = guyou(scale*long, scale*long_line)
+    for longitude in range(-limit, limit+1, step):
+        x, y = guyou(scale*longitude, scale*long_line)
         plt.plot(x, y, 'b')
     #plt.xlabel('longitude')
     plt.ylabel('latitude')
+    plt.title('forward transform')
 
     plt.subplot(212)
@@ -211,6 +212,8 @@
     plt.xlabel('longitude')
     plt.ylabel('latitude')
+    plt.title('inverse transform')
 
 def main():
+    """Show the Guyou transformation"""
     plot_grid()
     import matplotlib.pyplot as plt
Index: sasmodels/jitter.py
===================================================================
--- sasmodels/jitter.py	(revision 802c41208e9f4bf5067b84d7ec5a9b78ceafaa69)
+++ sasmodels/jitter.py	(revision b3703f51bafe67afb8acefe986ad7466720f755b)
@@ -99,5 +99,5 @@
     _draw_crystal(axes, size, view, jitter, atoms=atoms)
 
-def _draw_crystal(axes, size, view, jitter, steps=None, alpha=1, atoms=None):
+def _draw_crystal(axes, size, view, jitter, atoms=None):
     atoms, size = np.asarray(atoms, 'd').T, np.asarray(size, 'd')
     x, y, z = atoms*size[:, None]
@@ -298,12 +298,12 @@
     # TODO: try Kent distribution instead of a gaussian warped by projection
 
-    t = np.linspace(-1, 1, n)
-    weights = np.ones_like(t)
+    dist_x = np.linspace(-1, 1, n)
+    weights = np.ones_like(dist_x)
     if dist == 'gaussian':
-        t *= 3
-        weights = exp(-0.5*t**2)
+        dist_x *= 3
+        weights = exp(-0.5*dist_x**2)
     elif dist == 'rectangle':
         # Note: uses sasmodels ridiculous definition of rectangle width
-        t *= sqrt(3)
+        dist_x *= sqrt(3)
     elif dist == 'uniform':
         pass
@@ -392,13 +392,13 @@
     z = np.matrix([[0], [0], [radius]])
     points = np.hstack([_rotate(theta_i, phi_j)*z
-                        for theta_i in dtheta*t
-                        for phi_j in dphi*t])
-    w = np.array([_weight(theta_i, phi_j, w_i, w_j)
-                  for w_i, theta_i in zip(weights, dtheta*t)
-                  for w_j, phi_j in zip(weights, dphi*t)])
-    #print(max(w), min(w), min(w[w>0]))
-    points = points[:, w > 0]
-    w = w[w > 0]
-    w /= max(w)
+                        for theta_i in dtheta*dist_x
+                        for phi_j in dphi*dist_x])
+    dist_w = np.array([_weight(theta_i, phi_j, w_i, w_j)
+                       for w_i, theta_i in zip(weights, dtheta*dist_x)
+                       for w_j, phi_j in zip(weights, dphi*dist_x)])
+    #print(max(dist_w), min(dist_w), min(dist_w[dist_w > 0]))
+    points = points[:, dist_w > 0]
+    dist_w = dist_w[dist_w > 0]
+    dist_w /= max(dist_w)
 
     # rotate relative to beam
@@ -407,5 +407,5 @@
     x, y, z = [np.array(v).flatten() for v in points]
     #plt.figure(2); plt.clf(); plt.hist(z, bins=np.linspace(-1, 1, 51))
-    axes.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.)
+    axes.scatter(x, y, z, c=dist_w, marker='o', vmin=0., vmax=1.)
 
 def draw_labels(axes, view, jitter, text):
@@ -430,24 +430,24 @@
     """Construct a matrix to rotate points about *x* by *angle* degrees."""
     angle = radians(angle)
-    R = [[1, 0, 0],
-         [0, +cos(angle), -sin(angle)],
-         [0, +sin(angle), +cos(angle)]]
-    return np.matrix(R)
+    rot = [[1, 0, 0],
+           [0, +cos(angle), -sin(angle)],
+           [0, +sin(angle), +cos(angle)]]
+    return np.matrix(rot)
 
 def Ry(angle):
     """Construct a matrix to rotate points about *y* by *angle* degrees."""
     angle = radians(angle)
-    R = [[+cos(angle), 0, +sin(angle)],
-         [0, 1, 0],
-         [-sin(angle), 0, +cos(angle)]]
-    return np.matrix(R)
+    rot = [[+cos(angle), 0, +sin(angle)],
+           [0, 1, 0],
+           [-sin(angle), 0, +cos(angle)]]
+    return np.matrix(rot)
 
 def Rz(angle):
     """Construct a matrix to rotate points about *z* by *angle* degrees."""
     angle = radians(angle)
-    R = [[+cos(angle), -sin(angle), 0],
-         [+sin(angle), +cos(angle), 0],
-         [0, 0, 1]]
-    return np.matrix(R)
+    rot = [[+cos(angle), -sin(angle), 0],
+           [+sin(angle), +cos(angle), 0],
+           [0, 0, 1]]
+    return np.matrix(rot)
 
 def transform_xyz(view, jitter, x, y, z):
@@ -564,5 +564,5 @@
     elif 1:
         axes.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1,
-                    levels=np.linspace(vmin, vmax, 24))
+                      levels=np.linspace(vmin, vmax, 24))
     else:
         axes.pcolormesh(qx, qy, Iqxy)
@@ -583,10 +583,10 @@
     for details.
     """
-    from sasmodels.core import load_model_info, build_model
+    from sasmodels.core import load_model_info, build_model as build_sasmodel
     from sasmodels.data import empty_data2D
     from sasmodels.direct_model import DirectModel
 
     model_info = load_model_info(model_name)
-    model = build_model(model_info) #, dtype='double!')
+    model = build_sasmodel(model_info) #, dtype='double!')
     q = np.linspace(-qmax, qmax, n)
     data = empty_data2D(q, q)
Index: sasmodels/multiscat.py
===================================================================
--- sasmodels/multiscat.py	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
+++ sasmodels/multiscat.py	(revision b3703f51bafe67afb8acefe986ad7466720f755b)
@@ -73,5 +73,4 @@
 import argparse
 import time
-import os.path
 
 import numpy as np
@@ -81,7 +80,5 @@
 from sasmodels import core
 from sasmodels import compare
-from sasmodels import resolution2d
 from sasmodels.resolution import Resolution, bin_edges
-from sasmodels.data import empty_data1D, empty_data2D, plot_data
 from sasmodels.direct_model import call_kernel
 import sasmodels.kernelcl
@@ -106,9 +103,38 @@
 USE_FAST = True  # OpenCL faster, less accurate math
 
-class NumpyCalculator:
+class ICalculator:
+    """
+    Multiple scattering calculator
+    """
+    def fft(self, Iq):
+        """
+        Compute the forward FFT for an image, real -> complex.
+        """
+        raise NotImplementedError()
+
+    def ifft(self, Iq):
+        """
+        Compute the inverse FFT for an image, complex -> complex.
+        """
+        raise NotImplementedError()
+
+    def mulitple_scattering(self, Iq):
+        r"""
+        Compute multiple scattering for I(q) given scattering probability p.
+
+        Given a probability p of scattering with the thickness, the expected
+        number of scattering events, $\lambda$ is $-\log(1 - p)$, giving a
+        Poisson weighted sum of single, double, triple, etc. scattering patterns.
+        The number of patterns used is based on coverage (default 99%).
+        """
+        raise NotImplementedError()
+
+class NumpyCalculator(ICalculator):
+    """
+    Multiple scattering calculator using numpy fft.
+    """
     def __init__(self, dims=None, dtype=PRECISION):
         self.dtype = dtype
         self.complex_dtype = np.dtype('F') if dtype == np.dtype('f') else np.dtype('D')
-        pass
 
     def fft(self, Iq):
@@ -127,12 +153,4 @@
 
     def multiple_scattering(self, Iq, p, coverage=0.99):
-        r"""
-        Compute multiple scattering for I(q) given scattering probability p.
-
-        Given a probability p of scattering with the thickness, the expected
-        number of scattering events, $\lambda$ is $-\log(1 - p)$, giving a
-        Poisson weighted sum of single, double, triple, etc. scattering patterns.
-        The number of patterns used is based on coverage (default 99%).
-        """
         #t0 = time.time()
         coeffs = scattering_coeffs(p, coverage)
@@ -140,7 +158,7 @@
         scale = np.sum(Iq)
         frame = _forward_shift(Iq/scale, dtype=self.dtype)
-        F = np.fft.fft2(frame)
-        F_convolved = F * np.polyval(poly, F)
-        frame = np.fft.ifft2(F_convolved)
+        fourier_frame = np.fft.fft2(frame)
+        convolved = fourier_frame * np.polyval(poly, fourier_frame)
+        frame = np.fft.ifft2(convolved)
         result = scale * _inverse_shift(frame.real, dtype=self.dtype)
         #print("numpy multiscat time", time.time()-t0)
@@ -173,5 +191,8 @@
 """
 
-class OpenclCalculator(NumpyCalculator):
+class OpenclCalculator(ICalculator):
+    """
+    Multiple scattering calculator using OpenCL via pyfft.
+    """
     polyval1f = None
     polyval1d = None
@@ -180,5 +201,5 @@
         context = env.get_context(dtype)
         if dtype == np.dtype('f'):
-            if self.polyval1f is None:
+            if OpenclCalculator.polyval1f is None:
                 program = sasmodels.kernelcl.compile_model(
                     context, POLYVAL1_KERNEL, dtype, fast=USE_FAST)
@@ -187,7 +208,7 @@
             self.dtype = dtype
             self.complex_dtype = np.dtype('F')
-            self.polyval1 = self.polyval1f
+            self.polyval1 = OpenclCalculator.polyval1f
         else:
-            if self.polyval1d is None:
+            if OpenclCalculator.polyval1d is None:
                 program = sasmodels.kernelcl.compile_model(
                     context, POLYVAL1_KERNEL, dtype, fast=False)
@@ -196,5 +217,5 @@
             self.dtype = dtype
             self.complex_type = np.dtype('D')
-            self.polyval1 = self.polyval1d
+            self.polyval1 = OpenclCalculator.polyval1d
         self.queue = env.get_queue(dtype)
         self.plan = pyfft.cl.Plan(dims, queue=self.queue)
@@ -229,8 +250,8 @@
         gpu_poly = cl_array.to_device(self.queue, poly)
         self.plan.execute(gpu_data.data)
-        degree, n = poly.shape[0], frame.shape[0]*frame.shape[1]
+        degree, data_size= poly.shape[0], frame.shape[0]*frame.shape[1]
         self.polyval1(
-            self.queue, [n], None,
-            np.int32(degree), gpu_poly.data, np.int32(n), gpu_data.data)
+            self.queue, [data_size], None,
+            np.int32(degree), gpu_poly.data, np.int32(data_size), gpu_data.data)
         self.plan.execute(gpu_data.data, inverse=True)
         frame = gpu_data.get()
@@ -251,6 +272,6 @@
     """
     if transform is None:
-        nx, ny = Iq.shape
-        transform = Calculator(dims=(nx*2, ny*2), dtype=dtype)
+        n_x, n_y = Iq.shape
+        transform = Calculator(dims=(n_x*2, n_y*2), dtype=dtype)
     scale = np.sum(Iq)
     frame = _forward_shift(Iq/scale, dtype=dtype)
@@ -528,4 +549,7 @@
 def parse_pars(model, opts):
     # type: (ModelInfo, argparse.Namespace) -> Dict[str, float]
+    """
+    Parse par=val arguments from the command line.
+    """
 
     seed = np.random.randint(1000000) if opts.random and opts.seed < 0 else opts.seed
@@ -541,5 +565,6 @@
         'is2d': opts.is2d,
     }
-    pars, pars2 = compare.parse_pars(compare_opts)
+    # Note: sascomp allows comparison on a pair of models, so ignore the second.
+    pars, _ = compare.parse_pars(compare_opts)
     return pars
 
@@ -550,14 +575,24 @@
         formatter_class=argparse.ArgumentDefaultsHelpFormatter,
         )
-    parser.add_argument('-p', '--probability', type=float, default=0.1, help="scattering probability")
-    parser.add_argument('-n', '--nq', type=int, default=1024, help='number of mesh points')
-    parser.add_argument('-q', '--qmax', type=float, default=0.5, help='max q')
-    parser.add_argument('-w', '--window', type=float, default=2.0, help='q calc = q max * window')
-    parser.add_argument('-2', '--2d', dest='is2d', action='store_true', help='oriented sample')
-    parser.add_argument('-s', '--seed', default=-1, help='random pars with given seed')
-    parser.add_argument('-r', '--random', action='store_true', help='random pars with random seed')
-    parser.add_argument('-o', '--outfile', type=str, default="", help='random pars with random seed')
-    parser.add_argument('model', type=str, help='sas model name such as cylinder')
-    parser.add_argument('pars', type=str, nargs='*', help='model parameters such as radius=30')
+    parser.add_argument('-p', '--probability', type=float, default=0.1,
+                        help="scattering probability")
+    parser.add_argument('-n', '--nq', type=int, default=1024,
+                        help='number of mesh points')
+    parser.add_argument('-q', '--qmax', type=float, default=0.5,
+                        help='max q')
+    parser.add_argument('-w', '--window', type=float, default=2.0,
+                        help='q calc = q max * window')
+    parser.add_argument('-2', '--2d', dest='is2d', action='store_true',
+                        help='oriented sample')
+    parser.add_argument('-s', '--seed', default=-1,
+                        help='random pars with given seed')
+    parser.add_argument('-r', '--random', action='store_true',
+                        help='random pars with random seed')
+    parser.add_argument('-o', '--outfile', type=str, default="",
+                        help='random pars with random seed')
+    parser.add_argument('model', type=str,
+                        help='sas model name such as cylinder')
+    parser.add_argument('pars', type=str, nargs='*',
+                        help='model parameters such as radius=30')
     opts = parser.parse_args()
     assert opts.nq%2 == 0, "require even # points"
@@ -607,4 +642,8 @@
             plotxy((res._q_steps, res._q_steps), res.Iqxy+background)
             pylab.title("total scattering for p=%g" % probability)
+            if res.resolution is not None:
+                pylab.figure()
+                plotxy((res._q_steps, res._q_steps), result)
+                pylab.title("total scattering with resolution")
     else:
         q = res._q
@@ -624,4 +663,6 @@
             # Plot 1D pattern for partial scattering
             pylab.loglog(q, res.Iq+background, label="total for p=%g"%probability)
+            if res.resolution is not None:
+                pylab.loglog(q, result, label="total with dQ")
             #new_annulus = annular_average(res._radius, res.Iqxy, res._edges)
             #pylab.loglog(q, new_annulus+background, label="new total for p=%g"%probability)
