Index: sasmodels/models/core_shell_parallelepiped.c
===================================================================
--- sasmodels/models/core_shell_parallelepiped.c	(revision e0772316ffb087990c2e1b341c744594e3087097)
+++ sasmodels/models/core_shell_parallelepiped.c	(revision dbf1a60059f2daf9ef240db7a7fbf91cf5de9b9e)
@@ -59,17 +59,19 @@
 
     // outer integral (with gauss points), integration limits = 0, 1
+    // substitute d_cos_alpha for sin_alpha d_alpha
     double outer_sum = 0; //initialize integral
     for( int i=0; i<GAUSS_N; i++) {
         const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 );
         const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha);
-
-        // inner integral (with gauss points), integration limits = 0, pi/2
         const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q);
         const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q);
+
+        // inner integral (with gauss points), integration limits = 0, 1
+        // substitute beta = PI/2 u (so 2/PI * d_(PI/2 * beta) = d_beta)
         double inner_sum = 0.0;
         for(int j=0; j<GAUSS_N; j++) {
-            const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 );
+            const double u = 0.5 * ( GAUSS_Z[j] + 1.0 );
             double sin_beta, cos_beta;
-            SINCOS(M_PI_2*beta, sin_beta, cos_beta);
+            SINCOS(M_PI_2*u, sin_beta, cos_beta);
             const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta);
             const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta);
@@ -91,8 +93,10 @@
             inner_sum += GAUSS_W[j] * f * f;
         }
+        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5
         inner_sum *= 0.5;
         // now sum up the outer integral
         outer_sum += GAUSS_W[i] * inner_sum;
     }
+    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5
     outer_sum *= 0.5;
 
Index: sasmodels/models/core_shell_parallelepiped.py
===================================================================
--- sasmodels/models/core_shell_parallelepiped.py	(revision 97be87700f140bc9317b791ad4691cc31fa30981)
+++ sasmodels/models/core_shell_parallelepiped.py	(revision 5bc6d21cb8000cd9bb3d14f07d5f5de178c9b063)
@@ -11,47 +11,60 @@
 .. math::
 
-    I(q) = \text{scale}\frac{\langle f^2 \rangle}{V} + \text{background}
+    I(q) = \frac{\text{scale}}{V} \langle P(q,\alpha,\beta) \rangle 
+    + \text{background}
 
 where $\langle \ldots \rangle$ is an average over all possible orientations
-of the rectangular solid.
-
-The function calculated is the form factor of the rectangular solid below.
+of the rectangular solid, and the usual $\Delta \rho^2 \ V^2$ term cannot be
+pulled out of the form factor term due to the multiple slds in the model.
+
 The core of the solid is defined by the dimensions $A$, $B$, $C$ such that
 $A < B < C$.
 
-.. image:: img/core_shell_parallelepiped_geometry.jpg
+.. figure:: img/parallelepiped_geometry.jpg
+
+   Core of the core shell parallelepiped with the corresponding definition
+   of sides.
+
 
 There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension
 (on the $BC$ faces). There are similar slabs on the $AC$ $(=t_B)$ and $AB$
-$(=t_C)$ faces. The projection in the $AB$ plane is then
-
-.. image:: img/core_shell_parallelepiped_projection.jpg
-
-The volume of the solid is
+$(=t_C)$ faces. The projection in the $AB$ plane is
+
+.. figure:: img/core_shell_parallelepiped_projection.jpg
+
+   AB cut through the core-shell parallelipiped showing the cross secion of
+   four of the six shell slabs. As can be seen, this model leaves **"gaps"**
+   at the corners of the solid.
+
+
+The total volume of the solid is thus given as
 
 .. math::
 
     V = ABC + 2t_ABC + 2t_BAC + 2t_CAB
-
-**meaning that there are "gaps" at the corners of the solid.**
 
 The intensity calculated follows the :ref:`parallelepiped` model, with the
 core-shell intensity being calculated as the square of the sum of the
