Index: sasmodels/convert.py
===================================================================
--- sasmodels/convert.py	(revision 34d6cabda48bc12af0423cded9c76c0e4bf2572b)
+++ sasmodels/convert.py	(revision 44bd2be361a937362169e016a821f398577e4661)
@@ -132,4 +132,8 @@
         _remove_pd(oldpars, 'num_pearls', name)
         _remove_pd(oldpars, 'thick_string', name)
+    elif name == 'core_shell_parallelepiped':
+        _remove_pd(oldpars, 'rimA', name)
+        _remove_pd(oldpars, 'rimB', name)
+        _remove_pd(oldpars, 'rimC', name)
     elif name == 'rpa':
         # convert scattering lengths from femtometers to centimeters
Index: sasmodels/models/core_shell_parallelepiped.c
===================================================================
--- sasmodels/models/core_shell_parallelepiped.c	(revision 44bd2be361a937362169e016a821f398577e4661)
+++ sasmodels/models/core_shell_parallelepiped.c	(revision 44bd2be361a937362169e016a821f398577e4661)
@@ -0,0 +1,281 @@
+double form_volume(double a_side, double b_side, double c_side, 
+                   double arim_thickness, double brim_thickness, double crim_thickness);
+double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld,
+          double solvent_sld, double a_side, double b_side, double c_side,
+          double arim_thickness, double brim_thickness, double crim_thickness);
+double Iqxy(double qx, double qy, double core_sld, double arim_sld, double brim_sld,
+            double crim_sld, double solvent_sld, double a_side, double b_side,
+            double c_side, double arim_thickness, double brim_thickness,
+            double crim_thickness, double theta, double phi, double psi);
+
+double form_volume(double a_side, double b_side, double c_side, 
+                   double arim_thickness, double brim_thickness, double crim_thickness)
+{
+    //return a_side * b_side * c_side;
+    return a_side * b_side * c_side + 
+           2.0 * arim_thickness * b_side * c_side + 
+           2.0 * brim_thickness * a_side * c_side +
+           2.0 * crim_thickness * a_side * b_side;
+}
+
+double Iq(double q,
+    double core_sld,
+    double arim_sld,
+    double brim_sld,
+    double crim_sld,
+    double solvent_sld,
+    double a_side,
+    double b_side,
+    double c_side,
+    double arim_thickness,
+    double brim_thickness,
+    double crim_thickness)
+{
+    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled
+    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
+    
+    double t1, t2, t3, t4, tmp, answer;   
+    double mu = q * b_side;
+    
+    //calculate volume before rescaling (in original code, but not used)
+    //double vol = form_volume(a_side, b_side, c_side, arim_thickness, brim_thickness, crim_thickness);		
+    //double vol = a_side * b_side * c_side + 
+    //       2.0 * arim_thickness * b_side * c_side + 
+    //       2.0 * brim_thickness * a_side * c_side +
+    //       2.0 * crim_thickness * a_side * b_side;
+    
+    // Scale sides by B
+    double a_scaled = a_side / b_side;
+    double c_scaled = c_side / b_side;
+    double arim_scaled = arim_thickness / b_side;
+    double brim_scaled = brim_thickness / b_side;
+       
+    // DelRho values (note that drC is not used later)       
+	double dr0 = core_sld-solvent_sld;
+	double drA = arim_sld-solvent_sld;
+	double drB = brim_sld-solvent_sld;
+	//double drC = crim_sld-solvent_sld;
+
+    //Order of integration
+    int nordi=76;			        
+    int nordj=76;
+    
+    // outer integral (with nordi gauss points), integration limits = 0, 1
+    double summ = 0; //initialize integral
+
+    for( int i=0; i<nordi; i++) {
+		
+        // inner integral (with nordj gauss points), integration limits = 0, 1
+	
+        double summj = 0.0;
+	    double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );		
+		
+	    for(int j=0; j<nordj; j++) {
+
+            double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
+            double mudum = mu * sqrt(1.0-sigma*sigma);
+
+	        double Vin = a_side * b_side * c_side;
+	        //double Vot = (a_side * b_side * c_side +
+            //            2.0 * arim_thickness * b_side * c_side +
+            //            2.0 * a_side * brim_thickness * c_side +
+            //            2.0 * a_side * b_side * crim_thickness);
+	        double V1 = (2.0 * arim_thickness * b_side * c_side);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc)
+	        double V2 = (2.0 * a_side * brim_thickness * c_side);    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
+	
+            // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG)
+            // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c,
+            // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions
+            // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!)            
+            double ta = (a_scaled+2.0*arim_thickness)/b_side; 
+            double tb = (a_scaled+2.0*brim_thickness)/b_side;
+    
+	        double arg1 = (0.5*mudum*a_scaled) * sin(0.5*M_PI*uu);
+	        double arg2 = (0.5*mudum) * cos(0.5*M_PI*uu);
+	        double arg3=  (0.5*mudum*ta) * sin(0.5*M_PI*uu);
+	        double arg4=  (0.5*mudum*tb) * cos(0.5*M_PI*uu);
+
+	        if(arg1==0.0){
+		        t1 = 1.0;
+	        } else {
+		        t1 = (sin(arg1)/arg1);                //defn for CSPP model sin(arg1)/arg1    test:  (sin(arg1)/arg1)*(sin(arg1)/arg1)   
+	        }
+	        if(arg2==0.0){
+		        t2 = 1.0;
+	        } else {
+		        t2 = (sin(arg2)/arg2);           //defn for CSPP model sin(arg2)/arg2   test: (sin(arg2)/arg2)*(sin(arg2)/arg2)    
+	        }	
+	        if(arg3==0.0){
+		        t3 = 1.0;
+	        } else {
+		        t3 = sin(arg3)/arg3;
+	        }
+	        if(arg4==0.0){
+		        t4 = 1.0;
+	        } else {
+		        t4 = sin(arg4)/arg4;
+	        }
+            
+            // Expression in libCylinder.c (neither drC nor Vot are used)
+	        tmp =( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 );   //  correct FF : square of sum of phase factors            
+            
+            // To note also that in csparallelepiped.cpp, there is a function called
+            // cspkernel, which seems to make more sense and has the following comment:
+            //   This expression is different from NIST/IGOR package (I strongly believe the IGOR is wrong!!!). 10/15/2010.
+            //   tmp =( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3)*
+            //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors
+            // This is the function called by csparallelepiped_analytical_2D_scaled,
+            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c        
+            
+            summj += Gauss76Wt[j] * tmp;
+
+        }
+		
+        // value of the inner integral
+        answer = 0.5 * summj;
+
+		// finish the outer integral 
+		double arg = 0.5 * mu* c_scaled *sigma;
+		if ( arg == 0.0 ) {
+			answer *= 1.0;
+		} else {
+	        answer *= sin(arg)*sin(arg)/arg/arg;
+		}
+		
+		// now sum up the outer integral
+		summ += Gauss76Wt[i] * answer;
+
+    }
+    
+	answer = 0.5 * summ;
+
+	//convert from [1e-12 A-1] to [cm-1]
+	answer *= 1.0e-4;
+	
+	return answer;
+}
+
+double Iqxy(double qx, double qy,
+    double core_sld,
+    double arim_sld,
+    double brim_sld,
+    double crim_sld,
+    double solvent_sld,
+    double a_side,
+    double b_side,
+    double c_side,
+    double arim_thickness,
+    double brim_thickness,
+    double crim_thickness,
+    double theta,
+    double phi,
+    double psi)
+{
+    double tmp1, tmp2, tmp3, tmpt1, tmpt2, tmpt3;   
+
+    double q = sqrt(qx*qx+qy*qy);
+    double qx_scaled = qx/q;
+    double qy_scaled = qy/q;
+
+    // Convert angles given in degrees to radians
+    theta *= M_PI_180;
+    phi   *= M_PI_180;
+    psi   *= M_PI_180;
+    
+    // Parallelepiped c axis orientation
+    double cparallel_x = cos(theta) * cos(phi);
+    double cparallel_y = sin(theta);
+    
+    // Compute angle between q and parallelepiped axis
+    double cos_val_c = cparallel_x*qx_scaled + cparallel_y*qy_scaled;// + cparallel_z*qz;
+
+    // Parallelepiped a axis orientation
+    double parallel_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
+    double parallel_y = sin(psi)*cos(theta);
+    double cos_val_a = parallel_x*qx_scaled + parallel_y*qy_scaled;
+
+    // Parallelepiped b axis orientation
+    double bparallel_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
+    double bparallel_y = cos(theta)*cos(psi);
+    double cos_val_b = bparallel_x*qx_scaled + bparallel_y*qy_scaled;
+
+    // The following tests should always pass
+    if (fabs(cos_val_c)>1.0) {
+      //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+      cos_val_c = 1.0;
+    }
+    if (fabs(cos_val_a)>1.0) {
+      //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+      cos_val_a = 1.0;
+    }
+    if (fabs(cos_val_b)>1.0) {
+      //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+      cos_val_b = 1.0;
+    }
+    
+    // cspkernel in csparallelepiped recoded here 
+    double dr0 = core_sld-solvent_sld;
+	double drA = arim_sld-solvent_sld;
+	double drB = brim_sld-solvent_sld;
+	double drC = crim_sld-solvent_sld;
+	double Vin = a_side * b_side * c_side;
+    // As for the 1D case, Vot is not used
+	//double Vot = (a_side * b_side * c_side +
+    //              2.0 * arim_thickness * b_side * c_side +
+    //              2.0 * a_side * brim_thickness * c_side +
+    //              2.0 * a_side * b_side * crim_thickness);
+	double V1 = (2.0 * arim_thickness * b_side * c_side);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc)
+	double V2 = (2.0 * a_side * brim_thickness * c_side);    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
+    double V3 = (2.0 * a_side * b_side * crim_thickness);
+    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
+    // the scaling by B. The use of a_side for the 3 of them seems clearly a mistake to me,
+    // but for the moment I let it like this until understanding better the code.
+	double ta = a_side + 2.0*arim_thickness;
+    double tb = a_side + 2.0*brim_thickness;
+    double tc = a_side + 2.0*crim_thickness;
+    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
+    double argA = 0.5*q*a_side*cos_val_a;
+    double argB = 0.5*q*b_side*cos_val_b;
+    double argC = 0.5*q*c_side*cos_val_c;
+    double argtA = 0.5*q*ta*cos_val_a;
+    double argtB = 0.5*q*tb*cos_val_b;
+    double argtC = 0.5*q*tc*cos_val_c;
+    
+    if(argA==0.0) {
+        tmp1 = 1.0;
+    } else {
+        tmp1 = sin(argA)/argA;
+    }
+    if (argB==0.0) {
+        tmp2 = 1.0;
+    } else {
+        tmp2 = sin(argB)/argB;
+    }
+    if (argC==0.0) {
+        tmp3 = 1.0;
+    } else {
+        tmp3 = sin(argC)/argC;
+    }
+    if(argtA==0.0) {
+        tmpt1 = 1.0;
+    } else {
+        tmpt1 = sin(argtA)/argtA;
+    }
+    if (argtB==0.0) {
+        tmpt2 = 1.0;
+    } else {
+        tmpt2 = sin(argtB)/argtB;
+    }
+    if (argtC==0.0) {
+        tmpt3 = 1.0;
+    } else {
+        tmpt3 = sin(argtC)*sin(argtC)/argtC/argtC;
+    }
+
+    // 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*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1 + 
+               drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);
+   
+    return 1.0e-4 * f * f;
+}
Index: sasmodels/models/core_shell_parallelepiped.py
===================================================================
--- sasmodels/models/core_shell_parallelepiped.py	(revision 44bd2be361a937362169e016a821f398577e4661)
+++ sasmodels/models/core_shell_parallelepiped.py	(revision 44bd2be361a937362169e016a821f398577e4661)
@@ -0,0 +1,193 @@
+# core_shell_parallelepiped model
+# Note: model title and parameter table are inserted automatically
+r"""
+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 all three (pairs) of faces.**
+
+The form factor is normalized by the particle volume *V* such that
+
+*I(q)* = *scale* \* <*f*\ :sup:`2`> / *V* + *background*
+
+where < > is an average over all possible orientations of the rectangular solid.
+
+An instrument resolution smeared version of the model is also provided.
+
+
+Definition
+----------
+
+The function calculated is the form factor of the rectangular solid below.
+The core of the solid is defined by the dimensions *A*, *B*, *C* such that
+*A* < *B* < *C*.
+
+.. image:: img/core_shell_parallelepiped.jpg
+
+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
+
+.. 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 shell, in the same manner as a core-shell model.
+
+**For the calculation of the form factor to be valid, the sides of the solid
+MUST be chosen such that** *A* < *B* < *C*.
+**If this inequality is not satisfied, the model will not report an error,
+and the calculation will not be correct.**
+
+FITTING NOTES
+If the scale is set equal to the particle volume fraction, |phi|, the returned
+value is the scattered intensity per unit volume; ie, *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, and used as the effective radius
+for *S(Q)* when *P(Q)* \* *S(Q)* is applied.
+
+.. Comment by Miguel Gonzalez:
+   The later seems to contradict the previous statement that interparticle interference
+   effects are not included.
+
+To provide easy access to the orientation of the parallelepiped, we define the
+axis of the cylinder using three angles |theta|, |phi| and |bigpsi|.
+These angles are defined on Figure 2 of the :ref:`cylinder` model.
+The angle |bigpsi| is the rotational angle around the *long_c* axis against the
+*q* plane. For example, |bigpsi| = 0 when the *short_b* axis is parallel to the
+*x*-axis of the detector.
+
+.. figure:: img/parallelepiped_angles_definition.jpg
+
+    Definition of the angles for oriented core-shell parallelepipeds.
+
+.. figure:: img/parallelepiped_angles_examples.jpg
+
+    Examples of the angles for oriented core-shell parallelepipeds against the
+    detector plane.
+
+Validation
+----------
+
+The model uses the form factor calculations implemented in a c-library provided
+by the NIST Center for Neutron Research (Kline, 2006).
+
+REFERENCE
+
+P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
+Equations (1), (13-14). (in German)
+
+"""
+
+import numpy as np
+from numpy import pi, inf, sqrt
+
+name = "core_shell_parallelepiped"
+title = "Rectangular solid with a core-shell structure."
+description = """
+     P(q)= 
+"""
+category = "shape:parallelepiped"
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["core_sld", "1e-6/Ang^2", 1, [-inf, inf], "",
+               "Parallelepiped core scattering length density"],
+              ["arim_sld", "1e-6/Ang^2", 2, [-inf, inf], "",
+               "Parallelepiped A rim scattering length density"],
+              ["brim_sld", "1e-6/Ang^2", 4, [-inf, inf], "",
+               "Parallelepiped B rim scattering length density"],
+              ["crim_sld", "1e-6/Ang^2", 2, [-inf, inf], "",
+               "Parallelepiped C rim scattering length density"],
+              ["solvent_sld", "1e-6/Ang^2", 6, [-inf, inf], "",
+               "Solvent scattering length density"],
+              ["a_side", "Ang", 35, [0, inf], "volume",
+               "Shorter side of the parallelepiped"],
+              ["b_side", "Ang", 75, [0, inf], "volume",
+               "Second side of the parallelepiped"],
+              ["c_side", "Ang", 400, [0, inf], "volume",
+               "Larger side of the parallelepiped"],
+              ["arim_thickness", "Ang", 10, [0, inf], "volume",
+               "Thickness of A rim"],
+              ["brim_thickness", "Ang", 10, [0, inf], "volume",
+               "Thickness of B rim"],
+              ["crim_thickness", "Ang", 10, [0, inf], "volume",
+               "Thickness of C rim"],
+              ["theta", "degrees", 0, [-inf, inf], "orientation",
+               "In plane angle"],
+              ["phi", "degrees", 0, [-inf, inf], "orientation",
+               "Out of plane angle"],
+              ["psi", "degrees", 0, [-inf, inf], "orientation",
+               "Rotation angle around its own c axis against q plane"],
+             ]
+
+source = ["lib/J1.c", "lib/gauss76.c", "core_shell_parallelepiped.c"]
+
+
+def ER(a_side, b_side, c_side, arim_thickness, brim_thickness, crim_thickness):
+    """
+        Return equivalent radius (ER)
+    """
+
+    # surface average radius (rough approximation)
+    surf_rad = sqrt((a_side + 2.0*arim_thickness) * (b_side + 2.0*brim_thickness) / pi)
+
+    height = c_side + 2.0*crim_thickness
+
+    ddd = 0.75 * surf_rad * (2 * surf_rad * height + (height + surf_rad) * (height + pi * surf_rad))
+    return 0.5 * (ddd) ** (1. / 3.)
+
+# VR defaults to 1.0
+
+# parameters for demo
+demo = dict(scale=1, background=0.0,
+            core_sld=1e-6, arim_sld=2e-6, brim_sld=4e-6,
+            crim_sld=2e-6, solvent_sld=6e-6,
+            a_side=35, b_side=75, c_side=400,
+            arim_thickness=10, brim_thickness=10, crim_thickness=10,
+            theta=0, phi=0, psi=0,
+            a_side_pd=0.1, a_side_pd_n=1,
+            b_side_pd=0.1, b_side_pd_n=1,
+            c_side_pd=0.1, c_side_pd_n=1,
+            arim_thickness_pd=0.1, arim_thickness_pd_n=1,
+            brim_thickness_pd=0.1, brim_thickness_pd_n=1,
+            crim_thickness_pd=0.1, crim_thickness_pd_n=1,
+            theta_pd=10, theta_pd_n=1,
+            phi_pd=10, phi_pd_n=1,
+            psi_pd=10, psi_pd_n=10)
+
+
+# For testing against the old sasview models, include the converted parameter
+# names and the target sasview model name.
+oldname = 'CSParallelepipedModel'
+oldpars = dict(theta='parallel_theta', phi='parallel_phi', psi='parallel_psi',
+               core_sld='sld_pcore', arim_sld='sld_rimA', brim_sld='sld_rimB',
+               crim_sld='sld_rimC', solvent_sld='sld_solv',
+               a_side='shortA', b_side='midB', c_side='longC',
+               arim_thickness='rimA', brim_thickness='rimB', crim_thickness='rimC')
+
+qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5)
+tests = [[{}, 0.2, 0.532149288477],
+         [{}, [0.2], [0.532149288477]],
+         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.031102135569],
+         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.031102135569]],
+        ]
+del qx, qy  # not necessary to delete, but cleaner
Index: sasmodels/models/hollow_rectangular_prism.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism.c	(revision deb7ee06caa9916d8c0504bdc15ffbc16df623f6)
+++ sasmodels/models/hollow_rectangular_prism.c	(revision deb7ee06caa9916d8c0504bdc15ffbc16df623f6)
@@ -0,0 +1,115 @@
+double form_volume(double a_side, double b2a_ratio, double c2a_ratio, double thickness);
+double Iq(double q, double sld, double solvent_sld, double a_side, 
+          double b2a_ratio, double c2a_ratio, double thickness);
+double Iqxy(double qx, double qy, double sld, double solvent_sld, 
+            double a_side, double b2a_ratio, double c2a_ratio, double thickness);
+
+double form_volume(double a_side, double b2a_ratio, double c2a_ratio, double thickness)
+{
+    double b_side = a_side * b2a_ratio;
+    double c_side = a_side * c2a_ratio;
+    double a_core = a_side - 2.0*thickness;
+    double b_core = b_side - 2.0*thickness;
+    double c_core = c_side - 2.0*thickness;
+    double vol_core = a_core * b_core * c_core;
+    double vol_total = a_side * b_side * c_side;
+    double vol_shell = vol_total - vol_core;
+    return vol_shell;
+}
+
+double Iq(double q,
+    double sld,
+    double solvent_sld,
+    double a_side,
+    double b2a_ratio,
+    double c2a_ratio,
+    double thickness)
+{
+    double termA1, termA2, termB1, termB2, termC1, termC2;
+    
+    double b_side = a_side * b2a_ratio;
+    double c_side = a_side * c2a_ratio;
+    double a_half = 0.5 * a_side;
+    double b_half = 0.5 * b_side;
+    double c_half = 0.5 * c_side;
+
+   //Integration limits to use in Gaussian quadrature
+    double v1a = 0.0;
+    double v1b = 0.5 * M_PI;  //theta integration limits
+    double v2a = 0.0;
+    double v2b = 0.5 * M_PI;  //phi integration limits
+    
+    //Order of integration
+    int nordi=76;			        
+    int nordj=76;
+
+    double sumi = 0.0;
+    
+    for(int i=0; i<nordi; i++) {
+
+	    double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );	
+
+	    double arg = q * c_half * cos(theta);
+	    if (fabs(arg) > 1.e-16) {termC1 = sin(arg)/arg;} else {termC1 = 1.0;}
+	    arg = q * (c_half-thickness)*cos(theta);
+	    if (fabs(arg) > 1.e-16) {termC2 = sin(arg)/arg;} else {termC2 = 1.0;}
+
+	    double sumj = 0.0;
+        
+	    for(int j=0; j<nordj; j++) {
+
+            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
+
+            // Amplitude AP from eqn. (13), rewritten to avoid round-off effects when arg=0
+
+	        arg = q * a_half * sin(theta) * sin(phi);
+	        if (fabs(arg) > 1.e-16) {termA1 = sin(arg)/arg;} else {termA1 = 1.0;}
+	        arg = q * (a_half-thickness) * sin(theta) * sin(phi);
+	        if (fabs(arg) > 1.e-16) {termA2 = sin(arg)/arg;} else {termA2 = 1.0;}
+
+	        arg = q * b_half * sin(theta) * cos(phi);
+	        if (fabs(arg) > 1.e-16) {termB1 = sin(arg)/arg;} else {termB1 = 1.0;}
+	        arg = q * (b_half-thickness) * sin(theta) * cos(phi);
+	        if (fabs(arg) > 1.e-16) {termB2 = sin(arg)/arg;} else {termB2 = 1.0;}
+
+            double AP1 = (a_side*b_side*c_side) * termA1 * termB1 * termC1;
+            double AP2 = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness) * termA2 * termB2 * termC2;
+            double AP = AP1 - AP2;
+
+	        sumj += Gauss76Wt[j] * (AP*AP);
+
+	    }
+
+	    sumj = 0.5 * (v2b-v2a) * sumj;
+	    sumi += Gauss76Wt[i] * sumj * sin(theta);
+
+    }
+
+    double answer = 0.5*(v1b-v1a)*sumi;
+
+    // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization)
+    // The factor 2 is due to the different theta integration limit (pi/2 instead of pi)
+    answer *= (2.0/M_PI);
+
+    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization.
+    answer *= (sld-solvent_sld)*(sld-solvent_sld);
+
+    // Convert from [1e-12 A-1] to [cm-1]
+    answer *= 1.0e-4;
+
+    return answer;
+    
+}
+
+double Iqxy(double qx, double qy,
+    double sld,
+    double solvent_sld,
+    double a_side,
+    double b2a_ratio,
+    double c2a_ratio,
+    double thickness)
+{
+    double q = sqrt(qx*qx + qy*qy);
+    double intensity = Iq(q, sld, solvent_sld, a_side, b2a_ratio, c2a_ratio, thickness); 
+    return intensity;    
+}
Index: sasmodels/models/hollow_rectangular_prism.py
===================================================================
--- sasmodels/models/hollow_rectangular_prism.py	(revision deb7ee06caa9916d8c0504bdc15ffbc16df623f6)
+++ sasmodels/models/hollow_rectangular_prism.py	(revision deb7ee06caa9916d8c0504bdc15ffbc16df623f6)
@@ -0,0 +1,161 @@
+# rectangular_prism model
+# Note: model title and parameter table are inserted automatically
+r"""
+
+This model provides the form factor, *P(q)*, for a hollow rectangular
+parallelepiped with a wall of thickness |bigdelta|.
+It computes only the 1D scattering, not the 2D.
+
+Definition
+----------
+
+The 1D scattering intensity for this model is calculated by forming
+the difference of the amplitudes of two massive parallelepipeds
+differing in their outermost dimensions in each direction by the
+same length increment :math:`2\Delta` (Nayuk, 2012).
+
+As in the case of the massive parallelepiped model (:ref:`rectangular-prism`),
+the scattering amplitude is computed for a particular orientation of the
+parallelepiped with respect to the scattering vector and then averaged over all
+possible orientations, giving
+
+.. math::
+  P(q) =  \frac{1}{V^2} \frac{2}{\pi} \times \, \int_0^{\frac{\pi}{2}} \,
+  \int_0^{\frac{\pi}{2}} A_{P\Delta}^2(q) \, \sin\theta \, d\theta \, d\phi
+
+where |theta| is the angle between the *z* axis and the longest axis
+of the parallelepiped, |phi| is the angle between the scattering vector
+(lying in the *xy* plane) and the *y* axis, and
+
+.. math::
+  \begin{align}
+  A_{P\Delta}(q) & =  A B C
+    \left[\frac{\sin \bigl( q \frac{C}{2} \cos\theta \bigr)}
+    {\left( q \frac{C}{2} \cos\theta \right)} \right]
+    \left[\frac{\sin \bigl( q \frac{A}{2} \sin\theta \sin\phi \bigr)}
+    {\left( q \frac{A}{2} \sin\theta \sin\phi \right)}\right]
+    \left[\frac{\sin \bigl( q \frac{B}{2} \sin\theta \cos\phi \bigr)}
+    {\left( q \frac{B}{2} \sin\theta \cos\phi \right)}\right] \\
+    & - 8
+    \left(\frac{A}{2}-\Delta\right) \left(\frac{B}{2}-\Delta\right) \left(\frac{C}{2}-\Delta\right)
+    \left[ \frac{\sin \bigl[ q \bigl(\frac{C}{2}-\Delta\bigr) \cos\theta \bigr]}
+    {q \bigl(\frac{C}{2}-\Delta\bigr) \cos\theta} \right]
+    \left[ \frac{\sin \bigl[ q \bigl(\frac{A}{2}-\Delta\bigr) \sin\theta \sin\phi \bigr]}
+    {q \bigl(\frac{A}{2}-\Delta\bigr) \sin\theta \sin\phi} \right]
+    \left[ \frac{\sin \bigl[ q \bigl(\frac{B}{2}-\Delta\bigr) \sin\theta \cos\phi \bigr]}
+    {q \bigl(\frac{B}{2}-\Delta\bigr) \sin\theta \cos\phi} \right]
+  \end{align}
+
+where *A*, *B* and *C* are the external sides of the parallelepiped fulfilling
+:math:`A \le B \le C`, and the volume *V* of the parallelepiped is
+
+.. math::
+  V = A B C \, - \, (A - 2\Delta) (B - 2\Delta) (C - 2\Delta)
+
+The 1D scattering intensity is then calculated as
+
+.. math::
+  I(q) = \mbox{scale} \times V \times (\rho_{\mbox{p}} -
+  \rho_{\mbox{solvent}})^2 \times P(q)
+
+where :math:`\rho_{\mbox{p}}` is the scattering length of the parallelepiped,
+:math:`\rho_{\mbox{solvent}}` is the scattering length of the solvent,
+and (if the data are in absolute units) *scale* represents the volume fraction
+(which is unitless).
+
+**The 2D scattering intensity is not computed by this model.**
+
+
+Validation
+----------
+
+Validation of the code was conducted by qualitatively comparing the output
+of the 1D model to the curves shown in (Nayuk, 2012).
+
+REFERENCES
+
+R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
+
+"""
+
+from numpy import pi, inf, sqrt
+
+name = "hollow_rectangular_prism"
+title = "Hollow rectangular parallelepiped with uniform scattering length density."
+description = """
+    I(q)= scale*V*(sld - solvent_sld)^2*P(q,theta,phi)+background
+        P(q,theta,phi) = (2/pi/V^2) * double integral from 0 to pi/2 of ...
+           (AP1-AP2)^2(q)*sin(theta)*dtheta*dphi
+        AP1 = S(q*C*cos(theta)/2) * S(q*A*sin(theta)*sin(phi)/2) * S(q*B*sin(theta)*cos(phi)/2)
+        AP2 = S(q*C'*cos(theta)) * S(q*A'*sin(theta)*sin(phi)) * S(q*B'*sin(theta)*cos(phi))
+        C' = (C/2-thickness)
+        B' = (B/2-thickness)
+        A' = (A/2-thickness)
+        S(x) = sin(x)/x
+"""
+category = "shape:parallelepiped"
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["sld", "1e-6/Ang^2", 6.3, [-inf, inf], "",
+               "Parallelepiped scattering length density"],
+              ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "",
+               "Solvent scattering length density"],
+              ["a_side", "Ang", 35, [0, inf], "volume",
+               "Shorter side of the parallelepiped"],
+              ["b2a_ratio", "Ang", 1, [0, inf], "volume",
+               "Ratio sides b/a"],
+              ["c2a_ratio", "Ang", 1, [0, inf], "volume",
+               "Ratio sides c/a"],
+              ["thickness", "Ang", 1, [0, inf], "volume",
+               "Thickness of parallelepiped"],
+             ]
+
+source = ["lib/J1.c", "lib/gauss76.c", "hollow_rectangular_prism.c"]
+
+def ER(a_side, b2a_ratio, c2a_ratio, thickness):
+    """
+        Return equivalent radius (ER)
+        thickness parameter not used
+    """
+    b_side = a_side * b2a_ratio
+    c_side = a_side * c2a_ratio
+
+    # surface average radius (rough approximation)
+    surf_rad = sqrt(a_side * b_side / pi)
+
+    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
+    return 0.5 * (ddd) ** (1. / 3.)
+
+def VR(a_side, b2a_ratio, c2a_ratio, thickness):
+    """
+        Return shell volume and total volume
+    """
+    b_side = a_side * b2a_ratio
+    c_side = a_side * c2a_ratio
+    a_core = a_side - 2.0*thickness
+    b_core = b_side - 2.0*thickness
+    c_core = c_side - 2.0*thickness
+    vol_core = a_core * b_core * c_core
+    vol_total = a_side * b_side * c_side
+    vol_shell = vol_total - vol_core
+    return vol_total, vol_shell
+
+
+# parameters for demo
+demo = dict(scale=1, background=0,
+            sld=6.3e-6, solvent_sld=1.0e-6,
+            a_side=35, b2a_ratio=1, c2a_ratio=1, thickness=1,
+            a_side_pd=0.1, a_side_pd_n=10,
+            b2a_ratio_pd=0.1, b2a_ratio_pd_n=1,
+            c2a_ratio_pd=0.1, c2a_ratio_pd_n=1)
+
+# For testing against the old sasview models, include the converted parameter
+# names and the target sasview model name.
+oldname = 'RectangularHollowPrismModel'
+oldpars = dict(a_side='short_side', b2a_ratio='b2a_ratio', c_side='c2a_ratio',
+               thickness='thickness', sld='sldPipe', solvent_sld='sldSolv')
+
+tests = [[{}, 0.2, 0.76587283098],
+         [{}, [0.2], [0.76587283098]],
+        ]
+
Index: sasmodels/models/hollow_rectangular_prism_infinitely_thin_walls.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism_infinitely_thin_walls.c	(revision deb7ee06caa9916d8c0504bdc15ffbc16df623f6)
+++ sasmodels/models/hollow_rectangular_prism_infinitely_thin_walls.c	(revision deb7ee06caa9916d8c0504bdc15ffbc16df623f6)
@@ -0,0 +1,97 @@
+double form_volume(double a_side, double b2a_ratio, double c2a_ratio);
+double Iq(double q, double sld, double solvent_sld, double a_side, 
+          double b2a_ratio, double c2a_ratio);
+double Iqxy(double qx, double qy, double sld, double solvent_sld, 
+            double a_side, double b2a_ratio, double c2a_ratio);
+
+double form_volume(double a_side, double b2a_ratio, double c2a_ratio)
+{
+    double b_side = a_side * b2a_ratio;
+    double c_side = a_side * c2a_ratio;
+    double vol_shell = 2.0 * (a_side*b_side + a_side*c_side + b_side*c_side);
+    return vol_shell;
+}
+
+double Iq(double q,
+    double sld,
+    double solvent_sld,
+    double a_side,
+    double b2a_ratio,
+    double c2a_ratio)
+{
+    double b_side = a_side * b2a_ratio;
+    double c_side = a_side * c2a_ratio;
+    double a_half = 0.5 * a_side;
+    double b_half = 0.5 * b_side;
+    double c_half = 0.5 * c_side;
+
+   //Integration limits to use in Gaussian quadrature
+    double v1a = 0.0;
+    double v1b = 0.5 * M_PI;  //theta integration limits
+    double v2a = 0.0;
+    double v2b = 0.5 * M_PI;  //phi integration limits
+    
+    //Order of integration
+    int nordi=76;			        
+    int nordj=76;
+
+    double sumi = 0.0;
+    
+    for(int i=0; i<nordi; i++) {
+
+	    double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );	
+        
+        // To check potential problems if denominator goes to zero here !!!
+        double termAL_theta = 8.0*cos(q*c_half*cos(theta)) / (q*q*sin(theta)*sin(theta));
+        double termAT_theta = 8.0*sin(q*c_half*cos(theta)) / (q*q*sin(theta)*cos(theta));
+
+	    double sumj = 0.0;
+        
+	    for(int j=0; j<nordj; j++) {
+
+            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
+            
+            // Amplitude AL from eqn. (7c)
+            double AL = termAL_theta * sin(q*a_half*sin(theta)*sin(phi)) * 
+                sin(q*b_half*sin(theta)*cos(phi)) / (sin(phi)*cos(phi));
+
+            // Amplitude AT from eqn. (9)
+            double AT = termAT_theta * (  (cos(q*a_half*sin(theta)*sin(phi))*sin(q*b_half*sin(theta)*cos(phi))/cos(phi)) 
+                + (cos(q*b_half*sin(theta)*cos(phi))*sin(q*a_half*sin(theta)*sin(phi))/sin(phi)) );
+
+            sumj += Gauss76Wt[j] * (AL+AT)*(AL+AT);
+
+	    }
+
+	    sumj = 0.5 * (v2b-v2a) * sumj;
+	    sumi += Gauss76Wt[i] * sumj * sin(theta);
+
+    }
+
+    double answer = 0.5*(v1b-v1a)*sumi;
+
+    // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization)
+    // The factor 2 is due to the different theta integration limit (pi/2 instead of pi)
+    answer *= (2.0/M_PI);
+
+    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization.
+    answer *= (sld-solvent_sld)*(sld-solvent_sld);
+
+    // Convert from [1e-12 A-1] to [cm-1]
+    answer *= 1.0e-4;
+
+    return answer;
+    
+}
+
+double Iqxy(double qx, double qy,
+    double sld,
+    double solvent_sld,
+    double a_side,
+    double b2a_ratio,
+    double c2a_ratio)
+{
+    double q = sqrt(qx*qx + qy*qy);
+    double intensity = Iq(q, sld, solvent_sld, a_side, b2a_ratio, c2a_ratio); 
+    return intensity;    
+}
Index: sasmodels/models/hollow_rectangular_prism_infinitely_thin_walls.py
===================================================================
--- sasmodels/models/hollow_rectangular_prism_infinitely_thin_walls.py	(revision deb7ee06caa9916d8c0504bdc15ffbc16df623f6)
+++ sasmodels/models/hollow_rectangular_prism_infinitely_thin_walls.py	(revision deb7ee06caa9916d8c0504bdc15ffbc16df623f6)
@@ -0,0 +1,143 @@
+# rectangular_prism model
+# Note: model title and parameter table are inserted automatically
+r"""
+
+This model provides the form factor, *P(q)*, for a hollow rectangular
+prism with infinitely thin walls. It computes only the 1D scattering, not the 2D.
+
+
+Definition
+----------
+
+The 1D scattering intensity for this model is calculated according to the
+equations given by Nayuk and Huber (Nayuk, 2012).
+
+Assuming a hollow parallelepiped with infinitely thin walls, edge lengths
+:math:`A \le B \le C` and presenting an orientation with respect to the
+scattering vector given by |theta| and |phi|, where |theta| is the angle
+between the *z* axis and the longest axis of the parallelepiped *C*, and
+|phi| is the angle between the scattering vector (lying in the *xy* plane)
+and the *y* axis, the form factor is given by
+
+.. math::
+  P(q) =  \frac{1}{V^2} \frac{2}{\pi} \int_0^{\frac{\pi}{2}}
+  \int_0^{\frac{\pi}{2}} [A_L(q)+A_T(q)]^2 \sin\theta d\theta d\phi
+
+where
+
+.. math::
+  V = 2AB + 2AC + 2BC
+
+.. math::
+  A_L(q) =  8 \times \frac{ \sin \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr)
+                              \sin \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr)
+                              \cos \bigl( q \frac{C}{2} \cos\theta \bigr) }
+                            {q^2 \, \sin^2\theta \, \sin\phi \cos\phi}
+
+.. math::
+  A_T(q) =  A_F(q) \times \frac{2 \, \sin \bigl( q \frac{C}{2} \cos\theta \bigr)}{q \, \cos\theta}
+
+and
+
+.. math::
+  A_F(q) =  4 \frac{ \cos \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr)
+                       \sin \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr) }
+                     {q \, \cos\phi \, \sin\theta} +
+              4 \frac{ \sin \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr)
+                       \cos \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr) }
+                     {q \, \sin\phi \, \sin\theta}
+
+The 1D scattering intensity is then calculated as
+
+.. math::
+  I(q) = \mbox{scale} \times V \times (\rho_{\mbox{p}} - \rho_{\mbox{solvent}})^2 \times P(q)
+
+where *V* is the volume of the rectangular prism, :math:`\rho_{\mbox{p}}`
+is the scattering length of the parallelepiped, :math:`\rho_{\mbox{solvent}}`
+is the scattering length of the solvent, and (if the data are in absolute
+units) *scale* represents the volume fraction (which is unitless).
+
+**The 2D scattering intensity is not computed by this model.**
+
+
+Validation
+----------
+
+Validation of the code was conducted  by qualitatively comparing the output
+of the 1D model to the curves shown in (Nayuk, 2012).
+
+REFERENCES
+
+R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
+
+"""
+
+from numpy import pi, inf, sqrt
+
+name = "hollow_rectangular_prism_infinitely_thin_walls"
+title = "Hollow rectangular parallelepiped with infinitely thin walls."
+description = """
+    I(q)= scale*V*(sld - solvent_sld)^2*P(q)+background
+        with P(q) being the form factor corresponding to a hollow rectangular
+        parallelepiped with infinitely thin walls.
+"""
+category = "shape:parallelepiped"
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["sld", "1e-6/Ang^2", 6.3, [-inf, inf], "",
+               "Parallelepiped scattering length density"],
+              ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "",
+               "Solvent scattering length density"],
+              ["a_side", "Ang", 35, [0, inf], "volume",
+               "Shorter side of the parallelepiped"],
+              ["b2a_ratio", "Ang", 1, [0, inf], "volume",
+               "Ratio sides b/a"],
+              ["c2a_ratio", "Ang", 1, [0, inf], "volume",
+               "Ratio sides c/a"],
+             ]
+
+source = ["lib/J1.c", "lib/gauss76.c", "hollow_rectangular_prism_infinitely_thin_walls.c"]
+
+def ER(a_side, b2a_ratio, c2a_ratio):
+    """
+        Return equivalent radius (ER)
+    """
+    b_side = a_side * b2a_ratio
+    c_side = a_side * c2a_ratio
+
+    # surface average radius (rough approximation)
+    surf_rad = sqrt(a_side * b_side / pi)
+
+    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
+    return 0.5 * (ddd) ** (1. / 3.)
+
+def VR(a_side, b2a_ratio, c2a_ratio):
+    """
+        Return shell volume and total volume
+    """
+    b_side = a_side * b2a_ratio
+    c_side = a_side * c2a_ratio
+    vol_total = a_side * b_side * c_side
+    vol_shell = 2.0 * (a_side*b_side + a_side*c_side + b_side*c_side)
+    return vol_shell, vol_total
+
+
+# parameters for demo
+demo = dict(scale=1, background=0,
+            sld=6.3e-6, solvent_sld=1.0e-6,
+            a_side=35, b2a_ratio=1, c2a_ratio=1,
+            a_side_pd=0.1, a_side_pd_n=10,
+            b2a_ratio_pd=0.1, b2a_ratio_pd_n=1,
+            c2a_ratio_pd=0.1, c2a_ratio_pd_n=1)
+
+# For testing against the old sasview models, include the converted parameter
+# names and the target sasview model name.
+oldname = 'RectangularHollowPrismInfThinWallsModel'
+oldpars = dict(a_side='short_side', b2a_ratio='b2a_ratio', c_side='c2a_ratio',
+               sld='sldPipe', solvent_sld='sldSolv')
+
+tests = [[{}, 0.2, 0.836719188592],
+         [{}, [0.2], [0.836719188592]],
+        ]
+
+
Index: sasmodels/models/parallelepiped.c
===================================================================
--- sasmodels/models/parallelepiped.c	(revision c5b7d07435739547137368520d072f6e87ac3208)
+++ sasmodels/models/parallelepiped.c	(revision deb7ee06caa9916d8c0504bdc15ffbc16df623f6)
@@ -9,7 +9,7 @@
     double argA,argB,argC,tmp1,tmp2,tmp3;
     //handle arg=0 separately, as sin(t)/t -> 1 as t->0
