Index: sasmodels/kernel_template.c
===================================================================
--- sasmodels/kernel_template.c	(revision 6a0d6aac6a8b68592f9fd9c4b56139af8869b275)
+++ sasmodels/kernel_template.c	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -264,6 +264,6 @@
       #if defined(IQXY_HAS_THETA)
         // Force a nominal value for the spherical correction even when
-        // theta is +90/-90 so that there are no divide by zero problems.
-        // For cos(theta) fixed at 90, we effectively multiply top and bottom
+        // theta is +0/180 so that there are no divide by zero problems.
+        // For sin(theta) fixed at 0 and 180, we effectively multiply top and bottom
         // by 1e-6, so the effect cancels.
         const double spherical_correction = fmax(fabs(cos(M_PI_180*theta)), 1.e-6);
Index: sasmodels/kernelpy.py
===================================================================
--- sasmodels/kernelpy.py	(revision 14a15a34e512ac50d2cc517fc8d4243ec0bcd9ce)
+++ sasmodels/kernelpy.py	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -10,5 +10,5 @@
 
 import numpy as np  # type: ignore
-from numpy import pi, cos  #type: ignore
+from numpy import pi, sin, cos  #type: ignore
 
 from . import details
@@ -220,5 +220,5 @@
             partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:])
             if call_details.theta_par >= 0:
-                cor = cos(pi / 180 * parameters[call_details.theta_par])
+                cor = sin(pi / 180 * parameters[call_details.theta_par])
                 spherical_correction = max(abs(cor), 1e-6)
             p0_index = loop_index%p0_length
Index: sasmodels/models/barbell.c
===================================================================
--- sasmodels/models/barbell.c	(revision 22221346d45e39a667fe231ce7ca5ecfc79b0ac5)
+++ sasmodels/models/barbell.c	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -95,7 +95,7 @@
     // Compute angle alpha between q and the cylinder axis
     double sn, cn; // slots to hold sincos function output
-    SINCOS(theta*M_PI_180, sn, cn);
-    const double q = sqrt(qx*qx+qy*qy);
-    const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);
+    SINCOS(phi*M_PI_180, sn, cn);
+    const double q = sqrt(qx*qx + qy*qy);
+    const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q);
     const double alpha = acos(cos_val); // rod angle relative to q
 
Index: sasmodels/models/barbell.py
===================================================================
--- sasmodels/models/barbell.py	(revision b0c42719ac9096c35137bf644db59e1c5915b06b)
+++ sasmodels/models/barbell.py	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -72,8 +72,4 @@
     Definition of the angles for oriented 2D barbells.
 
-.. figure:: img/cylinder_angle_projection.jpg
-    :width: 600px
-
-    Examples of the angles for oriented pp against the detector plane.
 
 References
Index: sasmodels/models/capped_cylinder.c
===================================================================
--- sasmodels/models/capped_cylinder.c	(revision 22221346d45e39a667fe231ce7ca5ecfc79b0ac5)
+++ sasmodels/models/capped_cylinder.c	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -116,7 +116,7 @@
     // Compute angle alpha between q and the cylinder axis
     double sn, cn;
-    SINCOS(theta*M_PI_180, sn, cn);
+    SINCOS(phi*M_PI_180, sn, cn);
     const double q = sqrt(qx*qx+qy*qy);
-    const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);
+    const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q);
     const double alpha = acos(cos_val); // rod angle relative to q
 
Index: sasmodels/models/capped_cylinder.py
===================================================================
--- sasmodels/models/capped_cylinder.py	(revision b0c42719ac9096c35137bf644db59e1c5915b06b)
+++ sasmodels/models/capped_cylinder.py	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -75,8 +75,4 @@
     Definition of the angles for oriented 2D cylinders.
 
-.. figure:: img/cylinder_angle_projection.jpg
-    :width: 600px
-
-    Examples of the angles for oriented 2D cylinders against the detector plane.
 
 References
@@ -132,6 +128,6 @@
               ["radius_cap", "Ang",     20, [0, inf],    "volume", "Cap radius"],
               ["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, [-inf, inf], "orientation", "inclination angle"],
+              ["phi",        "degrees", 60, [-inf, inf], "orientation", "deflection angle"],
              ]
 # pylint: enable=bad-whitespace, line-too-long
