Index: explore/asymint.py
===================================================================
--- explore/asymint.py	(revision 1820208ff3a9fba87725a7037153e1cf974a72a9)
+++ explore/asymint.py	(revision 49eb251955404dc1006cd29ad5a1b9f69c53343e)
@@ -86,11 +86,37 @@
     a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
     def Fq(qa, qb, qc):
-        siA = env.sas_sinx_x(0.5*a*qa/2)
-        siB = env.sas_sinx_x(0.5*b*qb/2)
-        siC = env.sas_sinx_x(0.5*c*qc/2)
+        siA = env.sas_sinx_x(a*qa/2)
+        siB = env.sas_sinx_x(b*qb/2)
+        siC = env.sas_sinx_x(c*qc/2)
         return siA * siB * siC
     Fq.__doc__ = "parallelepiped a=%g, b=%g c=%g"%(a, b, c)
     volume = a*b*c
     norm = CONTRAST**2*volume/10000
+    return norm, Fq
+
+def make_core_shell_parallelepiped(a, b, c, da, db, dc, slda, sldb, sldc, env=NPenv):
+    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
+    da, db, dc = env.mpf(da), env.mpf(db), env.mpf(dc)
+    slda, sldb, sldc = env.mpf(slda), env.mpf(sldb), env.mpf(sldc)
+    drV0 = CONTRAST*a*b*c
+    dra, drb, drc = slda-SLD_SOLVENT, sldb-SLD_SOLVENT, sldc-SLD_SOLVENT
+    Aa, Ab, Ac = b*c, a*c, a*b
+    Ta, Tb, Tc = a + 2*da, b + 2*db, c + 2*dc
+    drVa, drVb, drVc = dra*a*Aa, drb*b*Ab, drc*c*Ac
+    drVTa, drVTb, drVTc = dra*Ta*Aa, drb*Tb*Ab, drc*Tc*Ac
+    def Fq(qa, qb, qc):
+        siA = env.sas_sinx_x(a*qa/2)
+        siB = env.sas_sinx_x(b*qb/2)
+        siC = env.sas_sinx_x(c*qc/2)
+        siAt = env.sas_sinx_x(Ta*qa/2)
+        siBt = env.sas_sinx_x(Tb*qb/2)
+        siCt = env.sas_sinx_x(Tc*qc/2)
+        return (drV0*siA*siB*siC
+            + (drVTa*siAt-drVa*siA)*siB*siC
+            + siA*(drVTb*siBt-drVb*siB)*siC
+            + siA*siB*(drVTc*siCt-drVc*siC))
+    Fq.__doc__ = "core-shell parallelepiped a=%g, b=%g c=%g"%(a, b, c)
+    volume = a*b*c + 2*da*Aa + 2*db*Ab + 2*dc*Ac
+    norm = 1/(volume*10000)
     return norm, Fq
 
@@ -184,4 +210,20 @@
     NORM, KERNEL = make_parallelepiped(A, B, C)
     NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv)
+elif shape == 'core_shell_parallelepiped':
+    #A, B, C = 4450, 14000, 47
+    #A, B, C = 445, 140, 47  # integer for the sake of mpf
+    A, B, C = 6800, 114, 1380
+    DA, DB, DC = 2300, 21, 58
+    SLDA, SLDB, SLDC = "5", "-0.3", "11.5"
+    if 1: # C shortest
+        B, C = C, B
+        DB, DC = DC, DB
+        SLDB, SLDC = SLDC, SLDB
+    elif 0: # C longest
+        A, C = C, A
+        DA, DC = DC, DA
+        SLDA, SLDC = SLDC, SLDA
+    NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC)
+    NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv)
 elif shape == 'paracrystal':
     LATTICE = 'bcc'
@@ -342,9 +384,12 @@
     print("gauss-150", *gauss_quad_2d(Q, n=150))
     print("gauss-500", *gauss_quad_2d(Q, n=500))