-    argA = a*ala/2.0;
-    argB = b*alb/2.0;
-    argC = c*alc/2.0;
+    argA = 0.5*a*ala;
+    argB = 0.5*b*alb;
+    argC = 0.5*c*alc;
     if(argA==0.0) {
         tmp1 = 1.0;
@@ -31,8 +31,10 @@
 }
 
+
 double form_volume(double a_side, double b_side, double c_side)
 {
     return a_side * b_side * c_side;
 }
+
 
 double Iq(double q,
@@ -43,5 +45,5 @@
     double c_side)
 {
-    double tmp1, tmp2, yyy;
+    double tmp1, tmp2;
     
     double mu = q * b_side;
@@ -50,21 +52,25 @@
     double a_scaled = a_side / b_side;
     double c_scaled = c_side / b_side;
-    
-    // outer integral (with 76 gauss points), integration limits = 0, 1
+        
+    //Order of integration
+    int nordi=76;			        
+    int nordj=76;
+
+    // outer integral (with nordi gauss points), integration limits = 0, 1
     double summ = 0; //initialize integral
 
-    for( int i=0; i<76; i++) {
+    for( int i=0; i<nordi; i++) {
 		
-        // inner integral (with 76 gauss points), integration limits = 0, 1
+        // inner integral (with nordj gauss points), integration limits = 0, 1
 	
         double summj = 0.0;
-	double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );		
+	    double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );		
 		
-	for(int j=0; j<76; j++) {
+	    for(int j=0; j<nordj; j++) {
 
             double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
             double mudum = mu * sqrt(1.0-sigma*sigma);
-	    double arg1 = 0.5 * mudum * cos(0.5*M_PI*uu);
-	    double arg2 = 0.5 * mudum * a_scaled * sin(0.5*M_PI*uu);
+	        double arg1 = 0.5 * mudum * cos(0.5*M_PI*uu);
+	        double arg2 = 0.5 * mudum * a_scaled * sin(0.5*M_PI*uu);
             if(arg1==0.0) {
 	        tmp1 = 1.0;
@@ -77,7 +83,6 @@
                 tmp2 = sin(arg2)*sin(arg2)/arg2/arg2;
             }
-            yyy = Gauss76Wt[j] * tmp1 * tmp2;
 
-            summj += yyy;
+            summj += Gauss76Wt[j] * tmp1 * tmp2;
         }
 		
@@ -92,11 +97,12 @@
         }
 		
-	// sum of outer integral
-	yyy = Gauss76Wt[i] * answer;
-        summ += yyy;
+	    // sum of outer integral
+        summ += Gauss76Wt[i] * answer;
         
     }	
    
     const double vd = (sld-solvent_sld) * form_volume(a_side, b_side, c_side);