Index: sasmodels/models/core_shell_bicelle.c
===================================================================
--- sasmodels/models/core_shell_bicelle.c	(revision 22221346d45e39a667fe231ce7ca5ecfc79b0ac5)
+++ sasmodels/models/core_shell_bicelle.c	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -117,6 +117,6 @@
 
     // Cylinder orientation
-    const double cyl_x = cos(theta) * cos(phi);
-    const double cyl_y = sin(theta);
+    const double cyl_x = sin(theta) * cos(phi);
+    const double cyl_y = sin(theta) * sin(phi);
 
     // Compute the angle btw vector q and the axis of the cylinder
Index: sasmodels/models/core_shell_bicelle.py
===================================================================
--- sasmodels/models/core_shell_bicelle.py	(revision 416f5c7f5af38d5e5d926a5b6b140a94e7e60502)
+++ sasmodels/models/core_shell_bicelle.py	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -71,8 +71,4 @@
     Definition of the angles for the oriented core shell bicelle tmodel.
 
-.. figure:: img/cylinder_angle_projection.jpg
-    :width: 600px
-
-    Examples of the angles for oriented pp against the detector plane.
 
 References
@@ -161,39 +157,2 @@
 
 qx, qy = 0.4 * cos(90), 0.5 * sin(0)
-tests = [
-    # Accuracy tests based on content in test/utest_other_models.py
-    [{'radius': 20.0,
-      'thick_rim': 10.0,
-      'thick_face': 10.0,
-      'length': 400.0,
-      'sld_core': 1.0,
-      'sld_face': 4.0,
-      'sld_rim': 4.0,
-      'sld_solvent': 1.0,
-      'background': 0.0,
-     }, 0.001, 353.550],
-
-    [{'radius': 20.0,
-      'thick_rim': 10.0,
-      'thick_face': 10.0,
-      'length': 400.0,
-      'sld_core': 1.0,
-      'sld_face': 4.0,
-      'sld_rim': 4.0,
-      'sld_solvent': 1.0,
-      'theta': 90.0,
-      'phi': 0.0,
-      'background': 0.00,
-     }, (qx, qy), 24.9167],
-
-    # Additional tests with larger range of parameters
-    [{'radius': 3.0,
-      'thick_rim': 100.0,
-      'thick_face': 100.0,
-      'length': 1200.0,
-      'sld_core': 5.0,
-      'sld_face': 41.0,
-      'sld_rim': 42.0,
-      'sld_solvent': 21.0,
-     }, 0.05, 1670.1828],
-    ]
Index: sasmodels/models/core_shell_cylinder.c
===================================================================
--- sasmodels/models/core_shell_cylinder.c	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
+++ sasmodels/models/core_shell_cylinder.c	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -69,10 +69,10 @@
 
     // Compute angle alpha between q and the cylinder axis
-    SINCOS(theta*M_PI_180, sn, cn);
+    SINCOS(phi*M_PI_180, sn, cn);
     // # The following correction factor exists in sasview, but it can't be
     // # right, so we are leaving it out for now.
     // const double correction = fabs(cn)*M_PI_2;
     const double q = sqrt(qx*qx+qy*qy);
-    const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);
+    const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q);
     const double alpha = acos(cos_val);
 
Index: sasmodels/models/core_shell_ellipsoid.c
===================================================================
--- sasmodels/models/core_shell_ellipsoid.c	(revision 5031ca39a761c97db704b9b8e3a0f211847bb0db)
+++ sasmodels/models/core_shell_ellipsoid.c	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -96,6 +96,6 @@
 
     // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis.
-    const double cyl_x = cos(theta) * cos(phi);
-    const double cyl_y = sin(theta);
+    const double cyl_x = sin(theta) * cos(phi);
+    const double cyl_y = sin(theta) * sin(phi);
 
     const double sldcs = core_sld - shell_sld;
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision 03cac087353e872f40154bf57ef8161eb5367fd4)
+++ sasmodels/models/cylinder.c	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -1,3 +1,5 @@
 double form_volume(double radius, double length);
