Index: sasmodels/models/core_shell_parallelepiped.c
===================================================================
--- sasmodels/models/core_shell_parallelepiped.c	(revision fc0b7aa68c08867563a9e1bbe169ec028bd4ed78)
+++ sasmodels/models/core_shell_parallelepiped.c	(revision 4493288f50c1e858772d5058d6ca027bfa72f50e)
@@ -43,34 +43,14 @@
     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c
     // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
-    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)
+    // Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)
+    // Code rewritten (PAK)
 
-    const double mu = 0.5 * q * length_b;
+    const double half_q = 0.5*q;
 
-    // Scale sides by B
-    const double a_over_b = length_a / length_b;
-    const double c_over_b = length_c / length_b;
+    const double tA = length_a + 2.0*thick_rim_a;
+    const double tB = length_b + 2.0*thick_rim_b;
+    const double tC = length_c + 2.0*thick_rim_c;
 
-    double tA_over_b = a_over_b + 2.0*thick_rim_a/length_b;
-    double tB_over_b = 1+ 2.0*thick_rim_b/length_b;
-    double tC_over_b = c_over_b + 2.0*thick_rim_c/length_b;
-
-    double Vin = length_a * length_b * length_c;
-#if OVERLAPPING
-    const double capA_area = length_b*length_c;
-    const double capB_area = (length_a+2.*thick_rim_a)*length_c;
-    const double capC_area = (length_a+2.*thick_rim_a)*(length_b+2.*thick_rim_b);
-#else
-    const double capA_area = length_b*length_c;
-    const double capB_area = length_a*length_c;
-    const double capC_area = length_a*length_b;
-#endif
-    const double Va = length_a * capA_area;
-    const double Vb = length_b * capB_area;
-    const double Vc = length_c * capC_area;
-    const double Vat = Va + 2.0 * thick_rim_a * capA_area;
-    const double Vbt = Vb + 2.0 * thick_rim_b * capB_area;
-    const double Vct = Vc + 2.0 * thick_rim_c * capC_area;
-
-    // Scale factors (note that drC is not used later)
+    // Scale factors
     const double dr0 = (core_sld-solvent_sld);
     const double drA = (arim_sld-solvent_sld);
