Index: explore/angular_pd.py
===================================================================
--- explore/angular_pd.py	(revision 12eb36b023c96fc01e95ceccb2c48bf9a52055ae)
+++ explore/angular_pd.py	(revision 8267e0b2278b86ad580015ef3fe90b1d7cf31109)
@@ -47,5 +47,5 @@
 
 def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'):
-    theta_center = radians(theta)
+    theta_center = radians(90-theta)
     phi_center = radians(phi)
     flow_center = radians(flow)
@@ -137,5 +137,5 @@
                               radius=11., dist=dist)
         if not axis.startswith('d'):
-            ax.view_init(elev=theta, azim=phi)
+            ax.view_init(elev=90-theta if use_new else theta, azim=phi)
         plt.gcf().canvas.draw()
 
Index: explore/shimmy.py
===================================================================
--- explore/shimmy.py	(revision 1b693baedc46a8c69820ad2850cc105c701a1905)
+++ explore/shimmy.py	(revision 1b693baedc46a8c69820ad2850cc105c701a1905)
@@ -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: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision c4e321550a7eecc695933ac3c012bddd93f534d7)
+++ sasmodels/compare.py	(revision c4e321550a7eecc695933ac3c012bddd93f534d7)
@@ -84,5 +84,5 @@
     -edit starts the parameter explorer
     -default/-demo* use demo vs default parameters
-    -html shows the model docs instead of running the model
+    -help/-html shows the model docs instead of running the model
     -title="note" adds note to the plot title, after the model name
     -data="path" uses q, dq from the data file
@@ -832,5 +832,5 @@
     'linear', 'log', 'q4',
     'hist', 'nohist',
-    'edit', 'html',
+    'edit', 'html', 'help',
     'demo', 'default',
     ])
@@ -1001,4 +1001,5 @@
         elif arg == '-default':    opts['use_demo'] = False
         elif arg == '-html':    opts['html'] = True
+        elif arg == '-help':    opts['html'] = True
     # pylint: enable=bad-whitespace
 
Index: sasmodels/kernel_header.c
===================================================================
--- sasmodels/kernel_header.c	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
+++ sasmodels/kernel_header.c	(revision b00a646c6ab83c01056acef9b46d1e392433af0b)
@@ -148,4 +148,15 @@
 inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; }
 
+// To rotate from the canonical position to theta, phi, psi, first rotate by
+// psi about the major axis, oriented along z, which is a rotation in the
+// detector plane xy. Next rotate by theta about the y axis, aligning the major
+// axis in the xz plane. Finally, rotate by phi in the detector plane xy.
+// To compute the scattering, undo these rotations in reverse order:
+//     rotate in xy by -phi, rotate in xz by -theta, rotate in xy by -psi
+// The returned q is the length of the q vector and (xhat, yhat, zhat) is a unit
+// vector in the q direction.
+// To change between counterclockwise and clockwise rotation, change the
+// sign of phi and psi.
+
 #if 1
 //think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017
@@ -166,4 +177,24 @@
 #endif
 
+#if 1
+#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat) do { \
+    q = sqrt(qx*qx + qy*qy); \
+    const double qxhat = qx/q; \
+    const double qyhat = qy/q; \
+    double sin_theta, cos_theta; \
+    double sin_phi, cos_phi; \
+    double sin_psi, cos_psi; \
+    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \
+    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \
+    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \
+    xhat = qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi) \
+         + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi); \
+    yhat = qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi) \
+         + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi); \
+    zhat = qxhat*(-sin_theta*cos_phi) \
+         + qyhat*(-sin_theta*sin_phi); \
+    } while (0)
+#else
+// SasView 3.x definition of orientation
 #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \
     q = sqrt(qx*qx + qy*qy); \
@@ -180,2 +211,3 @@
     cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \
     } while (0)
+#endif
Index: sasmodels/models/bcc_paracrystal.c
===================================================================
--- sasmodels/models/bcc_paracrystal.c	(revision 4962519cc1910422a467e4f643dba223623e0e09)
+++ sasmodels/models/bcc_paracrystal.c	(revision 50beefeac18f2f952bf0863526ec9bbd3009e29c)
@@ -90,10 +90,10 @@
     double theta, double phi, double psi)
 {
-    double q, cos_a1, cos_a2, cos_a3;
-    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1);
+    double q, zhat, yhat, xhat;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
 