-amplitudes of the core and the slabs on the edges.
-
-the scattering amplitude is computed for a particular orientation of the
-core-shell parallelepiped with respect to the scattering vector and then
-averaged over all possible orientations, where $\alpha$ is the angle between
-the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is
-the angle between projection of the particle in the $xy$ detector plane
-and the $y$ axis.
-
-.. math::
-
-    F(Q)
+amplitudes of the core and the slabs on the edges. The scattering amplitude is
+computed for a particular orientation of the core-shell parallelepiped with
+respect to the scattering vector and then averaged over all possible
+orientations, where $\alpha$ is the angle between the $z$ axis and the $C$ axis
+of the parallelepiped, and $\beta$ is the angle between the projection of the
+particle in the $xy$ detector plane and the $y$ axis.
+
+.. math::
+
+    P(q)=\frac {\int_{0}^{\pi/2}\int_{0}^{\pi/2}F^2(q,\alpha,\beta) \ sin\alpha
+    \ d\alpha \ d\beta} {\int_{0}^{\pi/2} \ sin\alpha \ d\alpha \ d\beta}
+
+and
+
+.. math::
+
+    F(q,\alpha,\beta)
     &= (\rho_\text{core}-\rho_\text{solvent})
        S(Q_A, A) S(Q_B, B) S(Q_C, C) \\
     &+ (\rho_\text{A}-\rho_\text{solvent})
-        \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\
+        \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\
     &+ (\rho_\text{B}-\rho_\text{solvent})
         S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\
@@ -63,5 +76,5 @@
 .. math::
 
-    S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L}
+    S(Q_X, L) = L \frac{\sin (\tfrac{1}{2} Q_X L)}{\tfrac{1}{2} Q_X L}
 
 and
@@ -69,13 +82,22 @@
 .. math::
 
-    Q_A &= \sin\alpha \sin\beta \\
-    Q_B &= \sin\alpha \cos\beta \\
-    Q_C &= \cos\alpha
+    Q_A &= q \sin\alpha \sin\beta \\
+    Q_B &= q \sin\alpha \cos\beta \\
+    Q_C &= q \cos\alpha
 
 
 where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$
-are the scattering length of the parallelepiped core, and the rectangular
+are the scattering lengths of the parallelepiped core, and the rectangular
 slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$
 is the scattering length of the solvent.
+
+.. note:: 
+
+   the code actually implements two substitutions: $d(cos\alpha)$ is
+   substituted for -$sin\alpha \ d\alpha$ (note that in the
+   :ref:`parallelepiped` code this is explicitly implemented with
+   $\sigma = cos\alpha$), and $\beta$ is set to $\beta = u \pi/2$ so that
+   $du = \pi/2 \ d\beta$.  Thus both integrals go from 0 to 1 rather than 0
+   to $\pi/2$.
 
 FITTING NOTES
@@ -94,17 +116,17 @@
 based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$
 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 give an oblate or prolate particle, to give an effective radius
+for $S(q)$ when $P(q) * S(q)$ is applied.
 
 For 2d data the orientation of the particle is required, described using
-angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further
+angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below. For further
 details of the calculation and angular dispersions see :ref:`orientation`.
 The angle $\Psi$ is the rotational angle around the *long_c* axis. For example,
 $\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector.
 
-For 2d, constraints must be applied during fitting to ensure that the
-inequality $A < B < C$ is not violated, and hence the correct definition
-of angles is preserved. The calculation will not report an error,
-but the results may be not correct.
+.. note:: For 2d, constraints must be applied during fitting to ensure that the
+   inequality $A < B < C$ is not violated, and hence the correct definition
+   of angles is preserved. The calculation will not report an error,
+   but the results may be not correct.
 
 .. figure:: img/parallelepiped_angle_definition.png
@@ -113,5 +135,5 @@
     Note that rotation $\theta$, initially in the $xz$ plane, is carried
     out first, then rotation $\phi$ about the $z$ axis, finally rotation
-    $\Psi$ is now around the axis of the cylinder. The neutron or X-ray
+    $\Psi$ is now around the axis of the particle. The neutron or X-ray
     beam is along the $z$ axis.
 
