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