-    const double a1 = +cos_a3 - cos_a1 + cos_a2;
-    const double a2 = +cos_a3 + cos_a1 - cos_a2;
-    const double a3 = -cos_a3 + cos_a1 + cos_a2;
+    const double a1 = +xhat - zhat + yhat;
+    const double a2 = +xhat + zhat - yhat;
+    const double a3 = -xhat + zhat + yhat;
 
     const double qd = 0.5*q*dnn;
Index: sasmodels/models/core_shell_bicelle.c
===================================================================
--- sasmodels/models/core_shell_bicelle.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
+++ sasmodels/models/core_shell_bicelle.c	(revision b260926c42c2ae859a32c20b98935378d51fab4d)
@@ -30,9 +30,9 @@
 
 static double
-bicelle_kernel(double qq,
+bicelle_kernel(double q,
               double rad,
               double radthick,
               double facthick,
-              double length,
+              double halflength,
               double rhoc,
               double rhoh,
@@ -42,21 +42,15 @@
               double cos_alpha)
 {
-    double si1,si2,be1,be2;
-
     const double dr1 = rhoc-rhoh;
     const double dr2 = rhor-rhosolv;
     const double dr3 = rhoh-rhor;
-    const double vol1 = M_PI*rad*rad*(2.0*length);
-    const double vol2 = M_PI*(rad+radthick)*(rad+radthick)*2.0*(length+facthick);
-    const double vol3 = M_PI*rad*rad*2.0*(length+facthick);
-    double besarg1 = qq*rad*sin_alpha;
-    double besarg2 = qq*(rad+radthick)*sin_alpha;
-    double sinarg1 = qq*length*cos_alpha;
-    double sinarg2 = qq*(length+facthick)*cos_alpha;
+    const double vol1 = M_PI*square(rad)*2.0*(halflength);
+    const double vol2 = M_PI*square(rad+radthick)*2.0*(halflength+facthick);
+    const double vol3 = M_PI*square(rad)*2.0*(halflength+facthick);
 
-    be1 = sas_2J1x_x(besarg1);
-    be2 = sas_2J1x_x(besarg2);
-    si1 = sas_sinx_x(sinarg1);
-    si2 = sas_sinx_x(sinarg2);
+    const double be1 = sas_2J1x_x(q*(rad)*sin_alpha);
+    const double be2 = sas_2J1x_x(q*(rad+radthick)*sin_alpha);
+    const double si1 = sas_sinx_x(q*(halflength)*cos_alpha);
+    const double si2 = sas_sinx_x(q*(halflength+facthick)*cos_alpha);
 
     const double t = vol1*dr1*si1*be1 +
@@ -64,5 +58,5 @@
                      vol3*dr3*si2*be1;
 
-    const double retval = t*t*sin_alpha;
+    const double retval = t*t;
 
     return retval;
@@ -71,5 +65,5 @@
 
 static double
-bicelle_integration(double qq,
+bicelle_integration(double q,
                    double rad,
                    double radthick,
@@ -83,5 +77,5 @@
     // set up the integration end points
     const double uplim = M_PI_4;
-    const double halfheight = 0.5*length;
+    const double halflength = 0.5*length;
 
     double summ = 0.0;
@@ -90,8 +84,8 @@
         double sin_alpha, cos_alpha; // slots to hold sincos function output
         SINCOS(alpha, sin_alpha, cos_alpha);
-        double yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick,
-                             halfheight, rhoc, rhoh, rhor, rhosolv,
+        double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick,
+                             halflength, rhoc, rhoh, rhor, rhosolv,
                              sin_alpha, cos_alpha);
-        summ += yyy;
+        summ += yyy*sin_alpha;
     }
 
@@ -119,9 +113,6 @@
     double answer = bicelle_kernel(q, radius, thick_rim, thick_face,
                            0.5*length, core_sld, face_sld, rim_sld,
-                           solvent_sld, sin_alpha, cos_alpha) / fabs(sin_alpha);
-
-    answer *= 1.0e-4;
-
-    return answer;
+                           solvent_sld, sin_alpha, cos_alpha);
+    return 1.0e-4*answer;
 }
 