+    print("gauss-1025", *gauss_quad_2d(Q, n=1025))
+    print("gauss-2049", *gauss_quad_2d(Q, n=2049))
     #gridded_2d(Q, n=2**8+1)
     gridded_2d(Q, n=2**10+1)
-    #gridded_2d(Q, n=2**13+1)
+    #gridded_2d(Q, n=2**12+1)
     #gridded_2d(Q, n=2**15+1)
-    if shape != 'paracrystal':  # adaptive forms are too slow!
+    if shape not in ('paracrystal', 'core_shell_parallelepiped'):
+        # adaptive forms on models for which the calculations are fast enough
         print("dblquad", *scipy_dblquad_2d(Q))
         print("semi-romberg-100", *semi_romberg_2d(Q, n=100))
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/compare.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
@@ -693,5 +693,5 @@
         data = empty_data2D(q, resolution=res)
         data.accuracy = opts['accuracy']
-        set_beam_stop(data, 0.0004)
+        set_beam_stop(data, qmin)
         index = ~data.mask
     else:
@@ -1274,5 +1274,5 @@
         if model_info != model_info2:
             pars2 = randomize_pars(model_info2, pars2)
-            limit_dimensions(model_info, pars2, maxdim)
+            limit_dimensions(model_info2, pars2, maxdim)
             # Share values for parameters with the same name
             for k, v in pars.items():
Index: sasmodels/models/core_shell_parallelepiped.c
===================================================================
--- sasmodels/models/core_shell_parallelepiped.c	(revision c69d6d6b7d1b4a8e1acc9c734f6ec3b26084cdd7)
+++ sasmodels/models/core_shell_parallelepiped.c	(revision 4493288f50c1e858772d5058d6ca027bfa72f50e)
@@ -1,11 +1,28 @@
+// Set OVERLAPPING to 1 in order to fill in the edges of the box, with
+// c endcaps and b overlapping a.  With the proper choice of parameters,
+// (setting rim slds to sld, core sld to solvent, rim thickness to thickness
+// and subtracting 2*thickness from length, this should match the hollow
+// rectangular prism.)  Set it to 0 for the documented behaviour.
+#define OVERLAPPING 0
 static double
 form_volume(double length_a, double length_b, double length_c,
     double thick_rim_a, double thick_rim_b, double thick_rim_c)
 {
-    //return length_a * length_b * length_c;
-    return length_a * length_b * length_c +
-           2.0 * thick_rim_a * length_b * length_c +
-           2.0 * thick_rim_b * length_a * length_c +
-           2.0 * thick_rim_c * length_a * length_b;
+    return
+#if OVERLAPPING
+        // Hollow rectangular prism only includes the volume of the shell
+        // so uncomment the next line when comparing.  Solid rectangular
+        // prism, or parallelepiped want filled cores, so comment when
+        // comparing.
+        //-length_a * length_b * length_c +
+        (length_a + 2.0*thick_rim_a) *
+        (length_b + 2.0*thick_rim_b) *
+        (length_c + 2.0*thick_rim_c);
+#else
+        length_a * length_b * length_c +
+        2.0 * thick_rim_a * length_b * length_c +
+        2.0 * length_a * thick_rim_b * length_c +
+        2.0 * length_a * length_b * thick_rim_c;
+#endif
 }
 