+double fq(double q, double sn, double cn,double radius, double length);
+double orient_avg_1D(double q, double radius, double length);
 double Iq(double q, double sld, double solvent_sld, double radius, double length);
 double Iqxy(double qx, double qy, double sld, double solvent_sld,
@@ -11,4 +13,30 @@
 }
 
+double fq(double q, double sn, double cn, double radius, double length)
+{
+    // precompute qr and qh to save time in the loop
+    const double qr = q*radius;
+    const double qh = q*0.5*length; 
+    return  sas_J1c(qr*sn) * sinc(qh*cn) ;
+}
+
+double orient_avg_1D(double q, double radius, double length)
+{
+    // 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 total = 0.0;
+    for (int i=0; i<76 ;i++) {
+        const double alpha = Gauss76Z[i]*zm + zb;
+        double sn, cn; // slots to hold sincos function output
+        // alpha(theta,phi) the projection of the cylinder on the detector plane
+        SINCOS(alpha, sn, cn);
+        total += Gauss76Wt[i] * fq(q, sn, cn, radius, length)* fq(q, sn, cn, radius, length) * sn;
+    }
+    // translate dx in [-1,1] to dx in [lower,upper]
+    return total*zm;
+}
+
 double Iq(double q,
     double sld,
@@ -17,24 +45,6 @@
     double length)
 {
-    // precompute qr and qh to save time in the loop
-    const double qr = q*radius;
-    const double qh = q*0.5*length;
-
-    // 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 total = 0.0;
-    for (int i=0; i<76 ;i++) {
-        const double alpha = Gauss76Z[i]*zm + zb;
-        double sn, cn;
-        SINCOS(alpha, sn, cn);
-        const double fq = sinc(qh*cn) * sas_J1c(qr*sn);
-        total += Gauss76Wt[i] * fq*fq * sn;
-    }
-    // translate dx in [-1,1] to dx in [lower,upper]
-    const double form = total*zm;
     const double s = (sld - solvent_sld) * form_volume(radius, length);
-    return 1.0e-4 * s * s * form;
+    return 1.0e-4 * s * s * orient_avg_1D(q, radius, length);
 }
 
@@ -51,12 +61,12 @@
 
     // Compute angle alpha between q and the cylinder axis
-    SINCOS(theta*M_PI_180, sn, cn);
+    SINCOS(phi*M_PI_180, sn, cn);
     const double q = sqrt(qx*qx + qy*qy);
-    const double cos_val = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q);
+    const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q);
+
     const double alpha = acos(cos_val);
 
     SINCOS(alpha, sn, cn);
-    const double fq = sinc(q*0.5*length*cn) * sas_J1c(q*radius*sn);
     const double s = (sld-solvent_sld) * form_volume(radius, length);
-    return 1.0e-4 * square(s * fq);
+    return 1.0e-4 * square(s * fq(q, sn, cn, radius, length));
 }
Index: sasmodels/models/cylinder.py
===================================================================
--- sasmodels/models/cylinder.py	(revision 416f5c7f5af38d5e5d926a5b6b140a94e7e60502)
+++ sasmodels/models/cylinder.py	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -35,5 +35,5 @@
 .. math::
 
-    F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}
+    F^2(q)=\int_{0}^{\pi/2}{F^2(q,\theta)\sin(\theta)d\theta}
 
 
@@ -48,7 +48,4 @@
     Definition of the angles for oriented cylinders.
 
-.. figure:: img/cylinder_angle_projection.jpg
-
-    Examples of the angles for oriented cylinders against the detector plane.
 
 NB: The 2nd virial coefficient of the cylinder is calculated based on the
@@ -77,8 +74,8 @@
 
     P(q) = \int_0^{\pi/2} d\phi
-        \int_0^\pi p(\theta, \phi) P_0(q,\alpha) \sin \theta\ d\theta
+        \int_0^\pi p(\alpha) P_0(q,\alpha) \sin \alpha\ d\alpha
 
 
-where $p(\theta,\phi)$ is the probability distribution for the orientation
+where $p(\theta,\phi) = 1$ is the probability distribution for the orientation
 and $P_0(q,\alpha)$ is the scattering intensity for the fully oriented
 system, and then comparing to the 1D result.
@@ -87,6 +84,6 @@
 ----------
 