Index: sasmodels/models/core_shell_bicelle_elliptical.c
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
+++ sasmodels/models/core_shell_bicelle_elliptical.c	(revision f4f85b3a43d7e158d487593ba7eee6f3d8e98fc8)
@@ -32,15 +32,14 @@
 }
 
-double 
-                Iq(double qq,
-                   double rad,
-                   double x_core,
-                   double radthick,
-                   double facthick,
-                   double length,
-                   double rhoc,
-                   double rhoh,
-                   double rhor,
-                   double rhosolv)
+double Iq(double q,
+          double rad,
+          double x_core,
+          double radthick,
+          double facthick,
+          double length,
+          double rhoc,
+          double rhoh,
+          double rhor,
+          double rhosolv)
 {
     double si1,si2,be1,be2;
@@ -74,6 +73,6 @@
         const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
         double inner_sum=0;
-        double sinarg1 = qq*halfheight*cos_alpha;
-        double sinarg2 = qq*(halfheight+facthick)*cos_alpha;
+        double sinarg1 = q*halfheight*cos_alpha;
+        double sinarg2 = q*(halfheight+facthick)*cos_alpha;
         si1 = sas_sinx_x(sinarg1);
         si2 = sas_sinx_x(sinarg2);
@@ -83,6 +82,6 @@
             const double beta = ( Gauss76Z[j] +1.0)*M_PI_2;
             const double rr = sqrt(rA - rB*cos(beta));
-            double besarg1 = qq*rr*sin_alpha;
-            double besarg2 = qq*(rr+radthick)*sin_alpha;
+            double besarg1 = q*rr*sin_alpha;
+            double besarg2 = q*(rr+radthick)*sin_alpha;
             be1 = sas_2J1x_x(besarg1);
             be2 = sas_2J1x_x(besarg2);
@@ -114,6 +113,6 @@
 {
        // THIS NEEDS TESTING
-    double qq, cos_val, cos_mu, cos_nu;
-    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, qq, cos_val, cos_mu, cos_nu);
+    double q, xhat, yhat, zhat;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
     const double dr1 = rhoc-rhoh;
     const double dr2 = rhor-rhosolv;
@@ -125,12 +124,16 @@
     const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick);
 
-    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2)
+    // 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 r = rad*sqrt(square(x_core*cos_nu) + cos_mu*cos_mu);
-    const double be1 = sas_2J1x_x( qq*r );
-    const double be2 = sas_2J1x_x( qq*(r + radthick ) );
-    const double si1 = sas_sinx_x( qq*halfheight*cos_val );
-    const double si2 = sas_sinx_x( qq*(halfheight + facthick)*cos_val );
+    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));
+    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 Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1);
     //const double vol = form_volume(radius_minor, r_ratio, length);
Index: sasmodels/models/core_shell_parallelepiped.c
===================================================================
--- sasmodels/models/core_shell_parallelepiped.c	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
+++ sasmodels/models/core_shell_parallelepiped.c	(revision 92dfe0c540274c8b800a9929c3f1e1a7d4cba9dc)
@@ -134,6 +134,6 @@
     double psi)
 {
-    double q, cos_val_a, cos_val_b, cos_val_c;
-    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a);
+    double q, zhat, yhat, xhat;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
 
     // cspkernel in csparallelepiped recoded here
@@ -160,10 +160,10 @@
     double tc = length_a + 2.0*thick_rim_c;
     //handle arg=0 separately, as sin(t)/t -> 1 as t->0
