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/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: 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/barbell.py
===================================================================
--- sasmodels/models/barbell.py	(revision fcb33e4672814b4a8fe1ee865d430090b5fc03e6)
+++ sasmodels/models/barbell.py	(revision 9802ab3945a154b871264715dcb463d4ec644afa)
@@ -68,5 +68,5 @@
 The 2D scattering intensity is calculated similar to the 2D cylinder model.
 
-.. figure:: img/cylinder_angle_definition.jpg
+.. figure:: img/cylinder_angle_definition.png
 
     Definition of the angles for oriented 2D barbells.
@@ -87,5 +87,5 @@
 * **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017
 """
-from numpy import inf
+from numpy import inf, sin, cos, pi
 
 name = "barbell"
@@ -108,6 +108,6 @@
               ["radius",      "Ang",         20, [0, inf],    "volume",      "Cylindrical bar radius"],
               ["length",      "Ang",        400, [0, inf],    "volume",      "Cylinder bar length"],
-              ["theta",       "degrees",     60, [-inf, inf], "orientation", "In plane angle"],
-              ["phi",         "degrees",     60, [-inf, inf], "orientation", "Out of plane angle"],
+              ["theta",       "degrees",     60, [-360, 360], "orientation", "Barbell axis to beam angle"],
+              ["phi",         "degrees",     60, [-360, 360], "orientation", "Rotation about beam"],
              ]
 # pylint: enable=bad-whitespace, line-too-long
@@ -125,2 +125,10 @@
             phi_pd=15, phi_pd_n=0,
            )
+q = 0.1
+# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
+qx = q*cos(pi/6.0)
+qy = q*sin(pi/6.0)
+tests = [[{}, 0.075, 25.5691260532],
+        [{'theta':80., 'phi':10.}, (qx, qy), 3.04233067789],
+        ]
+del qx, qy  # not necessary to delete, but cleaner
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/bcc_paracrystal.py
===================================================================
--- sasmodels/models/bcc_paracrystal.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
+++ sasmodels/models/bcc_paracrystal.py	(revision 69e1afc57270f67b2a0d7157a6d263fc2db98cc0)
@@ -79,7 +79,8 @@
 be accurate.
 
-.. figure:: img/bcc_angle_definition.png
+.. figure:: img/parallelepiped_angle_definition.png
 
-    Orientation of the crystal with respect to the scattering plane.
+    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).
 
 References
@@ -99,5 +100,5 @@
 """
 
-from numpy import inf
+from numpy import inf, pi
 
 name = "bcc_paracrystal"
@@ -122,7 +123,7 @@
               ["sld",         "1e-6/Ang^2",  4,    [-inf, inf], "sld",         "Particle scattering length density"],
               ["sld_solvent", "1e-6/Ang^2",  1,    [-inf, inf], "sld",         "Solvent scattering length density"],
-              ["theta",       "degrees",    60,    [-inf, inf], "orientation", "In plane angle"],
-              ["phi",         "degrees",    60,    [-inf, inf], "orientation", "Out of plane angle"],
-              ["psi",         "degrees",    60,    [-inf, inf], "orientation", "Out of plane angle"]
+              ["theta",       "degrees",    60,    [-360, 360], "orientation", "c axis to beam angle"],
+              ["phi",         "degrees",    60,    [-360, 360], "orientation", "rotation about beam"],
+              ["psi",         "degrees",    60,    [-360, 360], "orientation", "rotation about c axis"]
              ]
 # pylint: enable=bad-whitespace, line-too-long
@@ -141,2 +142,11 @@
     psi_pd=15, psi_pd_n=0,
     )
+# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
+# add 2d test later
+q =4.*pi/220.
+tests = [
+    [{ },
+     [0.001, q, 0.215268], [1.46601394721, 2.85851284174, 0.00866710287078]],
+    [{'theta':20.0,'phi':30,'psi':40.0},(-0.017,0.035),2082.20264399 ],
+    [{'theta':20.0,'phi':30,'psi':40.0},(-0.081,0.011),0.436323144781 ]
+    ]
Index: sasmodels/models/capped_cylinder.py
===================================================================
--- sasmodels/models/capped_cylinder.py	(revision fcb33e4672814b4a8fe1ee865d430090b5fc03e6)
+++ sasmodels/models/capped_cylinder.py	(revision 9802ab3945a154b871264715dcb463d4ec644afa)
@@ -71,5 +71,5 @@
 The 2D scattering intensity is calculated similar to the 2D cylinder model.
 
-.. figure:: img/cylinder_angle_definition.jpg
+.. figure:: img/cylinder_angle_definition.png
 
     Definition of the angles for oriented 2D cylinders.
@@ -91,5 +91,5 @@
 
 """
-from numpy import inf
+from numpy import inf, sin, cos, pi
 
 name = "capped_cylinder"
@@ -129,6 +129,6 @@
               ["radius_cap", "Ang",     20, [0, inf],    "volume", "Cap radius"],
               ["length",     "Ang",    400, [0, inf],    "volume", "Cylinder length"],
-              ["theta",      "degrees", 60, [-inf, inf], "orientation", "inclination angle"],
-              ["phi",        "degrees", 60, [-inf, inf], "orientation", "deflection angle"],
+              ["theta",      "degrees", 60, [-360, 360], "orientation", "cylinder axis to beam angle"],
+              ["phi",        "degrees", 60, [-360, 360], "orientation", "rotation about beam"],
              ]
 # pylint: enable=bad-whitespace, line-too-long
@@ -145,2 +145,10 @@
             theta_pd=15, theta_pd_n=45,
             phi_pd=15, phi_pd_n=1)
+q = 0.1
+# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
+qx = q*cos(pi/6.0)
+qy = q*sin(pi/6.0)
+tests = [[{}, 0.075, 26.0698570695],
+        [{'theta':80., 'phi':10.}, (qx, qy), 0.561811990502],
+        ]
+del qx, qy  # not necessary to delete, but cleaner
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.py
===================================================================
--- sasmodels/models/core_shell_bicelle.py	(revision 8afefaee3d363dde220aa55fe83181a73b621f98)
+++ sasmodels/models/core_shell_bicelle.py	(revision 9802ab3945a154b871264715dcb463d4ec644afa)
@@ -47,5 +47,5 @@
 .. math::
 
-        \begin{align}    
+    \begin{align}    
     F(Q,\alpha) = &\bigg[ 
     (\rho_c - \rho_f) V_c \frac{2J_1(QRsin \alpha)}{QRsin\alpha}\frac{sin(QLcos\alpha/2)}{Q(L/2)cos\alpha} \\
@@ -63,11 +63,14 @@
 cylinders is then given by integrating over all possible $\theta$ and $\phi$.
 
-The *theta* and *phi* parameters are not used for the 1D output.
+For oriented bicelles the *theta*, and *phi* orientation parameters will appear when fitting 2D data, 
+see the :ref:`cylinder` model for further information.
 Our implementation of the scattering kernel and the 1D scattering intensity
 use the c-library from NIST.
 
-.. figure:: img/cylinder_angle_definition.jpg
+.. figure:: img/cylinder_angle_definition.png
 
-    Definition of the angles for the oriented core shell bicelle tmodel.
+    Definition of the angles for the oriented core shell bicelle model,
+    note that the cylinder axis of the bicelle starts along the beam direction
+    when $\theta  = \phi = 0$.
 
 
@@ -88,5 +91,5 @@
 """
 
