Index: explore/jitter.py
===================================================================
--- explore/jitter.py	(revision 85190c2070ee71c354d018a96403e7d1e6045492)
+++ explore/jitter.py	(revision 85190c2070ee71c354d018a96403e7d1e6045492)
@@ -0,0 +1,243 @@
+"""
+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],
+    ]
+    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, shimmy, steps=100, alpha=0.8)
+    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
+    dtheta, dphi, dpsi = 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', -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)
+    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()
Index: plore/shimmy.py
===================================================================
--- explore/shimmy.py	(revision 1b693baedc46a8c69820ad2850cc105c701a1905)
+++ 	(revision )
@@ -1,243 +1,0 @@
-"""
-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],
-    ]
-    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, shimmy, steps=100, alpha=0.8)
-    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
-    dtheta, dphi, dpsi = 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', -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)
-    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()
Index: sasmodels/models/core_shell_bicelle_elliptical.c
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical.c	(revision f4f85b3a43d7e158d487593ba7eee6f3d8e98fc8)
+++ sasmodels/models/core_shell_bicelle_elliptical.c	(revision 44e8a93e1773cb4d0897a436b8770e1d4e5ee343)
@@ -1,45 +1,24 @@
-double form_volume(double radius, double x_core, double thick_rim, double thick_face, double length);
-double Iq(double q,
-          double radius,
-          double x_core,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld);
-
-
-double Iqxy(double qx, double qy,
-          double radius,
-          double x_core,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld,
-          double theta,
-          double phi,
-          double psi);
-
 // NOTE that "length" here is the full height of the core!
-double form_volume(double radius, double x_core, double thick_rim, double thick_face, double length)
+static double
+form_volume(double r_minor,
+        double x_core,
+        double thick_rim,
+        double thick_face,
+        double length)
 {
-    return M_PI*(radius+thick_rim)*(radius*x_core+thick_rim)*(length+2.0*thick_face);
+    return M_PI*(r_minor+thick_rim)*(r_minor*x_core+thick_rim)*(length+2.0*thick_face);
 }
 
-double Iq(double q,
-          double rad,
-          double x_core,
-          double radthick,
-          double facthick,
-          double length,
-          double rhoc,
-          double rhoh,
-          double rhor,
-          double rhosolv)
+static double
+Iq(double q,
+        double r_minor,
+        double x_core,
+        double thick_rim,
+        double thick_face,
+        double length,
+        double rhoc,
+        double rhoh,
+        double rhor,
+        double rhosolv)
 {
     double si1,si2,be1,be2;
@@ -54,13 +33,13 @@
     //const double vbj=M_PI;
 
-    const double radius_major = rad * x_core;
-    const double rA = 0.5*(square(radius_major) + square(rad));
-    const double rB = 0.5*(square(radius_major) - square(rad));
-    const double dr1 = (rhoc-rhoh)   *M_PI*rad*radius_major*(2.0*halfheight);;
-    const double dr2 = (rhor-rhosolv)*M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick);
-    const double dr3 = (rhoh-rhor)   *M_PI*rad*radius_major*2.0*(halfheight+facthick);
-    //const double vol1 = M_PI*rad*radius_major*(2.0*halfheight);
-    //const double vol2 = M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick);
-    //const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick);
+    const double r_major = r_minor * x_core;
+    const double rA = 0.5*(square(r_major) + square(r_minor));
+    const double rB = 0.5*(square(r_major) - square(r_minor));
+    const double dr1 = (rhoc-rhoh)   *M_PI*r_minor*r_major*(2.0*halfheight);;
+    const double dr2 = (rhor-rhosolv)*M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);
+    const double dr3 = (rhoh-rhor)   *M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
+    //const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight);
+    //const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);
+    //const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
 
     //initialize integral
@@ -74,5 +53,5 @@
         double inner_sum=0;
         double sinarg1 = q*halfheight*cos_alpha;
-        double sinarg2 = q*(halfheight+facthick)*cos_alpha;
+        double sinarg2 = q*(halfheight+thick_face)*cos_alpha;
         si1 = sas_sinx_x(sinarg1);
         si2 = sas_sinx_x(sinarg2);
@@ -83,5 +62,5 @@
             const double rr = sqrt(rA - rB*cos(beta));
             double besarg1 = q*rr*sin_alpha;
-            double besarg2 = q*(rr+radthick)*sin_alpha;
+            double besarg2 = q*(rr+thick_rim)*sin_alpha;
             be1 = sas_2J1x_x(besarg1);
             be2 = sas_2J1x_x(besarg2);
@@ -97,10 +76,10 @@
 }
 