-None
-
+J. S. Pedersen, Adv. Colloid Interface Sci. 70, 171-210 (1997).
+G. Fournet, Bull. Soc. Fr. Mineral. Cristallogr. 74, 39-113 (1951).
 """
 
@@ -123,10 +120,10 @@
                "Cylinder length"],
               ["theta", "degrees", 60, [-inf, inf], "orientation",
-               "In plane angle"],
+               "latitude"],
               ["phi", "degrees", 60, [-inf, inf], "orientation",
-               "Out of plane angle"],
+               "longitude"],
              ]
 
-source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"]
+source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c",  "cylinder.c"]
 
 def ER(radius, length):
@@ -148,9 +145,10 @@
 
 qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5)
-tests = [[{}, 0.2, 0.042761386790780453],
-         [{}, [0.2], [0.042761386790780453]],
-         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852],
-         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]],
-        ]
+# After redefinition of angles, find new tests values 
+#tests = [[{}, 0.2, 0.042761386790780453],
+#         [{}, [0.2], [0.042761386790780453]],
+#         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852],
+#         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]],
+#        ]
 del qx, qy  # not necessary to delete, but cleaner
 # ADDED by:  RKH  ON: 18Mar2016 renamed sld's etc
Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/ellipsoid.c	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -52,8 +52,9 @@
 
     const double q = sqrt(qx*qx + qy*qy);
-    SINCOS(theta*M_PI_180, sn, cn);
-    // TODO: check if this is actually sin(alpha), not cos(alpha)
-    const double cos_alpha = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);
-    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha);
+    SINCOS(phi*M_PI_180, sn, cn);
+    const double cos_alpha = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q);
+    const double alpha = acos(cos_alpha);
+    SINCOS(alpha, sn, cn);
+    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, sn);
     const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial);
 
Index: sasmodels/models/ellipsoid.py
===================================================================
--- sasmodels/models/ellipsoid.py	(revision 416f5c7f5af38d5e5d926a5b6b140a94e7e60502)
+++ sasmodels/models/ellipsoid.py	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -57,9 +57,5 @@
 The $\theta$ and $\phi$ parameters are not used for the 1D output.
 
-.. _ellipsoid-geometry:
 
-.. figure:: img/ellipsoid_angle_projection.jpg
-
-    The angles for oriented ellipsoid, shown here as oblate, $a$ = $R_p$ and $b$ = $R_e$
 
 Validation
Index: sasmodels/models/hollow_cylinder.c
===================================================================
--- sasmodels/models/hollow_cylinder.c	(revision aea2e2ab3f11c42b7fc1a72a02b6784076c41d3a)
+++ sasmodels/models/hollow_cylinder.c	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -54,6 +54,6 @@
 
     // Cylinder orientation
-    cyl_x = cos(theta) * cos(phi);
-    cyl_y = sin(theta);
+    cyl_x = sin(theta) * cos(phi);
+    cyl_y = sin(theta) * sin(phi);
     //cyl_z = -cos(theta) * sin(phi);
 
Index: sasmodels/models/stacked_disks.c
===================================================================
--- sasmodels/models/stacked_disks.c	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/stacked_disks.c	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -142,6 +142,6 @@
 
     // parallelepiped orientation
-    const double cyl_x = ct * cp;
-    const double cyl_y = st;
+    const double cyl_x = st * cp;
+    const double cyl_y = st * sp;
 
     // Compute the angle btw vector q and the
Index: sasmodels/models/stacked_disks.py
===================================================================
--- sasmodels/models/stacked_disks.py	(revision 7c578619a7733b31ffab2ce827e4701af4643dea)
+++ sasmodels/models/stacked_disks.py	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
@@ -77,12 +77,9 @@
 the axis of the cylinder using two angles $\theta$ and $\varphi$.
 
-.. figure:: img/stacked_disks_angle_definition.jpg
-
-    Examples of the angles for oriented stacked disks against
+.. figure:: img/cylinder_angle_definition.jpg
+
+    Examples of the angles against
     the detector plane.
 
-.. figure:: img/stacked_disks_angle_projection.jpg
-
-    Examples of the angles for oriented pp against the detector plane.
 
 Our model uses the form factor calculations implemented in a c-library provided