-from numpy import inf, sin, cos
+from numpy import inf, sin, cos, pi
 
 name = "core_shell_bicelle"
@@ -135,6 +138,6 @@
     ["sld_rim",        "1e-6/Ang^2", 4, [-inf, inf], "sld",         "Cylinder rim scattering length density"],
     ["sld_solvent",    "1e-6/Ang^2", 1, [-inf, inf], "sld",         "Solvent scattering length density"],
-    ["theta",          "degrees",   90, [-inf, inf], "orientation", "In plane angle"],
-    ["phi",            "degrees",    0, [-inf, inf], "orientation", "Out of plane angle"],
+    ["theta",          "degrees",   90, [-360, 360], "orientation", "cylinder axis to beam angle"],
+    ["phi",            "degrees",    0, [-360, 360], "orientation", "rotation about beam"]
     ]
 
@@ -155,4 +158,11 @@
             theta=90,
             phi=0)
+q = 0.1
+# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
+qx = q*cos(pi/6.0)
+qy = q*sin(pi/6.0)
+tests = [[{}, 0.05, 7.4883545957],
+        [{'theta':80., 'phi':10.}, (qx, qy), 2.81048892474 ]
+        ]
+del qx, qy  # not necessary to delete, but cleaner
 
-#qx, qy = 0.4 * cos(pi/2.0), 0.5 * sin(0)
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 dedcf34f547ba57a5e8ca3c5e9892ede11d1fcaf)
@@ -1,50 +1,28 @@
-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 qq,
-                   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;
      // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle
-     // tested against limiting cases of cylinder, elliptical_cylinder and core_shell_bicelle
+     // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle
      //    const double uplim = M_PI_4;
     const double halfheight = 0.5*length;
@@ -55,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 r2A = 0.5*(square(r_major) + square(r_minor));
+    const double r2B = 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,6 +52,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+thick_face)*cos_alpha;
         si1 = sas_sinx_x(sinarg1);
         si2 = sas_sinx_x(sinarg2);
@@ -82,7 +60,7 @@
             //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
             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;
+            const double rr = sqrt(r2A - r2B*cos(beta));
+            double besarg1 = q*rr*sin_alpha;
+            double besarg2 = q*(rr+thick_rim)*sin_alpha;
             be1 = sas_2J1x_x(besarg1);
             be2 = sas_2J1x_x(besarg2);
@@ -98,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,
@@ -114,25 +92,24 @@
 {
        // 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;
     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*cos_nu)^2 + (radius_minor*cos_mu)^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 );
+    // 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 + 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/core_shell_bicelle_elliptical.py
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical.py	(revision 8afefaee3d363dde220aa55fe83181a73b621f98)
+++ sasmodels/models/core_shell_bicelle_elliptical.py	(revision 9802ab3945a154b871264715dcb463d4ec644afa)
@@ -76,12 +76,12 @@
 bicelles is then given by integrating over all possible $\alpha$ and $\psi$.
 
-For oriented bicellles the *theta*, *phi* and *psi* orientation parameters only appear when fitting 2D data, 
+For oriented bicelles the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data, 
 see the :ref:`elliptical-cylinder` model for further information.
 
 
-.. figure:: img/elliptical_cylinder_angle_definition.jpg
+.. figure:: img/elliptical_cylinder_angle_definition.png
 
-    Definition of the angles for the oriented core_shell_bicelle_elliptical model.
-    Note that *theta* and *phi* are currently defined differently to those for the core_shell_bicelle model.
+    Definition of the angles for the oriented core_shell_bicelle_elliptical particles.   
+
 
 
@@ -99,5 +99,5 @@
 """
 
-from numpy import inf, sin, cos
+from numpy import inf, sin, cos, pi
 
 name = "core_shell_bicelle_elliptical"
@@ -119,14 +119,14 @@
     ["radius",         "Ang",       30, [0, inf],    "volume",      "Cylinder core radius"],
     ["x_core",        "None",       3,  [0, inf],    "volume",      "axial ratio of core, X = r_polar/r_equatorial"],
-    ["thick_rim",  "Ang",        8, [0, inf],    "volume",      "Rim shell thickness"],
-    ["thick_face", "Ang",       14, [0, inf],    "volume",      "Cylinder face thickness"],
-    ["length",         "Ang",      50, [0, inf],    "volume",      "Cylinder length"],
+    ["thick_rim",  "Ang",            8, [0, inf],    "volume",      "Rim shell thickness"],
+    ["thick_face", "Ang",           14, [0, inf],    "volume",      "Cylinder face thickness"],
+    ["length",         "Ang",       50, [0, inf],    "volume",      "Cylinder length"],
     ["sld_core",       "1e-6/Ang^2", 4, [-inf, inf], "sld",         "Cylinder core scattering length density"],
     ["sld_face",       "1e-6/Ang^2", 7, [-inf, inf], "sld",         "Cylinder face scattering length density"],
     ["sld_rim",        "1e-6/Ang^2", 1, [-inf, inf], "sld",         "Cylinder rim scattering length density"],
     ["sld_solvent",    "1e-6/Ang^2", 6, [-inf, inf], "sld",         "Solvent scattering length density"],
-    ["theta",          "degrees",   90, [-360, 360], "orientation", "In plane angle"],
-    ["phi",            "degrees",    0, [-360, 360], "orientation", "Out of plane angle"],
-    ["psi",            "degrees",    0, [-360, 360], "orientation", "Major axis angle relative to Q"],
+    ["theta",       "degrees",    90.0, [-360, 360], "orientation", "cylinder axis to beam angle"],
+    ["phi",         "degrees",    0,    [-360, 360], "orientation", "rotation about beam"],
+    ["psi",         "degrees",    0,    [-360, 360], "orientation", "rotation about cylinder axis"]
     ]
 
@@ -150,5 +150,8 @@
             psi=0)
 
-#qx, qy = 0.4 * cos(pi/2.0), 0.5 * sin(0)
+q = 0.1
+# april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct!
+qx = q*cos(pi/6.0)
+qy = q*sin(pi/6.0)
 
 tests = [
@@ -159,3 +162,6 @@
     'sld_core':4.0, 'sld_face':7.0, 'sld_rim':1.0, 'sld_solvent':6.0, 'background':0.0},
     0.015, 286.540286],
-]
+#    [{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ],
+        ]
+
+del qx, qy  # not necessary to delete, but cleaner
Index: sasmodels/models/core_shell_cylinder.py
===================================================================
--- sasmodels/models/core_shell_cylinder.py	(revision 8e68ea06a7246c4bd95ce82c755719ced6fa5a60)
+++ sasmodels/models/core_shell_cylinder.py	(revision 9b79f2966c4083a48b4894b6bf0587bef7a5be09)
@@ -73,5 +73,5 @@
 """
 