-double 
+static double
 Iqxy(double qx, double qy,
-          double rad,
+          double r_minor,
           double x_core,
-          double radthick,
-          double facthick,
+          double thick_rim,
+          double thick_face,
           double length,
           double rhoc,
@@ -118,24 +97,19 @@
     const double dr2 = rhor-rhosolv;
     const double dr3 = rhoh-rhor;
-    const double radius_major = rad*x_core;
+    const double r_major = r_minor*x_core;
     const double halfheight = 0.5*length;
-    const double vol1 = M_PI*rad*radius_major*length;
-    const double vol2 = M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick);
-    const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick);
+    const double vol1 = M_PI*r_minor*r_major*length;
+    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);
+    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
 
-    // Compute:  r = sqrt((radius_major*zhat)^2 + (radius_minor*yhat)^2)
-    // Given:    radius_major = r_ratio * radius_minor  
-    // ASSUME the sin_alpha is included in the separate integration over orientation of rod angle
-    const double rad_minor = rad;
-    const double rad_major = rad*x_core;
-    const double r_hat = sqrt(square(rad_major*xhat) + square(rad_minor*yhat));
-    const double rshell_hat = sqrt(square((rad_major+radthick)*xhat)
-                                   + square((rad_minor+radthick)*yhat));
+    // Compute effective radius in rotated coordinates
+    const double r_hat = sqrt(square(r_major*xhat) + square(r_minor*yhat));
+    const double rshell_hat = sqrt(square((r_major+thick_rim)*xhat)
+                                   + square((r_minor+thick_rim)*yhat));
     const double be1 = sas_2J1x_x( q*r_hat );
     const double be2 = sas_2J1x_x( q*rshell_hat );
     const double si1 = sas_sinx_x( q*halfheight*zhat );
-    const double si2 = sas_sinx_x( q*(halfheight + facthick)*zhat );
+    const double si2 = sas_sinx_x( q*(halfheight + thick_face)*zhat );
     const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1);
-    //const double vol = form_volume(radius_minor, r_ratio, length);
     return 1.0e-4 * Aq;
 }
Index: sasmodels/models/parallelepiped.py
===================================================================
--- sasmodels/models/parallelepiped.py	(revision 9b79f2966c4083a48b4894b6bf0587bef7a5be09)
+++ sasmodels/models/parallelepiped.py	(revision 9b79f2966c4083a48b4894b6bf0587bef7a5be09)
@@ -23,6 +23,6 @@
 
    The edge of the solid used to have to satisfy the condition that $A < B < C$.
-   After some improvements to the effective radius calculation, used with an S(Q),
-   it is beleived that this is no longer the case.
+   After some improvements to the effective radius calculation, used with
+   an S(Q), it is beleived that this is no longer the case.
 
 The 1D scattering intensity $I(q)$ is calculated as:
@@ -72,6 +72,6 @@
 
 NB: The 2nd virial coefficient of the parallelepiped is calculated based on
-the averaged effective radius, after appropriately
-sorting the three dimensions, to give an oblate or prolate particle, $(=\sqrt{A B / \pi})$ and
+the averaged effective radius, after appropriately sorting the three
+dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and
 length $(= C)$ values, and used as the effective radius for
 $S(q)$ when $P(q) \cdot S(q)$ is applied.
@@ -106,5 +106,5 @@
 .. figure:: img/parallelepiped_angle_definition.png
 
-    Definition of the angles for oriented parallelepiped, shown with $A < B < C$.
+    Definition of the angles for oriented parallelepiped, shown with $A<B<C$.
 
 .. figure:: img/parallelepiped_angle_projection.png
@@ -167,6 +167,6 @@
 ----------------------------
 
-* **Author:** This model is based on form factor calculations implemented in a c-library
-provided by the NIST Center for Neutron Research (Kline, 2006).
+* **Author:** This model is based on form factor calculations implemented
+in a c-library provided by the NIST Center for Neutron Research (Kline, 2006).
 * **Last Modified by:**  Paul Kienzle **Date:** April 05, 2017
 * **Last Reviewed by:**  Richard Heenan **Date:** April 06, 2017
Index: sasmodels/models/sc_paracrystal.py
===================================================================
--- sasmodels/models/sc_paracrystal.py	(revision 9b79f2966c4083a48b4894b6bf0587bef7a5be09)
+++ sasmodels/models/sc_paracrystal.py	(revision 9b79f2966c4083a48b4894b6bf0587bef7a5be09)
@@ -85,5 +85,5 @@
 .. figure:: img/parallelepiped_angle_definition.png
 
-    Orientation of the crystal with respect to the scattering plane, when 
+    Orientation of the crystal with respect to the scattering plane, when
     $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis).
 
