Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision 3221de0ed9b3480163246ceb0503432126bd5134)
+++ sasmodels/compare.py	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
@@ -718,5 +718,5 @@
     model = core.build_model(model_info, dtype=dtype, platform="ocl")
     calculator = DirectModel(data, model, cutoff=cutoff)
-    engine_type = calculator._model.__class__.__name__.replace('Model','').upper()
+    engine_type = calculator._model.__class__.__name__.replace('Model', '').upper()
     bits = calculator._model.dtype.itemsize*8
     precision = "fast" if getattr(calculator._model, 'fast', False) else str(bits)
@@ -1491,5 +1491,6 @@
             vmin, vmax = limits
             self.limits = vmax*1e-7, 1.3*vmax
-            import pylab; pylab.clf()
+            import pylab
+            pylab.clf()
             plot_models(self.opts, result, limits=self.limits)
 
Index: sasmodels/data.py
===================================================================
--- sasmodels/data.py	(revision f549e37e13404ccb975dec4ce01b1bd567aa14b4)
+++ sasmodels/data.py	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
@@ -706,5 +706,5 @@
     else:
         vmin_scaled, vmax_scaled = vmin, vmax
-    nx, ny = len(data.x_bins), len(data.y_bins)
+    #nx, ny = len(data.x_bins), len(data.y_bins)
     x_bins, y_bins, image = _build_matrix(data, plottable)
     plt.imshow(image,
@@ -772,5 +772,5 @@
         if loop >= max_loop:  # this protects never-ending loop
             break
-        image = _fillup_pixels(self, image=image, weights=weights)
+        image = _fillup_pixels(image=image, weights=weights)
         loop += 1
 
@@ -817,5 +817,5 @@
     return x_bins, y_bins
 
-def _fillup_pixels(self, image=None, weights=None):
+def _fillup_pixels(image=None, weights=None):
     """
     Fill z values of the empty cells of 2d image matrix
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision f9c8474501a646f958c5fcf104d9d611d53d05e9)
+++ sasmodels/generate.py	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
@@ -290,12 +290,11 @@
     import loops.
     """
-    if (info.source and any(lib.startswith('lib/gauss') for lib in info.source)):
-        import os.path
+    if info.source and any(lib.startswith('lib/gauss') for lib in info.source):
         from .gengauss import gengauss
-        path = os.path.join(MODEL_PATH, "lib", "gauss%d.c"%n)
-        if not os.path.exists(path):
+        path = joinpath(MODEL_PATH, "lib", "gauss%d.c"%n)
+        if not exists(path):
             gengauss(n, path)
         info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss')
-                        else lib for lib in info.source]
+                       else lib for lib in info.source]
 
 def format_units(units):
@@ -321,10 +320,10 @@
                      for w, h in zip(column_widths, PARTABLE_HEADERS)]
 
-    sep = " ".join("="*w for w in column_widths)
+    underbar = " ".join("="*w for w in column_widths)
     lines = [
-        sep,
+        underbar,
         " ".join("%-*s" % (w, h)
                  for w, h in zip(column_widths, PARTABLE_HEADERS)),
-        sep,
+        underbar,
         ]
     for p in pars:
@@ -335,5 +334,5 @@
             "%*g" % (column_widths[3], p.default),
             ]))
-    lines.append(sep)
+    lines.append(underbar)
     return "\n".join(lines)
 