-from numpy import pi, inf
+from numpy import pi, inf, sin, cos
 
 name = "core_shell_cylinder"
@@ -117,8 +117,8 @@
               ["length", "Ang", 400, [0, inf], "volume",
                "Cylinder length"],
-              ["theta", "degrees", 60, [-inf, inf], "orientation",
-               "In plane angle"],
-              ["phi", "degrees", 60, [-inf, inf], "orientation",
-               "Out of plane angle"],
+              ["theta", "degrees", 60, [-360, 360], "orientation",
+               "cylinder axis to beam angle"],
+              ["phi", "degrees",   60, [-360, 360], "orientation",
+               "rotation about beam"],
              ]
 
@@ -151,3 +151,10 @@
             theta_pd=15, theta_pd_n=45,
             phi_pd=15, phi_pd_n=1)
-
+q = 0.1
+# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
+qx = q*cos(pi/6.0)
+qy = q*sin(pi/6.0)
+tests = [[{}, 0.075, 10.8552692237],
+        [{}, (qx, qy), 0.444618752741 ],
+        ]
+del qx, qy  # not necessary to delete, but cleaner
Index: sasmodels/models/core_shell_ellipsoid.py
===================================================================
--- sasmodels/models/core_shell_ellipsoid.py	(revision 8e68ea06a7246c4bd95ce82c755719ced6fa5a60)
+++ sasmodels/models/core_shell_ellipsoid.py	(revision 9802ab3945a154b871264715dcb463d4ec644afa)
@@ -77,4 +77,6 @@
    F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}
 
+For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data, 
+see the :ref:`elliptical-cylinder` model for further information.
 
 References
@@ -132,6 +134,6 @@
     ["sld_shell",     "1e-6/Ang^2", 1,   [-inf, inf], "sld",         "Shell scattering length density"],
     ["sld_solvent",   "1e-6/Ang^2", 6.3, [-inf, inf], "sld",         "Solvent scattering length density"],
-    ["theta",         "degrees",    0,   [-inf, inf], "orientation", "Oblate orientation wrt incoming beam"],
-    ["phi",           "degrees",    0,   [-inf, inf], "orientation", "Oblate orientation in the plane of the detector"],
+    ["theta",         "degrees",    0,   [-360, 360], "orientation", "elipsoid axis to beam angle"],
+    ["phi",           "degrees",    0,   [-360, 360], "orientation", "rotation about beam"],
     ]
 # pylint: enable=bad-whitespace, line-too-long
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 9b79f2966c4083a48b4894b6bf0587bef7a5be09)
@@ -6,10 +6,10 @@
 The thickness and the scattering length density of the shell or 
 "rim" can be different on each (pair) of faces. However at this time
-the model does **NOT** actually calculate a c face rim despite the presence of
-the parameter.
+the 1D calculation does **NOT** actually calculate a c face rim despite the presence of
+the parameter. Some other aspects of the 1D calculation may be wrong.
 
 .. note::
-   This model was originally ported from NIST IGOR macros. However,t is not
-   yet fully understood by the SasView developers and is currently review.
+   This model was originally ported from NIST IGOR macros. However, it is not
+   yet fully understood by the SasView developers and is currently under review.
 
 The form factor is normalized by the particle volume $V$ such that
@@ -48,14 +48,22 @@
 amplitudes of the core and shell, in the same manner as a core-shell model.
 
+.. math::
+
+    F_{a}(Q,\alpha,\beta)=
+    \left[\frac{\sin(\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha \sin\beta)}{\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha\sin\beta}
+    - \frac{\sin(\tfrac{1}{2}QL_A\sin\alpha \sin\beta)}{\tfrac{1}{2}QL_A\sin\alpha \sin\beta} \right]
+    \left[\frac{\sin(\tfrac{1}{2}QL_B\sin\alpha \sin\beta)}{\tfrac{1}{2}QL_B\sin\alpha \sin\beta} \right]
+    \left[\frac{\sin(\tfrac{1}{2}QL_C\sin\alpha \sin\beta)}{\tfrac{1}{2}QL_C\sin\alpha \sin\beta} \right]
 
 .. note::
 
+    Why does t_B not appear in the above equation?
     For the calculation of the form factor to be valid, the sides of the solid
-    MUST be chosen such that** $A < B < C$.
+    MUST (perhaps not any more?) be chosen such that** $A < B < C$.
     If this inequality is not satisfied, the model will not report an error,
     but the calculation will not be correct and thus the result wrong.
 
 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
@@ -73,6 +81,7 @@
 NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated
 based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$