@@ -24,86 +41,62 @@
     double thick_rim_c)
 {
-    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled
+    // 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;
 
-    //calculate volume before rescaling (in original code, but not used)
-    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
-    //double vol = length_a * length_b * length_c +
-    //       2.0 * thick_rim_a * length_b * length_c +
-    //       2.0 * thick_rim_b * length_a * length_c +
-    //       2.0 * thick_rim_c * length_a * 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;
 
-    // Scale sides by B
-    const double a_scaled = length_a / length_b;
-    const double c_scaled = length_c / length_b;
-
-    double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b;
-    double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b;
-    double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present
-
-    double Vin = length_a * length_b * length_c;
-    //double Vot = (length_a * length_b * length_c +
-    //            2.0 * thick_rim_a * length_b * length_c +
-    //            2.0 * length_a * thick_rim_b * length_c +
-    //            2.0 * length_a * length_b * thick_rim_c);
-    double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc)
-    double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
-    double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present
-
-    // Scale factors (note that drC is not used later)
-    const double drho0 = (core_sld-solvent_sld);
-    const double drhoA = (arim_sld-solvent_sld);
-    const double drhoB = (brim_sld-solvent_sld);
-    const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot;
-
-
-    // Precompute scale factors for combining cross terms from the shape
-    const double scale23 = drhoA*V1;
-    const double scale14 = drhoB*V2;
-    const double scale24 = drhoC*V3;
-    const double scale11 = drho0*Vin;
-    const double scale12 = drho0*Vin - scale23 - scale14 - scale24;
+    // Scale factors
+    const double dr0 = (core_sld-solvent_sld);
+    const double drA = (arim_sld-solvent_sld);
+    const double drB = (brim_sld-solvent_sld);
+    const double drC = (crim_sld-solvent_sld);
 
     // outer integral (with gauss points), integration limits = 0, 1
-    double outer_total = 0; //initialize integral
+    double outer_sum = 0; //initialize integral
+    for( int i=0; i<76; i++) {
+        const double cos_alpha = 0.5 * ( Gauss76Z[i] + 1.0 );
+        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha);
 
-    for( int i=0; i<76; i++) {
-        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
-        double mu_proj = mu * sqrt(1.0-sigma*sigma);
+        // 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);
+        double inner_sum = 0.0;
+        for(int j=0; j<76; j++) {
+            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);
 
-        // inner integral (with gauss points), integration limits = 0, 1
-        double inner_total = 0.0;
-        double inner_total_crim = 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 si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
-            const double si2 = sas_sinx_x(mu_proj * cos_uu );
-            const double si3 = sas_sinx_x(mu_proj * sin_uu * ta);
-            const double si4 = sas_sinx_x(mu_proj * cos_uu * tb);
+#if OVERLAPPING
+            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*siA*siB*siC
+                + drA*(siAt-siA)*siB*siC
+                + drB*siA*(siBt-siB)*siC
+                + drC*siA*siB*(siCt-siC);
+#endif
 
-            // Expression in libCylinder.c (neither drC nor Vot are used)
-            const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4;
-            const double form_crim = scale11*si1*si2;
-
-            //  correct FF : sum of square of phase factors
-            inner_total += Gauss76Wt[j] * form * form;
-            inner_total_crim += Gauss76Wt[j] * form_crim * form_crim;
+            inner_sum += Gauss76Wt[j] * f * f;
         }
-        inner_total *= 0.5;
-        inner_total_crim *= 0.5;
+        inner_sum *= 0.5;
         // now sum up the outer integral
-        const double si = sas_sinx_x(mu * c_scaled * sigma);
-        const double si_crim = sas_sinx_x(mu * tc * sigma);
-        outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim);
+        outer_sum += Gauss76Wt[i] * inner_sum;
     }
-    outer_total *= 0.5;
+    outer_sum *= 0.5;
 
     //convert from [1e-12 A-1] to [cm-1]
-    return 1.0e-4 * outer_total;
+    return 1.0e-4 * outer_sum;
 }
 
@@ -128,34 +121,27 @@
     const double drC = crim_sld-solvent_sld;
 
-    double Vin = length_a * length_b * length_c;
-    double V1 = 2.0 * thick_rim_a * length_b * length_c;    // incorrect V1(aa*bb*cc+2*ta*bb*cc)
-    double V2 = 2.0 * length_a * thick_rim_b * length_c;    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
-    double V3 = 2.0 * length_a * length_b * thick_rim_c;
-    // As for the 1D case, Vot is not used
-    //double Vot = (length_a * length_b * length_c +
-    //              2.0 * thick_rim_a * length_b * length_c +
-    //              2.0 * length_a * thick_rim_b * length_c +
-    //              2.0 * length_a * length_b * thick_rim_c);
-
     // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
     // the scaling by B.