@@ -120,4 +142,13 @@
     Examples of the angles for oriented core-shell parallelepipeds against the
     detector plane.
+
+
+Validation
+----------
+
+Cross-checked against hollow rectangular prism and rectangular prism for equal
+thickness overlapping sides, and by Monte Carlo sampling of points within the
+shape for non-uniform, non-overlapping sides.
+
 
 References
@@ -135,9 +166,6 @@
 
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
-* **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016
+* **Converted to sasmodels by:** Miguel Gonzalez **Date:** February 26, 2016
 * **Last Modified by:** Paul Kienzle **Date:** October 17, 2017
-* Cross-checked against hollow rectangular prism and rectangular prism for
-  equal thickness overlapping sides, and by Monte Carlo sampling of points
-  within the shape for non-uniform, non-overlapping sides.
 """
 
Index: sasmodels/models/parallelepiped.c
===================================================================
--- sasmodels/models/parallelepiped.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
+++ sasmodels/models/parallelepiped.c	(revision dbf1a60059f2daf9ef240db7a7fbf91cf5de9b9e)
@@ -38,4 +38,5 @@
             inner_total += GAUSS_W[j] * square(si1 * si2);
         }
+        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5
         inner_total *= 0.5;
 
@@ -43,4 +44,5 @@
         outer_total += GAUSS_W[i] * inner_total * si * si;
     }
+    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5
     outer_total *= 0.5;
 
Index: sasmodels/models/parallelepiped.py
===================================================================
--- sasmodels/models/parallelepiped.py	(revision ef07e950965796e2a3a8156fcade98ce5e403001)
+++ sasmodels/models/parallelepiped.py	(revision b343226d68c6e3af776c858761eacbaf49e0dd8f)
@@ -2,14 +2,12 @@
 # Note: model title and parameter table are inserted automatically
 r"""
-The form factor is normalized by the particle volume.
-For information about polarised and magnetic scattering, see
-the :ref:`magnetism` documentation.
-
 Definition
 ----------
 
- 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`. For
+information about polarised and magnetic scattering, see
+the :ref:`magnetism` documentation.
 
 .. _parallelepiped-image:
@@ -26,5 +24,6 @@
 error, or fixing of some dimensions at expected values, may help.
 
-The 1D scattering intensity $I(q)$ is calculated as:
+The form factor is normalized by the particle volume and the 1D scattering
+intensity $I(q)$ is then calculated as:
 
 .. Comment by Miguel Gonzalez:
@@ -39,14 +38,16 @@
 
     I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2
-           \left< P(q, \alpha) \right> + \text{background}
+           \left< P(q, \alpha, \beta) \right> + \text{background}
 
 where the volume $V = A B C$, the contrast is defined as
-$\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$,
-$P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented
-at an angle $\alpha$ (angle between the long axis C and $\vec q$),
-and the averaging $\left<\ldots\right>$ is applied over all orientations.
+$\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, $P(q, \alpha, \beta)$
+is the form factor corresponding to a parallelepiped oriented
+at an angle $\alpha$ (angle between the long axis C and $\vec q$), and $\beta$
+(the angle between the projection of the particle in the $xy$ detector plane
+and the $y$ axis) and the averaging $\left<\ldots\right>$ is applied over all
+orientations.
 
 Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the
-form factor is given by (Mittelbach and Porod, 1961)
+form factor is given by (Mittelbach and Porod, 1961 [#Mittelbach]_)
 
 .. math::
@@ -66,5 +67,6 @@
     \mu &= qB
 
-The scattering intensity per unit volume is returned in units of |cm^-1|.
+where substitution of $\sigma = cos\alpha$ and $\beta = \pi/2 \ u$ have been
+applied.
 
 NB: The 2nd virial coefficient of the parallelepiped is calculated based on
@@ -120,7 +122,10 @@
 .. math::
 
-    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
+    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
@@ -160,7 +165,7 @@
 ----------
 
-P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
-
-R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
+.. [#Mittelbach] P Mittelbach and G Porod, *Acta Physica Austriaca*,
+   14 (1961) 185-211
+.. [#] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
 
 Authorship and Verification