-    double siA = sas_sinx_x(0.5*q*length_a*cos_val_a);
-    double siB = sas_sinx_x(0.5*q*length_b*cos_val_b);
-    double siC = sas_sinx_x(0.5*q*length_c*cos_val_c);
-    double siAt = sas_sinx_x(0.5*q*ta*cos_val_a);
-    double siBt = sas_sinx_x(0.5*q*tb*cos_val_b);
-    double siCt = sas_sinx_x(0.5*q*tc*cos_val_c);
+    double siA = sas_sinx_x(0.5*q*length_a*xhat);
+    double siB = sas_sinx_x(0.5*q*length_b*yhat);
+    double siC = sas_sinx_x(0.5*q*length_c*zhat);
+    double siAt = sas_sinx_x(0.5*q*ta*xhat);
+    double siBt = sas_sinx_x(0.5*q*tb*yhat);
+    double siCt = sas_sinx_x(0.5*q*tc*zhat);
     
 
Index: sasmodels/models/core_shell_parallelepiped.py
===================================================================
--- sasmodels/models/core_shell_parallelepiped.py	(revision cb0dc22aeaff717fb8130fd48c3134fa935d7ceb)
+++ sasmodels/models/core_shell_parallelepiped.py	(revision cb0dc22aeaff717fb8130fd48c3134fa935d7ceb)
@@ -48,4 +48,11 @@
 amplitudes of the core and shell, in the same manner as a core-shell model.
 
+.. math::
+
+    F_{a}(Q,\alpha,\beta)=
+    \left[\frac{\sin(Q(L_A+2t_A)/2\sin\alpha \sin\beta)}{Q(L_A+2t_A)/2\sin\alpha\sin\beta}
+    - \frac{\sin(QL_A/2\sin\alpha \sin\beta)}{QL_A/2\sin\alpha \sin\beta} \right]
+    \left[\frac{\sin(QL_B/2\sin\alpha \sin\beta)}{QL_B/2\sin\alpha \sin\beta} \right]
+    \left[\frac{\sin(QL_C/2\sin\alpha \sin\beta)}{QL_C/2\sin\alpha \sin\beta} \right]
 
 .. note::
@@ -57,5 +64,5 @@
 
 FITTING NOTES
-If the scale is set equal to the particle volume fraction, |phi|, the returned
+If the scale is set equal to the particle volume fraction, $\phi$, the returned
 value is the scattered intensity per unit volume, $I(q) = \phi P(q)$.
 However, **no interparticle interference effects are included in this
Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision 130d4c75420726aaf82e363091daa51b41ec0468)
+++ sasmodels/models/ellipsoid.c	(revision 3b571ae4950c5a99059d2ce8c74e80d0eb908a8f)
@@ -3,26 +3,4 @@
 double Iqxy(double qx, double qy, double sld, double sld_solvent,
     double radius_polar, double radius_equatorial, double theta, double phi);
-
-static double
-_ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double cos_alpha)
-{
-    double ratio = radius_polar/radius_equatorial;
-    // Using ratio v = Rp/Re, we can expand the following to match the
-    // form given in Guinier (1955)
-    //     r = Re * sqrt(1 + cos^2(T) (v^2 - 1))
-    //       = Re * sqrt( (1 - cos^2(T)) + v^2 cos^2(T) )
-    //       = Re * sqrt( sin^2(T) + v^2 cos^2(T) )
-    //       = sqrt( Re^2 sin^2(T) + Rp^2 cos^2(T) )
-    //
-    // Instead of using pythagoras we could pass in sin and cos; this may be
-    // slightly better for 2D which has already computed it, but it introduces
-    // an extra sqrt and square for 1-D not required by the current form, so
-    // leave it as is.
-    const double r = radius_equatorial
-                     * sqrt(1.0 + cos_alpha*cos_alpha*(ratio*ratio - 1.0));
-    const double f = sas_3j1x_x(q*r);
-
-    return f*f;
-}
 
 double form_volume(double radius_polar, double radius_equatorial)