+    
+    // convert from [1e-12 A-1] to [cm-1] and 0.5 factor for outer integral
     return 1.0e-4 * 0.5 * vd * vd * summ;
     
Index: sasmodels/models/parallelepiped.py
===================================================================
--- sasmodels/models/parallelepiped.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/parallelepiped.py	(revision deb7ee06caa9916d8c0504bdc15ffbc16df623f6)
@@ -7,17 +7,8 @@
 ----------
 
-This model provides the form factor, $P(q)$, for a rectangular parallelepiped
-(below) where the form factor is normalized by the volume of the
-parallelepiped. If you need to apply polydispersity, see also
-rectangular_prism_.
-
-The calculated form factor is:
-
-.. math::
-
-    P(q) = \frac{\text{scale}}{V} F^2(q) + \text{background}
-
-where the volume $V = A B C$ and the averaging $\left<\ldots\right>$ is
-applied over all orientations for 1D.
+| This model calculates the scattering from a rectangular parallelepiped (:num:`Figure #parallelepiped-image`).
+| If you need to apply polydispersity, see also :ref:`rectangular-prism`.
+
+.. _parallelepiped-image:
 
 .. figure:: img/parallelepiped.jpg
@@ -25,12 +16,37 @@
    Parallelepiped with the corresponding definition of sides.
 
