Index: sasmodels/jitter.py
===================================================================
--- sasmodels/jitter.py	(revision d86f0fc0fff835c5ddf398e64c046518fa4536fc)
+++ sasmodels/jitter.py	(revision 802c41208e9f4bf5067b84d7ec5a9b78ceafaa69)
@@ -16,15 +16,14 @@
 
 import matplotlib.pyplot as plt
-from matplotlib.widgets import Slider, CheckButtons
-from matplotlib import cm
+from matplotlib.widgets import Slider
 import numpy as np
 from numpy import pi, cos, sin, sqrt, exp, degrees, radians
 
-def draw_beam(ax, view=(0, 0)):
+def draw_beam(axes, view=(0, 0)):
     """
     Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1)
     """
-    #ax.plot([0,0],[0,0],[1,-1])
-    #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8)
+    #axes.plot([0,0],[0,0],[1,-1])
+    #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8)
 
     steps = 25
@@ -43,7 +42,7 @@
     x, y, z = [v.reshape(shape) for v in points]
 
-    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5)
-
-def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1):
+    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5)
+
+def draw_ellipsoid(axes, size, view, jitter, steps=25, alpha=1):
     """Draw an ellipsoid."""
     a, b, c = size
@@ -55,7 +54,7 @@
     x, y, z = transform_xyz(view, jitter, x, y, z)
 
-    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha)
-
-    draw_labels(ax, view, jitter, [
+    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha)
+
+    draw_labels(axes, view, jitter, [
         ('c+', [+0, +0, +c], [+1, +0, +0]),
         ('c-', [+0, +0, -c], [+0, +0, -1]),
@@ -66,10 +65,10 @@
     ])
 
-def draw_sc(ax, size, view, jitter, steps=None, alpha=1):
+def draw_sc(axes, 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_crystal(axes, size, view, jitter, atoms=atoms)
+
+def draw_fcc(axes, size, view, jitter, steps=None, alpha=1):
     """Draw points for face-centered cubic paracrystal"""
     # Build the simple cubic crystal
@@ -84,7 +83,7 @@
     # y and z planes can be generated by substituting x for y and z respectively
     atoms.extend(zip(x+y+z, y+z+x, z+x+y))
-    _draw_crystal(ax, size, view, jitter, atoms=atoms)
-
-def draw_bcc(ax, size, view, jitter, steps=None, alpha=1):
+    _draw_crystal(axes, size, view, jitter, atoms=atoms)
+
+def draw_bcc(axes, size, view, jitter, steps=None, alpha=1):
     """Draw points for body-centered cubic paracrystal"""
     # Build the simple cubic crystal
@@ -98,12 +97,12 @@
     )
     atoms.extend(zip(x, y, z))
-    _draw_crystal(ax, size, view, jitter, atoms=atoms)
-
-def _draw_crystal(ax, size, view, jitter, steps=None, alpha=1, atoms=None):
+    _draw_crystal(axes, size, view, jitter, atoms=atoms)
+
+def _draw_crystal(axes, size, view, jitter, steps=None, alpha=1, atoms=None):
     atoms, size = np.asarray(atoms, 'd').T, np.asarray(size, 'd')
     x, y, z = atoms*size[:, None]
     x, y, z = transform_xyz(view, jitter, x, y, z)
-    ax.scatter([x[0]], [y[0]], [z[0]], c='yellow', marker='o')
-    ax.scatter(x[1:], y[1:], z[1:], c='r', marker='o')
+    axes.scatter([x[0]], [y[0]], [z[0]], c='yellow', marker='o')
+    axes.scatter(x[1:], y[1:], z[1:], c='r', marker='o')
 
 def _build_sc():
@@ -124,5 +123,5 @@
     return atoms
 
-def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1):
+def draw_parallelepiped(axes, size, view, jitter, steps=None, alpha=1):
     """Draw a parallelepiped."""
     a, b, c = size
@@ -142,5 +141,5 @@
 
     x, y, z = transform_xyz(view, jitter, x, y, z)
-    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha)
+    axes.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha)
 
     # Draw pink face on box.
@@ -153,7 +152,7 @@
         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)