@@ -37,12 +15,23 @@
     double radius_equatorial)
 {
+    // Using ratio v = Rp/Re, we can implement the form given in Guinier (1955)
+    //     i(h) = int_0^pi/2 Phi^2(h a sqrt(cos^2 + v^2 sin^2) cos dT
+    //          = int_0^pi/2 Phi^2(h a sqrt((1-sin^2) + v^2 sin^2) cos dT
+    //          = int_0^pi/2 Phi^2(h a sqrt(1 + sin^2(v^2-1)) cos dT
+    // u-substitution of
+    //     u = sin, du = cos dT
+    //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du
+    const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0;
+
     // translate a point in [-1,1] to a point in [0, 1]
+    // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2;
     const double zm = 0.5;
     const double zb = 0.5;
     double total = 0.0;
     for (int i=0;i<76;i++) {
-        //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2;
-        const double cos_alpha = Gauss76Z[i]*zm + zb;
-        total += Gauss76Wt[i] * _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha);
+        const double u = Gauss76Z[i]*zm + zb;
+        const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one);
+        const double f = sas_3j1x_x(q*r);
+        total += Gauss76Wt[i] * f * f;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -62,8 +51,10 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
-    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha);
+    const double r = sqrt(square(radius_equatorial*sin_alpha)
+                          + square(radius_polar*cos_alpha));
+    const double f = sas_3j1x_x(q*r);
     const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial);
 
-    return 1.0e-4 * form * s * s;
+    return 1.0e-4 * square(f * s);
 }
 
Index: sasmodels/models/ellipsoid.py
===================================================================
--- sasmodels/models/ellipsoid.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
+++ sasmodels/models/ellipsoid.py	(revision 3b571ae4950c5a99059d2ce8c74e80d0eb908a8f)
@@ -18,28 +18,38 @@
 .. math::
 
-    F(q,\alpha) = \frac{3 \Delta \rho V (\sin[qr(R_p,R_e,\alpha)]
-                - \cos[qr(R_p,R_e,\alpha)])}
-                {[qr(R_p,R_e,\alpha)]^3}
+    F(q,\alpha) = \Delta \rho V \frac{3(\sin qr  - qr \cos qr)}{(qr)^3}
 
-and
+for
 
 .. math::
 
-    r(R_p,R_e,\alpha) = \left[ R_e^2 \sin^2 \alpha
-        + R_p^2 \cos^2 \alpha \right]^{1/2}
+    r = \left[ R_e^2 \sin^2 \alpha + R_p^2 \cos^2 \alpha \right]^{1/2}
 
 
 $\alpha$ is the angle between the axis of the ellipsoid and $\vec q$,
-$V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid , $R_p$ is the polar radius along the
-rotational axis of the ellipsoid, $R_e$ is the equatorial radius perpendicular
-to the rotational axis of the ellipsoid and $\Delta \rho$ (contrast) is the
-scattering length density difference between the scatterer and the solvent.
+$V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid, $R_p$ is the polar
+radius along the rotational axis of the ellipsoid, $R_e$ is the equatorial
+radius perpendicular to the rotational axis of the ellipsoid and
+$\Delta \rho$ (contrast) is the scattering length density difference between
+the scatterer and the solvent.
 
-For randomly oriented particles:
+For randomly oriented particles use the orientational average,
 
 .. math::
 
-   F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}
+   \langle F^2(q) \rangle = \int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)\,d\alpha}
 
+
+computed via substitution of $u=\sin(\alpha)$, $du=\cos(\alpha)\,d\alpha$ as
+
+.. math::
+
+    \langle F^2(q) \rangle = \int_0^1{F^2(q, u)\,du}
+
+with
+
+.. math::
+
+    r = R_e \left[ 1 + u^2\left(R_p^2/R_e^2 - 1\right)\right]^{1/2}
 
 To provide easy access to the orientation of the ellipsoid, we define
@@ -48,5 +58,6 @@
 :ref:`cylinder orientation figure <cylinder-angle-definition>`.
 For the ellipsoid, $\theta$ is the angle between the rotational axis
-and the $z$ -axis.
+and the $z$ -axis in the $xz$ plane followed by a rotation by $\phi$
+in the $xy$ plane.
 
 NB: The 2nd virial coefficient of the solid ellipsoid is calculated based
@@ -90,4 +101,8 @@
 than 500.
 
+Model was also tested against the triaxial ellipsoid model with equal major
+and minor equatorial radii.  It is also consistent with the cyclinder model
+with polar radius equal to length and equatorial radius equal to radius.
+
 References
 ----------
@@ -96,4 +111,11 @@
 *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*,
 Plenum Press, New York, 1987.