-    double ta = length_a + 2.0*thick_rim_a;
-    double tb = length_b + 2.0*thick_rim_b;
-    double tc = length_c + 2.0*thick_rim_c;
-    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
-    double siA = sas_sinx_x(0.5*length_a*qa);
-    double siB = sas_sinx_x(0.5*length_b*qb);
-    double siC = sas_sinx_x(0.5*length_c*qc);
-    double siAt = sas_sinx_x(0.5*ta*qa);
-    double siBt = sas_sinx_x(0.5*tb*qb);
-    double siCt = sas_sinx_x(0.5*tc*qc);
+    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;
+    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);
 
-
-    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed
-    // in the 1D code, but should be checked!
-    double f = ( dr0*siA*siB*siC*Vin
-               + drA*(siAt-siA)*siB*siC*V1
-               + drB*siA*(siBt-siB)*siC*V2
-               + drC*siA*siB*(siCt-siC)*V3);
+#if OVERLAPPING
+    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*siA*siB*siC
+        + drA*(siAt-siA)*siB*siC
+        + drB*siA*(siBt-siB)*siC
+        + drC*siA*siB*(siCt-siC);
+#endif
 
     return 1.0e-4 * f * f;
Index: sasmodels/models/core_shell_parallelepiped.py
===================================================================
--- sasmodels/models/core_shell_parallelepiped.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/core_shell_parallelepiped.py	(revision 10ee8389f2d5511546ae78ffddc6d16e36d2a6a8)
@@ -5,11 +5,5 @@
 Calculates the form factor for a rectangular solid with a core-shell structure.
 The thickness and the scattering length density of the shell or
-"rim" can be different on each (pair) of faces. However at this time the 1D
-calculation does **NOT** actually calculate a c face rim despite the presence
-of the parameter. Some other aspects of the 1D calculation may be wrong.
-
-.. note::
-   This model was originally ported from NIST IGOR macros. However, it is not
-   yet fully understood by the SasView developers and is currently under review.
+"rim" can be different on each (pair) of faces.
 
 The form factor is normalized by the particle volume $V$ such that
@@ -21,5 +15,4 @@
 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.
@@ -41,57 +34,77 @@
     V = ABC + 2t_ABC + 2t_BAC + 2t_CAB
 