@@ -613,6 +612,6 @@
     """
     spaces = " "*depth
-    sep = "\n" + spaces
-    return spaces + sep.join(s.split("\n"))
+    interline_separator = "\n" + spaces
+    return spaces + interline_separator.join(s.split("\n"))
 
 
@@ -620,4 +619,7 @@
 def load_template(filename):
     # type: (str) -> str
+    """
+    Load template file from sasmodels resource directory.
+    """
     path = joinpath(DATA_PATH, filename)
     mtime = getmtime(path)
Index: sasmodels/guyou.py
===================================================================
--- sasmodels/guyou.py	(revision 0d5a655d6ddebe6096fad050347df4796b956767)
+++ sasmodels/guyou.py	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
@@ -31,9 +31,17 @@
 #
 # 2017-11-01 Paul Kienzle
-# * converted to python, using degrees rather than radians
+# * converted to python, with degrees rather than radians
+"""
+Convert between latitude-longitude and Guyou map coordinates.
+"""
+
 from __future__ import division, print_function
 
 import numpy as np
-from numpy import sqrt, pi, tan, cos, sin, log, exp, arctan2 as atan2, sign, radians, degrees
+from numpy import sqrt, pi, tan, cos, sin, sign, radians, degrees
+from numpy import sinh, arctan as atan
+
+# scipy version of special functions
+from scipy.special import ellipj as ellipticJ, ellipkinc as ellipticF
 
 _ = """
@@ -59,19 +67,15 @@
 """
 
-# scipy version of special functions
-from scipy.special import ellipj as ellipticJ, ellipkinc as ellipticF
-from numpy import sinh, sign, arctan as atan
-
 def ellipticJi(u, v, m):
     scalar = np.isscalar(u) and np.isscalar(v) and np.isscalar(m)
     u, v, m = np.broadcast_arrays(u, v, m)
-    result = np.empty_like([u,u,u], 'D')
-    real = v==0
-    imag = u==0
+    result = np.empty_like([u, u, u], 'D')
+    real = (v == 0)
+    imag = (u == 0)
     mixed = ~(real|imag)
     result[:, real] = _ellipticJi_real(u[real], m[real])
     result[:, imag] = _ellipticJi_imag(v[imag], m[imag])
     result[:, mixed] = _ellipticJi(u[mixed], v[mixed], m[mixed])
-    return result[0,:] if scalar else result
+    return result[0, :] if scalar else result
 
 def _ellipticJi_real(u, m):
@@ -104,5 +108,6 @@
         result = np.empty_like(phi, 'D')
         index = (phi == 0)
-        result[index] = ellipticF(atan(sinh(abs(phi[index]))), 1-m[index]) * sign(psi[index])
+        result[index] = ellipticF(atan(sinh(abs(phi[index]))),
+                                  1-m[index]) * sign(psi[index])
         result[~index] = ellipticFi(phi[~index], psi[~index], m[~index])
         return result.reshape(1)[0] if scalar else result
@@ -117,8 +122,9 @@
     cotlambda2 = (-b + sqrt(b * b - 4 * c)) / 2
     re = ellipticF(atan(1 / sqrt(cotlambda2)), m) * sign(phi)
-    im = ellipticF(atan(sqrt(np.maximum(0,(cotlambda2 / cotphi2 - 1) / m))), 1 - m) * sign(psi)
+    im = ellipticF(atan(sqrt(np.maximum(0, (cotlambda2 / cotphi2 - 1) / m))),
+                   1 - m) * sign(psi)
     return re + 1j*im
 
-sqrt2 = sqrt(2)
+SQRT2 = sqrt(2)
 
 # [PAK] renamed k_ => cos_u, k => sin_u, k*k => sinsq_u to avoid k,K confusion
@@ -127,4 +133,5 @@
 # K = 3.165103454447431823666270142140819753058976299237578486994...
 def guyou(lam, phi):
+    """Transform from (latitude, longitude) to point (x, y)"""
     # [PAK] wrap into [-pi/2, pi/2] radians
     x, y = np.asarray(lam), np.asarray(phi)
@@ -135,5 +142,5 @@
 
     # Compute constant K
-    cos_u = (sqrt2 - 1) / (sqrt2 + 1)
+    cos_u = (SQRT2 - 1) / (SQRT2 + 1)
     sinsq_u = 1 - cos_u**2
     K = ellipticF(pi/2, sinsq_u)
@@ -144,5 +151,5 @@
     at = atan(r * (cos(lam) - 1j*sin(lam)))
     t = ellipticFi(at.real, at.imag, sinsq_u)
-    x, y = (-t.imag, sign(phi + (phi==0))*(0.5 * K - t.real))
+    x, y = (-t.imag, sign(phi + (phi == 0))*(0.5 * K - t.real))
 
     # [PAK] convert to degrees, and return to original tile
@@ -150,4 +157,5 @@
 
 def guyou_invert(x, y):
+    """Transform from point (x, y) on plot to (latitude, longitude)"""
     # [PAK] wrap into [-pi/2, pi/2] radians
     x, y = np.asarray(x), np.asarray(y)
@@ -158,5 +166,5 @@
 
     # compute constant K
-    cos_u = (sqrt2 - 1) / (sqrt2 + 1)
+    cos_u = (SQRT2 - 1) / (SQRT2 + 1)
     sinsq_u = 1 - cos_u**2
     K = ellipticF(pi/2, sinsq_u)
@@ -174,4 +182,5 @@
 
 def plot_grid():
+    """Plot the latitude-longitude grid for Guyou transform"""
     import matplotlib.pyplot as plt
     from numpy import linspace
@@ -203,8 +212,11 @@
     plt.ylabel('latitude')
 
+def main():
+    plot_grid()
+    import matplotlib.pyplot as plt
+    plt.show()
+
 if __name__ == "__main__":
-    plot_grid()
-    import matplotlib.pyplot as plt; plt.show()
-
+    main()
 
 _ = """