+
+Authorship and Verification
+----------------------------
+
+* **Author:** NIST IGOR/DANSE **Date:** pre 2010
+* **Converted to sasmodels by:** Helen Park **Date:** July 9, 2014
+* **Last Modified by:** Paul Kienzle **Date:** March 22, 2017
 """
 
Index: sasmodels/models/elliptical_cylinder.c
===================================================================
--- sasmodels/models/elliptical_cylinder.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
+++ sasmodels/models/elliptical_cylinder.c	(revision 61104c829b100869b42c1afd75f30ed30a934cd1)
@@ -67,12 +67,12 @@
      double theta, double phi, double psi)
 {
-    double q, cos_val, cos_mu, cos_nu;
-    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val, cos_mu, cos_nu);
+    double q, xhat, yhat, zhat;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
 
     // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2)
     // Given:    radius_major = r_ratio * radius_minor
-    const double r = radius_minor*sqrt(square(r_ratio*cos_nu) + cos_mu*cos_mu);
+    const double r = radius_minor*sqrt(square(r_ratio*xhat) + square(yhat));
     const double be = sas_2J1x_x(q*r);
-    const double si = sas_sinx_x(q*0.5*length*cos_val);
+    const double si = sas_sinx_x(q*zhat*0.5*length);
     const double Aq = be * si;
     const double delrho = sld - solvent_sld;
Index: sasmodels/models/fcc_paracrystal.c
===================================================================
--- sasmodels/models/fcc_paracrystal.c	(revision 4962519cc1910422a467e4f643dba223623e0e09)
+++ sasmodels/models/fcc_paracrystal.c	(revision 50beefeac18f2f952bf0863526ec9bbd3009e29c)
@@ -90,10 +90,10 @@
     double theta, double phi, double psi)
 {
-    double q, cos_a1, cos_a2, cos_a3;
-    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1);
+    double q, zhat, yhat, xhat;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
 
-    const double a1 = cos_a2 + cos_a3;
-    const double a2 = cos_a3 + cos_a1;
-    const double a3 = cos_a2 + cos_a1;
+    const double a1 = yhat + xhat;
+    const double a2 = xhat + zhat;
+    const double a3 = yhat + zhat;
     const double qd = 0.5*q*dnn;
     const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3);
Index: sasmodels/models/parallelepiped.c
===================================================================
--- sasmodels/models/parallelepiped.c	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
+++ sasmodels/models/parallelepiped.c	(revision d6050806338b4bcae5a997acb08f4b1e7e14dd36)
@@ -67,10 +67,10 @@
     double psi)
 {
-    double q, cos_val_a, cos_val_b, cos_val_c;
-    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a);
+    double q, xhat, yhat, zhat;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
 
-    const double siA = sas_sinx_x(0.5*q*length_a*cos_val_a);
-    const double siB = sas_sinx_x(0.5*q*length_b*cos_val_b);
-    const double siC = sas_sinx_x(0.5*q*length_c*cos_val_c);
+    const double siA = sas_sinx_x(0.5*length_a*q*xhat);
+    const double siB = sas_sinx_x(0.5*length_b*q*yhat);
+    const double siC = sas_sinx_x(0.5*length_c*q*zhat);
     const double V = form_volume(length_a, length_b, length_c);
     const double drho = (sld - solvent_sld);
Index: sasmodels/models/sc_paracrystal.c
===================================================================
--- sasmodels/models/sc_paracrystal.c	(revision 4962519cc1910422a467e4f643dba223623e0e09)
+++ sasmodels/models/sc_paracrystal.c	(revision 50beefeac18f2f952bf0863526ec9bbd3009e29c)
@@ -111,6 +111,6 @@
           double psi)
 {
-    double q, cos_a1, cos_a2, cos_a3;
-    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1);
+    double q, zhat, yhat, xhat;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
 
     const double qd = q*dnn;
@@ -118,7 +118,7 @@
     const double tanh_qd = tanh(arg);
     const double cosh_qd = cosh(arg);
-    const double Zq = tanh_qd/(1. - cos(qd*cos_a1)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*cos_a2)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*cos_a3)/cosh_qd);
+    const double Zq = tanh_qd/(1. - cos(qd*zhat)/cosh_qd)
+                    * tanh_qd/(1. - cos(qd*yhat)/cosh_qd)
+                    * tanh_qd/(1. - cos(qd*xhat)/cosh_qd);
 
     const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq;
Index: sasmodels/models/triaxial_ellipsoid.c
===================================================================
--- sasmodels/models/triaxial_ellipsoid.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision 68dd6a95741269738dcbfb051d93fe9f26b89c43)
@@ -20,28 +20,29 @@
     double radius_polar)
 {
-    double sn, cn;
-    // translate a point in [-1,1] to a point in [0, 1]
-    const double zm = 0.5;
-    const double zb = 0.5;
+    const double pa = square(radius_equat_minor/radius_equat_major) - 1.0;
+    const double pc = square(radius_polar/radius_equat_major) - 1.0;
+    // translate a point in [-1,1] to a point in [0, pi/2]
+    const double zm = M_PI_4;
+    const double zb = M_PI_4;
     double outer = 0.0;
     for (int i=0;i<76;i++) {
-        //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2;
-        const double x = 0.5*(Gauss76Z[i] + 1.0);
-        SINCOS(M_PI_2*x, sn, cn);
-        const double acosx2 = radius_equat_minor*radius_equat_minor*cn*cn;
-        const double bsinx2 = radius_equat_major*radius_equat_major*sn*sn;
-        const double c2 = radius_polar*radius_polar;
+        //const double u = Gauss76Z[i]*(upper-lower)/2 + (upper + lower)/2;
+        const double phi = Gauss76Z[i]*zm + zb;
+        const double pa_sinsq_phi = pa*square(sin(phi));
 
         double inner = 0.0;
+        const double um = 0.5;
+        const double ub = 0.5;
         for (int j=0;j<76;j++) {
-            const double ysq = square(Gauss76Z[j]*zm + zb);
-            const double t = q*sqrt(acosx2 + bsinx2*(1.0-ysq) + c2*ysq);
-            const double fq = sas_3j1x_x(t);
-            inner += Gauss76Wt[j] * fq * fq ;
+            // translate a point in [-1,1] to a point in [0, 1]
+            const double usq = square(Gauss76Z[j]*um + ub);
+            const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0-usq) + 1.0 + pc*usq);
+            const double fq = sas_3j1x_x(q*r);
+            inner += Gauss76Wt[j] * fq * fq;
         }
-        outer += Gauss76Wt[i] * 0.5 * inner;
+        outer += Gauss76Wt[i] * inner;  // correcting for dx later
     }
-    // translate dx in [-1,1] to dx in [lower,upper]
-    const double fqsq = outer*zm;
+    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi
+    const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI);
     const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar);
     return 1.0e-4 * s * s * fqsq;
@@ -58,11 +59,11 @@
     double psi)
 {
-    double q, calpha, cmu, cnu;
-    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, calpha, cmu, cnu);
+    double q, xhat, yhat, zhat;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
 
-    const double t = q*sqrt(radius_equat_minor*radius_equat_minor*cnu*cnu
-                          + radius_equat_major*radius_equat_major*cmu*cmu
-                          + radius_polar*radius_polar*calpha*calpha);
-    const double fq = sas_3j1x_x(t);
+    const double r = sqrt(square(radius_equat_minor*xhat)
+                          + square(radius_equat_major*yhat)
+                          + square(radius_polar*zhat));
+    const double fq = sas_3j1x_x(q*r);
     const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar);
 
Index: sasmodels/models/triaxial_ellipsoid.py
===================================================================
--- sasmodels/models/triaxial_ellipsoid.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
+++ sasmodels/models/triaxial_ellipsoid.py	(revision 28d30676424390bfed1a13e70bcba3332cd86c75)
@@ -2,28 +2,23 @@
 # Note: model title and parameter table are inserted automatically
 r"""
