Index: explore/shimmy.py
===================================================================
--- explore/shimmy.py	(revision 782dd1f27298956f047ebe2379cedea7f07f4de0)
+++ explore/shimmy.py	(revision 782dd1f27298956f047ebe2379cedea7f07f4de0)
@@ -0,0 +1,237 @@
+"""
+Application to explore the difference between sasview 3.x orientation
+dispersity and possible replacement algorithms.
+"""
+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):
+    #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 = np.outer(np.ones_like(u), v)
+
+    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5)
+    
+def draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi):
+    size=[0.1, 0.4, 1.0] 
+    view=[theta, phi, psi]
+    shimmy=[0,0,0]
+    #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],
+    ]
+    draw_shape(ax, size, view, shimmy, steps=100, alpha=0.8)
+    cloud = []
+    for point in cloud:
+        shimmy=[dtheta*point[0], dphi*point[1], dpsi*point[2]]
+        draw_shape(ax, size, view, shimmy, 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, shimmy, steps=25, alpha=1):
+    a,b,c = size
+    theta, phi, psi = view
+    dtheta, dphi, dpsi = shimmy
+
+    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))
+
+    shape = x.shape
+    points = np.matrix([x.flatten(),y.flatten(),z.flatten()])
+    points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points
+    points = Rz(phi)*Ry(theta)*Rz(psi)*points
+    x,y,z = [v.reshape(shape) for v in points]
+
+    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha)
+
+def draw_parallelepiped(ax, size, view, shimmy, alpha=1):
+    a,b,c = size
+    theta, phi, psi = view
+    dtheta, dphi, dpsi = shimmy
+
+    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
+    ])
+
+    points = np.matrix([x,y,z])
+    points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points
+    points = Rz(phi)*Ry(theta)*Rz(psi)*points
+	
+    x,y,z = [np.array(v).flatten() for v in points]
+    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha)
+
+def draw_sphere(ax, radius=10., steps=100):
+    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')
+
+def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'):
+    theta_center = radians(theta)
+    phi_center = radians(phi)
+    flow_center = radians(flow)
+    dtheta = radians(dtheta)
+    dphi = radians(dphi)
+
+    # 10 point 3-sigma gaussian weights
+    t = np.linspace(-3., 3., 11)
+    if dist == 'gauss':
+        weights = exp(-0.5*t**2)
+    elif dist == 'rect':
+        weights = np.ones_like(t)
+    else:
+        raise ValueError("expected dist to be 'gauss' or 'rect'")
+    theta = dtheta*t
+    phi = dphi*t
+
+    x = radius * np.outer(cos(phi), cos(theta))
+    y = radius * np.outer(sin(phi), cos(theta))
+    z = radius * np.outer(np.ones_like(phi), sin(theta))
+    #w = np.outer(weights, weights*abs(cos(dtheta*t)))
+    w = np.outer(weights, weights*abs(cos(theta)))
+
+    x, y, z, w = [v.flatten() for v in (x,y,z,w)]
+    x, y, z = rotate(x, y, z, phi_center, theta_center, flow_center)
+
+    ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.)
+
+def rotate(x, y, z, phi, theta, psi):
+    R = Rz(phi)*Ry(theta)*Rz(psi)
+    p = np.vstack([x,y,z])
+    return R*p
+
+def Rx(angle):
+    a = radians(angle)
+    R = [[1., 0., 0.],
+         [0.,  cos(a), sin(a)],
+         [0., -sin(a), cos(a)]]
+    return np.matrix(R)
+
+def Ry(angle):
+    a = radians(angle)
+    R = [[cos(a), 0., -sin(a)],
+         [0., 1., 0.],
+         [sin(a), 0.,  cos(a)]]
+    return np.matrix(R)
+
+def Rz(angle):
+    a = radians(angle)
+    R = [[cos(a), -sin(a), 0.],
+         [sin(a),  cos(a), 0.],
+         [0., 0., 1.]]
+    return np.matrix(R)
+
+def main():
+    #plt.hold(True)
+    plt.set_cmap('gist_earth')
+    plt.clf()
+    #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')
+
+    theta, dtheta = 70., 10.
+    phi, dphi = -45., 3.
+    psi, dpsi = -45., 3.
+    theta, phi, psi = 0, 0, 0
+    #dist = 'rect'
+    dist = 'gauss'
+
+    axcolor = 'lightgoldenrodyellow'
+    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', -180, 180, 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)
+    sdtheta = Slider(axdtheta, 'dTheta', 0, 30, valinit=dtheta)
+    sdphi = Slider(axdphi, 'dPhi', 0, 30, valinit=dphi)
+    sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dphi)
+
+    def update(val, axis=None):
+        theta, phi, psi = stheta.val, sphi.val, spsi.val
+        dtheta, dphi, dpsi = sdtheta.val, sdphi.val, sdpsi.val
+        ax.cla()
+        draw_beam(ax)
+        draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi)
+        #if not axis.startswith('d'):
+        #    ax.view_init(elev=theta, azim=phi)
+        plt.gcf().canvas.draw()
+
+    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'))
+
+    update(None, 'phi')
+
+    plt.show()
+
+if __name__ == "__main__":
+    main()