-and length $(C+2t_C)$ values, and used as the effective radius
-for $S(Q)$ when $P(Q) * S(Q)$ is applied.
+and length $(C+2t_C)$ values, after appropriately
+sorting the three dimensions to give an oblate or prolate particle, to give an 
+effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied.
 
 To provide easy access to the orientation of the parallelepiped, we define the
@@ -83,9 +92,9 @@
 *x*-axis of the detector.
 
-.. figure:: img/parallelepiped_angle_definition.jpg
+.. figure:: img/parallelepiped_angle_definition.png
 
     Definition of the angles for oriented core-shell parallelepipeds.
 
-.. figure:: img/parallelepiped_angle_projection.jpg
+.. figure:: img/parallelepiped_angle_projection.png
 
     Examples of the angles for oriented core-shell parallelepipeds against the
@@ -112,5 +121,5 @@
 
 import numpy as np
-from numpy import pi, inf, sqrt
+from numpy import pi, inf, sqrt, cos, sin
 
 name = "core_shell_parallelepiped"
@@ -144,10 +153,10 @@
               ["thick_rim_c", "Ang", 10, [0, inf], "volume",
                "Thickness of C rim"],
-              ["theta", "degrees", 0, [-inf, inf], "orientation",
-               "In plane angle"],
-              ["phi", "degrees", 0, [-inf, inf], "orientation",
-               "Out of plane angle"],
-              ["psi", "degrees", 0, [-inf, inf], "orientation",
-               "Rotation angle around its own c axis against q plane"],
+              ["theta", "degrees", 0, [-360, 360], "orientation",
+               "c axis to beam angle"],
+              ["phi", "degrees", 0, [-360, 360], "orientation",
+               "rotation about beam"],
+              ["psi", "degrees", 0, [-360, 360], "orientation",
+               "rotation about c axis"],
              ]
 
@@ -186,9 +195,10 @@
             psi_pd=10, psi_pd_n=1)
 
-qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5)
+# rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models
+qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
 tests = [[{}, 0.2, 0.533149288477],
          [{}, [0.2], [0.533149288477]],
-         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.032102135569],
-         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.032102135569]],
+         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222],
+         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]],
         ]
 del qx, qy  # not necessary to delete, but cleaner
Index: sasmodels/models/cylinder.py
===================================================================
--- sasmodels/models/cylinder.py	(revision 3330bb44faf647c84e605c5a75bed3f385a25a52)
+++ sasmodels/models/cylinder.py	(revision 9802ab3945a154b871264715dcb463d4ec644afa)
@@ -38,6 +38,6 @@
 
 
-Numerical integration is simplified by a change of variable to $u = cos(\alpha)$ with 
-$sin(\alpha)=\sqrt{1-u^2}$. 
+Numerical integration is simplified by a change of variable to $u = cos(\alpha)$ with
+$sin(\alpha)=\sqrt{1-u^2}$.
 
 The output of the 1D scattering intensity function for randomly oriented
@@ -61,9 +61,21 @@
 .. _cylinder-angle-definition:
 
-.. figure:: img/cylinder_angle_definition.jpg
+.. figure:: img/cylinder_angle_definition.png
 
-    Definition of the angles for oriented cylinders.
+    Definition of the $\theta$ and $\phi$ orientation angles for a cylinder relative 
+    to the beam line coordinates, plus an indication of their orientation distributions 
+    which are described as rotations about each of the perpendicular axes $\delta_1$ and $\delta_2$ 
+    in the frame of the cylinder itself, which when $\theta = \phi = 0$ are parallel to the $Y$ and $X$ axes.
 
-The $\theta$ and $\phi$ parameters only appear in the model when fitting 2d data.
+.. figure:: img/cylinder_angle_projection.png
+
+    Examples for oriented cylinders.
+
+The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data. 
+On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will
+appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, which when $\theta = \phi = 0$ are parallel 
+to the $Y$ and $X$ axes of the instrument respectively. Some experimentation may be required to understand the 2d patterns fully.
+(Earlier implementations had numerical integration issues in some circumstances when orientation distributions passed through 90 degrees, such 
+situations, with very broad distributions, should still be approached with care.) 
 
 Validation
@@ -123,8 +135,8 @@
               ["length", "Ang", 400, [0, inf], "volume",
                "Cylinder length"],
-              ["theta", "degrees", 60, [-inf, inf], "orientation",
-               "latitude"],
-              ["phi", "degrees", 60, [-inf, inf], "orientation",
-               "longitude"],
+              ["theta", "degrees", 60, [-360, 360], "orientation",
+               "cylinder axis to beam angle"],
+              ["phi", "degrees",   60, [-360, 360], "orientation",
+               "rotation about beam"],
              ]
 