-All three axes are of different lengths with $R_a \leq R_b \leq R_c$
-**Users should maintain this inequality for all calculations**.
+Definition
+----------
+
+.. figure:: img/triaxial_ellipsoid_geometry.jpg
+
+    Ellipsoid with $R_a$ as *radius_equat_minor*, $R_b$ as *radius_equat_major*
+    and $R_c$ as *radius_polar*.
+
+Given an ellipsoid
 
 .. math::
 
-    P(q) = \text{scale} V \left< F^2(q) \right> + \text{background}
+    \frac{X^2}{R_a^2} + \frac{Y^2}{R_b^2} + \frac{Z^2}{R_c^2} = 1
 
-where the volume $V = 4/3 \pi R_a R_b R_c$, and the averaging
-$\left<\ldots\right>$ is applied over all orientations for 1D.
-
-.. figure:: img/triaxial_ellipsoid_geometry.jpg
-
-    Ellipsoid schematic.
-
-Definition
-----------
-
-The form factor calculated is
+the scattering is defined by the average over all orientations $\Omega$,
 
 .. math::
 
-    P(q) = \frac{\text{scale}}{V}\int_0^1\int_0^1
-        \Phi^2(qR_a^2\cos^2( \pi x/2) + qR_b^2\sin^2(\pi y/2)(1-y^2) + R_c^2y^2)
-        dx dy
+    P(q) = \text{scale}\frac{V}{4 \pi}\int_\Omega \Phi^2(qr) d\Omega + \text{background}
 
 where