-**meaning that there are "gaps" at the corners of the solid.**  Again note that
-$t_C = 0$ currently.
+**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 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 (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.
+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)
+    &= (\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
+
+.. math::
+
+    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
 known values, or you will certainly end up at a solution that is unphysical.
 
-Constraints must be applied during fitting to ensure that the inequality
-$A < B < C$ is not violated. The calculation will not report an error,
-but the results will not be correct.
-
 The returned value is in units of |cm^-1|, on absolute scale.
 
 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,
+but the results may be not correct.
 
 .. figure:: img/parallelepiped_angle_definition.png
@@ -114,5 +127,5 @@
     Equations (1), (13-14). (in German)
 .. [#] D Singh (2009). *Small angle scattering studies of self assembly in
-   lipid mixtures*, John's Hopkins University Thesis (2009) 223-225. `Available
+   lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available
    from Proquest <http://search.proquest.com/docview/304915826?accountid
    =26379>`_
@@ -175,12 +188,10 @@
         Return equivalent radius (ER)
     """
-
-    # surface average radius (rough approximation)
-    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))
-    return 0.5 * (ddd) ** (1. / 3.)
+    from .parallelepiped import ER as ER_p
+
+    a = length_a + 2*thick_rim_a
+    b = length_b + 2*thick_rim_b
+    c = length_c + 2*thick_rim_c
+    return ER_p(a, b, c)
 
 # VR defaults to 1.0
@@ -216,5 +227,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.)
Index: sasmodels/models/hollow_rectangular_prism.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism.c	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
+++ sasmodels/models/hollow_rectangular_prism.c	(revision 8de147789a271b990dd7ffc4244215e5def052c8)
@@ -1,4 +1,4 @@
 double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness);
-double Iq(double q, double sld, double solvent_sld, double length_a, 
+double Iq(double q, double sld, double solvent_sld, double length_a,
           double b2a_ratio, double c2a_ratio, double thickness);
 
@@ -37,5 +37,5 @@
     const double v2a = 0.0;
     const double v2b = M_PI_2;  //phi integration limits
-    
+
     double outer_sum = 0.0;
     for(int i=0; i<76; i++) {
@@ -84,2 +84,39 @@
     return 1.0e-4 * delrho * delrho * form;
 }
+
+double Iqxy(double qa, double qb, double qc,
+    double sld,
+    double solvent_sld,
+    double length_a,
+    double b2a_ratio,
+    double c2a_ratio,
+    double thickness)
+{
+    const double length_b = length_a * b2a_ratio;
+    const double length_c = length_a * c2a_ratio;
+    const double a_half = 0.5 * length_a;
+    const double b_half = 0.5 * length_b;
+    const double c_half = 0.5 * length_c;
+    const double vol_total = length_a * length_b * length_c;
+    const double vol_core = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness);
+
+    // Amplitude AP from eqn. (13)
+
+    const double termA1 = sas_sinx_x(qa * a_half);
+    const double termA2 = sas_sinx_x(qa * (a_half-thickness));
+
+    const double termB1 = sas_sinx_x(qb * b_half);
+    const double termB2 = sas_sinx_x(qb * (b_half-thickness));
+
+    const double termC1 = sas_sinx_x(qc * c_half);
+    const double termC2 = sas_sinx_x(qc * (c_half-thickness));
+
+    const double AP1 = vol_total * termA1 * termB1 * termC1;
+    const double AP2 = vol_core * termA2 * termB2 * termC2;
+
+    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization.
+    const double delrho = sld - solvent_sld;
+
+    // Convert from [1e-12 A-1] to [cm-1]
+    return 1.0e-4 * square(delrho * (AP1-AP2));
+}
Index: sasmodels/models/hollow_rectangular_prism.py
===================================================================
--- sasmodels/models/hollow_rectangular_prism.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/hollow_rectangular_prism.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
@@ -5,5 +5,5 @@
 This model provides the form factor, $P(q)$, for a hollow rectangular
 parallelepiped with a wall of thickness $\Delta$.
-It computes only the 1D scattering, not the 2D.
+
 
 Definition
@@ -66,5 +66,25 @@
 (which is unitless).
 
-**The 2D scattering intensity is not computed by this model.**
+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` .
+The angle $\Psi$ is the rotational angle around the long *C* axis. For example,
+$\Psi = 0$ when the *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.
+
+.. figure:: img/parallelepiped_angle_definition.png
+
+    Definition of the angles for oriented hollow rectangular prism.
+    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 prism.
+    The neutron or X-ray beam is along the $z$ axis.
+
+.. figure:: img/parallelepiped_angle_projection.png
+
+    Examples of the angles for oriented hollow rectangular prisms against the
+    detector plane.
 
 
@@ -113,4 +133,10 @@
               ["thickness", "Ang", 1, [0, inf], "volume",
                "Thickness of parallelepiped"],
+              ["theta", "degrees", 0, [-360, 360], "orientation",
+               "c axis to beam angle"],
+              ["phi", "degrees", 0, [-360, 360], "orientation",
+               "rotation about beam"],
+              ["psi", "degrees", 0, [-360, 360], "orientation",
+               "rotation about c axis"],
              ]
 
Index: sasmodels/models/rectangular_prism.c
===================================================================
--- sasmodels/models/rectangular_prism.c	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
+++ sasmodels/models/rectangular_prism.c	(revision 8de147789a271b990dd7ffc4244215e5def052c8)
@@ -1,4 +1,4 @@
 double form_volume(double length_a, double b2a_ratio, double c2a_ratio);
-double Iq(double q, double sld, double solvent_sld, double length_a, 
+double Iq(double q, double sld, double solvent_sld, double length_a,
           double b2a_ratio, double c2a_ratio);
 
@@ -26,5 +26,5 @@
     const double v2a = 0.0;
     const double v2b = M_PI_2;  //phi integration limits
-    
+
     double outer_sum = 0.0;
     for(int i=0; i<76; i++) {
@@ -53,8 +53,8 @@
     double answer = 0.5*(v1b-v1a)*outer_sum;
 
-    // Normalize by Pi (Eqn. 16). 
-    // The term (ABC)^2 does not appear because it was introduced before on 
+    // Normalize by Pi (Eqn. 16).
+    // The term (ABC)^2 does not appear because it was introduced before on
     // the definitions of termA, termB, termC.
-    // The factor 2 appears because the theta integral has been defined between 
+    // The factor 2 appears because the theta integral has been defined between
     // 0 and pi/2, instead of 0 to pi.
     answer /= M_PI_2; //Form factor P(q)
@@ -64,7 +64,37 @@
     answer *= square((sld-solvent_sld)*volume);
 
-    // Convert from [1e-12 A-1] to [cm-1] 
+    // Convert from [1e-12 A-1] to [cm-1]
     answer *= 1.0e-4;
 
     return answer;
 }
+
+
+double Iqxy(double qa, double qb, double qc,
+    double sld,
+    double solvent_sld,
+    double length_a,
+    double b2a_ratio,
+    double c2a_ratio)
+{
+    const double length_b = length_a * b2a_ratio;
+    const double length_c = length_a * c2a_ratio;
+    const double a_half = 0.5 * length_a;
+    const double b_half = 0.5 * length_b;
+    const double c_half = 0.5 * length_c;
+    const double volume = length_a * length_b * length_c;
+
+    // Amplitude AP from eqn. (13)
+
+    const double termA = sas_sinx_x(qa * a_half);
+    const double termB = sas_sinx_x(qb * b_half);
+    const double termC = sas_sinx_x(qc * c_half);
+
+    const double AP = termA * termB * termC;
+
+    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization.
+    const double delrho = sld - solvent_sld;
+
+    // Convert from [1e-12 A-1] to [cm-1]
+    return 1.0e-4 * square(volume * delrho * AP);
+}
Index: sasmodels/models/rectangular_prism.py
===================================================================
--- sasmodels/models/rectangular_prism.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/rectangular_prism.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
@@ -12,6 +12,4 @@
 the prism (e.g. setting $b/a = 1$ and $c/a = 1$ and applying polydispersity
 to *a* will generate a distribution of cubes of different sizes).
-Note also that, contrary to :ref:`parallelepiped`, it does not compute
-the 2D scattering.
 
 
@@ -26,6 +24,5 @@
 that reference), with $\theta$ corresponding to $\alpha$ in that paper,
 and not to the usual convention used for example in the
-:ref:`parallelepiped` model. As the present model does not compute
-the 2D scattering, this has no further consequences.
+:ref:`parallelepiped` model.
 
 In this model the scattering from a massive parallelepiped with an
@@ -65,5 +62,26 @@
 units) *scale* represents the volume fraction (which is unitless).
 
-**The 2D scattering intensity is not computed by this model.**
+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` .
+The angle $\Psi$ is the rotational angle around the long *C* axis. For example,
+$\Psi = 0$ when the *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.
+
+.. figure:: img/parallelepiped_angle_definition.png
+
+    Definition of the angles for oriented core-shell parallelepipeds.
+    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 beam is along the $z$ axis.
+
+.. figure:: img/parallelepiped_angle_projection.png
+
+    Examples of the angles for oriented rectangular prisms against the
+    detector plane.
+
 
 
@@ -108,4 +126,10 @@
               ["c2a_ratio", "", 1, [0, inf], "volume",
                "Ratio sides c/a"],
+              ["theta", "degrees", 0, [-360, 360], "orientation",
+               "c axis to beam angle"],
+              ["phi", "degrees", 0, [-360, 360], "orientation",
+               "rotation about beam"],
+              ["psi", "degrees", 0, [-360, 360], "orientation",
+               "rotation about c axis"],
              ]
 