@@ -81,30 +61,30 @@
     double outer_sum = 0; //initialize integral
     for( int i=0; i<76; i++) {
-        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
-        double mu_proj = mu * sqrt(1.0-sigma*sigma);
+        const double cos_alpha = 0.5 * ( Gauss76Z[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 = sas_sinx_x(mu * sigma * c_over_b);
-        const double siCt = sas_sinx_x(mu * sigma * tC_over_b);
+        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);
         double inner_sum = 0.0;
         for(int j=0; j<76; j++) {
-            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
-            double sin_uu, cos_uu;
-            SINCOS(M_PI_2*uu, sin_uu, cos_uu);
-            const double siA = sas_sinx_x(mu_proj * sin_uu * a_over_b);
-            const double siB = sas_sinx_x(mu_proj * cos_uu );
-            const double siAt = sas_sinx_x(mu_proj * sin_uu * tA_over_b);
-            const double siBt = sas_sinx_x(mu_proj * cos_uu * tB_over_b);
+            const double beta = 0.5 * ( Gauss76Z[j] + 1.0 );
+            double sin_beta, cos_beta;
+            SINCOS(M_PI_2*beta, 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);
+            const double siAt = tA * sas_sinx_x(tA * mu * sin_beta);
+            const double siBt = tB * sas_sinx_x(tB * mu * cos_beta);
 
 #if OVERLAPPING
-            const double f = dr0*Vin*siA*siB*siC
-                + drA*(Vat*siAt-Va*siA)*siB*siC
-                + drB*siAt*(Vbt*siBt-Vb*siB)*siC
-                + drC*siAt*siBt*(Vct*siCt-Vc*siC);
+            const double f = dr0*siA*siB*siC
+                + drA*(siAt-siA)*siB*siC
+                + drB*siAt*(siBt-siB)*siC
+                + drC*siAt*siBt*(siCt-siC);
 #else
-            const double f = dr0*Vin*siA*siB*siC
-                + drA*(Vat*siAt-Va*siA)*siB*siC
-                + drB*siA*(Vbt*siBt-Vb*siB)*siC
-                + drC*siA*siB*(Vct*siCt-Vc*siC);
+            const double f = dr0*siA*siB*siC
+                + drA*(siAt-siA)*siB*siC
+                + drB*siA*(siBt-siB)*siC
+                + drC*siA*siB*(siCt-siC);
 #endif
 
@@ -141,21 +121,4 @@
     const double drC = crim_sld-solvent_sld;
 
-    double Vin = length_a * length_b * length_c;
-#if OVERLAPPING
-    const double capA_area = length_b*length_c;
-    const double capB_area = (length_a+2.*thick_rim_a)*length_c;
-    const double capC_area = (length_a+2.*thick_rim_a)*(length_b+2.*thick_rim_b);
-#else
-    const double capA_area = length_b*length_c;
-    const double capB_area = length_a*length_c;
-    const double capC_area = length_a*length_b;
-#endif
-    const double Va = length_a * capA_area;
-    const double Vb = length_b * capB_area;
-    const double Vc = length_c * capC_area;
-    const double Vat = Va + 2.0 * thick_rim_a * capA_area;
-    const double Vbt = Vb + 2.0 * thick_rim_b * capB_area;
-    const double Vct = Vc + 2.0 * thick_rim_c * capC_area;
-
     // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
     // the scaling by B.
@@ -163,21 +126,21 @@
     const double tB = length_b + 2.0*thick_rim_b;
     const double tC = length_c + 2.0*thick_rim_c;
-    const double siA = sas_sinx_x(0.5*length_a*qa);
-    const double siB = sas_sinx_x(0.5*length_b*qb);
-    const double siC = sas_sinx_x(0.5*length_c*qc);
-    const double siAt = sas_sinx_x(0.5*tA*qa);
-    const double siBt = sas_sinx_x(0.5*tB*qb);
-    const double siCt = sas_sinx_x(0.5*tC*qc);
+    const double siA = length_a*sas_sinx_x(0.5*length_a*qa);
+    const double siB = length_b*sas_sinx_x(0.5*length_b*qb);
+    const double siC = length_c*sas_sinx_x(0.5*length_c*qc);
+    const double siAt = tA*sas_sinx_x(0.5*tA*qa);
+    const double siBt = tB*sas_sinx_x(0.5*tB*qb);
+    const double siCt = tC*sas_sinx_x(0.5*tC*qc);
 
 #if OVERLAPPING
-    const double f = dr0*Vin*siA*siB*siC
-        + drA*(Vat*siAt-Va*siA)*siB*siC
-        + drB*siAt*(Vbt*siBt-Vb*siB)*siC
-        + drC*siAt*siBt*(Vct*siCt-Vc*siC);
+    const double f = dr0*siA*siB*siC
+        + drA*(siAt-siA)*siB*siC
+        + drB*siAt*(siBt-siB)*siC
+        + drC*siAt*siBt*(siCt-siC);
 #else
-    const double f = dr0*Vin*siA*siB*siC
-        + drA*(Vat*siAt-Va*siA)*siB*siC
-        + drB*siA*(Vbt*siBt-Vb*siB)*siC
-        + drC*siA*siB*(Vct*siCt-Vc*siC);
+    const double f = dr0*siA*siB*siC
+        + drA*(siAt-siA)*siB*siC
+        + drB*siA*(siBt-siB)*siC
+        + drC*siA*siB*(siCt-siC);
 #endif
 
Index: sasmodels/models/core_shell_parallelepiped.py
===================================================================
--- sasmodels/models/core_shell_parallelepiped.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/core_shell_parallelepiped.py	(revision 4493288f50c1e858772d5058d6ca027bfa72f50e)
@@ -40,16 +40,22 @@
 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 longest axis $C$
-of the parallelepiped, $\beta$ is the angle between projection of the particle in the $xy$ detector plane and the $y$ axis.
-
-.. math::
-    \begin{align*}
-    F(Q)&=A B C (\rho_\text{core}-\rho_\text{solvent})  S(A \sin\alpha \sin\beta)S(B \sin\alpha \cos\beta)S(C \cos\alpha) \\
-    &+ 2t_A B C (\rho_\text{A}-\rho_\text{solvent})  \left[S((A+t_A) \sin\alpha \sin\beta)-S(A \sin\alpha \sin\beta)\right] S(B \sin\alpha \cos\beta) S(C \cos\alpha)\\
-    &+ 2 A t_B C (\rho_\text{B}-\rho_\text{solvent})  S(A \sin\alpha \sin\beta) \left[S((B+t_B) \sin\alpha \cos\beta)-S(B \sin\alpha \cos\beta)\right] S(C \cos\alpha)\\
-    &+ 2 A B t_C (\rho_\text{C}-\rho_\text{solvent}) S(A \sin\alpha \sin\beta) S(B \sin\alpha \cos\beta) \left[S((C+t_C) \cos\alpha)-S(C \cos\alpha)\right]
-    \end{align*}
+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)
+    &= (\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) \\
+    &+ (\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) \\
+    &+ (\rho_\text{C}-\rho_\text{solvent})
+        S(Q_A, A) S(Q_B, B) \left[S(Q_C, C+2t_C) - S(Q_C, C)\right]
 
 with
@@ -57,16 +63,26 @@
 .. math::
 
-    S(x) = \frac{\sin \tfrac{1}{2}Q x}{\tfrac{1}{2}Q x}
-
-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 slabs of
-thickness $t_A$, $t_B$ and $t_C$, respectively.
-$\rho_\text{solvent}$ is the scattering length of the solvent.
+    S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L}
+
+and
+
+.. math::
+
+    Q_A &= \sin\alpha \sin\beta \\
+    Q_B &= \sin\alpha \cos\beta \\
+    Q_C &= \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
+slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$
+is the scattering length of the solvent.
 
 FITTING NOTES
+~~~~~~~~~~~~~
+
 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
-calculation.**
+value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However,
+**no interparticle interference effects are included in this calculation.**
 
 There are many parameters in this model. Hold as many fixed as possible with
@@ -77,16 +93,17 @@
 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, 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.
+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.
 
 For 2d data the orientation of the particle is required, described using
-angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details
-of the calculation and angular dispersions see :ref:`orientation` .
+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,
+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.
 
@@ -173,9 +190,11 @@
 
     # surface average radius (rough approximation)
-    surf_rad = sqrt((length_a + 2.0*thick_rim_a) * (length_b + 2.0*thick_rim_b) / pi)
+    surf_rad = sqrt((length_a + 2.0*thick_rim_a)
+                    * (length_b + 2.0*thick_rim_b) / pi)
 
     height = length_c + 2.0*thick_rim_c
 
-    ddd = 0.75 * surf_rad * (2 * surf_rad * height + (height + surf_rad) * (height + pi * surf_rad))
+    ddd = (0.75 * surf_rad * (2 * surf_rad * height
+           + (height + surf_rad) * (height + pi * surf_rad)))
     return 0.5 * (ddd) ** (1. / 3.)
 
@@ -212,5 +231,6 @@
             psi_pd=10, psi_pd_n=1)
 
-# rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models
+# rkh 7/4/17 add random unit test for 2d, note make all params different,
+# 2d values not tested against other codes or models
 if 0:  # pak: model rewrite; need to update tests
     qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