Index: sasmodels/jitter.py
===================================================================
--- sasmodels/jitter.py	(revision 199bd0739ec5cc5d9325cc9e50fc75c67119b351)
+++ sasmodels/jitter.py	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
@@ -7,7 +7,4 @@
 """
 from __future__ import division, print_function
-
-import sys, os
-sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__))))
 
 import argparse
@@ -50,5 +47,5 @@
 def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1):
     """Draw an ellipsoid."""
-    a,b,c = size
+    a, b, c = size
     u = np.linspace(0, 2 * np.pi, steps)
     v = np.linspace(0, np.pi, steps)
@@ -61,17 +58,19 @@
 
     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]),
+        ('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_sc(ax, size, view, jitter, steps=None, alpha=1):
+    """Draw points for simple cubic paracrystal"""
     atoms = _build_sc()
     _draw_crystal(ax, size, view, jitter, atoms=atoms)
 
 def draw_fcc(ax, size, view, jitter, steps=None, alpha=1):
+    """Draw points for face-centered cubic paracrystal"""
     # Build the simple cubic crystal
     atoms = _build_sc()
@@ -88,4 +87,5 @@
 
 def draw_bcc(ax, size, view, jitter, steps=None, alpha=1):
+    """Draw points for body-centered cubic paracrystal"""
     # Build the simple cubic crystal
     atoms = _build_sc()
@@ -127,16 +127,16 @@
     """Draw a parallelepiped."""
     a, b, c = size
-    x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1])
-    y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1])
-    z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1])
+    x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1])
+    y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1])
+    z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1])
     tri = np.array([
         # counter clockwise triangles
         # z: up/down, x: right/left, y: front/back
-        [0,1,2], [3,2,1], # top face
-        [6,5,4], [5,6,7], # bottom face
-        [0,2,6], [6,4,0], # right face
-        [1,5,7], [7,3,1], # left face
-        [2,3,6], [7,6,3], # front face
-        [4,1,0], [5,1,4], # back face
+        [0, 1, 2], [3, 2, 1], # top face
+        [6, 5, 4], [5, 6, 7], # bottom face
+        [0, 2, 6], [6, 4, 0], # right face
+        [1, 5, 7], [7, 3, 1], # left face
+        [2, 3, 6], [7, 6, 3], # front face
+        [4, 1, 0], [5, 1, 4], # back face
     ])
 
@@ -149,17 +149,17 @@
     # rotate that face.
     if 1:
-        x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1])
-        y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1])
-        z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1])
+        x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1])
+        y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1])
+        z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1])
         x, y, z = transform_xyz(view, jitter, x, y, abs(z)+0.001)
-        ax.plot_trisurf(x, y, triangles=tri, Z=z, color=[1,0.6,0.6], alpha=alpha)
+        ax.plot_trisurf(x, y, triangles=tri, Z=z, color=[1, 0.6, 0.6], alpha=alpha)
 
     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]),
+        ('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]),
     ])
 
@@ -187,30 +187,30 @@
     cloud = [
         [-1, -1, -1],
-        [-1, -1,  0],
-        [-1, -1,  1],
-        [-1,  0, -1],
-        [-1,  0,  0],
-        [-1,  0,  1],
-        [-1,  1, -1],
-        [-1,  1,  0],
-        [-1,  1,  1],
-        [ 0, -1, -1],
-        [ 0, -1,  0],
-        [ 0, -1,  1],
-        [ 0,  0, -1],
-        [ 0,  0,  0],
-        [ 0,  0,  1],
-        [ 0,  1, -1],
-        [ 0,  1,  0],
-        [ 0,  1,  1],
-        [ 1, -1, -1],
-        [ 1, -1,  0],
-        [ 1, -1,  1],
-        [ 1,  0, -1],
-        [ 1,  0,  0],
-        [ 1,  0,  1],
-        [ 1,  1, -1],
-        [ 1,  1,  0],
-        [ 1,  1,  1],
+        [-1, -1, +0],
+        [-1, -1, +1],
+        [-1, +0, -1],
+        [-1, +0, +0],
+        [-1, +0, +1],
+        [-1, +1, -1],
+        [-1, +1, +0],
+        [-1, +1, +1],
+        [+0, -1, -1],
+        [+0, -1, +0],
+        [+0, -1, +1],
+        [+0, +0, -1],
+        [+0, +0, +0],
+        [+0, +0, +1],
+        [+0, +1, -1],
+        [+0, +1, +0],
+        [+0, +1, +1],
+        [+1, -1, -1],
+        [+1, -1, +0],
+        [+1, -1, +1],
+        [+1, +0, -1],
+        [+1, +0, +0],
+        [+1, +0, +1],
+        [+1, +1, -1],
+        [+1, +1, +0],
+        [+1, +1, +1],
     ]
     dtheta, dphi, dpsi = jitter
@@ -228,5 +228,5 @@
     for v in 'xyz':
         a, b, c = size
-        lim = np.sqrt(a**2+b**2+c**2)
+        lim = np.sqrt(a**2 + b**2 + c**2)
         getattr(ax, 'set_'+v+'lim')([-lim, lim])
         getattr(ax, v+'axis').label.set_text(v)
@@ -297,7 +297,4 @@
         Should allow free movement in theta, but phi is distorted.
     """
