Index: doc/developer/overview.rst
===================================================================
--- doc/developer/overview.rst	(revision 2a7e20e9c4ba72da5c6088e743a64f3345ee554d)
+++ doc/developer/overview.rst	(revision 0d5a655d6ddebe6096fad050347df4796b956767)
@@ -190,5 +190,5 @@
 *sasmodels.generate.PROJECTION*, with the default *PROJECTION=1* for
 equirectangular and *PROJECTION=2* for sinusoidal.  The more complicated
-Guyou and Postel projections are not implemented. See explore.jitter.draw_mesh
+Guyou and Postel projections are not implemented. See jitter.draw_mesh
 for details.
 
@@ -274,6 +274,6 @@
 code.
 
-*explore/jitter.py* is for exploring different options for handling
-orientation and orientation dispersity.  It uses *explore/guyou.py* to
+*sasmodels/jitter.py* is for exploring different options for handling
+orientation and orientation dispersity.  It uses *sasmodels/guyou.py* to
 generate the Guyou projection.
 
Index: doc/genapi.py
===================================================================
--- doc/genapi.py	(revision 706f46620c115bdbc319d902ba7678c9d247d6d0)
+++ doc/genapi.py	(revision 0d5a655d6ddebe6096fad050347df4796b956767)
@@ -69,4 +69,6 @@
     ('exception', 'Annotate exceptions'),
     ('generate', 'Model parser'),
+    ('jitter', 'Orientation explorer'),
+    ('guyou', 'Guyou map projection'),
     ('kernel', 'Evaluator type definitions'),
     ('kernelcl', 'OpenCL model evaluator'),
Index: doc/guide/orientation/orientation.rst
===================================================================
--- doc/guide/orientation/orientation.rst	(revision 5fb0634b45bba451ca824eecf40346b1ed5305a7)
+++ doc/guide/orientation/orientation.rst	(revision 0d5a655d6ddebe6096fad050347df4796b956767)
@@ -52,4 +52,11 @@
 yaw angle $\theta$ about the $b$-axis and pitch angle $\phi$ about the
 $a$-axis.
+
+You can explore the view and jitter angles interactively using
+:func:`sasmodels.jitter.run`.  Enter the following into the python
+interpreter::
+
+    from sasmodels import jitter
+    jitter.run()
 
 More formally, starting with axes $a$-$b$-$c$ of the particle aligned