-
-    draw_labels(ax, view, jitter, [
+        axes.plot_trisurf(x, y, triangles=tri, Z=z, color=[1, 0.6, 0.6], alpha=alpha)
+
+    draw_labels(axes, view, jitter, [
         ('c+', [+0, +0, +c], [+1, +0, +0]),
         ('c-', [+0, +0, -c], [+0, +0, -1]),
@@ -164,5 +163,5 @@
     ])
 
-def draw_sphere(ax, radius=10., steps=100):
+def draw_sphere(axes, radius=10., steps=100):
     """Draw a sphere"""
     u = np.linspace(0, 2 * np.pi, steps)
@@ -172,7 +171,7 @@
     y = radius * np.outer(np.sin(u), np.sin(v))
     z = radius * np.outer(np.ones(np.size(u)), np.cos(v))
-    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w')
-
-def draw_jitter(ax, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0),
+    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w')
+
+def draw_jitter(axes, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0),
                 draw_shape=draw_parallelepiped):
     """
@@ -221,14 +220,14 @@
     if dpsi == 0:
         cloud = [v for v in cloud if v[2] == 0]
-    draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8)
+    draw_shape(axes, size, view, [0, 0, 0], steps=100, alpha=0.8)
     scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist]
     for point in cloud:
         delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]]
-        draw_shape(ax, size, view, delta, alpha=0.8)
+        draw_shape(axes, size, view, delta, alpha=0.8)
     for v in 'xyz':
         a, b, c = size
         lim = np.sqrt(a**2 + b**2 + c**2)
-        getattr(ax, 'set_'+v+'lim')([-lim, lim])
-        getattr(ax, v+'axis').label.set_text(v)
+        getattr(axes, 'set_'+v+'lim')([-lim, lim])
+        getattr(axes, v+'axis').label.set_text(v)
 
 PROJECTIONS = [
@@ -238,5 +237,5 @@
     'azimuthal_equal_area',
 ]
-def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian',
+def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian',
               projection='equirectangular'):
     """
@@ -297,4 +296,6 @@
         Should allow free movement in theta, but phi is distorted.
     """
+    # TODO: try Kent distribution instead of a gaussian warped by projection
+
     t = np.linspace(-1, 1, n)
     weights = np.ones_like(t)
@@ -313,6 +314,6 @@
         def _rotate(theta_i, phi_j):
             return Rx(phi_j)*Ry(theta_i)
-        def _weight(theta_i, phi_j, wi, wj):
-            return wi*wj*abs(cos(radians(theta_i)))
+        def _weight(theta_i, phi_j, w_i, w_j):
+            return w_i*w_j*abs(cos(radians(theta_i)))
     elif projection == 'sinusoidal':  #define PROJECTION 2
         def _rotate(theta_i, phi_j):
@@ -322,17 +323,17 @@
             #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, w_i, w_j):
             latitude = theta_i
             scale = cos(radians(latitude))
-            w = 1 if abs(phi_j) < abs(scale)*180 else 0
-            return w*wi*wj
+            active = 1 if abs(phi_j) < abs(scale)*180 else 0
+            return active*w_i*w_j
     elif projection == 'guyou':  #define PROJECTION 3  (eventually?)
         def _rotate(theta_i, phi_j):
-            from guyou import guyou_invert
+            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):
-            return wi*wj
+        def _weight(theta_i, phi_j, w_i, w_j):
+            return w_i*w_j
     elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry
         def _rotate(theta_i, phi_j):
@@ -341,5 +342,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, w_i, w_j):
             # Weighting for each point comes from the integral:
             #     \int\int I(q, lat, log) sin(lat) dlat dlog
@@ -371,17 +372,17 @@
             # the entire sphere, and treats theta and phi identically.
             latitude = sqrt(theta_i**2 + phi_j**2)
-            w = sin(radians(latitude))/latitude if latitude != 0 else 1
-            return w*wi*wj if latitude < 180 else 0
+            weight = sin(radians(latitude))/latitude if latitude != 0 else 1
+            return weight*w_i*w_j if latitude < 180 else 0
     elif projection == 'azimuthal_equal_area':
         def _rotate(theta_i, phi_j):
-            R = min(1, sqrt(theta_i**2 + phi_j**2)/180)
-            latitude = 180-degrees(2*np.arccos(R))
+            radius = min(1, sqrt(theta_i**2 + phi_j**2)/180)
+            latitude = 180-degrees(2*np.arccos(radius))
             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, w_i, w_j):
             latitude = sqrt(theta_i**2 + phi_j**2)
-            w = sin(radians(latitude))/latitude if latitude != 0 else 1
-            return w*wi*wj if latitude < 180 else 0
+            weight = sin(radians(latitude))/latitude if latitude != 0 else 1
+            return weight*w_i*w_j if latitude < 180 else 0
     else:
         raise ValueError("unknown projection %r"%projection)
@@ -393,8 +394,7 @@
                         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)
-                  for wi, theta_i in zip(weights, dtheta*t)
-                  for wj, phi_j in zip(weights, 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]
@@ -402,14 +402,4 @@
     w /= max(w)
 
-    if 0: # Kent distribution
-        points = np.hstack([Rx(phi_j)*Ry(theta_i)*z for theta_i in 30*t for phi_j in 60*t])
-        xp, yp, zp = [np.array(v).flatten() for v in points]
-        kappa = max(1e6, radians(dtheta)/(2*pi))
-        beta = 1/max(1e-6, radians(dphi)/(2*pi))/kappa
-        w = exp(kappa*zp) #+ beta*(xp**2 + yp**2)
-        print(kappa, dtheta, radians(dtheta), min(w), max(w), sum(w))
-        #w /= abs(cos(radians(
-        #w /= sum(w)
-
     # rotate relative to beam
     points = orient_relative_to_beam(view, points)
@@ -417,7 +407,7 @@
     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))
-    ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.)
-
-def draw_labels(ax, view, jitter, text):
+    axes.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.)
+
+def draw_labels(axes, view, jitter, text):
     """
     Draw text at a particular location.
@@ -433,5 +423,5 @@
     for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)):
         zdir = np.asarray(zdir).flatten()