@@ -31,5 +26,40 @@
 .. math::
 
-    \Phi(u) = 3 u^{-3} (\sin u - u \cos u)
+    \Phi(qr) &= 3 j_1(qr)/qr = 3 (\sin qr - qr \cos qr)/(qr)^3 \\
+    r^2 &= R_a^2e^2 + R_b^2f^2 + R_c^2g^2 \\
+    V &= \tfrac{4}{3} \pi R_a R_b R_c
+
+The $e$, $f$ and $g$ terms are the projections of the orientation vector on $X$,
+$Y$ and $Z$ respectively.  Keeping the orientation fixed at the canonical
+axes, we can integrate over the incident direction using polar angle
+$-\pi/2 \le \gamma \le \pi/2$ and equatorial angle $0 \le \phi \le 2\pi$
+(as defined in ref [1]),
+
+ .. math::
+
+     \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$.
+A little algebra yields
+
+.. math::
+
+    r^2 = b^2(p_a \sin^2 \phi \cos^2 \gamma + 1 + p_c \sin^2 \gamma)
+
+for
+
+.. math::
+
+    p_a = \frac{a^2}{b^2} - 1 \text{ and } p_c = \frac{c^2}{b^2} - 1
+
+Due to symmetry, the ranges can be restricted to a single quadrant
+$0 \le \gamma \le \pi/2$ and $0 \le \phi \le \pi/2$, scaling the resulting
+integral by 8. The computation is done using the substitution $u = \sin\gamma$,
+$du = \cos\gamma\,d\gamma$, giving
+
+.. math::
+
+    \langle\Phi^2\rangle &= 8 \int_0^{\pi/2} \int_0^1 \Phi^2(qr) du d\phi \\
+    r^2 &= b^2(p_a \sin^2(\phi)(1 - u^2) + 1 + p_c u^2)
 
 To provide easy access to the orientation of the triaxial ellipsoid,
@@ -69,6 +99,8 @@
 ----------
 
-L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray
-and Neutron Scattering*, Plenum, New York, 1987.
+[1] Finnigan, J.A., Jacobs, D.J., 1971.
+*Light scattering by ellipsoidal particles in solution*,
+J. Phys. D: Appl. Phys. 4, 72-77. doi:10.1088/0022-3727/4/1/310
+
 """
 
@@ -91,9 +123,9 @@
                "Solvent scattering length density"],
               ["radius_equat_minor", "Ang", 20, [0, inf], "volume",
-               "Minor equatorial radius"],
+               "Minor equatorial radius, Ra"],
               ["radius_equat_major", "Ang", 400, [0, inf], "volume",
-               "Major equatorial radius"],
+               "Major equatorial radius, Rb"],
               ["radius_polar", "Ang", 10, [0, inf], "volume",
-               "Polar radius"],
+               "Polar radius, Rc"],
               ["theta", "degrees", 60, [-inf, inf], "orientation",
                "In plane angle"],