-The edge of the solid must satisfy the condition that $A < B < C$.
-Then, assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the
-form factor is
-
-.. math::
-
-    P(q) = \frac{\text{scale}}{V}\int_0^1
-        \phi\left(\mu \sqrt{1-\sigma^2},a\right)
+.. 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 1D scattering intensity $I(q)$ is calculated as:
+
+.. Comment by Miguel Gonzalez:
+   I am modifying the original text because I find the notation a little bit
+   confusing. I think that in most textbooks/papers, the notation P(Q) is
+   used for the form factor (adim, P(Q=0)=1), although F(q) seems also to
+   be used. But here (as for many other models), P(q) is used to represent
+   the scattering intensity (in cm-1 normally). It would be good to agree on
+   a common notation.
+
+.. math::
+
+    I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 \left< P(q, \alpha) \right>
+    + \text{background}
+
+where the volume $V = A B C$, the contrast is defined as
+:math:`\Delta\rho = \rho_{\textstyle p} - \rho_{\textstyle solvent}`,
+$P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented
+at an angle $\alpha$ (angle between the long axis C and :math:`\vec q`),
+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)
+
+.. math::
+
+    P(q, \alpha) = \int_0^1 \phi_Q\left(\mu \sqrt{1-\sigma^2},a\right)
         \left[S(\mu c \sigma/2)\right]^2 d\sigma
 