-    theta, phi, psi = view
-    dtheta, dphi, dpsi = jitter
-
     t = np.linspace(-1, 1, n)
     weights = np.ones_like(t)
@@ -314,10 +311,10 @@
 
     if projection == 'equirectangular':  #define PROJECTION 1
-        def rotate(theta_i, phi_j):
+        def _rotate(theta_i, phi_j):
             return Rx(phi_j)*Ry(theta_i)
-        def weight(theta_i, phi_j, wi, wj):
+        def _weight(theta_i, phi_j, wi, wj):
             return wi*wj*abs(cos(radians(theta_i)))
     elif projection == 'sinusoidal':  #define PROJECTION 2
-        def rotate(theta_i, phi_j):
+        def _rotate(theta_i, phi_j):
             latitude = theta_i
             scale = cos(radians(latitude))
@@ -325,5 +322,5 @@
             #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
             return Rx(longitude)*Ry(latitude)
-        def weight(theta_i, phi_j, wi, wj):
+        def _weight(theta_i, phi_j, wi, wj):
             latitude = theta_i
             scale = cos(radians(latitude))
@@ -331,18 +328,18 @@
             return w*wi*wj
     elif projection == 'guyou':  #define PROJECTION 3  (eventually?)
-        def rotate(theta_i, phi_j):
+        def _rotate(theta_i, phi_j):
             from guyou import guyou_invert
             #latitude, longitude = guyou_invert([theta_i], [phi_j])
             longitude, latitude = guyou_invert([phi_j], [theta_i])
             return Rx(longitude[0])*Ry(latitude[0])