@@ -152,5 +164,5 @@
 tests = [[{}, 0.2, 0.042761386790780453],
         [{}, [0.2], [0.042761386790780453]],
-#  new coords    
+#  new coords
         [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852],
         [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]],
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 9b79f2966c4083a48b4894b6bf0587bef7a5be09)
@@ -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,7 +111,14 @@
 *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
 """
 
-from numpy import inf
+from numpy import inf, sin, cos, pi
 
 name = "ellipsoid"
@@ -129,8 +151,8 @@
               ["radius_equatorial", "Ang", 400, [0, inf], "volume",
                "Equatorial radius"],
-              ["theta", "degrees", 60, [-inf, inf], "orientation",
-               "In plane angle"],
-              ["phi", "degrees", 60, [-inf, inf], "orientation",
-               "Out of plane angle"],
+              ["theta", "degrees", 60, [-360, 360], "orientation",
+               "ellipsoid axis to beam angle"],
+              ["phi", "degrees", 60, [-360, 360], "orientation",
+               "rotation about beam"],
              ]
 
@@ -139,5 +161,5 @@
 def ER(radius_polar, radius_equatorial):
     import numpy as np
-
+# see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
     ee = np.empty_like(radius_polar)
     idx = radius_polar > radius_equatorial
@@ -168,2 +190,10 @@
             theta_pd=15, theta_pd_n=45,
             phi_pd=15, phi_pd_n=1)
+q = 0.1
+# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
+qx = q*cos(pi/6.0)
+qy = q*sin(pi/6.0)
+tests = [[{}, 0.05, 54.8525847025],
+        [{'theta':80., 'phi':10.}, (qx, qy), 1.74134670026 ],
+        ]
+del qx, qy  # not necessary to delete, but cleaner
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/elliptical_cylinder.py
===================================================================
--- sasmodels/models/elliptical_cylinder.py	(revision fcb33e4672814b4a8fe1ee865d430090b5fc03e6)
+++ sasmodels/models/elliptical_cylinder.py	(revision 9802ab3945a154b871264715dcb463d4ec644afa)
@@ -57,19 +57,26 @@
 define the axis of the cylinder using two angles $\theta$, $\phi$ and $\Psi$
 (see :ref:`cylinder orientation <cylinder-angle-definition>`). The angle
-$\Psi$ is the rotational angle around its own long_c axis against the $q$ plane.
-For example, $\Psi = 0$ when the $r_\text{minor}$ axis is parallel to the
-$x$ axis of the detector.
+$\Psi$ is the rotational angle around its own long_c axis. 
 
 All angle parameters are valid and given only for 2D calculation; ie, an
 oriented system.
 
-.. figure:: img/elliptical_cylinder_angle_definition.jpg
+.. figure:: img/elliptical_cylinder_angle_definition.png
 
-    Definition of angles for 2D
+    Definition of angles for oriented elliptical cylinder, where axis_ratio is drawn >1,
+    and angle $\Psi$ is now a rotation around the axis of the cylinder.
 
-.. figure:: img/cylinder_angle_projection.jpg
+.. figure:: img/elliptical_cylinder_angle_projection.png
 
     Examples of the angles for oriented elliptical cylinders against the
-    detector plane.
+    detector plane, with $\Psi$ = 0.
+
+The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data. 
+On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will
+appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, the $b$ and $a$ axes of the 
+cylinder cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.) 
+The third orientation distribution, in $\psi$, is about the $c$ axis of the particle. Some experimentation may be required to 
+understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation 
+distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.) 
 
 NB: The 2nd virial coefficient of the cylinder is calculated based on the
@@ -108,5 +115,5 @@
 """
 
-from numpy import pi, inf, sqrt
+from numpy import pi, inf, sqrt, sin, cos
 
 name = "elliptical_cylinder"
@@ -125,7 +132,7 @@
               ["sld",         "1e-6/Ang^2", 4.0,   [-inf, inf], "sld",         "Cylinder scattering length density"],
               ["sld_solvent", "1e-6/Ang^2", 1.0,   [-inf, inf], "sld",         "Solvent scattering length density"],
-              ["theta",       "degrees",    90.0,  [-360, 360], "orientation", "In plane angle"],
-              ["phi",         "degrees",    0,     [-360, 360], "orientation", "Out of plane angle"],
-              ["psi",         "degrees",    0,     [-360, 360], "orientation", "Major axis angle relative to Q"]]
+              ["theta",       "degrees",    90.0,  [-360, 360], "orientation", "cylinder axis to beam angle"],
+              ["phi",         "degrees",    0,     [-360, 360], "orientation", "rotation about beam"],
+              ["psi",         "degrees",    0,     [-360, 360], "orientation", "rotation about cylinder axis"]]
 
 # pylint: enable=bad-whitespace, line-too-long
@@ -149,4 +156,8 @@
                            + (length + radius) * (length + pi * radius))
     return 0.5 * (ddd) ** (1. / 3.)
+q = 0.1
+# april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct!
+qx = q*cos(pi/6.0)
+qy = q*sin(pi/6.0)
 
 tests = [
@@ -158,3 +169,4 @@
       'sld_solvent':1.0, 'background':0.0},
      0.001, 675.504402],
+#    [{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ],
 ]
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/fcc_paracrystal.py
===================================================================
--- sasmodels/models/fcc_paracrystal.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
+++ sasmodels/models/fcc_paracrystal.py	(revision 69e1afc57270f67b2a0d7157a6d263fc2db98cc0)
@@ -76,7 +76,8 @@
 2D model computation.
 
-.. figure:: img/bcc_angle_definition.png
+.. figure:: img/parallelepiped_angle_definition.png
 
-    Orientation of the crystal with respect to the scattering plane.
+    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).
 
 References
@@ -90,5 +91,5 @@
 """
 
-from numpy import inf
+from numpy import inf, pi
 
 name = "fcc_paracrystal"
@@ -110,7 +111,7 @@
               ["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Particle scattering length density"],
               ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Solvent scattering length density"],
-              ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"],
-              ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"],
-              ["psi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"]
+              ["theta",       "degrees",    60,    [-360, 360], "orientation", "c axis to beam angle"],
+              ["phi",         "degrees",    60,    [-360, 360], "orientation", "rotation about beam"],
+              ["psi",         "degrees",    60,    [-360, 360], "orientation", "rotation about c axis"]
              ]
 # pylint: enable=bad-whitespace, line-too-long
@@ -128,2 +129,10 @@
             psi_pd=15, psi_pd_n=0,
            )
+# april 10 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
+q =4.*pi/220.
+tests = [
+    [{ },
+     [0.001, q, 0.215268], [0.275164706668, 5.7776842567, 0.00958167119232]],
+     [{}, (-0.047,-0.007), 238.103096286],
+     [{}, (0.053,0.063), 0.863609587796 ],
+]
Index: sasmodels/models/hollow_cylinder.py
===================================================================
--- sasmodels/models/hollow_cylinder.py	(revision aea2e2ab3f11c42b7fc1a72a02b6784076c41d3a)
+++ sasmodels/models/hollow_cylinder.py	(revision 9b79f2966c4083a48b4894b6bf0587bef7a5be09)
@@ -60,5 +60,5 @@
 """
 
-from numpy import pi, inf
+from numpy import pi, inf, sin, cos
 
 name = "hollow_cylinder"
@@ -82,6 +82,6 @@
     ["sld",         "1/Ang^2",  6.3, [-inf, inf], "sld",         "Cylinder sld"],
     ["sld_solvent", "1/Ang^2",  1,   [-inf, inf], "sld",         "Solvent sld"],
-    ["theta",       "degrees", 90,   [-360, 360], "orientation", "Theta angle"],
-    ["phi",         "degrees",  0,   [-360, 360], "orientation", "Phi angle"],
+    ["theta",       "degrees", 90,   [-360, 360], "orientation", "Cylinder axis to beam angle"],
+    ["phi",         "degrees",  0,   [-360, 360], "orientation", "Rotation about beam"],
     ]
 # pylint: enable=bad-whitespace, line-too-long