@@ -39,7 +55,7 @@
 .. math::
 
-    \phi(\mu,a) = \int_0^1
-        \left\{S\left[\frac{\mu}{2}\cos(\frac{\pi}{2}u)\right]
-               S\left[\frac{\mu a}{2}\sin(\frac{\pi}{2}u)\right]
+    \phi_Q(\mu,a) = \int_0^1
+        \left\{S\left[\frac{\mu}{2}\cos\left(\frac{\pi}{2}u\right)\right]
+               S\left[\frac{\mu a}{2}\sin\left(\frac{\pi}{2}u\right)\right]
                \right\}^2 du
 
@@ -48,14 +64,8 @@
     \mu = qB
 
-and the contrast is defined as
-
-.. math::
-
-    \Delta\rho = \rho_{\textstyle p} - \rho_{\textstyle solvent}
-
-The scattering intensity per unit volume is returned in units of |cm^-1|;
-i.e., $I(q) = \phi P(q)$.
-
-NB: The 2nd virial coefficient of the parallelpiped is calculated based on
+
+The scattering intensity per unit volume is returned in units of |cm^-1|.
+
+NB: The 2nd virial coefficient of the parallelepiped is calculated based on
 the averaged effective radius $(=\sqrt{A B / \pi})$ and
 length $(= C)$ values, and used as the effective radius for
@@ -65,18 +75,65 @@
 three angles $\theta$, $\phi$ and $\Psi$. The definition of $\theta$ and
 $\phi$ is the same as for the cylinder model (see also figures below).
-The angle $\Psi$ is the rotational angle around the $C$ axis against
-the $q$ plane. For example, $\Psi = 0$ when the $B$ axis is parallel
-to the $x$-axis of the detector.
-
+
+.. Comment by Miguel Gonzalez:
+   The following text has been commented because I think there are two
+   mistakes. Psi is the rotational angle around C (but I cannot understand
+   what it means against the q plane) and psi=0 corresponds to a||x and b||y.
+
+   The angle $\Psi$ is the rotational angle around the $C$ axis against
+   the $q$ plane. For example, $\Psi = 0$ when the $B$ axis is parallel
+   to the $x$-axis of the detector.
+
+The angle $\Psi$ is the rotational angle around the $C$ axis.
+For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the $B$ axis
+oriented parallel to the y-axis of the detector with $A$ along the z-axis.
+For other $\theta$, $\phi$ values, the parallelepiped has to be first rotated
+$\theta$ degrees around $z$ and $\phi$ degrees around $y$,
+before doing a final rotation of $\Psi$ degrees around the resulting $C$ to
+obtain the final orientation of the parallelepiped.
+For example, for $\theta = 0$ and $\phi = 90$, we have that $\Psi = 0$
+corresponds to $A$ along $x$ and $B$ along $y$,
+while for $\theta = 90$ and $\phi = 0$, $\Psi = 0$ corresponds to
+$A$ along $z$ and $B$ along $x$.
 
 .. _parallelepiped-orientation:
 
-.. figure:: img/orientation.jpg
+.. figure:: img/parallelepiped_angles_definition.jpg
 
     Definition of the angles for oriented parallelepipeds.
 
-.. figure:: img/orientation2.jpg
+.. figure:: img/parallelepiped_angles_examples.jpg
 
     Examples of the angles for oriented parallelepipeds against the detector plane.
+
+For a given orientation of the parallelepiped, the 2D form factor is calculated as
+
+.. 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
+
+with
+
+.. math::
+
+    \cos\alpha = \hat A \cdot \hat q, \;
+    \cos\beta  = \hat B \cdot \hat q, \;
+    \cos\gamma = \hat C \cdot \hat q
+
+and the scattering intensity as:
+
+.. math::
+
+    I(q_x, q_y) = \frac{\text{scale}}{V} V^2 \Delta\rho^2 P(q_x, q_y) + \text{background}
+
+.. Comment by Miguel Gonzalez:
+   This reflects the logic of the code, as in parallelepiped.c the call
+   to _pkernel returns P(q_x, q_y) and then this is multiplied by V^2 * (delta rho)^2.
+   And finally outside parallelepiped.c it will be multiplied by scale, normalized by
+   V and the background added. But mathematically it makes more sense to write
+   I(q_x, q_y) = \text{scale} V \Delta\rho^2 P(q_x, q_y) + \text{background},
+   with scale being the volume fraction.
 
 
@@ -97,9 +154,10 @@
    Comparison between 1D and averaged 2D.
 
-This model reimplements the form factor calculations implemented in a c-library
+This model is based on form factor calculations implemented in a c-library
 provided by the NIST Center for Neutron Research (Kline, 2006).
 
 """
 
+import numpy as np
 from numpy import pi, inf, sqrt
 
@@ -107,13 +165,17 @@
 title = "Rectangular parallelepiped with uniform scattering length density."
 description = """
-     P(q)= scale/V*integral from 0 to 1 of ...
+    I(q)= scale*V*(sld - solvent_sld)^2*P(q,alpha)+background
+        P(q,alpha) = integral from 0 to 1 of ...
            phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma
-
+        with
             phi(mu,a) = integral from 0 to 1 of ..
-        (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du
+            (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du
             S(x) = sin(x)/x
-        mu = q*B
+            mu = q*B
+        V: Volume of the rectangular parallelepiped
+        alpha: angle between the long axis of the 
+            parallelepiped and the q-vector for 1D
 """
-category = "shape:parallelpiped"
+category = "shape:parallelepiped"
 
 #             ["name", "units", default, [lower, upper], "type","description"],
@@ -139,16 +201,15 @@
 
 def ER(a_side, b_side, c_side):
+    """
+        Return effective radius (ER) for P(q)*S(q)
+    """
 
     # surface average radius (rough approximation)
     surf_rad = sqrt(a_side * b_side / pi)
 
-    # DiamCyl recoded here (to check and possibly put in a library?)
-    a = surf_rad
-    b = 0.5 * c_side
-    t1 = a * a * b
-    t2 = 1.0 + (b / a) * (1.0 + a / b / 2.0) * (1.0 + pi * a / b / 2.0)
-    ddd = 3.0 * t1 * t2
-
+    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
     return 0.5 * (ddd) ** (1. / 3.)
+
+# VR defaults to 1.0
 
 # parameters for demo
@@ -171,2 +232,10 @@
                sld='sldPipe', solvent_sld='sldSolv')
 
+
+qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5)
+tests = [[{}, 0.2, 0.17658004974],
+         [{}, [0.2], [0.17658004974]],
+         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.00460296014],
+         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.00460296014]],
+        ]
+del qx, qy  # not necessary to delete, but cleaner
Index: sasmodels/models/rectangular_prism.c
===================================================================
--- sasmodels/models/rectangular_prism.c	(revision deb7ee06caa9916d8c0504bdc15ffbc16df623f6)
+++ sasmodels/models/rectangular_prism.c	(revision deb7ee06caa9916d8c0504bdc15ffbc16df623f6)
@@ -0,0 +1,101 @@
+double form_volume(double a_side, double b2a_ratio, double c2a_ratio);
+double Iq(double q, double sld, double solvent_sld, double a_side, 
+          double b2a_ratio, double c2a_ratio);
+double Iqxy(double qx, double qy, double sld, double solvent_sld, 
+            double a_side, double b2a_ratio, double c2a_ratio);
+
+double form_volume(double a_side, double b2a_ratio, double c2a_ratio)
+{
+    return a_side * (a_side*b2a_ratio) * (a_side*c2a_ratio);
+}
+
+double Iq(double q,
+    double sld,
+    double solvent_sld,
+    double a_side,
+    double b2a_ratio,
+    double c2a_ratio)
+{
+    double termA, termB, termC;
+    
+    double b_side = a_side * b2a_ratio;
+    double c_side = a_side * c2a_ratio;
+    double volume = a_side * b_side * c_side;
+    double a_half = 0.5 * a_side;
+    double b_half = 0.5 * b_side;
+    double c_half = 0.5 * c_side;
+
+   //Integration limits to use in Gaussian quadrature
+    double v1a = 0.0;
+    double v1b = 0.5 * M_PI;  //theta integration limits
+    double v2a = 0.0;
+    double v2b = 0.5 * M_PI;  //phi integration limits
+    
+    //Order of integration
+    int nordi=76;			        
+    int nordj=76;
+
+    double sumi = 0.0;
+    
+    for(int i=0; i<nordi; i++) {
+
+	    double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );	
+
+	    double arg = q * c_half * cos(theta);
+	    if (fabs(arg) > 1.e-16) {termC = sin(arg)/arg;} else {termC = 1.0;}  
+
+	    double sumj = 0.0;
+        
+	    for(int j=0; j<nordj; j++) {
+
+            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
+
+	        // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0
+
+	        arg = q * a_half * sin(theta) * sin(phi); 
+	        if (fabs(arg) > 1.e-16) {termA = sin(arg)/arg;} else {termA = 1.0;}
+	       
+	        arg = q * b_half * sin(theta) * cos(phi); 
+	        if (fabs(arg) > 1.e-16) {termB = sin(arg)/arg;} else {termB = 1.0;}	  
+               
+	        double AP = termA * termB * termC;  
+
+	        sumj += Gauss76Wt[j] * (AP*AP);
+
+	    }
+
+	    sumj = 0.5 * (v2b-v2a) * sumj;
+	    sumi += Gauss76Wt[i] * sumj * sin(theta);
+
+    }
+
+    double answer = 0.5*(v1b-v1a)*sumi;
+
+    // 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 
+    // 0 and pi/2, instead of 0 to pi.
+    answer *= (2.0/M_PI); //Form factor P(q)
+
+    // Multiply by contrast^2 and volume^2
+    answer *= (sld-solvent_sld)*(sld-solvent_sld)*volume*volume;
+
+    // Convert from [1e-12 A-1] to [cm-1] 
+    answer *= 1.0e-4;
+
+    return answer;
+    
+}
+
+double Iqxy(double qx, double qy,
+    double sld,
+    double solvent_sld,
+    double a_side,
+    double b2a_ratio,
+    double c2a_ratio)
+{
+    double q = sqrt(qx*qx + qy*qy);
+    double intensity = Iq(q, sld, solvent_sld, a_side, b2a_ratio, c2a_ratio); 
+    return intensity;    
+}
Index: sasmodels/models/rectangular_prism.py
===================================================================
--- sasmodels/models/rectangular_prism.py	(revision deb7ee06caa9916d8c0504bdc15ffbc16df623f6)
+++ sasmodels/models/rectangular_prism.py	(revision deb7ee06caa9916d8c0504bdc15ffbc16df623f6)
@@ -0,0 +1,139 @@
+# rectangular_prism model
+# Note: model title and parameter table are inserted automatically
+r"""
+
+This model provides the form factor, *P(q)*, for a rectangular prism.
+
+Note that this model is almost totally equivalent to the existing
+:ref:`parallelepiped` model.
+The only difference is that the way the relevant
+parameters are defined here (*a*, *b/a*, *c/a* instead of *a*, *b*, *c*)
+allows to use polydispersity with this model while keeping the shape of
+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.
+
+
+Definition
+----------
+
+The 1D scattering intensity for this model was calculated by Mittelbach and
+Porod (Mittelbach, 1961), but the implementation here is closer to the
+equations given by Nayuk and Huber (Nayuk, 2012).
+Note also that the angle definitions used in the code and the present
+documentation correspond to those used in (Nayuk, 2012) (see Fig. 1 of
+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.
+
+In this model the scattering from a massive parallelepiped with an
+orientation with respect to the scattering vector given by |theta|
+and |phi|
+
+.. math::
+  A_P\,(q) =  \frac{\sin \bigl( q \frac{C}{2} \cos\theta \bigr)}{\left( q \frac{C}{2}
+  \cos\theta \right)} \, \times \, \frac{\sin \bigl( q \frac{A}{2} \sin\theta \sin\phi
+  \bigr)}{\left( q \frac{A}{2} \sin\theta \sin\phi \right)} \, \times \, \frac{\sin \bigl(
+  q \frac{B}{2} \sin\theta \cos\phi \bigr)}{\left( q \frac{B}{2} \sin\theta \cos\phi \right)}
+
+where *A*, *B* and *C* are the sides of the parallelepiped and must fulfill
+:math:`A \le B \le C`, |theta| is the angle between the *z* axis and the
+longest axis of the parallelepiped *C*, and |phi| is the angle between the
+scattering vector (lying in the *xy* plane) and the *y* axis.
+
+The normalized form factor in 1D is obtained averaging over all possible
+orientations
+
+.. math::
+  P(q) =  \frac{2}{\pi} \times \, \int_0^{\frac{\pi}{2}} \,
+  \int_0^{\frac{\pi}{2}} A_P^2(q) \, \sin\theta \, d\theta \, d\phi
+
+And the 1D scattering intensity is calculated as
+
+.. math::
+  I(q) = \mbox{scale} \times V \times (\rho_{\mbox{p}} -
+  \rho_{\mbox{solvent}})^2 \times P(q)
+
+where *V* is the volume of the rectangular prism, :math:`\rho_{\mbox{p}}`
+is the scattering length of the parallelepiped, :math:`\rho_{\mbox{solvent}}`
+is the scattering length of the solvent, and (if the data are in absolute
+units) *scale* represents the volume fraction (which is unitless).
+
+**The 2D scattering intensity is not computed by this model.**
+
+
+Validation
+----------
+
+Validation of the code was conducted by comparing the output of the 1D model
+to the output of the existing :ref:`parallelepiped` model.
+
+REFERENCES
+
+P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
+
+R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
+
+"""
+
+from numpy import pi, inf, sqrt
+
+name = "rectangular_prism"
+title = "Rectangular parallelepiped with uniform scattering length density."
+description = """
+    I(q)= scale*V*(sld - solvent_sld)^2*P(q,theta,phi)+background
+        P(q,theta,phi) = (2/pi) * double integral from 0 to pi/2 of ...
+           AP^2(q)*sin(theta)*dtheta*dphi
+        AP = S(q*C*cos(theta)/2) * S(q*A*sin(theta)*sin(phi)/2) * S(q*B*sin(theta)*cos(phi)/2)
+        S(x) = sin(x)/x
+"""
+category = "shape:parallelepiped"
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["sld", "1e-6/Ang^2", 6.3, [-inf, inf], "",
+               "Parallelepiped scattering length density"],
+              ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "",
+               "Solvent scattering length density"],
+              ["a_side", "Ang", 35, [0, inf], "volume",
+               "Shorter side of the parallelepiped"],
+              ["b2a_ratio", "Ang", 1, [0, inf], "volume",
+               "Ratio sides b/a"],
+              ["c2a_ratio", "Ang", 1, [0, inf], "volume",
+               "Ratio sides c/a"],
+             ]
+
+source = ["lib/J1.c", "lib/gauss76.c", "rectangular_prism.c"]
+
+def ER(a_side, b2a_ratio, c2a_ratio):
+    """
+        Return equivalent radius (ER)
+    """
+    b_side = a_side * b2a_ratio
+    c_side = a_side * c2a_ratio
+
+    # surface average radius (rough approximation)
+    surf_rad = sqrt(a_side * b_side / pi)
+
+    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
+    return 0.5 * (ddd) ** (1. / 3.)
+
+
+# parameters for demo
+demo = dict(scale=1, background=0,
+            sld=6.3e-6, solvent_sld=1.0e-6,
+            a_side=35, b2a_ratio=1, c2a_ratio=1,
+            a_side_pd=0.1, a_side_pd_n=10,
+            b2a_ratio_pd=0.1, b2a_ratio_pd_n=1,
+            c2a_ratio_pd=0.1, c2a_ratio_pd_n=1)
+
+# For testing against the old sasview models, include the converted parameter
+# names and the target sasview model name.
+oldname = 'RectangularPrismModel'
+oldpars = dict(a_side='short_side', b2a_ratio='b2a_ratio', c_side='c2a_ratio',
+               sld='sldPipe', solvent_sld='sldSolv')
+
+tests = [[{}, 0.2, 0.374248406825],
+         [{}, [0.2], [0.374248406825]],
+        ]