Index: plore/guyou.py
===================================================================
--- explore/guyou.py	(revision 9d9fcbd8b087f9dd571a6eefd725819968dfea54)
+++ 	(revision )
@@ -1,375 +1,0 @@
-# Copyright 2013-2016 Mike Bostock
-# All rights reserved.
-#
-# Redistribution and use in source and binary forms, with or without modification,
-# are permitted provided that the following conditions are met:
-#
-# * Redistributions of source code must retain the above copyright notice, this
-#   list of conditions and the following disclaimer.
-#
-# * Redistributions in binary form must reproduce the above copyright notice,
-#   this list of conditions and the following disclaimer in the documentation
-#   and/or other materials provided with the distribution.
-#
-# * Neither the name of the author nor the names of contributors may be used to
-#   endorse or promote products derived from this software without specific prior
-#   written permission.
-#
-# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
-# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
-# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
-# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
-# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
-# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
-# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
-# ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
-# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-#
-# https://github.com/d3/d3-geo-projection
-# commit fd2886555e46b35163c7b898d43c7d1bcebbba7c 2016-07-02
-#
-# 2017-11-01 Paul Kienzle
-# * converted to python, using degrees rather than radians
-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
-
-_ = """
-# mpmath version of special functions
-import mpmath as mp
-sn, cn, dn = [mp.ellipfun(v) for v in 'sn', 'cn', 'dn']
-
-def ellipticJi(u, v, m):
-    z = u+v*1j
-    return sn(z, m), cn(z, m), dn(z, m)
-
-def ellipticJ(u, m):
-    s, c, d = sn(u, m), cn(u, m), dn(u, m)
-    phi = mp.asin(s)
-    return s, c, d, phi
-
-def ellipticFi(phi, psi, m):
-    z = phi + psi*1j
-    return mp.ellipf(z, m)
-
-def ellipticF(phi, m):
-    return mp.ellipf(phi, m)
-"""
-
-# 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
-    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
-
-def _ellipticJi_real(u, m):
-    sn, cn, dn, phi = ellipticJ(u, m)
-    return sn, cn, dn
-
-def _ellipticJi_imag(v, m):
-    sn, cn, dn, phi = ellipticJ(v, 1-m)
-    return 1j*sn/cn, 1/cn, dn/cn
-
-def _ellipticJi(u, v, m):
-    # Ignoring special cases for now.
-    #   u=0: (1j*b[0]/b[1], 1/b[1], b[2]/b[1])
-    #   v=0: (a[0], a[1], a[2])
-    a = ellipticJ(u, m)
-    b = ellipticJ(v, 1 - m)
-    c = b[1]**2 + m * (a[0] * b[0])**2
-    return [
-        (a[0] * b[2] / c) + 1j*(a[1] * a[2] * b[0] * b[1] / c),
-        (a[1] * b[1] / c) + 1j*(-a[0] * a[2] * b[0] * b[2] / c),
-        (a[2] * b[1] * b[2] / c) + 1j*(-m * a[0] * a[1] * b[0] / c),
-    ]
-
-# calculate F(phi+ipsi|m).
-# see Abramowitz and Stegun, 17.4.11.
-def ellipticFi(phi, psi, m):
-    if np.any(phi == 0):
-        scalar = np.isscalar(phi) and np.isscalar(psi) and np.isscalar(m)
-        phi, psi, m = np.broadcast_arrays(phi, psi, m)
-        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] = ellipticFi(phi[~index], psi[~index], m[~index])
-        return result.reshape(1)[0] if scalar else result
-
-    r = abs(phi)
-    i = abs(psi)
-    sinhpsi2 = sinh(i)**2
-    cscphi2 = 1 / sin(r)**2
-    cotphi2 = 1 / tan(r)**2
-    b = -(cotphi2 + m * (sinhpsi2 * cscphi2) - 1 + m)
-    c = (m - 1) * cotphi2
-    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)
-    return re + 1j*im
-
-sqrt2 = sqrt(2)
-
-# [PAK] renamed k_ => cos_u, k => sin_u, k*k => sinsq_u to avoid k,K confusion
-# cos_u = 0.171572875253809902396622551580603842860656249246103853646...
-# sinsq_u = 0.970562748477140585620264690516376942836062504523376878120...
-# K = 3.165103454447431823666270142140819753058976299237578486994...
-def guyou(lam, phi):
-    # [PAK] wrap into [-pi/2, pi/2] radians
-    x, y = np.asarray(lam), np.asarray(phi)
-    xn, x = divmod(x+90, 180)
-    yn, y = divmod(y+90, 180)
-    xn, lam = xn*180, radians(x-90)
-    yn, phi = yn*180, radians(y-90)
-
-    # Compute constant K
-    cos_u = (sqrt2 - 1) / (sqrt2 + 1)
-    sinsq_u = 1 - cos_u**2
-    K = ellipticF(pi/2, sinsq_u)
-
-    # [PAK] simplify expressions, using the fact that f = -1
-    # Note: exp(f log(x)) => 1/x,  cos(f x) => cos(x), sin(f x) => -sin(x)
-    r = 1/(tan(pi/4 + abs(phi)/2) * sqrt(cos_u))
-    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))
-
-    # [PAK] convert to degrees, and return to original tile
-    return degrees(x)+xn, degrees(y)+yn
-
-def guyou_invert(x, y):
-    # [PAK] wrap into [-pi/2, pi/2] radians
-    x, y = np.asarray(x), np.asarray(y)
-    xn, x = divmod(x+90, 180)
-    yn, y = divmod(y+90, 180)
-    xn, x = xn*180, radians(x-90)
-    yn, y = yn*180, radians(y-90)
-
-    # compute constant K
-    cos_u = (sqrt2 - 1) / (sqrt2 + 1)
-    sinsq_u = 1 - cos_u**2
-    K = ellipticF(pi/2, sinsq_u)
-
-    # [PAK] simplify expressions, using the fact that f = -1
-    j = ellipticJi(K/2 - y, -x, sinsq_u)
-    tn = j[0]/j[1]  # j[0], j[1] are complex
-    # Note: -atan2(im(x)/re(x)) => angle(x)
-    lam = -np.angle(tn)
-    # Note: exp(0.5/f log(a re(x)^2 + a im(x)^2)) => 1/(sqrt(a) |x|)
-    phi = 2*atan(1/sqrt(cos_u)/abs(tn)) - pi/2
-
-    # [PAK] convert to degrees, and return to original tile
-    return degrees(lam)+xn, degrees(phi)+yn
-
-def plot_grid():
-    import matplotlib.pyplot as plt
-    from numpy import linspace
-    lat_line = linspace(-90, 90, 400)
-    long_line = linspace(-90, 90, 400)
-
-    #scale = 1
-    limit, step, scale = 90, 10, 2
-    plt.subplot(211)
-    for lat in range(-limit, limit+1, step):
-        x, y = guyou(scale*lat_line, scale*lat)
-        plt.plot(x, y, 'g')
-
-    for long in range(-limit, limit+1, step):
-        x, y = guyou(scale*long, scale*long_line)
-        plt.plot(x, y, 'b')
-    #plt.xlabel('longitude')
-    plt.ylabel('latitude')
-
-    plt.subplot(212)
-    for lat in range(-limit, limit+1, step):
-        x, y = guyou_invert(scale*lat_line, scale*lat)
-        plt.plot(x, y, 'g')
-
-    for long in range(-limit, limit+1, step):
-        x, y = guyou_invert(scale*long, scale*long_line)
-        plt.plot(x, y, 'b')
-    plt.xlabel('longitude')
-    plt.ylabel('latitude')
-
-if __name__ == "__main__":
-    plot_grid()
-    import matplotlib.pyplot as plt; plt.show()
-
-
-_ = """
-// Javascript source for elliptic functions
-//
-// Returns [sn, cn, dn](u + iv|m).
-export function ellipticJi(u, v, m) {
-  var a, b, c;
-  if (!u) {
-    b = ellipticJ(v, 1 - m);
-    return [
-      [0, b[0] / b[1]],
-      [1 / b[1], 0],
-      [b[2] / b[1], 0]
-    ];
-  }
-  a = ellipticJ(u, m);
-  if (!v) return [[a[0], 0], [a[1], 0], [a[2], 0]];
-  b = ellipticJ(v, 1 - m);
-  c = b[1] * b[1] + m * a[0] * a[0] * b[0] * b[0];
-  return [
-    [a[0] * b[2] / c, a[1] * a[2] * b[0] * b[1] / c],
-    [a[1] * b[1] / c, -a[0] * a[2] * b[0] * b[2] / c],
-    [a[2] * b[1] * b[2] / c, -m * a[0] * a[1] * b[0] / c]
-  ];
-}
-
-// Returns [sn, cn, dn, ph](u|m).
-export function ellipticJ(u, m) {
-  var ai, b, phi, t, twon;
-  if (m < epsilon) {
-    t = sin(u);
-    b = cos(u);
-    ai = m * (u - t * b) / 4;
-    return [
-      t - ai * b,
-      b + ai * t,
-      1 - m * t * t / 2,
-      u - ai
-    ];
-  }
-  if (m >= 1 - epsilon) {
-    ai = (1 - m) / 4;
-    b = cosh(u);
-    t = tanh(u);
-    phi = 1 / b;
-    twon = b * sinh(u);
-    return [
-      t + ai * (twon - u) / (b * b),
-      phi - ai * t * phi * (twon - u),
-      phi + ai * t * phi * (twon + u),
-      2 * atan(exp(u)) - halfPi + ai * (twon - u) / b
-    ];
-  }
-
-  var a = [1, 0, 0, 0, 0, 0, 0, 0, 0],
-      c = [sqrt(m), 0, 0, 0, 0, 0, 0, 0, 0],
-      i = 0;
-  b = sqrt(1 - m);
-  twon = 1;
-
-  while (abs(c[i] / a[i]) > epsilon && i < 8) {
-    ai = a[i++];
-    c[i] = (ai - b) / 2;
-    a[i] = (ai + b) / 2;
-    b = sqrt(ai * b);
-    twon *= 2;
-  }
-
-  phi = twon * a[i] * u;
-  do {
-    t = c[i] * sin(b = phi) / a[i];
-    phi = (asin(t) + phi) / 2;
-  } while (--i);
-
-  // [PAK] Cephes uses dn = sqrt(1 - m*sin^2 phi) rather than cos(phi)/cos(phi-b)
-  // DLMF says the second version is unstable near x = K.
-  return [sin(phi), t = cos(phi), t / cos(phi - b), phi];
-}
-
-// calculate F(phi+ipsi|m).
-// see Abramowitz and Stegun, 17.4.11.
-export function ellipticFi(phi, psi, m) {
-  var r = abs(phi),
-      i = abs(psi),
-      sinhpsi = sinh(i);
-  if (r) {
-    var cscphi = 1 / sin(r),
-        cotphi2 = 1 / (tan(r) * tan(r)),
-        b = -(cotphi2 + m * (sinhpsi * sinhpsi * cscphi * cscphi) - 1 + m),
-        c = (m - 1) * cotphi2,
-        cotlambda2 = (-b + sqrt(b * b - 4 * c)) / 2;
-    return [
-      ellipticF(atan(1 / sqrt(cotlambda2)), m) * sign(phi),
-      ellipticF(atan(sqrt((cotlambda2 / cotphi2 - 1) / m)), 1 - m) * sign(psi)
-    ];
-  }
-  return [
-    0,
-    ellipticF(atan(sinhpsi), 1 - m) * sign(psi)
-  ];
-}
-
-// Calculate F(phi|m) where m = k^2 = sin(alpha)^2.
-// See Abramowitz and Stegun, 17.6.7.
-export function ellipticF(phi, m) {
-  if (!m) return phi;
-  if (m === 1) return log(tan(phi / 2 + quarterPi));
-  var a = 1,
-      b = sqrt(1 - m),
-      c = sqrt(m);
-  for (var i = 0; abs(c) > epsilon; i++) {
-    if (phi % pi) {
-      var dPhi = atan(b * tan(phi) / a);
-      if (dPhi < 0) dPhi += pi;
-      phi += dPhi + ~~(phi / pi) * pi;
-    } else phi += phi;
-    c = (a + b) / 2;
-    b = sqrt(a * b);
-    c = ((a = c) - b) / 2;
-  }
-  return phi / (pow(2, i) * a);
-
-
-export function guyouRaw(lambda, phi) {
-  var k_ = (sqrt2 - 1) / (sqrt2 + 1),
-      k = sqrt(1 - k_ * k_),
-      K = ellipticF(halfPi, k * k),
-      f = -1,
-      psi = log(tan(pi / 4 + abs(phi) / 2)),
-      r = exp(f * psi) / sqrt(k_),
-      at = guyouComplexAtan(r * cos(f * lambda), r * sin(f * lambda)),
-      t = ellipticFi(at[0], at[1], k * k);
-  return [-t[1], (phi >= 0 ? 1 : -1) * (0.5 * K - t[0])];
-}
-
-function guyouComplexAtan(x, y) {
-  var x2 = x * x,
-      y_1 = y + 1,
-      t = 1 - x2 - y * y;
-  return [
-   0.5 * ((x >= 0 ? halfPi : -halfPi) - atan2(t, 2 * x)),
-    -0.25 * log(t * t + 4 * x2) +0.5 * log(y_1 * y_1 + x2)
-  ];
-}
-
-function guyouComplexDivide(a, b) {
-  var denominator = b[0] * b[0] + b[1] * b[1];
-  return [
-    (a[0] * b[0] + a[1] * b[1]) / denominator,
-    (a[1] * b[0] - a[0] * b[1]) / denominator
-  ];
-}
-
-guyouRaw.invert = function(x, y) {
-  var k_ = (sqrt2 - 1) / (sqrt2 + 1),
-      k = sqrt(1 - k_ * k_),
-      K = ellipticF(halfPi, k * k),
-      f = -1,
-      j = ellipticJi(0.5 * K - y, -x, k * k),
-      tn = guyouComplexDivide(j[0], j[1]),
-      lambda = atan2(tn[1], tn[0]) / f;
-  return [
-    lambda,
-    2 * atan(exp(0.5 / f * log(k_ * tn[0] * tn[0] + k_ * tn[1] * tn[1]))) - halfPi
-  ];
-};
-"""
Index: plore/jitter.py
===================================================================
--- explore/jitter.py	(revision 31a02f1520275e9438e6d60aa7e24838f67d3746)
+++ 	(revision )
@@ -1,726 +1,0 @@
-#!/usr/bin/env python
-"""
-Application to explore the difference between sasview 3.x orientation
-dispersity and possible replacement algorithms.
-"""
-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
-
-import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot
-import matplotlib.pyplot as plt
-from matplotlib.widgets import Slider, CheckButtons
-from matplotlib import cm
-import numpy as np
-from numpy import pi, cos, sin, sqrt, exp, degrees, radians
-
-def draw_beam(ax, 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)
-
-    steps = 25
-    u = np.linspace(0, 2 * np.pi, steps)
-    v = np.linspace(-1, 1, steps)
-
-    r = 0.02
-    x = r*np.outer(np.cos(u), np.ones_like(v))
-    y = r*np.outer(np.sin(u), np.ones_like(v))
-    z = 1.3*np.outer(np.ones_like(u), v)
-
-    theta, phi = view
-    shape = x.shape
-    points = np.matrix([x.flatten(), y.flatten(), z.flatten()])
-    points = Rz(phi)*Ry(theta)*points
-    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_jitter(ax, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0)):
-    """
-    Represent jitter as a set of shapes at different orientations.
-    """
-    # set max diagonal to 0.95
-    scale = 0.95/sqrt(sum(v**2 for v in size))
-    size = tuple(scale*v for v in size)
-    draw_shape = draw_parallelepiped
-    #draw_shape = draw_ellipsoid
-
-    #np.random.seed(10)
-    #cloud = np.random.randn(10,3)
-    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],
-    ]
-    dtheta, dphi, dpsi = jitter
-    if dtheta == 0:
-        cloud = [v for v in cloud if v[0] == 0]
-    if dphi == 0:
-        cloud = [v for v in cloud if v[1] == 0]
-    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)
-    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)
-    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)
-
-def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1):
-    """Draw an ellipsoid."""
-    a,b,c = size
-    u = np.linspace(0, 2 * np.pi, steps)
-    v = np.linspace(0, np.pi, steps)
-    x = a*np.outer(np.cos(u), np.sin(v))
-    y = b*np.outer(np.sin(u), np.sin(v))
-    z = c*np.outer(np.ones_like(u), np.cos(v))
-    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, [
-         ('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_parallelepiped(ax, size, view, jitter, steps=None, alpha=1):
-    """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])
-    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
-    ])
-
-    x, y, z = transform_xyz(view, jitter, x, y, z)
-    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha)
-
-    # Draw pink face on box.
-    # Since I can't control face color, instead draw a thin box situated just
-    # in front of the "a+" 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, y, z = transform_xyz(view, jitter, abs(x)*1.05, y, z)
-        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]),
-    ])
-
-def draw_sphere(ax, radius=10., steps=100):
-    """Draw a sphere"""
-    u = np.linspace(0, 2 * np.pi, steps)
-    v = np.linspace(0, np.pi, steps)
-
-    x = radius * np.outer(np.cos(u), np.sin(v))
-    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')
-
-PROJECTIONS = [
-    # in order of PROJECTION number; do not change without updating the
-    # constants in kernel_iq.c
-    'equirectangular', 'sinusoidal', 'guyou', 'azimuthal_equidistance',
-    'azimuthal_equal_area',
-]
-def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian',
-              projection='equirectangular'):
-    """
-    Draw the dispersion mesh showing the theta-phi orientations at which
-    the model will be evaluated.
-
-    jitter projections
-    <https://en.wikipedia.org/wiki/List_of_map_projections>
-
-    equirectangular (standard latitude-longitude mesh)
-        <https://en.wikipedia.org/wiki/Equirectangular_projection>
-        Allows free movement in phi (around the equator), but theta is
-        limited to +/- 90, and points are cos-weighted. Jitter in phi is
-        uniform in weight along a line of latitude.  With small theta and
-        phi ranging over +/- 180 this forms a wobbling disk.  With small
-        phi and theta ranging over +/- 90 this forms a wedge like a slice
-        of an orange.
-    azimuthal_equidistance (Postel)
-        <https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection>
-        Preserves distance from center, and so is an excellent map for
-        representing a bivariate gaussian on the surface.  Theta and phi
-        operate identically, cutting wegdes from the antipode of the viewing
-        angle.  This unfortunately does not allow free movement in either
-        theta or phi since the orthogonal wobble decreases to 0 as the body
-        rotates through 180 degrees.
-    sinusoidal (Sanson-Flamsteed, Mercator equal-area)
-        <https://en.wikipedia.org/wiki/Sinusoidal_projection>
-        Preserves arc length with latitude, giving bad behaviour at
-        theta near +/- 90.  Theta and phi operate somewhat differently,
-        so a system with a-b-c dtheta-dphi-dpsi will not give the same
-        value as one with b-a-c dphi-dtheta-dpsi, as would be the case
-        for azimuthal equidistance.  Free movement using theta or phi
-        uniform over +/- 180 will work, but not as well as equirectangular
-        phi, with theta being slightly worse.  Computationally it is much
-        cheaper for wide theta-phi meshes since it excludes points which
-        lie outside the sinusoid near theta +/- 90 rather than packing
-        them close together as in equirectangle.  Note that the poles
-        will be slightly overweighted for theta > 90 with the circle
-        from theta at 90+dt winding backwards around the pole, overlapping
-        the circle from theta at 90-dt.
-    Guyou (hemisphere-in-a-square) **not weighted**
-        <https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection>
-        With tiling, allows rotation in phi or theta through +/- 180, with
-        uniform spacing.  Both theta and phi allow free rotation, with wobble
-        in the orthogonal direction reasonably well behaved (though not as
-        good as equirectangular phi). The forward/reverse transformations
-        relies on elliptic integrals that are somewhat expensive, so the
-        behaviour has to be very good to justify the cost and complexity.
-        The weighting function for each point has not yet been computed.
-        Note: run the module *guyou.py* directly and it will show the forward
-        and reverse mappings.
-    azimuthal_equal_area  **incomplete**
-        <https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>
-        Preserves the relative density of the surface patches.  Not that
-        useful and not completely implemented
-    Gauss-Kreuger **not implemented**
-        <https://en.wikipedia.org/wiki/Transverse_Mercator_projection#Ellipsoidal_transverse_Mercator>
-        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)
-    if dist == 'gaussian':
-        t *= 3
-        weights = exp(-0.5*t**2)
-    elif dist == 'rectangle':
-        # Note: uses sasmodels ridiculous definition of rectangle width
-        t *= sqrt(3)
-    elif dist == 'uniform':
-        pass
-    else:
-        raise ValueError("expected dist to be gaussian, rectangle or uniform")
-
-    if projection == 'equirectangular':  #define PROJECTION 1
-        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)))
-    elif projection == 'sinusoidal':  #define PROJECTION 2
-        def rotate(theta_i, phi_j):
-            latitude = theta_i
-            scale = cos(radians(latitude))
-            longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0
-            #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):
-            latitude = theta_i
-            scale = cos(radians(latitude))
-            w = 1 if abs(phi_j) < abs(scale)*180 else 0
-            return w*wi*wj
-    elif projection == 'guyou':  #define PROJECTION 3  (eventually?)
-        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):
-            return wi*wj
-    elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry
-        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):
-            # Weighting for each point comes from the integral:
-            #     \int\int I(q, lat, log) sin(lat) dlat dlog
-            # We are doing a conformal mapping from disk to sphere, so we need
-            # a change of variables g(theta, phi) -> (lat, long):
-            #     lat, long = sqrt(theta^2 + phi^2), arctan(phi/theta)
-            # giving:
-            #     dtheta dphi = det(J) dlat dlong
-            # where J is the jacobian from the partials of g. Using
-            #     R = sqrt(theta^2 + phi^2),
-            # then
-            #     J = [[x/R, Y/R], -y/R^2, x/R^2]]
-            # and
-            #     det(J) = 1/R
-            # with the final integral being:
-            #    \int\int I(q, theta, phi) sin(R)/R dtheta dphi
-            #
-            # This does approximately the right thing, decreasing the weight
-            # of each point as you go farther out on the disk, but it hasn't
-            # yet been checked against the 1D integral results. Prior
-            # to declaring this "good enough" and checking that integrals
-            # work in practice, we will examine alternative mappings.
-            #
-            # The issue is that the mapping does not support the case of free
-            # rotation about a single axis correctly, with a small deviation
-            # in the orthogonal axis independent of the first axis.  Like the
-            # usual polar coordiates integration, the integrated sections
-            # form wedges, though at least in this case the wedge cuts through
-            # 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
-    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))
-            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):
-            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
-    else:
-        raise ValueError("unknown projection %r"%projection)
-
-    # mesh in theta, phi formed by rotating z
-    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])
-    # 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)])
-    #print(max(w), min(w), min(w[w>0]))
-    points = points[:, w>0]
-    w = w[w>0]
-    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)
-
-    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):
-    """
-    Draw text at a particular location.
-    """
-    labels, locations, orientations = zip(*text)
-    px, py, pz = zip(*locations)
-    dx, dy, dz = zip(*orientations)
-
-    px, py, pz = transform_xyz(view, jitter, px, py, pz)
-    dx, dy, dz = transform_xyz(view, jitter, dx, dy, dz)
-
-    # TODO: zdir for labels is broken, and labels aren't appearing.
-    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)
-
-# Definition of rotation matrices comes from wikipedia:
-#    https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations
-def Rx(angle):
-    """Construct a matrix to rotate points about *x* by *angle* degrees."""
-    a = radians(angle)
-    R = [[1, 0, 0],
-         [0, +cos(a), -sin(a)],
-         [0, +sin(a), +cos(a)]]
-    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)],
-         [0, 1, 0],
-         [-sin(a), 0, +cos(a)]]
-    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],
-         [0, 0, 1]]
-    return np.matrix(R)
-
-def transform_xyz(view, jitter, x, y, z):
-    """
-    Send a set of (x,y,z) points through the jitter and view transforms.
-    """
-    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 = apply_jitter(jitter, points)
-    points = orient_relative_to_beam(view, points)
-    x, y, z = [np.array(v).reshape(shape) for v in points]
-    return x, y, z
-
-def apply_jitter(jitter, points):
-    """
-    Apply the jitter transform to a set of points.
-
-    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple.
-    """
-    dtheta, dphi, dpsi = jitter
-    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points
-    return points
-
-def orient_relative_to_beam(view, points):
-    """
-    Apply the view transform to a set of points.
-
-    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple.
-    """
-    theta, phi, psi = view
-    points = Rz(phi)*Ry(theta)*Rz(psi)*points
-    return points
-
-# translate between number of dimension of dispersity and the number of
-# points along each dimension.
-PD_N_TABLE = {
-    (0, 0, 0): (0, 0, 0),     # 0
-    (1, 0, 0): (100, 0, 0),   # 100
-    (0, 1, 0): (0, 100, 0),
-    (0, 0, 1): (0, 0, 100),
-    (1, 1, 0): (30, 30, 0),   # 900
-    (1, 0, 1): (30, 0, 30),
-    (0, 1, 1): (0, 30, 30),
-    (1, 1, 1): (15, 15, 15),  # 3375
-}
-
-def clipped_range(data, portion=1.0, mode='central'):
-    """
-    Determine range from data.
-
-    If *portion* is 1, use full range, otherwise use the center of the range
-    or the top of the range, depending on whether *mode* is 'central' or 'top'.
-    """
-    if portion == 1.0:
-        return data.min(), data.max()
-    elif mode == 'central':
-        data = np.sort(data.flatten())
-        offset = int(portion*len(data)/2 + 0.5)
-        return data[offset], data[-offset]
-    elif mode == 'top':
-        data = np.sort(data.flatten())
-        offset = int(portion*len(data) + 0.5)
-        return data[offset], data[-1]
-
-def draw_scattering(calculator, ax, view, jitter, dist='gaussian'):
-    """
-    Plot the scattering for the particular view.
-
-    *calculator* is returned from :func:`build_model`.  *ax* 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
-    weight distributions.
-    """
-    if dist == 'uniform':  # uniform is not yet in this branch
-        dist, scale = 'rectangle', 1/sqrt(3)
-    else:
-        scale = 1
-
-    # add the orientation parameters to the model parameters
-    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)]
-    ## 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)]
-
-    pars = dict(
-        theta=theta, theta_pd=theta_pd, theta_pd_type=dist, theta_pd_n=theta_pd_n,
-        phi=phi, phi_pd=phi_pd, phi_pd_type=dist, phi_pd_n=phi_pd_n,
-        psi=psi, psi_pd=psi_pd, psi_pd_type=dist, psi_pd_n=psi_pd_n,
-    )
-    pars.update(calculator.pars)
-
-    # compute the pattern
-    qx, qy = calculator._data.x_bins, calculator._data.y_bins
-    Iqxy = calculator(**pars).reshape(len(qx), len(qy))
-
-    # scale it and draw it
-    Iqxy = np.log(Iqxy)
-    if calculator.limits:
-        # use limits from orientation (0,0,0)
-        vmin, vmax = calculator.limits
-    else:
-        vmax = Iqxy.max()
-        vmin = vmax*10**-7
-        #vmin, vmax = clipped_range(Iqxy, portion=portion, mode='top')
-    #print("range",(vmin,vmax))
-    #qx, qy = np.meshgrid(qx, qy)
-    if 0:
-        level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i')
-        level[level<0] = 0
-        colors = plt.get_cmap()(level)
-        ax.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,
-                    levels=np.linspace(vmin, vmax, 24))
-    else:
-        ax.pcolormesh(qx, qy, Iqxy)
-
-def build_model(model_name, n=150, qmax=0.5, **pars):
-    """
-    Build a calculator for the given shape.
-
-    *model_name* is any sasmodels model.  *n* and *qmax* define an n x n mesh
-    on which to evaluate the model.  The remaining parameters are stored in
-    the returned calculator as *calculator.pars*.  They are used by
-    :func:`draw_scattering` to set the non-orientation parameters in the
-    calculation.
-
-    Returns a *calculator* function which takes a dictionary or parameters and
-    produces Iqxy.  The Iqxy value needs to be reshaped to an n x n matrix
-    for plotting.  See the :class:`sasmodels.direct_model.DirectModel` class
-    for details.
-    """
-    from sasmodels.core import load_model_info, build_model
-    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!')
-    q = np.linspace(-qmax, qmax, n)
-    data = empty_data2D(q, q)
-    calculator = DirectModel(data, model)
-
-    # stuff the values for non-orientation parameters into the calculator
-    calculator.pars = pars.copy()
-    calculator.pars.setdefault('backgound', 1e-3)
-
-    # fix the data limits so that we can see if the pattern fades
-    # under rotation or angular dispersion
-    Iqxy = calculator(theta=0, phi=0, psi=0, **calculator.pars)
-    Iqxy = np.log(Iqxy)
-    vmin, vmax = clipped_range(Iqxy, 0.95, mode='top')
-    calculator.limits = vmin, vmax+1
-
-    return calculator
-
-def select_calculator(model_name, n=150, size=(10,40,100)):
-    """
-    Create a model calculator for the given shape.
-
-    *model_name* is one of sphere, cylinder, ellipsoid, triaxial_ellipsoid,
-    parallelepiped or bcc_paracrystal. *n* is the number of points to use
-    in the q range.  *qmax* is chosen based on model parameters for the
-    given model to show something intersting.
-
-    Returns *calculator* and tuple *size* (a,b,c) giving minor and major
-    equitorial axes and polar axis respectively.  See :func:`build_model`
-    for details on the returned calculator.
-    """
-    a, b, c = size
-    if model_name == 'sphere':
-        calculator = build_model('sphere', n=n, radius=c)
-        a = b = c
-    elif model_name == 'bcc_paracrystal':
-        calculator = build_model('bcc_paracrystal', n=n, dnn=c,
-                                  d_factor=0.06, radius=40)
-        a = b = c
-    elif model_name == 'cylinder':
-        calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c)
-        a = b
-    elif model_name == 'ellipsoid':
-        calculator = build_model('ellipsoid', n=n, qmax=1.0,
-                                 radius_polar=c, radius_equatorial=b)
-        a = b
-    elif model_name == 'triaxial_ellipsoid':
-        calculator = build_model('triaxial_ellipsoid', n=n, qmax=0.5,
-                                 radius_equat_minor=a,
-                                 radius_equat_major=b,
-                                 radius_polar=c)
-    elif model_name == 'parallelepiped':
-        calculator = build_model('parallelepiped', n=n, a=a, b=b, c=c)
-    else:
-        raise ValueError("unknown model %s"%model_name)
-
-    return calculator, (a, b, c)
-
-SHAPES = [
-    'parallelepiped', 'triaxial_ellipsoid', 'bcc_paracrystal',
-    'cylinder', 'ellipsoid',
-    'sphere',
- ]
-
-DISTRIBUTIONS = [
-    'gaussian', 'rectangle', 'uniform',
-]
-DIST_LIMITS = {
-    'gaussian': 30,
-    'rectangle': 90/sqrt(3),
-    'uniform': 90,
-}
-
-def run(model_name='parallelepiped', size=(10, 40, 100),
-        dist='gaussian', mesh=30,
-        projection='equirectangular'):
-    """
-    Show an interactive orientation and jitter demo.
-
-    *model_name* is one of the models available in :func:`select_model`.
-    """
-    # projection number according to 1-order position in list, but
-    # only 1 and 2 are implemented so far.
-    from sasmodels import generate
-    generate.PROJECTION = PROJECTIONS.index(projection) + 1
-    if generate.PROJECTION > 2:
-        print("*** PROJECTION %s not implemented in scattering function ***"%projection)
-        generate.PROJECTION = 2
-
-    # set up calculator
-    calculator, size = select_calculator(model_name, n=150, size=size)
-
-    ## uncomment to set an independent the colour range for every view
-    ## If left commented, the colour range is fixed for all views
-    calculator.limits = None
-
-    ## initial view
-    #theta, dtheta = 70., 10.
-    #phi, dphi = -45., 3.
-    #psi, dpsi = -45., 3.
-    theta, phi, psi = 0, 0, 0
-    dtheta, dphi, dpsi = 0, 0, 0
-
-    ## create the plot window
-    #plt.hold(True)
-    plt.subplots(num=None, figsize=(5.5, 5.5))
-    plt.set_cmap('gist_earth')
-    plt.clf()
-    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')
-    try:  # CRUFT: not all versions of matplotlib accept 'square' 3d projection
-        ax.axis('square')
-    except Exception:
-        pass
-
-    axcolor = 'lightgoldenrodyellow'
-
-    ## 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)
-    # 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)
-
-    ## callback to draw the new view
-    def update(val, axis=None):
-        view = stheta.val, sphi.val, spsi.val
-        jitter = sdtheta.val, sdphi.val, sdpsi.val
-        # set small jitter as 0 if multiple pd dims
-        dims = sum(v > 0 for v in jitter)
-        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_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)
-        plt.gcf().canvas.draw()
-
-    ## bind control widgets to view updater
-    stheta.on_changed(lambda v: update(v,'theta'))
-    sphi.on_changed(lambda v: update(v, 'phi'))
-    spsi.on_changed(lambda v: update(v, 'psi'))
-    sdtheta.on_changed(lambda v: update(v, 'dtheta'))
-    sdphi.on_changed(lambda v: update(v, 'dphi'))
-    sdpsi.on_changed(lambda v: update(v, 'dpsi'))
-
-    ## initialize view
-    update(None, 'phi')
-
-    ## go interactive
-    plt.show()
-
-def main():
-    parser = argparse.ArgumentParser(
-        description="Display jitter",
-        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')
-    opts = parser.parse_args()
-    size = tuple(int(v) for v in opts.size.split(','))
-    run(opts.shape, size=size,
-        mesh=opts.mesh, dist=opts.distribution,
-        projection=opts.projection)
-
-if __name__ == "__main__":
-    main()
Index: sasmodels/guyou.py
===================================================================
--- sasmodels/guyou.py	(revision 0d5a655d6ddebe6096fad050347df4796b956767)
+++ sasmodels/guyou.py	(revision 0d5a655d6ddebe6096fad050347df4796b956767)
@@ -0,0 +1,375 @@
+# Copyright 2013-2016 Mike Bostock
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without modification,
+# are permitted provided that the following conditions are met:
+#
+# * Redistributions of source code must retain the above copyright notice, this
+#   list of conditions and the following disclaimer.
+#
+# * Redistributions in binary form must reproduce the above copyright notice,
+#   this list of conditions and the following disclaimer in the documentation
+#   and/or other materials provided with the distribution.
+#
+# * Neither the name of the author nor the names of contributors may be used to
+#   endorse or promote products derived from this software without specific prior
+#   written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+# ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+#
+# https://github.com/d3/d3-geo-projection
+# commit fd2886555e46b35163c7b898d43c7d1bcebbba7c 2016-07-02
+#
+# 2017-11-01 Paul Kienzle
+# * converted to python, using degrees rather than radians
+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
+
+_ = """
+# mpmath version of special functions
+import mpmath as mp
+sn, cn, dn = [mp.ellipfun(v) for v in 'sn', 'cn', 'dn']
+
+def ellipticJi(u, v, m):
+    z = u+v*1j
+    return sn(z, m), cn(z, m), dn(z, m)
+
+def ellipticJ(u, m):
+    s, c, d = sn(u, m), cn(u, m), dn(u, m)
+    phi = mp.asin(s)
+    return s, c, d, phi
+
+def ellipticFi(phi, psi, m):
+    z = phi + psi*1j
+    return mp.ellipf(z, m)
+
+def ellipticF(phi, m):
+    return mp.ellipf(phi, m)
+"""
+
+# 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
+    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
+
+def _ellipticJi_real(u, m):
+    sn, cn, dn, phi = ellipticJ(u, m)
+    return sn, cn, dn
+
+def _ellipticJi_imag(v, m):
+    sn, cn, dn, phi = ellipticJ(v, 1-m)
+    return 1j*sn/cn, 1/cn, dn/cn
+
+def _ellipticJi(u, v, m):
+    # Ignoring special cases for now.
+    #   u=0: (1j*b[0]/b[1], 1/b[1], b[2]/b[1])
+    #   v=0: (a[0], a[1], a[2])
+    a = ellipticJ(u, m)
+    b = ellipticJ(v, 1 - m)
+    c = b[1]**2 + m * (a[0] * b[0])**2
+    return [
+        (a[0] * b[2] / c) + 1j*(a[1] * a[2] * b[0] * b[1] / c),
+        (a[1] * b[1] / c) + 1j*(-a[0] * a[2] * b[0] * b[2] / c),
+        (a[2] * b[1] * b[2] / c) + 1j*(-m * a[0] * a[1] * b[0] / c),
+    ]
+
+# calculate F(phi+ipsi|m).
+# see Abramowitz and Stegun, 17.4.11.
+def ellipticFi(phi, psi, m):
+    if np.any(phi == 0):
+        scalar = np.isscalar(phi) and np.isscalar(psi) and np.isscalar(m)
+        phi, psi, m = np.broadcast_arrays(phi, psi, m)
+        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] = ellipticFi(phi[~index], psi[~index], m[~index])
+        return result.reshape(1)[0] if scalar else result
+
+    r = abs(phi)
+    i = abs(psi)
+    sinhpsi2 = sinh(i)**2
+    cscphi2 = 1 / sin(r)**2
+    cotphi2 = 1 / tan(r)**2
+    b = -(cotphi2 + m * (sinhpsi2 * cscphi2) - 1 + m)
+    c = (m - 1) * cotphi2
+    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)
+    return re + 1j*im
+
+sqrt2 = sqrt(2)
+
+# [PAK] renamed k_ => cos_u, k => sin_u, k*k => sinsq_u to avoid k,K confusion
+# cos_u = 0.171572875253809902396622551580603842860656249246103853646...
+# sinsq_u = 0.970562748477140585620264690516376942836062504523376878120...
+# K = 3.165103454447431823666270142140819753058976299237578486994...
+def guyou(lam, phi):
+    # [PAK] wrap into [-pi/2, pi/2] radians
+    x, y = np.asarray(lam), np.asarray(phi)
+    xn, x = divmod(x+90, 180)
+    yn, y = divmod(y+90, 180)
+    xn, lam = xn*180, radians(x-90)
+    yn, phi = yn*180, radians(y-90)
+
+    # Compute constant K
+    cos_u = (sqrt2 - 1) / (sqrt2 + 1)
+    sinsq_u = 1 - cos_u**2
+    K = ellipticF(pi/2, sinsq_u)
+
+    # [PAK] simplify expressions, using the fact that f = -1
+    # Note: exp(f log(x)) => 1/x,  cos(f x) => cos(x), sin(f x) => -sin(x)
+    r = 1/(tan(pi/4 + abs(phi)/2) * sqrt(cos_u))
+    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))
+
+    # [PAK] convert to degrees, and return to original tile
+    return degrees(x)+xn, degrees(y)+yn
+
+def guyou_invert(x, y):
+    # [PAK] wrap into [-pi/2, pi/2] radians
+    x, y = np.asarray(x), np.asarray(y)
+    xn, x = divmod(x+90, 180)
+    yn, y = divmod(y+90, 180)
+    xn, x = xn*180, radians(x-90)
+    yn, y = yn*180, radians(y-90)
+
+    # compute constant K
+    cos_u = (sqrt2 - 1) / (sqrt2 + 1)
+    sinsq_u = 1 - cos_u**2
+    K = ellipticF(pi/2, sinsq_u)
+
+    # [PAK] simplify expressions, using the fact that f = -1
+    j = ellipticJi(K/2 - y, -x, sinsq_u)
+    tn = j[0]/j[1]  # j[0], j[1] are complex
+    # Note: -atan2(im(x)/re(x)) => angle(x)
+    lam = -np.angle(tn)
+    # Note: exp(0.5/f log(a re(x)^2 + a im(x)^2)) => 1/(sqrt(a) |x|)
+    phi = 2*atan(1/sqrt(cos_u)/abs(tn)) - pi/2
+
+    # [PAK] convert to degrees, and return to original tile
+    return degrees(lam)+xn, degrees(phi)+yn
+
+def plot_grid():
+    import matplotlib.pyplot as plt
+    from numpy import linspace
+    lat_line = linspace(-90, 90, 400)
+    long_line = linspace(-90, 90, 400)
+
+    #scale = 1
+    limit, step, scale = 90, 10, 2
+    plt.subplot(211)
+    for lat in range(-limit, limit+1, step):
+        x, y = guyou(scale*lat_line, scale*lat)
+        plt.plot(x, y, 'g')
+
+    for long in range(-limit, limit+1, step):
+        x, y = guyou(scale*long, scale*long_line)
+        plt.plot(x, y, 'b')
+    #plt.xlabel('longitude')
+    plt.ylabel('latitude')
+
+    plt.subplot(212)
+    for lat in range(-limit, limit+1, step):
+        x, y = guyou_invert(scale*lat_line, scale*lat)
+        plt.plot(x, y, 'g')
+
+    for long in range(-limit, limit+1, step):
+        x, y = guyou_invert(scale*long, scale*long_line)
+        plt.plot(x, y, 'b')
+    plt.xlabel('longitude')
+    plt.ylabel('latitude')
+
+if __name__ == "__main__":
+    plot_grid()
+    import matplotlib.pyplot as plt; plt.show()
+
+
+_ = """
+// Javascript source for elliptic functions
+//
+// Returns [sn, cn, dn](u + iv|m).
+export function ellipticJi(u, v, m) {
+  var a, b, c;
+  if (!u) {
+    b = ellipticJ(v, 1 - m);
+    return [
+      [0, b[0] / b[1]],
+      [1 / b[1], 0],
+      [b[2] / b[1], 0]
+    ];
+  }
+  a = ellipticJ(u, m);
+  if (!v) return [[a[0], 0], [a[1], 0], [a[2], 0]];
+  b = ellipticJ(v, 1 - m);
+  c = b[1] * b[1] + m * a[0] * a[0] * b[0] * b[0];
+  return [
+    [a[0] * b[2] / c, a[1] * a[2] * b[0] * b[1] / c],
+    [a[1] * b[1] / c, -a[0] * a[2] * b[0] * b[2] / c],
+    [a[2] * b[1] * b[2] / c, -m * a[0] * a[1] * b[0] / c]
+  ];
+}
+
+// Returns [sn, cn, dn, ph](u|m).
+export function ellipticJ(u, m) {
+  var ai, b, phi, t, twon;
+  if (m < epsilon) {
+    t = sin(u);
+    b = cos(u);
+    ai = m * (u - t * b) / 4;
+    return [
+      t - ai * b,
+      b + ai * t,
+      1 - m * t * t / 2,
+      u - ai
+    ];
+  }
+  if (m >= 1 - epsilon) {
+    ai = (1 - m) / 4;
+    b = cosh(u);
+    t = tanh(u);
+    phi = 1 / b;
+    twon = b * sinh(u);
+    return [
+      t + ai * (twon - u) / (b * b),
+      phi - ai * t * phi * (twon - u),
+      phi + ai * t * phi * (twon + u),
+      2 * atan(exp(u)) - halfPi + ai * (twon - u) / b
+    ];
+  }
+
+  var a = [1, 0, 0, 0, 0, 0, 0, 0, 0],
+      c = [sqrt(m), 0, 0, 0, 0, 0, 0, 0, 0],
+      i = 0;
+  b = sqrt(1 - m);
+  twon = 1;
+
+  while (abs(c[i] / a[i]) > epsilon && i < 8) {
+    ai = a[i++];
+    c[i] = (ai - b) / 2;
+    a[i] = (ai + b) / 2;
+    b = sqrt(ai * b);
+    twon *= 2;
+  }
+
+  phi = twon * a[i] * u;
+  do {
+    t = c[i] * sin(b = phi) / a[i];
+    phi = (asin(t) + phi) / 2;
+  } while (--i);
+
+  // [PAK] Cephes uses dn = sqrt(1 - m*sin^2 phi) rather than cos(phi)/cos(phi-b)
+  // DLMF says the second version is unstable near x = K.
+  return [sin(phi), t = cos(phi), t / cos(phi - b), phi];
+}
+
+// calculate F(phi+ipsi|m).
+// see Abramowitz and Stegun, 17.4.11.
+export function ellipticFi(phi, psi, m) {
+  var r = abs(phi),
+      i = abs(psi),
+      sinhpsi = sinh(i);
+  if (r) {
+    var cscphi = 1 / sin(r),
+        cotphi2 = 1 / (tan(r) * tan(r)),
+        b = -(cotphi2 + m * (sinhpsi * sinhpsi * cscphi * cscphi) - 1 + m),
+        c = (m - 1) * cotphi2,
+        cotlambda2 = (-b + sqrt(b * b - 4 * c)) / 2;
+    return [
+      ellipticF(atan(1 / sqrt(cotlambda2)), m) * sign(phi),
+      ellipticF(atan(sqrt((cotlambda2 / cotphi2 - 1) / m)), 1 - m) * sign(psi)
+    ];
+  }
+  return [
+    0,
+    ellipticF(atan(sinhpsi), 1 - m) * sign(psi)
+  ];
+}
+
+// Calculate F(phi|m) where m = k^2 = sin(alpha)^2.
+// See Abramowitz and Stegun, 17.6.7.
+export function ellipticF(phi, m) {
+  if (!m) return phi;
+  if (m === 1) return log(tan(phi / 2 + quarterPi));
+  var a = 1,
+      b = sqrt(1 - m),
+      c = sqrt(m);
+  for (var i = 0; abs(c) > epsilon; i++) {
+    if (phi % pi) {
+      var dPhi = atan(b * tan(phi) / a);
+      if (dPhi < 0) dPhi += pi;
+      phi += dPhi + ~~(phi / pi) * pi;
+    } else phi += phi;
+    c = (a + b) / 2;
+    b = sqrt(a * b);
+    c = ((a = c) - b) / 2;
+  }
+  return phi / (pow(2, i) * a);
+
+
+export function guyouRaw(lambda, phi) {
+  var k_ = (sqrt2 - 1) / (sqrt2 + 1),
+      k = sqrt(1 - k_ * k_),
+      K = ellipticF(halfPi, k * k),
+      f = -1,
+      psi = log(tan(pi / 4 + abs(phi) / 2)),
+      r = exp(f * psi) / sqrt(k_),
+      at = guyouComplexAtan(r * cos(f * lambda), r * sin(f * lambda)),
+      t = ellipticFi(at[0], at[1], k * k);
+  return [-t[1], (phi >= 0 ? 1 : -1) * (0.5 * K - t[0])];
+}
+
+function guyouComplexAtan(x, y) {
+  var x2 = x * x,
+      y_1 = y + 1,
+      t = 1 - x2 - y * y;
+  return [
+   0.5 * ((x >= 0 ? halfPi : -halfPi) - atan2(t, 2 * x)),
+    -0.25 * log(t * t + 4 * x2) +0.5 * log(y_1 * y_1 + x2)
+  ];
+}
+
+function guyouComplexDivide(a, b) {
+  var denominator = b[0] * b[0] + b[1] * b[1];
+  return [
+    (a[0] * b[0] + a[1] * b[1]) / denominator,
+    (a[1] * b[0] - a[0] * b[1]) / denominator
+  ];
+}
+
+guyouRaw.invert = function(x, y) {
+  var k_ = (sqrt2 - 1) / (sqrt2 + 1),
+      k = sqrt(1 - k_ * k_),
+      K = ellipticF(halfPi, k * k),
+      f = -1,
+      j = ellipticJi(0.5 * K - y, -x, k * k),
+      tn = guyouComplexDivide(j[0], j[1]),
+      lambda = atan2(tn[1], tn[0]) / f;
+  return [
+    lambda,
+    2 * atan(exp(0.5 / f * log(k_ * tn[0] * tn[0] + k_ * tn[1] * tn[1]))) - halfPi
+  ];
+};
+"""
Index: sasmodels/jitter.py
===================================================================
--- sasmodels/jitter.py	(revision 0d5a655d6ddebe6096fad050347df4796b956767)
+++ sasmodels/jitter.py	(revision 0d5a655d6ddebe6096fad050347df4796b956767)
@@ -0,0 +1,738 @@
+#!/usr/bin/env python
+"""
+Jitter Explorer
+===============
+
+Application to explore orientation angle and angular dispersity.
+"""
+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
+
+import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider, CheckButtons
+from matplotlib import cm
+import numpy as np
+from numpy import pi, cos, sin, sqrt, exp, degrees, radians
+
+def draw_beam(ax, 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)
+
+    steps = 25
+    u = np.linspace(0, 2 * np.pi, steps)
+    v = np.linspace(-1, 1, steps)
+
+    r = 0.02
+    x = r*np.outer(np.cos(u), np.ones_like(v))
+    y = r*np.outer(np.sin(u), np.ones_like(v))
+    z = 1.3*np.outer(np.ones_like(u), v)
+
+    theta, phi = view
+    shape = x.shape
+    points = np.matrix([x.flatten(), y.flatten(), z.flatten()])
+    points = Rz(phi)*Ry(theta)*points
+    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_jitter(ax, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0)):
+    """
+    Represent jitter as a set of shapes at different orientations.
+    """
+    # set max diagonal to 0.95
+    scale = 0.95/sqrt(sum(v**2 for v in size))
+    size = tuple(scale*v for v in size)
+    draw_shape = draw_parallelepiped
+    #draw_shape = draw_ellipsoid
+
+    #np.random.seed(10)
+    #cloud = np.random.randn(10,3)
+    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],
+    ]
+    dtheta, dphi, dpsi = jitter
+    if dtheta == 0:
+        cloud = [v for v in cloud if v[0] == 0]
+    if dphi == 0:
+        cloud = [v for v in cloud if v[1] == 0]
+    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)
+    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)
+    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)
+
+def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1):
+    """Draw an ellipsoid."""
+    a,b,c = size
+    u = np.linspace(0, 2 * np.pi, steps)
+    v = np.linspace(0, np.pi, steps)
+    x = a*np.outer(np.cos(u), np.sin(v))
+    y = b*np.outer(np.sin(u), np.sin(v))
+    z = c*np.outer(np.ones_like(u), np.cos(v))
+    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, [
+         ('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_parallelepiped(ax, size, view, jitter, steps=None, alpha=1):
+    """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])
+    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
+    ])
+
+    x, y, z = transform_xyz(view, jitter, x, y, z)
+    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha)
+
+    # Draw pink face on box.
+    # Since I can't control face color, instead draw a thin box situated just
+    # in front of the "a+" 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, y, z = transform_xyz(view, jitter, abs(x)*1.05, y, z)
+        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]),
+    ])
+
+def draw_sphere(ax, radius=10., steps=100):
+    """Draw a sphere"""
+    u = np.linspace(0, 2 * np.pi, steps)
+    v = np.linspace(0, np.pi, steps)
+
+    x = radius * np.outer(np.cos(u), np.sin(v))
+    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')
+
+PROJECTIONS = [
+    # in order of PROJECTION number; do not change without updating the
+    # constants in kernel_iq.c
+    'equirectangular', 'sinusoidal', 'guyou', 'azimuthal_equidistance',
+    'azimuthal_equal_area',
+]
+def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian',
+              projection='equirectangular'):
+    """
+    Draw the dispersion mesh showing the theta-phi orientations at which
+    the model will be evaluated.
+
+    jitter projections
+    <https://en.wikipedia.org/wiki/List_of_map_projections>
+
+    equirectangular (standard latitude-longitude mesh)
+        <https://en.wikipedia.org/wiki/Equirectangular_projection>
+        Allows free movement in phi (around the equator), but theta is
+        limited to +/- 90, and points are cos-weighted. Jitter in phi is
+        uniform in weight along a line of latitude.  With small theta and
+        phi ranging over +/- 180 this forms a wobbling disk.  With small
+        phi and theta ranging over +/- 90 this forms a wedge like a slice
+        of an orange.
+    azimuthal_equidistance (Postel)
+        <https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection>
+        Preserves distance from center, and so is an excellent map for
+        representing a bivariate gaussian on the surface.  Theta and phi
+        operate identically, cutting wegdes from the antipode of the viewing
+        angle.  This unfortunately does not allow free movement in either
+        theta or phi since the orthogonal wobble decreases to 0 as the body
+        rotates through 180 degrees.
+    sinusoidal (Sanson-Flamsteed, Mercator equal-area)
+        <https://en.wikipedia.org/wiki/Sinusoidal_projection>
+        Preserves arc length with latitude, giving bad behaviour at
+        theta near +/- 90.  Theta and phi operate somewhat differently,
+        so a system with a-b-c dtheta-dphi-dpsi will not give the same
+        value as one with b-a-c dphi-dtheta-dpsi, as would be the case
+        for azimuthal equidistance.  Free movement using theta or phi
+        uniform over +/- 180 will work, but not as well as equirectangular
+        phi, with theta being slightly worse.  Computationally it is much
+        cheaper for wide theta-phi meshes since it excludes points which
+        lie outside the sinusoid near theta +/- 90 rather than packing
+        them close together as in equirectangle.  Note that the poles
+        will be slightly overweighted for theta > 90 with the circle
+        from theta at 90+dt winding backwards around the pole, overlapping
+        the circle from theta at 90-dt.
+    Guyou (hemisphere-in-a-square) **not weighted**
+        <https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection>
+        With tiling, allows rotation in phi or theta through +/- 180, with
+        uniform spacing.  Both theta and phi allow free rotation, with wobble
+        in the orthogonal direction reasonably well behaved (though not as
+        good as equirectangular phi). The forward/reverse transformations
+        relies on elliptic integrals that are somewhat expensive, so the
+        behaviour has to be very good to justify the cost and complexity.
+        The weighting function for each point has not yet been computed.
+        Note: run the module *guyou.py* directly and it will show the forward
+        and reverse mappings.
+    azimuthal_equal_area  **incomplete**
+        <https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>
+        Preserves the relative density of the surface patches.  Not that
+        useful and not completely implemented
+    Gauss-Kreuger **not implemented**
+        <https://en.wikipedia.org/wiki/Transverse_Mercator_projection#Ellipsoidal_transverse_Mercator>
+        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)
+    if dist == 'gaussian':
+        t *= 3
+        weights = exp(-0.5*t**2)
+    elif dist == 'rectangle':
+        # Note: uses sasmodels ridiculous definition of rectangle width
+        t *= sqrt(3)
+    elif dist == 'uniform':
+        pass
+    else:
+        raise ValueError("expected dist to be gaussian, rectangle or uniform")
+
+    if projection == 'equirectangular':  #define PROJECTION 1
+        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)))
+    elif projection == 'sinusoidal':  #define PROJECTION 2
+        def rotate(theta_i, phi_j):
+            latitude = theta_i
+            scale = cos(radians(latitude))
+            longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0
+            #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):
+            latitude = theta_i
+            scale = cos(radians(latitude))
+            w = 1 if abs(phi_j) < abs(scale)*180 else 0
+            return w*wi*wj
+    elif projection == 'guyou':  #define PROJECTION 3  (eventually?)
+        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):
+            return wi*wj
+    elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry
+        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):
+            # Weighting for each point comes from the integral:
+            #     \int\int I(q, lat, log) sin(lat) dlat dlog
+            # We are doing a conformal mapping from disk to sphere, so we need
+            # a change of variables g(theta, phi) -> (lat, long):
+            #     lat, long = sqrt(theta^2 + phi^2), arctan(phi/theta)
+            # giving:
+            #     dtheta dphi = det(J) dlat dlong
+            # where J is the jacobian from the partials of g. Using
+            #     R = sqrt(theta^2 + phi^2),
+            # then
+            #     J = [[x/R, Y/R], -y/R^2, x/R^2]]
+            # and
+            #     det(J) = 1/R
+            # with the final integral being:
+            #    \int\int I(q, theta, phi) sin(R)/R dtheta dphi
+            #
+            # This does approximately the right thing, decreasing the weight
+            # of each point as you go farther out on the disk, but it hasn't
+            # yet been checked against the 1D integral results. Prior
+            # to declaring this "good enough" and checking that integrals
+            # work in practice, we will examine alternative mappings.
+            #
+            # The issue is that the mapping does not support the case of free
+            # rotation about a single axis correctly, with a small deviation
+            # in the orthogonal axis independent of the first axis.  Like the
+            # usual polar coordiates integration, the integrated sections
+            # form wedges, though at least in this case the wedge cuts through
+            # 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
+    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))
+            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):
+            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
+    else:
+        raise ValueError("unknown projection %r"%projection)
+
+    # mesh in theta, phi formed by rotating z
+    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])
+    # 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)])
+    #print(max(w), min(w), min(w[w>0]))
+    points = points[:, w>0]
+    w = w[w>0]
+    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)
+
+    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):
+    """
+    Draw text at a particular location.
+    """
+    labels, locations, orientations = zip(*text)
+    px, py, pz = zip(*locations)
+    dx, dy, dz = zip(*orientations)
+
+    px, py, pz = transform_xyz(view, jitter, px, py, pz)
+    dx, dy, dz = transform_xyz(view, jitter, dx, dy, dz)
+
+    # TODO: zdir for labels is broken, and labels aren't appearing.
+    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)
+
+# Definition of rotation matrices comes from wikipedia:
+#    https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations
+def Rx(angle):
+    """Construct a matrix to rotate points about *x* by *angle* degrees."""
+    a = radians(angle)
+    R = [[1, 0, 0],
+         [0, +cos(a), -sin(a)],
+         [0, +sin(a), +cos(a)]]
+    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)],
+         [0, 1, 0],
+         [-sin(a), 0, +cos(a)]]
+    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],
+         [0, 0, 1]]
+    return np.matrix(R)
+
+def transform_xyz(view, jitter, x, y, z):
+    """
+    Send a set of (x,y,z) points through the jitter and view transforms.
+    """
+    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 = apply_jitter(jitter, points)
+    points = orient_relative_to_beam(view, points)
+    x, y, z = [np.array(v).reshape(shape) for v in points]
+    return x, y, z
+
+def apply_jitter(jitter, points):
+    """
+    Apply the jitter transform to a set of points.
+
+    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple.
+    """
+    dtheta, dphi, dpsi = jitter
+    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points
+    return points
+
+def orient_relative_to_beam(view, points):
+    """
+    Apply the view transform to a set of points.
+
+    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple.
+    """
+    theta, phi, psi = view
+    points = Rz(phi)*Ry(theta)*Rz(psi)*points
+    return points
+
+# translate between number of dimension of dispersity and the number of
+# points along each dimension.
+PD_N_TABLE = {
+    (0, 0, 0): (0, 0, 0),     # 0
+    (1, 0, 0): (100, 0, 0),   # 100
+    (0, 1, 0): (0, 100, 0),
+    (0, 0, 1): (0, 0, 100),
+    (1, 1, 0): (30, 30, 0),   # 900
+    (1, 0, 1): (30, 0, 30),
+    (0, 1, 1): (0, 30, 30),
+    (1, 1, 1): (15, 15, 15),  # 3375
+}
+
+def clipped_range(data, portion=1.0, mode='central'):
+    """
+    Determine range from data.
+
+    If *portion* is 1, use full range, otherwise use the center of the range
+    or the top of the range, depending on whether *mode* is 'central' or 'top'.
+    """
+    if portion == 1.0:
+        return data.min(), data.max()
+    elif mode == 'central':
+        data = np.sort(data.flatten())
+        offset = int(portion*len(data)/2 + 0.5)
+        return data[offset], data[-offset]
+    elif mode == 'top':
+        data = np.sort(data.flatten())
+        offset = int(portion*len(data) + 0.5)
+        return data[offset], data[-1]
+
+def draw_scattering(calculator, ax, view, jitter, dist='gaussian'):
+    """
+    Plot the scattering for the particular view.
+
+    *calculator* is returned from :func:`build_model`.  *ax* 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
+    weight distributions.
+    """
+    if dist == 'uniform':  # uniform is not yet in this branch
+        dist, scale = 'rectangle', 1/sqrt(3)
+    else:
+        scale = 1
+
+    # add the orientation parameters to the model parameters
+    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)]
+    ## 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)]
+
+    pars = dict(
+        theta=theta, theta_pd=theta_pd, theta_pd_type=dist, theta_pd_n=theta_pd_n,
+        phi=phi, phi_pd=phi_pd, phi_pd_type=dist, phi_pd_n=phi_pd_n,
+        psi=psi, psi_pd=psi_pd, psi_pd_type=dist, psi_pd_n=psi_pd_n,
+    )
+    pars.update(calculator.pars)
+
+    # compute the pattern
+    qx, qy = calculator._data.x_bins, calculator._data.y_bins
+    Iqxy = calculator(**pars).reshape(len(qx), len(qy))
+
+    # scale it and draw it
+    Iqxy = np.log(Iqxy)
+    if calculator.limits:
+        # use limits from orientation (0,0,0)
+        vmin, vmax = calculator.limits
+    else:
+        vmax = Iqxy.max()
+        vmin = vmax*10**-7
+        #vmin, vmax = clipped_range(Iqxy, portion=portion, mode='top')
+    #print("range",(vmin,vmax))
+    #qx, qy = np.meshgrid(qx, qy)
+    if 0:
+        level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i')
+        level[level<0] = 0
+        colors = plt.get_cmap()(level)
+        ax.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,
+                    levels=np.linspace(vmin, vmax, 24))
+    else:
+        ax.pcolormesh(qx, qy, Iqxy)
+
+def build_model(model_name, n=150, qmax=0.5, **pars):
+    """
+    Build a calculator for the given shape.
+
+    *model_name* is any sasmodels model.  *n* and *qmax* define an n x n mesh
+    on which to evaluate the model.  The remaining parameters are stored in
+    the returned calculator as *calculator.pars*.  They are used by
+    :func:`draw_scattering` to set the non-orientation parameters in the
+    calculation.
+
+    Returns a *calculator* function which takes a dictionary or parameters and
+    produces Iqxy.  The Iqxy value needs to be reshaped to an n x n matrix
+    for plotting.  See the :class:`sasmodels.direct_model.DirectModel` class
+    for details.
+    """
+    from sasmodels.core import load_model_info, build_model
+    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!')
+    q = np.linspace(-qmax, qmax, n)
+    data = empty_data2D(q, q)
+    calculator = DirectModel(data, model)
+
+    # stuff the values for non-orientation parameters into the calculator
+    calculator.pars = pars.copy()
+    calculator.pars.setdefault('backgound', 1e-3)
+
+    # fix the data limits so that we can see if the pattern fades
+    # under rotation or angular dispersion
+    Iqxy = calculator(theta=0, phi=0, psi=0, **calculator.pars)
+    Iqxy = np.log(Iqxy)
+    vmin, vmax = clipped_range(Iqxy, 0.95, mode='top')
+    calculator.limits = vmin, vmax+1
+
+    return calculator
+
+def select_calculator(model_name, n=150, size=(10,40,100)):
+    """
+    Create a model calculator for the given shape.
+
+    *model_name* is one of sphere, cylinder, ellipsoid, triaxial_ellipsoid,
+    parallelepiped or bcc_paracrystal. *n* is the number of points to use
+    in the q range.  *qmax* is chosen based on model parameters for the
+    given model to show something intersting.
+
+    Returns *calculator* and tuple *size* (a,b,c) giving minor and major
+    equitorial axes and polar axis respectively.  See :func:`build_model`
+    for details on the returned calculator.
+    """
+    a, b, c = size
+    if model_name == 'sphere':
+        calculator = build_model('sphere', n=n, radius=c)
+        a = b = c
+    elif model_name == 'bcc_paracrystal':
+        calculator = build_model('bcc_paracrystal', n=n, dnn=c,
+                                  d_factor=0.06, radius=40)
+        a = b = c
+    elif model_name == 'cylinder':
+        calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c)
+        a = b
+    elif model_name == 'ellipsoid':
+        calculator = build_model('ellipsoid', n=n, qmax=1.0,
+                                 radius_polar=c, radius_equatorial=b)
+        a = b
+    elif model_name == 'triaxial_ellipsoid':
+        calculator = build_model('triaxial_ellipsoid', n=n, qmax=0.5,
+                                 radius_equat_minor=a,
+                                 radius_equat_major=b,
+                                 radius_polar=c)
+    elif model_name == 'parallelepiped':
+        calculator = build_model('parallelepiped', n=n, a=a, b=b, c=c)
+    else:
+        raise ValueError("unknown model %s"%model_name)
+
+    return calculator, (a, b, c)
+
+SHAPES = [
+    'parallelepiped', 'triaxial_ellipsoid', 'bcc_paracrystal',
+    'cylinder', 'ellipsoid',
+    'sphere',
+ ]
+
+DISTRIBUTIONS = [
+    'gaussian', 'rectangle', 'uniform',
+]
+DIST_LIMITS = {
+    'gaussian': 30,
+    'rectangle': 90/sqrt(3),
+    'uniform': 90,
+}
+
+def run(model_name='parallelepiped', size=(10, 40, 100),
+        dist='gaussian', mesh=30,
+        projection='equirectangular'):
+    """
+    Show an interactive orientation and jitter demo.
+
+    *model_name* is one of: sphere, cylinder, ellipsoid, parallelepiped,
+    triaxial_ellipsoid, or bcc_paracrystal.
+
+    *size* gives the dimensions (a, b, c) of the shape.
+
+    *dist* is the type of dispersition: gaussian, rectangle, or uniform.
+
+    *mesh* is the number of points in the dispersion mesh.
+
+    *projection* is the map projection to use for the mesh: equirectangular,
+    sinusoidal, guyou, azimuthal_equidistance, or azimuthal_equal_area.
+    """
+    # projection number according to 1-order position in list, but
+    # only 1 and 2 are implemented so far.
+    from sasmodels import generate
+    generate.PROJECTION = PROJECTIONS.index(projection) + 1
+    if generate.PROJECTION > 2:
+        print("*** PROJECTION %s not implemented in scattering function ***"%projection)
+        generate.PROJECTION = 2
+
+    # set up calculator
+    calculator, size = select_calculator(model_name, n=150, size=size)
+
+    ## uncomment to set an independent the colour range for every view
+    ## If left commented, the colour range is fixed for all views
+    calculator.limits = None
+
+    ## initial view
+    #theta, dtheta = 70., 10.
+    #phi, dphi = -45., 3.
+    #psi, dpsi = -45., 3.
+    theta, phi, psi = 0, 0, 0
+    dtheta, dphi, dpsi = 0, 0, 0
+
+    ## create the plot window
+    #plt.hold(True)
+    plt.subplots(num=None, figsize=(5.5, 5.5))
+    plt.set_cmap('gist_earth')
+    plt.clf()
+    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')
+    try:  # CRUFT: not all versions of matplotlib accept 'square' 3d projection
+        ax.axis('square')
+    except Exception:
+        pass
+
+    axcolor = 'lightgoldenrodyellow'
+
+    ## 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)
+    # 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)
+
+    ## callback to draw the new view
+    def update(val, axis=None):
+        view = stheta.val, sphi.val, spsi.val
+        jitter = sdtheta.val, sdphi.val, sdpsi.val
+        # set small jitter as 0 if multiple pd dims
+        dims = sum(v > 0 for v in jitter)
+        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_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)
+        plt.gcf().canvas.draw()
+
+    ## bind control widgets to view updater
+    stheta.on_changed(lambda v: update(v,'theta'))
+    sphi.on_changed(lambda v: update(v, 'phi'))
+    spsi.on_changed(lambda v: update(v, 'psi'))
+    sdtheta.on_changed(lambda v: update(v, 'dtheta'))
+    sdphi.on_changed(lambda v: update(v, 'dphi'))
+    sdpsi.on_changed(lambda v: update(v, 'dpsi'))
+
+    ## initialize view
+    update(None, 'phi')
+
+    ## go interactive
+    plt.show()
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Display jitter",
+        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')
+    opts = parser.parse_args()
+    size = tuple(int(v) for v in opts.size.split(','))
+    run(opts.shape, size=size,
+        mesh=opts.mesh, dist=opts.distribution,
+        projection=opts.projection)
+
+if __name__ == "__main__":
+    main()