-        def weight(theta_i, phi_j, wi, wj):
+        def _weight(theta_i, phi_j, wi, wj):
             return wi*wj
     elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry
-        def rotate(theta_i, phi_j):
+        def _rotate(theta_i, phi_j):
             latitude = sqrt(theta_i**2 + phi_j**2)
             longitude = degrees(np.arctan2(phi_j, theta_i))
             #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
             return Rz(longitude)*Ry(latitude)
-        def weight(theta_i, phi_j, wi, wj):
+        def _weight(theta_i, phi_j, wi, wj):
             # Weighting for each point comes from the integral:
             #     \int\int I(q, lat, log) sin(lat) dlat dlog
@@ -377,5 +374,5 @@
             return w*wi*wj if latitude < 180 else 0
     elif projection == 'azimuthal_equal_area':
-        def rotate(theta_i, phi_j):
+        def _rotate(theta_i, phi_j):
             R = min(1, sqrt(theta_i**2 + phi_j**2)/180)
             latitude = 180-degrees(2*np.arccos(R))
@@ -383,5 +380,5 @@
             #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
             return Rz(longitude)*Ry(latitude)
-        def weight(theta_i, phi_j, wi, wj):
+        def _weight(theta_i, phi_j, wi, wj):
             latitude = sqrt(theta_i**2 + phi_j**2)
             w = sin(radians(latitude))/latitude if latitude != 0 else 1
@@ -391,15 +388,16 @@
 
     # mesh in theta, phi formed by rotating z
+    dtheta, dphi, dpsi = jitter
     z = np.matrix([[0], [0], [radius]])