@@ -129,9 +129,14 @@
             theta_pd=10, theta_pd_n=5,
            )
-
+q = 0.1
+# april 6 2017, rkh added a 2d unit test, assume correct!
+qx = q*cos(pi/6.0)
+qy = q*sin(pi/6.0)
 # Parameters for unit tests
 tests = [
     [{}, 0.00005, 1764.926],
     [{}, 'VR', 1.8],
-    [{}, 0.001, 1756.76]
-    ]
+    [{}, 0.001, 1756.76],
+    [{}, (qx, qy), 2.36885476192  ],
+        ]
+del qx, qy  # not necessary to delete, but cleaner
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/parallelepiped.py
===================================================================
--- sasmodels/models/parallelepiped.py	(revision 3330bb44faf647c84e605c5a75bed3f385a25a52)
+++ sasmodels/models/parallelepiped.py	(revision 34a9e4eed11f1952a0e9680cbfbdbfa79e7cffe3)
@@ -9,10 +9,11 @@
 ----------
 
-| This model calculates the scattering from a rectangular parallelepiped
-| (\:numref:`parallelepiped-image`\).
-| If you need to apply polydispersity, see also :ref:`rectangular-prism`.
+ This model calculates the scattering from a rectangular parallelepiped
+ (\:numref:`parallelepiped-image`\).
+ If you need to apply polydispersity, see also :ref:`rectangular-prism`.
 
 .. _parallelepiped-image:
 
+
 .. figure:: img/parallelepiped_geometry.jpg
 
@@ -21,7 +22,9 @@
 .. note::
 
-   The edge of the solid must satisfy the condition that $A < B < C$.
-   This requirement is not enforced in the model, so it is up to the
-   user to check this during the analysis.
+The three dimensions of the parallelepiped (strictly here a cuboid) may be given in 
+$any$ size order. To avoid multiple fit solutions, especially
+with Monte-Carlo fit methods, it may be advisable to restrict their ranges. There may 
+be a number of closely similar "best fits", so some trial and error, or fixing of some 
+dimensions at expected values, may help.
 
 The 1D scattering intensity $I(q)$ is calculated as:
@@ -67,9 +70,9 @@
     \mu &= qB
 
-
 The scattering intensity per unit volume is returned in units of |cm^-1|.
 
 NB: The 2nd virial coefficient of the parallelepiped is calculated based on
-the averaged effective radius $(=\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.
@@ -102,13 +105,21 @@
 .. _parallelepiped-orientation:
 
-.. figure:: img/parallelepiped_angle_definition.jpg
-
-    Definition of the angles for oriented parallelepipeds.
-
-.. figure:: img/parallelepiped_angle_projection.jpg
-
-    Examples of the angles for oriented parallelepipeds against the
+.. figure:: img/parallelepiped_angle_definition.png
+
+    Definition of the angles for oriented parallelepiped, shown with $A<B<C$.
+
+.. figure:: img/parallelepiped_angle_projection.png
+
+    Examples of the angles for an oriented parallelepiped against the
     detector plane.
 
+On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will
+appear. These are actually rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped, perpendicular to the $a$ x $c$ and $b$ x $c$ faces. 
+(When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.) The third orientation distribution, in $\psi$, is 
+about the $c$ axis of the particle, perpendicular to the $a$ x $b$ face. Some experimentation may be required to 
+understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation 
+distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.) 
+
+    
 For a given orientation of the parallelepiped, the 2D form factor is
 calculated as
@@ -116,7 +127,7 @@
 .. math::
 
-    P(q_x, q_y) = \left[\frac{\sin(qA\cos\alpha/2)}{(qA\cos\alpha/2)}\right]^2
-                  \left[\frac{\sin(qB\cos\beta/2)}{(qB\cos\beta/2)}\right]^2
-                  \left[\frac{\sin(qC\cos\gamma/2)}{(qC\cos\gamma/2)}\right]^2
+    P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}{2}qA\cos\alpha)}\right]^2
+                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}{2}qB\cos\beta)}\right]^2
+                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}{2}qC\cos\gamma)}\right]^2
 
 with
@@ -154,6 +165,4 @@
 angles.
 
-This model is based on form factor calculations implemented in a c-library
-provided by the NIST Center for Neutron Research (Kline, 2006).
 
 References
@@ -163,8 +172,17 @@
 
 R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
+
+Authorship and Verification
+----------------------------
+
+* **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
+
 """
 
 import numpy as np
-from numpy import pi, inf, sqrt
+from numpy import pi, inf, sqrt, sin, cos
 
 name = "parallelepiped"
@@ -180,5 +198,5 @@
             mu = q*B
         V: Volume of the rectangular parallelepiped
-        alpha: angle between the long axis of the 
+        alpha: angle between the long axis of the
             parallelepiped and the q-vector for 1D
 """
@@ -196,10 +214,10 @@
               ["length_c", "Ang", 400, [0, inf], "volume",
                "Larger side of the parallelepiped"],
-              ["theta", "degrees", 60, [-inf, inf], "orientation",
-               "In plane angle"],
-              ["phi", "degrees", 60, [-inf, inf], "orientation",
-               "Out of plane angle"],
-              ["psi", "degrees", 60, [-inf, inf], "orientation",
-               "Rotation angle around its own c axis against q plane"],
+              ["theta", "degrees", 60, [-360, 360], "orientation",
+               "c axis to beam angle"],
+              ["phi", "degrees", 60, [-360, 360], "orientation",
+               "rotation about beam"],
+              ["psi", "degrees", 60, [-360, 360], "orientation",
+               "rotation about c axis"],
              ]
 
@@ -208,11 +226,16 @@
 def ER(length_a, length_b, length_c):
     """
-        Return effective radius (ER) for P(q)*S(q)
+    Return effective radius (ER) for P(q)*S(q)
     """
-
+    # now that axes can be in any size order, need to sort a,b,c where a~b and c is either much smaller
+    # or much larger
+    abc = np.vstack((length_a, length_b, length_c))
+    abc = np.sort(abc, axis=0)
+    selector = (abc[1] - abc[0]) > (abc[2] - abc[1])
+    length = np.where(selector, abc[0], abc[2])
     # surface average radius (rough approximation)