-        ax.text(p[0], p[1], p[2], label, zdir=zdir)
+        axes.text(p[0], p[1], p[2], label, zdir=zdir)
 
 # Definition of rotation matrices comes from wikipedia:
@@ -439,23 +429,23 @@
 def Rx(angle):
     """Construct a matrix to rotate points about *x* by *angle* degrees."""
-    a = radians(angle)
+    angle = radians(angle)
     R = [[1, 0, 0],
-         [0, +cos(a), -sin(a)],
-         [0, +sin(a), +cos(a)]]
+         [0, +cos(angle), -sin(angle)],
+         [0, +sin(angle), +cos(angle)]]
     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)],
+    angle = radians(angle)
+    R = [[+cos(angle), 0, +sin(angle)],
          [0, 1, 0],
-         [-sin(a), 0, +cos(a)]]
+         [-sin(angle), 0, +cos(angle)]]
     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],
+    angle = radians(angle)
+    R = [[+cos(angle), -sin(angle), 0],
+         [+sin(angle), +cos(angle), 0],
          [0, 0, 1]]
     return np.matrix(R)
@@ -524,9 +514,9 @@
         return data[offset], data[-1]
 
-def draw_scattering(calculator, ax, view, jitter, dist='gaussian'):
+def draw_scattering(calculator, axes, view, jitter, dist='gaussian'):
     """
     Plot the scattering for the particular view.
 
-    *calculator* is returned from :func:`build_model`.  *ax* are the 3D axes
+    *calculator* is returned from :func:`build_model`.  *axes* are the 3D axes
     on which the data will be plotted.  *view* and *jitter* are the current
     orientation and orientation dispersity.  *dist* is one of the sasmodels
@@ -571,10 +561,10 @@
         level[level < 0] = 0
         colors = plt.get_cmap()(level)
-        ax.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors)
+        axes.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors)
     elif 1:
-        ax.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1,
+        axes.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1,
                     levels=np.linspace(vmin, vmax, 24))
     else:
-        ax.pcolormesh(qx, qy, Iqxy)
+        axes.pcolormesh(qx, qy, Iqxy)
 
 def build_model(model_name, n=150, qmax=0.5, **pars):
@@ -749,8 +739,8 @@
     plt.gcf().canvas.set_window_title(projection)
     #gs = gridspec.GridSpec(2,1,height_ratios=[4,1])
-    #ax = plt.subplot(gs[0], projection='3d')
-    ax = plt.axes([0.0, 0.2, 1.0, 0.8], projection='3d')
+    #axes = plt.subplot(gs[0], projection='3d')
+    axes = plt.axes([0.0, 0.2, 1.0, 0.8], projection='3d')
     try:  # CRUFT: not all versions of matplotlib accept 'square' 3d projection
-        ax.axis('square')
+        axes.axis('square')
     except Exception:
         pass
@@ -759,21 +749,21 @@
 
     ## add control widgets to plot
-    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)
-    stheta = Slider(axtheta, 'Theta', -90, 90, valinit=theta)
-    sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi)
-    spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi)
-
-    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)
+    axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor)
+    axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor)
+    axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor)
+    stheta = Slider(axes_theta, 'Theta', -90, 90, valinit=theta)
+    sphi = Slider(axes_phi, 'Phi', -180, 180, valinit=phi)
+    spsi = Slider(axes_psi, 'Psi', -180, 180, valinit=psi)
+
+    axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor)
+    axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor)
+    axes_dpsi = 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
     # the maximum width to 90.
     dlimit = DIST_LIMITS[dist]
-    sdtheta = Slider(axdtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta)
-    sdphi = Slider(axdphi, 'dPhi', 0, 2*dlimit, valinit=dphi)
-    sdpsi = Slider(axdpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi)
+    sdtheta = Slider(axes_dtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta)
+    sdphi = Slider(axes_dphi, 'dPhi', 0, 2*dlimit, valinit=dphi)
+    sdpsi = Slider(axes_dpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi)
 
 
@@ -786,10 +776,10 @@
         limit = [0, 0, 0.5, 5][dims]
         jitter = [0 if v < limit else v for v in jitter]
-        ax.cla()
-        draw_beam(ax, (0, 0))
-        draw_jitter(ax, view, jitter, dist=dist, size=size, draw_shape=draw_shape)
-        #draw_jitter(ax, view, (0,0,0))
-        draw_mesh(ax, view, jitter, dist=dist, n=mesh, projection=projection)
-        draw_scattering(calculator, ax, view, jitter, dist=dist)
+        axes.cla()
+        draw_beam(axes, (0, 0))
+        draw_jitter(axes, view, jitter, dist=dist, size=size, draw_shape=draw_shape)
+        #draw_jitter(axes, view, (0,0,0))
+        draw_mesh(axes, view, jitter, dist=dist, n=mesh, projection=projection)
+        draw_scattering(calculator, axes, view, jitter, dist=dist)
         plt.gcf().canvas.draw()
 