-    points = np.hstack([rotate(theta_i, phi_j)*z
+    points = np.hstack([_rotate(theta_i, phi_j)*z
                         for theta_i in dtheta*t
                         for phi_j in dphi*t])
     # select just the active points (i.e., those with phi < 180
-    w = np.array([weight(theta_i, phi_j, wi, wj)
+    w = np.array([_weight(theta_i, phi_j, wi, wj)
                   for wi, theta_i in zip(weights, dtheta*t)
                   for wj, phi_j in zip(weights, dphi*t)])
     #print(max(w), min(w), min(w[w>0]))
-    points = points[:, w>0]
-    w = w[w>0]
+    points = points[:, w > 0]
+    w = w[w > 0]
     w /= max(w)
 
@@ -469,5 +467,5 @@
     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 = np.matrix([x.flatten(), y.flatten(), z.flatten()])
     points = apply_jitter(jitter, points)
     points = orient_relative_to_beam(view, points)
@@ -543,5 +541,5 @@
     theta, phi, psi = view
     theta_pd, phi_pd, psi_pd = [scale*v for v in jitter]
-    theta_pd_n, phi_pd_n, psi_pd_n = PD_N_TABLE[(theta_pd>0, phi_pd>0, psi_pd>0)]
+    theta_pd_n, phi_pd_n, psi_pd_n = PD_N_TABLE[(theta_pd > 0, phi_pd > 0, psi_pd > 0)]
     ## increase pd_n for testing jitter integration rather than simple viz
     #theta_pd_n, phi_pd_n, psi_pd_n = [5*v for v in (theta_pd_n, phi_pd_n, psi_pd_n)]
@@ -571,5 +569,5 @@
     if 0:
         level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i')
-        level[level<0] = 0
+        level[level < 0] = 0
         colors = plt.get_cmap()(level)
         ax.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors)
@@ -618,5 +616,5 @@
     return calculator
 
-def select_calculator(model_name, n=150, size=(10,40,100)):
+def select_calculator(model_name, n=150, size=(10, 40, 100)):
     """
     Create a model calculator for the given shape.
@@ -641,6 +639,6 @@
         radius = 0.5*c
         calculator = build_model('sc_paracrystal', n=n, dnn=dnn,
-                                  d_factor=d_factor, radius=(1-d_factor)*radius,
-                                  background=0)
+                                 d_factor=d_factor, radius=(1-d_factor)*radius,
+                                 background=0)
     elif model_name == 'fcc_paracrystal':
         a = b = c
@@ -650,6 +648,6 @@
         radius = sqrt(2)/4 * c
         calculator = build_model('fcc_paracrystal', n=n, dnn=dnn,
-                                  d_factor=d_factor, radius=(1-d_factor)*radius,
-                                  background=0)
+                                 d_factor=d_factor, radius=(1-d_factor)*radius,
+                                 background=0)
     elif model_name == 'bcc_paracrystal':
         a = b = c
@@ -659,6 +657,6 @@
         radius = sqrt(3)/2 * c
         calculator = build_model('bcc_paracrystal', n=n, dnn=dnn,
-                                  d_factor=d_factor, radius=(1-d_factor)*radius,
-                                  background=0)
+                                 d_factor=d_factor, radius=(1-d_factor)*radius,
+                                 background=0)
     elif model_name == 'cylinder':
         calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c)
@@ -685,5 +683,5 @@
     'cylinder',
     'fcc_paracrystal', 'bcc_paracrystal', 'sc_paracrystal',
- ]
+]
 
 DRAW_SHAPES = {
@@ -761,5 +759,5 @@
 
     ## add control widgets to plot
-    axtheta  = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor)
+    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)
     axpsi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor)
@@ -768,7 +766,7 @@
     spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi)
 
-    axdtheta  = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor)
+    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)
-    axdpsi= plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor)
+    axdpsi = plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor)
     # Note: using ridiculous definition of rectangle distribution, whose width
     # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep
@@ -797,5 +795,5 @@
 
     ## bind control widgets to view updater
-    stheta.on_changed(lambda v: update(v,'theta'))
+    stheta.on_changed(lambda v: update(v, 'theta'))
     sphi.on_changed(lambda v: update(v, 'phi'))
     spsi.on_changed(lambda v: update(v, 'psi'))
@@ -815,9 +813,16 @@
         formatter_class=argparse.ArgumentDefaultsHelpFormatter,
         )
-    parser.add_argument('-p', '--projection', choices=PROJECTIONS, default=PROJECTIONS[0], help='coordinate projection')
-    parser.add_argument('-s', '--size', type=str, default='10,40,100', help='a,b,c lengths')
-    parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, default=DISTRIBUTIONS[0], help='jitter distribution')
-    parser.add_argument('-m', '--mesh', type=int, default=30, help='#points in theta-phi mesh')
-    parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], help='oriented shape')
+    parser.add_argument('-p', '--projection', choices=PROJECTIONS,
+                        default=PROJECTIONS[0],
+                        help='coordinate projection')
+    parser.add_argument('-s', '--size', type=str, default='10,40,100',
+                        help='a,b,c lengths')
+    parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS,
+                        default=DISTRIBUTIONS[0],
+                        help='jitter distribution')
+    parser.add_argument('-m', '--mesh', type=int, default=30,
+                        help='#points in theta-phi mesh')
+    parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0],
+                        help='oriented shape')
     opts = parser.parse_args()
     size = tuple(int(v) for v in opts.size.split(','))
Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision aadec17f5a4db1c04d3dbefd23a3ac36b2dd5b96)
+++ sasmodels/kernel_iq.c	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
@@ -173,5 +173,5 @@
 // Apply the rotation matrix returned from qac_rotation to the point (qx,qy),
 // returning R*[qx,qy]' = [qa,qc]'
-static double
+static void
 qac_apply(
     QACRotation *rotation,
@@ -246,5 +246,5 @@
 // Apply the rotation matrix returned from qabc_rotation to the point (qx,qy),
 // returning R*[qx,qy]' = [qa,qb,qc]'
-static double
+static void
 qabc_apply(
     QABCRotation *rotation,
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision 6cbdcd4693191a58b26ef70549c1263ab01cd716)
+++ sasmodels/kernelcl.py	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
@@ -151,5 +151,5 @@
         if not HAVE_OPENCL:
             raise RuntimeError("OpenCL startup failed with ***"
-                            + OPENCL_ERROR + "***; using C compiler instead")
+                               + OPENCL_ERROR + "***; using C compiler instead")
         reset_environment()
         if ENV is None:
Index: sasmodels/models/core_shell_cylinder.c
===================================================================
--- sasmodels/models/core_shell_cylinder.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
+++ sasmodels/models/core_shell_cylinder.c	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
@@ -48,5 +48,6 @@
 
 
-double Iqac(double qab, double qc,
+static double
+Iqac(double qab, double qc,
     double core_sld,
     double shell_sld,
Index: sasmodels/models/hollow_rectangular_prism.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
+++ sasmodels/models/hollow_rectangular_prism.c	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
@@ -1,7 +1,4 @@
-double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness);
-double Iq(double q, double sld, double solvent_sld, double length_a,
-          double b2a_ratio, double c2a_ratio, double thickness);
-
-double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness)
+static double
+form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness)
 {
     double length_b = length_a * b2a_ratio;
@@ -16,5 +13,6 @@
 }
 
-double Iq(double q,
+static double
+Iq(double q,
     double sld,
     double solvent_sld,
@@ -85,5 +83,6 @@
 }
 
-double Iqabc(double qa, double qb, double qc,
+static double
+Iqabc(double qa, double qb, double qc,
     double sld,
     double solvent_sld,
Index: sasmodels/models/hollow_rectangular_prism_thin_walls.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism_thin_walls.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/hollow_rectangular_prism_thin_walls.c	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
@@ -1,7 +1,4 @@
-double form_volume(double length_a, double b2a_ratio, double c2a_ratio);
-double Iq(double q, double sld, double solvent_sld, double length_a,
-          double b2a_ratio, double c2a_ratio);
-
-double form_volume(double length_a, double b2a_ratio, double c2a_ratio)
+static double
+form_volume(double length_a, double b2a_ratio, double c2a_ratio)
 {
     double length_b = length_a * b2a_ratio;
@@ -11,5 +8,6 @@
 }
 
-double Iq(double q,
+static double
+Iq(double q,
     double sld,
     double solvent_sld,
Index: sasmodels/models/rectangular_prism.c
===================================================================
--- sasmodels/models/rectangular_prism.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
+++ sasmodels/models/rectangular_prism.c	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
@@ -1,12 +1,10 @@
-double form_volume(double length_a, double b2a_ratio, double c2a_ratio);
-double Iq(double q, double sld, double solvent_sld, double length_a,
-          double b2a_ratio, double c2a_ratio);
-
-double form_volume(double length_a, double b2a_ratio, double c2a_ratio)
+static double
+form_volume(double length_a, double b2a_ratio, double c2a_ratio)
 {
     return length_a * (length_a*b2a_ratio) * (length_a*c2a_ratio);
 }
 
-double Iq(double q,
+static double
+Iq(double q,
     double sld,
     double solvent_sld,
@@ -71,5 +69,6 @@
 
 
-double Iqabc(double qa, double qb, double qc,
+static double
+Iqabc(double qa, double qb, double qc,
     double sld,
     double solvent_sld,
Index: sasmodels/multiscat.py
===================================================================
--- sasmodels/multiscat.py	(revision 49d1f8b8b50df60b60d3cc8b9dfe3c9c22f24db1)
+++ sasmodels/multiscat.py	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
@@ -261,4 +261,10 @@
 
 def scattering_coeffs(p, coverage=0.99):
+    r"""
+    Return the coefficients of the scattering powers for transmission
+    probability *p*.  This is just the corresponding values for the
+    Poisson distribution for $\lambda = -\ln(1-p)$ such that
+    $\sum_{k = 0 \ldots n} P(k; \lambda)$ is larger than *coverage*.
+    """
     L = -np.log(1-p)
     num_scatter = truncated_poisson_invcdf(coverage, L)
@@ -266,7 +272,7 @@
     return coeffs
 
-def truncated_poisson_invcdf(p, L):
+def truncated_poisson_invcdf(coverage, L):
     r"""
-    Return smallest k such that cdf(k; L) > p from the truncated Poisson
+    Return smallest k such that cdf(k; L) > coverage from the truncated Poisson
     probability excluding k=0
     """
@@ -275,5 +281,5 @@
     pmf = -np.exp(-L) / np.expm1(-L)
     k = 0
-    while cdf < p:
+    while cdf < coverage:
         k += 1
         pmf *= L/k
@@ -305,32 +311,32 @@
 
 class MultipleScattering(Resolution):
+    r"""
+    Compute multiple scattering using Fourier convolution.
+
+    The fourier steps are determined by *qmax*, the maximum $q$ value
+    desired, *nq* the number of $q$ steps and *window*, the amount
+    of padding around the circular convolution.  The $q$ spacing
+    will be $\Delta q = 2 q_\mathrm{max} w / n_q$.  If *nq* is not
+    given it will use $n_q = 2^k$ such that $\Delta q < q_\mathrm{min}$.
+
+    *probability* is related to the expected number of scattering
+    events in the sample $\lambda$ as $p = 1 = e^{-\lambda}$.  As a
+    hack to allow probability to be a fitted parameter, the "value"
+    can be a function that takes no parameters and returns the current
+    value of the probability.  *coverage* determines how many scattering
+    steps to consider.  The default is 0.99, which sets $n$ such that
+    $1 \ldots n$ covers 99% of the Poisson probability mass function.
+
+    *is2d* is True then 2D scattering is used, otherwise it accepts
+    and returns 1D scattering.
+
+    *resolution* is the resolution function to apply after multiple
+    scattering.  If present, then the resolution $q$ vectors will provide
+    default values for *qmin*, *qmax* and *nq*.
+    """
     def __init__(self, qmin=None, qmax=None, nq=None, window=2,
                  probability=None, coverage=0.99,
                  is2d=False, resolution=None,
                  dtype=PRECISION):
-        r"""
-        Compute multiple scattering using Fourier convolution.
-
-        The fourier steps are determined by *qmax*, the maximum $q$ value
-        desired, *nq* the number of $q$ steps and *window*, the amount
-        of padding around the circular convolution.  The $q$ spacing
-        will be $\Delta q = 2 q_\mathrm{max} w / n_q$.  If *nq* is not
-        given it will use $n_q = 2^k$ such that $\Delta q < q_\mathrm{min}$.
-
-        *probability* is related to the expected number of scattering
-        events in the sample $\lambda$ as $p = 1 = e^{-\lambda}$.  As a
-        hack to allow probability to be a fitted parameter, the "value"
-        can be a function that takes no parameters and returns the current
-        value of the probability.  *coverage* determines how many scattering
-        steps to consider.  The default is 0.99, which sets $n$ such that
-        $1 \ldots n$ covers 99% of the Poisson probability mass function.
-
-        *is2d* is True then 2D scattering is used, otherwise it accepts
-        and returns 1D scattering.
-
-        *resolution* is the resolution function to apply after multiple
-        scattering.  If present, then the resolution $q$ vectors will provide
-        default values for *qmin*, *qmax* and *nq*.
-        """
         # Infer qmin, qmax from instrument resolution calculator, if present
         if resolution is not None:
@@ -424,4 +430,8 @@
         # Prepare the multiple scattering calculator (either numpy or OpenCL)
         self.transform = Calculator((2*nq, 2*nq), dtype=dtype)
+
+        # Iq and Iqxy will be set during apply
+        self.Iq = None # type: np.ndarray
+        self.Iqxy = None # type: np.ndarray
 
     def apply(self, theory):
@@ -471,4 +481,8 @@
 
     def radial_profile(self, Iqxy):
+        """
+        Compute that radial profile for the given Iqxy grid.  The grid should
+        be defined as for
+        """
         # circular average, no anti-aliasing
         Iq = np.histogram(self._radius, bins=self._edges, weights=Iqxy)[0]/self._norm
@@ -478,5 +492,6 @@
 def annular_average(qxy, Iqxy, qbins):
     """
-    Compute annular average of points at
+    Compute annular average of points in *Iqxy* at *qbins*.  The $q_x$, $q_y$
+    coordinates for *Iqxy* are given in *qxy*.
     """
     qxy, Iqxy = qxy.flatten(), Iqxy.flatten()