-    surf_rad = sqrt(length_a * length_b / pi)
-
-    ddd = 0.75 * surf_rad * (2 * surf_rad * length_c + (length_c + surf_rad) * (length_c + pi * surf_rad))
+    radius = np.sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi)
+
+    ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius))
     return 0.5 * (ddd) ** (1. / 3.)
 
@@ -230,10 +253,10 @@
             phi_pd=10, phi_pd_n=1,
             psi_pd=10, psi_pd_n=10)
-
-qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5)
+# rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models
+qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
 tests = [[{}, 0.2, 0.17758004974],
          [{}, [0.2], [0.17758004974]],
-         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.00560296014],
-         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.00560296014]],
+         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0089517140475],
+         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0089517140475]],
         ]
 del qx, qy  # not necessary to delete, but cleaner
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/sc_paracrystal.py
===================================================================
--- sasmodels/models/sc_paracrystal.py	(revision 3330bb44faf647c84e605c5a75bed3f385a25a52)
+++ sasmodels/models/sc_paracrystal.py	(revision 69e1afc57270f67b2a0d7157a6d263fc2db98cc0)
@@ -83,5 +83,8 @@
 model computation.
 
-.. figure:: img/sc_crystal_angle_definition.jpg
+.. figure:: img/parallelepiped_angle_definition.png
+
+    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).
 
 Reference
@@ -127,7 +130,7 @@
               ["sld",  "1e-6/Ang^2",         3.0, [0.0, inf],  "sld",         "Sphere scattering length density"],
               ["sld_solvent", "1e-6/Ang^2",  6.3, [0.0, inf],  "sld",         "Solvent scattering length density"],
-              ["theta",       "degrees",     0.0, [-inf, inf], "orientation", "Orientation of the a1 axis w/respect incoming beam"],
-              ["phi",         "degrees",     0.0, [-inf, inf], "orientation", "Orientation of the a2 in the plane of the detector"],
-              ["psi",         "degrees",     0.0, [-inf, inf], "orientation", "Orientation of the a3 in the plane of the detector"],
+              ["theta",       "degrees",    0,    [-360, 360], "orientation", "c axis to beam angle"],
+              ["phi",         "degrees",    0,    [-360, 360], "orientation", "rotation about beam"],
+              ["psi",         "degrees",    0,    [-360, 360], "orientation", "rotation about c axis"]
              ]
 # pylint: enable=bad-whitespace, line-too-long
@@ -146,8 +149,10 @@
 
 tests = [
-    # Accuracy tests based on content in test/utest_extra_models.py
+    # Accuracy tests based on content in test/utest_extra_models.py, 2d tests added April 10, 2017
     [{}, 0.001, 10.3048],
     [{}, 0.215268, 0.00814889],
-    [{}, (0.414467), 0.001313289]
+    [{}, (0.414467), 0.001313289],
+    [{'theta':10.0,'phi':20,'psi':30.0},(0.045,-0.035),18.0397138402 ],
+    [{'theta':10.0,'phi':20,'psi':30.0},(0.023,0.045),0.0177333171285 ]
     ]
 
Index: sasmodels/models/stacked_disks.py
===================================================================
--- sasmodels/models/stacked_disks.py	(revision b7e8b94b477e6d9a7e4e2666bb218d3568ee66d2)
+++ sasmodels/models/stacked_disks.py	(revision 9802ab3945a154b871264715dcb463d4ec644afa)
@@ -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::
@@ -77,5 +77,5 @@
 the axis of the cylinder using two angles $\theta$ and $\varphi$.
 
-.. figure:: img/cylinder_angle_definition.jpg
+.. figure:: img/cylinder_angle_definition.png
 
     Examples of the angles against the detector plane.
