Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision 650c6d2665bcebbe07614117e6b9841191a18ea2)
+++ sasmodels/compare.py	(revision 650c6d2665bcebbe07614117e6b9841191a18ea2)
@@ -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
@@ -836,5 +836,5 @@
     'linear', 'log', 'q4',
     'hist', 'nohist',
-    'edit', 'html',
+    'edit', 'html', 'help',
     'demo', 'default',
     ])
@@ -1005,4 +1005,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/barbell.py
===================================================================
--- sasmodels/models/barbell.py	(revision fcb33e4672814b4a8fe1ee865d430090b5fc03e6)
+++ sasmodels/models/barbell.py	(revision 0b56f3871b1019e3bb99fb618e1342d5f26c8de4)
@@ -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"
@@ -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 e2d6e3b6cf0d284837fe50be3eb7b7954f57edef)
@@ -99,5 +99,5 @@
 """
 
-from numpy import inf
+from numpy import inf, pi
 
 name = "bcc_paracrystal"
@@ -141,2 +141,9 @@
     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]],
+]
Index: sasmodels/models/capped_cylinder.py
===================================================================
--- sasmodels/models/capped_cylinder.py	(revision fcb33e4672814b4a8fe1ee865d430090b5fc03e6)
+++ sasmodels/models/capped_cylinder.py	(revision 0b56f3871b1019e3bb99fb618e1342d5f26c8de4)
@@ -91,5 +91,5 @@
 
 """
-from numpy import inf
+from numpy import inf, sin, cos, pi
 
 name = "capped_cylinder"
@@ -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 0b56f3871b1019e3bb99fb618e1342d5f26c8de4)
@@ -88,5 +88,5 @@
 """
 
-from numpy import inf, sin, cos
+from numpy import inf, sin, cos, pi
 
 name = "core_shell_bicelle"
@@ -155,4 +155,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 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_bicelle_elliptical.py
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical.py	(revision 8afefaee3d363dde220aa55fe83181a73b621f98)
+++ sasmodels/models/core_shell_bicelle_elliptical.py	(revision 15a90c1aa6a67491d6170c2aa26907b5939da063)
@@ -80,8 +80,8 @@
 
 
-.. 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"
@@ -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 0b56f3871b1019e3bb99fb618e1342d5f26c8de4)
@@ -73,5 +73,5 @@
 """
 
-from numpy import pi, inf
+from numpy import pi, inf, sin, cos
 
 name = "core_shell_cylinder"
@@ -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_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 1916c5233063b68d94c979bea7935dfa475ffa95)
@@ -10,6 +10,6 @@
 
 .. 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"
@@ -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 3330bb44faf647c84e605c5a75bed3f385a25a52)
@@ -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
@@ -64,4 +64,8 @@
 
     Definition of the angles for oriented cylinders.
+
+.. figure:: img/cylinder_angle_projection.png
+
+    Examples for oriented cylinders.
 
 The $\theta$ and $\phi$ parameters only appear in the model when fitting 2d data.
@@ -152,5 +156,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 0b56f3871b1019e3bb99fb618e1342d5f26c8de4)
@@ -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"
@@ -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 15a90c1aa6a67491d6170c2aa26907b5939da063)
@@ -64,12 +64,13 @@
 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 >1,
+    and angle $\Psi$ is 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.
 
 NB: The 2nd virial coefficient of the cylinder is calculated based on the
@@ -108,5 +109,5 @@
 """
 
-from numpy import pi, inf, sqrt
+from numpy import pi, inf, sqrt, sin, cos
 
 name = "elliptical_cylinder"
@@ -149,4 +150,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 +163,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 e2d6e3b6cf0d284837fe50be3eb7b7954f57edef)
@@ -90,5 +90,5 @@
 """
 
-from numpy import inf
+from numpy import inf, pi
 
 name = "fcc_paracrystal"
@@ -128,2 +128,9 @@
             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], [0.275164706668, 5.7776842567, 0.00958167119232]],
+]
Index: sasmodels/models/hollow_cylinder.py
===================================================================
--- sasmodels/models/hollow_cylinder.py	(revision aea2e2ab3f11c42b7fc1a72a02b6784076c41d3a)
+++ sasmodels/models/hollow_cylinder.py	(revision 0b56f3871b1019e3bb99fb618e1342d5f26c8de4)
@@ -60,5 +60,5 @@
 """
 
-from numpy import pi, inf
+from numpy import pi, inf, sin, cos
 
 name = "hollow_cylinder"
@@ -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 3330bb44faf647c84e605c5a75bed3f385a25a52)
@@ -15,4 +15,5 @@
 .. _parallelepiped-image:
 
+
 .. figure:: img/parallelepiped_geometry.jpg
 
@@ -21,7 +22,7 @@
 .. 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 edge of the solid used to have to satisfy the condition that $A < B < C$.
+   After some improvements to the effective radius calculation, used with an S(Q),
+   it is beleived that this is no longer the case.
 
 The 1D scattering intensity $I(q)$ is calculated as:
@@ -71,5 +72,6 @@
 
 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{A B / \pi})$ and
 length $(= C)$ values, and used as the effective radius for
 $S(q)$ when $P(q) \cdot S(q)$ is applied.
@@ -102,11 +104,11 @@
 .. _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.
 
@@ -116,7 +118,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 +156,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 +163,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 +189,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
 """
@@ -208,11 +217,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 +244,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/stacked_disks.py
===================================================================
--- sasmodels/models/stacked_disks.py	(revision b7e8b94b477e6d9a7e4e2666bb218d3568ee66d2)
+++ sasmodels/models/stacked_disks.py	(revision 0b56f3871b1019e3bb99fb618e1342d5f26c8de4)
@@ -103,5 +103,5 @@
 """
 
-from numpy import inf
+from numpy import inf, sin, cos, pi
 
 name = "stacked_disks"
@@ -152,5 +152,8 @@
 # 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.
@@ -186,4 +189,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 +244,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 3330bb44faf647c84e605c5a75bed3f385a25a52)
@@ -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 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,10 +26,50 @@
 .. 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,
 we define the axis of the cylinder using the angles $\theta$, $\phi$
-and $\psi$. These angles are defined on
-:numref:`triaxial-ellipsoid-angles` .
+and $\psi$. These angles are defined analogously to the elliptical_cylinder below
+
+.. figure:: img/elliptical_cylinder_angle_definition.png
+
+    Definition of angles for oriented triaxial ellipsoid, where radii shown here are $a < b << c$
+    and angle $\Psi$ is a rotation around the axis of the particle.
+
 The angle $\psi$ is the rotational angle around its own $c$ axis
 against the $q$ plane. For example, $\psi = 0$ when the
@@ -43,11 +78,11 @@
 .. _triaxial-ellipsoid-angles:
 
-.. figure:: img/triaxial_ellipsoid_angle_projection.jpg
+.. figure:: img/triaxial_ellipsoid_angle_projection.png
 
-    The angles for oriented ellipsoid.
+    Some example angles for oriented 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.
@@ -69,9 +104,18 @@
 ----------
 
-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"
@@ -91,9 +135,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"],
@@ -108,9 +152,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 +176,12 @@
             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!
+# add 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), 9999. ],
+        ]
+del qx, qy  # not necessary to delete, but cleaner