Index: sasmodels/models/stacked_disks.py
===================================================================
--- sasmodels/models/stacked_disks.py	(revision 9b79f2966c4083a48b4894b6bf0587bef7a5be09)
+++ sasmodels/models/stacked_disks.py	(revision 9b79f2966c4083a48b4894b6bf0587bef7a5be09)
@@ -58,5 +58,5 @@
 and $\sigma_d$ = the Gaussian standard deviation of the d-spacing (*sigma_d*).
 Note that $D\cos(\alpha)$ is the component of $D$ parallel to $q$ and the last
-term in the equation above is effectively a Debye-Waller factor term. 
+term in the equation above is effectively a Debye-Waller factor term.
 
 .. note::
@@ -158,5 +158,7 @@
 tests = [
 # Accuracy tests based on content in test/utest_extra_models.py.
-# Added 2 tests with n_stacked = 5 using SasView 3.1.2 - PDB; for which alas q=0.001 values seem closer to n_stacked=1 not 5, changed assuming my 4.1 code OK, RKH
+# Added 2 tests with n_stacked = 5 using SasView 3.1.2 - PDB;
+# for which alas q=0.001 values seem closer to n_stacked=1 not 5,
+# changed assuming my 4.1 code OK, RKH
     [{'thick_core': 10.0,
       'thick_layer': 15.0,
Index: sasmodels/models/triaxial_ellipsoid.py
===================================================================
--- sasmodels/models/triaxial_ellipsoid.py	(revision 9b79f2966c4083a48b4894b6bf0587bef7a5be09)
+++ sasmodels/models/triaxial_ellipsoid.py	(revision 9b79f2966c4083a48b4894b6bf0587bef7a5be09)
@@ -16,9 +16,11 @@
     \frac{X^2}{R_a^2} + \frac{Y^2}{R_b^2} + \frac{Z^2}{R_c^2} = 1
 
-the scattering for randomly oriented particles is defined by the average over all orientations $\Omega$ of:
+the scattering for randomly oriented particles is defined by the average over
+all orientations $\Omega$ of:
 
 .. math::
 
-    P(q) = \text{scale}(\Delta\rho)^2\frac{V}{4 \pi}\int_\Omega \Phi^2(qr) d\Omega + \text{background}
+    P(q) = \text{scale}(\Delta\rho)^2\frac{V}{4 \pi}\int_\Omega\Phi^2(qr)\,d\Omega
+           + \text{background}
 
 where
@@ -38,5 +40,6 @@
  .. math::
 
-     \langle\Phi^2\rangle = \int_0^{2\pi} \int_{-\pi/2}^{\pi/2} \Phi^2(qr) \cos \gamma\,d\gamma d\phi
+     \langle\Phi^2\rangle = \int_0^{2\pi} \int_{-\pi/2}^{\pi/2} \Phi^2(qr)
+                                                \cos \gamma\,d\gamma d\phi
 
 with $e = \cos\gamma \sin\phi$, $f = \cos\gamma \cos\phi$ and $g = \sin\gamma$.
@@ -69,6 +72,7 @@
 .. figure:: img/elliptical_cylinder_angle_definition.png
 
-    Definition of angles for oriented triaxial ellipsoid, where radii shown here are $a < b << c$
-    and angle $\Psi$ is a rotation around the axis of the particle.
+    Definition of angles for oriented triaxial ellipsoid, where radii shown
+    here are $a < b << c$ and angle $\Psi$ is a rotation around the axis
+    of the particle.
 
 The angle $\psi$ is the rotational angle around its own $c$ axis
@@ -115,6 +119,4 @@
 * **Last Reviewed by:** Paul Kienzle & Richard Heenan **Date:**  April 4, 2017
 
-"""
-
 from numpy import inf, sin, cos, pi
 
@@ -122,5 +124,5 @@
 title = "Ellipsoid of uniform scattering length density with three independent axes."
 
-description = """\
+description = """
 Note: During fitting ensure that the inequality ra<rb<rc is not
 	violated. Otherwise the calculation will
@@ -157,6 +159,6 @@
     from .ellipsoid import ER as ellipsoid_ER
 
-    # now that radii can be in any size order, radii need sorting a,b,c where a~b and c is either much smaller
-    # or much larger
+    # now that radii can be in any size order, radii need sorting a,b,c
+    # where a~b and c is either much smaller or much larger
     radii = np.vstack((radius_equat_major, radius_equat_minor, radius_polar))
     radii = np.sort(radii, axis=0)
@@ -178,5 +180,6 @@
 
 q = 0.1
-# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
+# april 6 2017, rkh add unit tests
+#     NOT compared with any other calc method, assume correct!
 # add 2d test after pull #890
 qx = q*cos(pi/6.0)