@@ -103,5 +103,5 @@
 """
 
-from numpy import inf
+from numpy import inf, sin, cos, pi
 
 name = "stacked_disks"
@@ -131,6 +131,6 @@
     ["sld_layer",   "1e-6/Ang^2",  0.0, [-inf, inf], "sld",         "Layer scattering length density"],
     ["sld_solvent", "1e-6/Ang^2",  5.0, [-inf, inf], "sld",         "Solvent scattering length density"],
-    ["theta",       "degrees",     0,   [-inf, inf], "orientation", "Orientation of the stacked disk axis w/respect incoming beam"],
-    ["phi",         "degrees",     0,   [-inf, inf], "orientation", "Orientation of the stacked disk in the plane of the detector"],
+    ["theta",       "degrees",     0,   [-360, 360], "orientation", "Orientation of the stacked disk axis w/respect incoming beam"],
+    ["phi",         "degrees",     0,   [-360, 360], "orientation", "Rotation about beam"],
     ]
 # pylint: enable=bad-whitespace, line-too-long
@@ -152,8 +152,13 @@
 # After redefinition of spherical coordinates -
 # tests had in old coords theta=0, phi=0; new coords theta=90, phi=0
-# but should not matter here as so far all the tests are 1D not 2D
+q = 0.1
+# april 6 2017, rkh added a 2d unit test, assume correct!
+qx = q*cos(pi/6.0)
+qy = q*sin(pi/6.0)
 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,
@@ -186,4 +191,17 @@
     [{'thick_core': 10.0,
       'thick_layer': 15.0,
+      'radius': 100.0,
+      'n_stacking': 5,
+      'sigma_d': 0.0,
+      'sld_core': 4.0,
+      'sld_layer': -0.4,
+      'solvent_sd': 5.0,
+      'theta': 90.0,
+      'phi': 20.0,
+      'scale': 0.01,
+      'background': 0.001},
+      (qx, qy), 0.0491167089952  ],
+    [{'thick_core': 10.0,
+      'thick_layer': 15.0,
       'radius': 3000.0,
       'n_stacking': 5,
@@ -228,4 +246,17 @@
       'background': 0.001,
      }, ([0.4, 0.5]), [0.00105074, 0.00121761]],
+    [{'thick_core': 10.0,
+      'thick_layer': 15.0,
+      'radius': 3000.0,
+      'n_stacking': 1.0,
+      'sigma_d': 0.0,
+      'sld_core': 4.0,
+      'sld_layer': -0.4,
+      'solvent_sd': 5.0,
+      'theta': 90.0,
+      'phi': 20.0,
+      'scale': 0.01,
+      'background': 0.001,
+     }, (qx, qy), 0.0341738733124 ],
 
     [{'thick_core': 10.0,
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 3330bb44faf647c84e605c5a75bed3f385a25a52)
+++ sasmodels/models/triaxial_ellipsoid.py	(revision 34a9e4eed11f1952a0e9680cbfbdbfa79e7cffe3)
@@ -2,28 +2,25 @@
 # 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 for randomly oriented particles is defined by the average over
+all orientations $\Omega$ of:
 
 .. 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}(\Delta\rho)^2\frac{V}{4 \pi}\int_\Omega\Phi^2(qr)\,d\Omega
+           + \text{background}
 
 where
@@ -31,29 +28,77 @@
 .. 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)
+
+Though for convenience we describe the three radii of the ellipsoid as equatorial
+and polar, they may be given in $any$ size order. To avoid multiple solutions, especially
+with Monte-Carlo fit methods, it may be advisable to restrict their ranges. For typical
+small angle diffraction situations there may be a number of closely similar "best fits",
+so some trial and error, or fixing of some radii at expected values, may help.
+    
 To provide easy access to the orientation of the triaxial ellipsoid,
 we define the axis of the cylinder using the angles $\theta$, $\phi$
-and $\psi$. These angles are defined on
-:numref:`triaxial-ellipsoid-angles` .
-The angle $\psi$ is the rotational angle around its own $c$ axis
-against the $q$ plane. For example, $\psi = 0$ when the
-$a$ axis is parallel to the $x$ axis of the detector.
+and $\psi$. These angles are defined analogously to the elliptical_cylinder below, note that
+angle $\phi$ is now NOT the same as in the equations above.
+
+.. figure:: img/elliptical_cylinder_angle_definition.png
+
+    Definition of angles for oriented triaxial ellipsoid, where radii are for illustration here 
+    $a < b << c$ and angle $\Psi$ is a rotation around the axis of the particle.
+
+For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data, 
+see the :ref:`elliptical-cylinder` model for further information.
 
 .. _triaxial-ellipsoid-angles:
 
-.. figure:: img/triaxial_ellipsoid_angle_projection.jpg
+.. figure:: img/triaxial_ellipsoid_angle_projection.png
 
-    The angles for oriented ellipsoid.
+    Some examples for an oriented triaxial ellipsoid.
 
 The radius-of-gyration for this system is  $R_g^2 = (R_a R_b R_c)^2/5$.
 
-The contrast is defined as SLD(ellipsoid) - SLD(solvent).  In the
+The contrast $\Delta\rho$ is defined as SLD(ellipsoid) - SLD(solvent).  In the
 parameters, $R_a$ is the minor equatorial radius, $R_b$ is the major
 equatorial radius, and $R_c$ is the polar radius of the ellipsoid.
 
 NB: The 2nd virial coefficient of the triaxial solid ellipsoid is
-calculated based on the polar radius $R_p = R_c$ and equatorial
-radius $R_e = \sqrt{R_a R_b}$, and used as the effective radius for
+calculated after sorting the three radii to give the most appropriate
+prolate or oblate form, from the new polar radius $R_p = R_c$ and effective equatorial
+radius,  $R_e = \sqrt{R_a R_b}$, to then be used as the effective radius for
 $S(q)$ when $P(q) \cdot S(q)$ is applied.
 
@@ -69,17 +114,23 @@
 ----------
 
-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
+
+Authorship and Verification
+----------------------------
+
+* **Author:** NIST IGOR/DANSE **Date:** pre 2010
+* **Last Modified by:** Paul Kienzle (improved calculation) **Date:** April 4, 2017
+* **Last Reviewed by:** Paul Kienzle & Richard Heenan **Date:**  April 4, 2017
 """
 
-from numpy import inf
+from numpy import inf, sin, cos, pi
 
 name = "triaxial_ellipsoid"
 title = "Ellipsoid of uniform scattering length density with three independent axes."
 
-description = """\
-Note: During fitting ensure that the inequality ra<rb<rc is not
-	violated. Otherwise the calculation will
-	not be correct.
+description = """
+   Triaxial ellipsoid - see main documentation.
 """
 category = "shape:ellipsoid"
@@ -91,15 +142,15 @@
                "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"],
-              ["theta", "degrees", 60, [-inf, inf], "orientation",
-               "In plane angle"],
-              ["phi", "degrees", 60, [-inf, inf], "orientation",
-               "Out of plane angle"],
-              ["psi", "degrees", 60, [-inf, inf], "orientation",
-               "Out of plane angle"],
+               "Polar radius, Rc"],
+              ["theta", "degrees", 60, [-360, 360], "orientation",
+               "polar axis to beam angle"],
+              ["phi", "degrees", 60, [-360, 360], "orientation",
+               "rotation about beam"],
+              ["psi", "degrees", 60, [-360, 360], "orientation",
+               "rotation about polar axis"],
              ]
 
@@ -108,9 +159,17 @@
 def ER(radius_equat_minor, radius_equat_major, radius_polar):
     """
-        Returns the effective radius used in the S*P calculation
+    Returns the effective radius used in the S*P calculation
     """
     import numpy as np
     from .ellipsoid import ER as ellipsoid_ER
-    return ellipsoid_ER(radius_polar, np.sqrt(radius_equat_minor * radius_equat_major))
+
+    # 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)
+    selector = (radii[1] - radii[0]) > (radii[2] - radii[1])
+    polar = np.where(selector, radii[0], radii[2])
+    equatorial = np.sqrt(np.where(~selector, radii[0]*radii[1], radii[1]*radii[2]))
+    return ellipsoid_ER(polar, equatorial)
 
 demo = dict(scale=1, background=0,
@@ -124,2 +183,13 @@
             phi_pd=15, phi_pd_n=1,
             psi_pd=15, psi_pd_n=1)
+
+q = 0.1
+# april 6 2017, rkh add unit tests
+#     NOT compared with any other calc method, assume correct!
+# check 2d test after pull #890
+qx = q*cos(pi/6.0)
+qy = q*sin(pi/6.0)
+tests = [[{}, 0.05, 24.8839548033],
+        [{'theta':80., 'phi':10.}, (qx, qy), 166.712060266 ],
+        ]
+del qx, qy  # not necessary to delete, but cleaner
