Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision a503bfddd742ec35fe2b801f1d75a0a5f11eccfa)
+++ sasmodels/generate.py	(revision 3271e207340eb2fd5f01b1d8ba963ab4fe9ede8c)
@@ -722,4 +722,6 @@
         parameters = COMMON_PARAMETERS + kernel_module.parameters,
         source = getattr(kernel_module, 'source', []),
+        oldname = kernel_module.oldname,
+        oldpars = kernel_module.oldpars,
         )
     # Fill in attributes which default to None
Index: sasmodels/models/cylinder.py
===================================================================
--- sasmodels/models/cylinder.py	(revision 5428233710b752fa243b118f71ef7801850b581f)
+++ sasmodels/models/cylinder.py	(revision 5428233710b752fa243b118f71ef7801850b581f)
@@ -43,5 +43,5 @@
 .. figure:: img/orientation2.jpg
 
-    Examples of the angles for oriented pp against the detector plane.
+    Examples of the angles for oriented cylinders against the detector plane.
 
 NB: The 2nd virial coefficient of the cylinder is calculated based on the
@@ -113,8 +113,8 @@
 title = "Right circular cylinder with uniform scattering length density."
 description = """
-     P(q)= 2*(sld - solvent_sld)*V*sin(qLcos(alpha/2))
-            /[qLcos(alpha/2)]*J1(qRsin(alpha/2))/[qRsin(alpha)]
+     f(q,alpha) = 2*(sld - solvent_sld)*V*sin(qLcos(alpha/2))
+                /[qLcos(alpha/2)]*J1(qRsin(alpha/2))/[qRsin(alpha)]
 
-            P(q,alpha)= scale/V*f(q)^(2)+background
+            P(q,alpha)= scale/V*f(q,alpha)^(2)+background
             V: Volume of the cylinder
             R: Radius of the cylinder
@@ -125,5 +125,5 @@
             :the ouput is P(q)=scale/V*integral
             from pi/2 to zero of...
-            f(q)^(2)*sin(alpha)*dalpha + background
+            f(q,alpha)^(2)*sin(alpha)*dalpha + background
 """
 
Index: sasmodels/models/fcc.c
===================================================================
--- sasmodels/models/fcc.c	(revision 3271e207340eb2fd5f01b1d8ba963ab4fe9ede8c)
+++ sasmodels/models/fcc.c	(revision 3271e207340eb2fd5f01b1d8ba963ab4fe9ede8c)
@@ -0,0 +1,175 @@
+double form_volume(double radius);
+double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld);
+double Iqxy(double qx, double qy, double dnn,
+    double d_factor, double radius,double sld, double solvent_sld,
+    double theta, double phi, double psi);
+
+double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
+double _BCCeval(double Theta, double Phi, double temp1, double temp3);
+double _sphereform(double q, double radius, double sld, double solvent_sld);
+
+
+double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) {
+
+	const double Da = d_factor*dnn;
+	const double temp1 = q*q*Da*Da;
+	const double temp3 = q*dnn;
+
+	double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI);
+	return(retVal);
+}
+
+double _BCCeval(double Theta, double Phi, double temp1, double temp3) {
+
+	double temp6,temp7,temp8,temp9,temp10;
+	double result;
+
+	temp6 = sin(Theta);
+	temp7 = sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)+cos(Theta);
+	temp8 = -1.0*sin(Theta)*cos(Phi)-sin(Theta)*sin(Phi)+cos(Theta);
+	temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)-cos(Theta);
+	temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
+	result = pow(1.0-(temp10*temp10),3)*temp6/((1.0-2.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10)));
+
+	return (result);
+}
+
+double _sphereform(double q, double radius, double sld, double solvent_sld){
+    const double qr = q*radius;
+    double sn, cn;
+    SINCOS(qr, sn, cn);
+    const double bes = (qr == 0.0 ? 1.0 : 3.0*(sn-qr*cn)/(qr*qr*qr));
+    const double fq = bes * (sld - solvent_sld)*form_volume(radius);
+    return 1.0e-4*fq*fq;
+}
+
+double form_volume(double radius){
+    return 1.333333333333333*M_PI*radius*radius*radius;
+}
+
+
+double Iq(double q, double dnn,
+  double d_factor, double radius,
+  double sld, double solvent_sld){
+
+	//Volume fraction calculated from lattice symmetry and sphere radius
+	const double s1 = dnn/sqrt(0.75);
+	const double latticescale = 2.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1);
+
+    const double va = 0.0;
+    const double vb = 2.0*M_PI;
+    const double vaj = 0.0;
+    const double vbj = M_PI;
+
+    double summ = 0.0;
+    double answer = 0.0;
+	for(int i=0; i<150; i++) {
+		//setup inner integral over the ellipsoidal cross-section
+		double summj=0.0;
+		const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;		//the outer dummy is phi
+		for(int j=0;j<150;j++) {
+			//20 gauss points for the inner integral
+			double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;		//the inner dummy is theta
+			double yyy = Gauss150Wt[j] * _BCC_Integrand(q,dnn,d_factor,ztheta,zphi);
+			summj += yyy;
+		}
+		//now calculate the value of the inner integral
+		double answer = (vbj-vaj)/2.0*summj;
+
+		//now calculate outer integral
+		summ = summ+(Gauss150Wt[i] * answer);
+	}		//final scaling is done at the end of the function, after the NT_FP64 case
+
+	answer = (vb-va)/2.0*summ;
+	answer = answer*_sphereform(q,radius,sld,solvent_sld)*latticescale;
+
+    return answer;
+
+
+}
+
+
+double Iqxy(double qx, double qy, double dnn,
+    double d_factor, double radius,double sld, double solvent_sld,
+    double theta, double phi, double psi){
+
+  double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
+  double q_z;
+  double cos_val_b3, cos_val_b2, cos_val_b1;
+  double a1_dot_q, a2_dot_q,a3_dot_q;
+  double answer;
+  double Zq, Fkq, Fkq_2;
+
+  //convert to q and make scaled values
+  double q = sqrt(qx*qx+qy*qy);
+  double q_x = qx/q;
+  double q_y = qy/q;
+
+  //convert angle degree to radian
+  theta = theta * M_PI_180;
+  phi = phi * M_PI_180;
+  psi = psi * M_PI_180;
+
+  const double Da = d_factor*dnn;
+  const double qDa_2 = q*q*Da*Da;
+  const double s1 = dnn/sqrt(0.75);
+
+
+  //the occupied volume of the lattice
+  const double latticescale = 2.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1);
+  // q vector
+  q_z = 0.0; // for SANS; assuming qz is negligible
+  /// Angles here are respect to detector coordinate
+  ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
+    // b3 axis orientation
+    b3_x = cos(theta) * cos(phi);
+    b3_y = sin(theta);
+    //b3_z = -cos(theta) * sin(phi);
+    cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z;
+
+    //alpha = acos(cos_val_b3);
+    // b1 axis orientation
+    b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
+    b1_y = sin(psi)*cos(theta);
+    cos_val_b1 = b1_x*q_x + b1_y*q_y;
+    // b2 axis orientation
+    b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
+  	b2_y = cos(theta)*cos(psi);
+    cos_val_b2 = b2_x*q_x + b2_y*q_y;
+
+    // The following test should always pass
+    if (fabs(cos_val_b3)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b3 = 1.0;
+    }
+    if (fabs(cos_val_b2)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b2 = 1.0;
+    }
+    if (fabs(cos_val_b1)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b1 = 1.0;
+    }
+    // Compute the angle btw vector q and the a3 axis
+    a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3);
+
+    // a1 axis
+    a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1);
+
+    // a2 axis
+    a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2);
+
+
+    // Get Fkq and Fkq_2
+    Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
+    Fkq_2 = Fkq*Fkq;
+    // Call Zq=Z1*Z2*Z3
+    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2);
+
+  // Use SphereForm directly from libigor
+  answer = _sphereform(q,radius,sld,solvent_sld)*Zq*latticescale;
+
+  return answer;
+ }
Index: sasmodels/models/fcc.py
===================================================================
--- sasmodels/models/fcc.py	(revision 3271e207340eb2fd5f01b1d8ba963ab4fe9ede8c)
+++ sasmodels/models/fcc.py	(revision 3271e207340eb2fd5f01b1d8ba963ab4fe9ede8c)
@@ -0,0 +1,123 @@
+#bcc paracrystal model
+#note model title and parameter table are automatically inserted
+#note - calculation requires double precision
+"""
+Calculates the scattering from a **body-centered cubic lattice** with paracrystalline distortion. Thermal vibrations
+are considered to be negligible, and the size of the paracrystal is infinitely large. Paracrystalline distortion is
+assumed to be isotropic and characterized by a Gaussian distribution.
+
+The returned value is scaled to units of |cm^-1|\ |sr^-1|, absolute scale.
+
+Definition
+----------
+
+The scattering intensity *I(q)* is calculated as
+
+.. image:: img/image167.jpg
+
+where *scale* is the volume fraction of spheres, *Vp* is the volume of the primary particle, *V(lattice)* is a volume
+correction for the crystal structure, *P(q)* is the form factor of the sphere (normalized), and *Z(q)* is the
+paracrystalline structure factor for a body-centered cubic structure.
+
+Equation (1) of the 1990 reference is used to calculate *Z(q)*, using equations (29)-(31) from the 1987 paper for
+*Z1*\ , *Z2*\ , and *Z3*\ .
+
+The lattice correction (the occupied volume of the lattice) for a body-centered cubic structure of particles of radius
+*R* and nearest neighbor separation *D* is
+
+.. image:: img/image159.jpg
+
+The distortion factor (one standard deviation) of the paracrystal is included in the calculation of *Z(q)*
+
+.. image:: img/image160.jpg
+
+where *g* is a fractional distortion based on the nearest neighbor distance.
+
+The body-centered cubic lattice is
+
+.. image:: img/image168.jpg
+
+For a crystal, diffraction peaks appear at reduced q-values given by
+
+.. image:: img/image162.jpg
+
+where for a body-centered cubic lattice, only reflections where (\ *h* + *k* + *l*\ ) = even are allowed and
+reflections where (\ *h* + *k* + *l*\ ) = odd are forbidden. Thus the peak positions correspond to (just the first 5)
+
+.. image:: img/image169.jpg
+
+**NB: The calculation of** *Z(q)* **is a double numerical integral that must be carried out with a high density of**
+**points to properly capture the sharp peaks of the paracrystalline scattering.** So be warned that the calculation is
+SLOW. Go get some coffee. Fitting of any experimental data must be resolution smeared for any meaningful fit. This
+makes a triple integral. Very, very slow. Go get lunch!
+
+This example dataset is produced using 200 data points, *qmin* = 0.001 |Ang^-1|, *qmax* = 0.1 |Ang^-1| and the above
+default values.
+
+.. image:: img/image170.jpg
+
+*Figure. 1D plot in the linear scale using the default values (w/200 data point).*
+
+The 2D (Anisotropic model) is based on the reference below where *I(q)* is approximated for 1d scattering. Thus the
+scattering pattern for 2D may not be accurate. Note that we are not responsible for any incorrectness of the 2D model
+computation.
+
+.. image:: img/image165.gif
+
+.. image:: img/image171.jpg
+
+*Figure. 2D plot using the default values (w/200X200 pixels).*
+
+REFERENCE
+
+Hideki Matsuoka et. al. *Physical Review B*, 36 (1987) 1754-1765
+(Original Paper)
+
+Hideki Matsuoka et. al. *Physical Review B*, 41 (1990) 3854 -3856
+(Corrections to FCC and BCC lattice structure calculation)
+"""
+
+from numpy import pi, inf
+
+name = "BCCparacrystal"
+title = "Body-centred cubic lattic with paracrystalline distortion"
+description = """
+    Calculates the scattering from a **body-centered cubic lattice** with paracrystalline distortion. Thermal vibrations
+    are considered to be negligible, and the size of the paracrystal is infinitely large. Paracrystalline distortion is
+    assumed to be isotropic and characterized by a Gaussian distribution.
+    """
+
+parameters = [
+#   [ "name", "units", default, [lower, upper], "type","description" ],
+    [ "dnn", "Ang", 220, [-inf,inf],"","Nearest neighbour distance"],
+    [ "d_factor", "", 0.06,[-inf,inf],"","Paracrystal distortion factor" ],
+    [ "radius", "Ang",  40, [0, inf], "volume","Particle radius" ],
+    [ "sld", "1e-6/Ang^2", 4, [-inf,inf], "", "Particle scattering length density" ],
+    [ "solvent_sld", "1e-6/Ang^2", 1, [-inf,inf], "","Solvent scattering length density" ],
+    [ "theta", "degrees", 60, [-inf, inf], "orientation","In plane angle" ],
+    [ "phi", "degrees", 60, [-inf, inf], "orientation","Out of plane angle" ],
+    [ "psi", "degrees", 60, [-inf,inf], "orientation","Out of plane angle"]
+    ]
+
+source = [ "lib/J1.c", "lib/gauss150.c", "bcc.c" ]
+
+def ER(radius, length):
+    return 0
+
+# parameters for demo
+demo = dict(
+    scale=1, background=0,
+    dnn=220, d_factor=0.06, sld=4, solvent_sld=1,
+    radius=40,
+    theta=60, phi=60, psi=60,
+    radius_pd=.2, radius_pd_n=0.2,
+    theta_pd=15, theta_pd_n=0,
+    phi_pd=15, phi_pd_n=0,
+    psi_pd=15, psi_pd_n=0,
+    )
+
+# For testing against the old sasview models, include the converted parameter
+# names and the target sasview model name.
+oldname='BCCrystalModel'
+oldpars=dict(sld='sldSph',
+             solvent_sld='sldSolv')
Index: sasmodels/models/parallelepiped.c
===================================================================
--- sasmodels/models/parallelepiped.c	(revision c5b7d07435739547137368520d072f6e87ac3208)
+++ sasmodels/models/parallelepiped.c	(revision c5b7d07435739547137368520d072f6e87ac3208)
@@ -0,0 +1,162 @@
+double form_volume(double a_side, double b_side, double c_side);
+double Iq(double q, double sld, double solvent_sld, double a_side, double b_side, double c_side);
+double Iqxy(double qx, double qy, double sld, double solvent_sld,
+    double a_side, double b_side, double c_side, double theta, double phi, double psi);
+
+// From Igor library
+double _pkernel(double a, double b,double c, double ala, double alb, double alc);
+double _pkernel(double a, double b,double c, double ala, double alb, double alc){
+    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;
+    if(argA==0.0) {
+        tmp1 = 1.0;
+    } else {
+        tmp1 = sin(argA)*sin(argA)/argA/argA;
+    }
+    if (argB==0.0) {
+        tmp2 = 1.0;
+    } else {
+        tmp2 = sin(argB)*sin(argB)/argB/argB;
+    }
+    if (argC==0.0) {
+        tmp3 = 1.0;
+    } else {
+        tmp3 = sin(argC)*sin(argC)/argC/argC;
+    }
+    return (tmp1*tmp2*tmp3);
+
+}
+
+double form_volume(double a_side, double b_side, double c_side)
+{
+    return a_side * b_side * c_side;
+}
+
+double Iq(double q,
+    double sld,
+    double solvent_sld,
+    double a_side,
+    double b_side,
+    double c_side)
+{
+    double tmp1, tmp2, yyy;
+    
+    double mu = q * b_side;
+    
+    // Scale sides by B
+    double a_scaled = a_side / b_side;
+    double c_scaled = c_side / b_side;
+    
+    // outer integral (with 76 gauss points), integration limits = 0, 1
+    double summ = 0; //initialize integral
+
+    for( int i=0; i<76; i++) {
+		
+        // inner integral (with 76 gauss points), integration limits = 0, 1
+	
+        double summj = 0.0;
+	double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );		
+		
+	for(int j=0; j<76; 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);
+            if(arg1==0.0) {
+	        tmp1 = 1.0;
+            } else {
+                tmp1 = sin(arg1)*sin(arg1)/arg1/arg1;
+            }
+            if (arg2==0.0) {
+                tmp2 = 1.0;
+            } else {
+                tmp2 = sin(arg2)*sin(arg2)/arg2/arg2;
+            }
+            yyy = Gauss76Wt[j] * tmp1 * tmp2;
+
+            summj += yyy;
+        }
+		
+        // value of the inner integral
+        double answer = 0.5 * summj;
+
+        double arg = 0.5 * mu * c_scaled * sigma;
+        if ( arg == 0.0 ) {
+            answer *= 1.0;
+        } else {
+            answer *= sin(arg)*sin(arg)/arg/arg;
+        }
+		
+	// sum of outer integral
+	yyy = Gauss76Wt[i] * answer;
+        summ += yyy;
+        
+    }	
+   
+    const double vd = (sld-solvent_sld) * form_volume(a_side, b_side, c_side);
+    return 1.0e-4 * 0.5 * vd * vd * summ;
+    
+}
+
+
+double Iqxy(double qx, double qy,
+    double sld,
+    double solvent_sld,
+    double a_side,
+    double b_side,
+    double c_side,
+    double theta,
+    double phi,
+    double psi)
+{
+    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;
+    }
+    
+    // Call the IGOR library function to get the kernel
+    double form = _pkernel( q*a_side, q*b_side, q*c_side, cos_val_a, cos_val_b, cos_val_c);
+  
+    // Multiply by contrast^2
+    const double vd = (sld - solvent_sld) * form_volume(a_side, b_side, c_side);
+    return 1.0e-4 * vd * vd * form;
+}
Index: sasmodels/models/parallelepiped.py
===================================================================
--- sasmodels/models/parallelepiped.py	(revision c5b7d07435739547137368520d072f6e87ac3208)
+++ sasmodels/models/parallelepiped.py	(revision c5b7d07435739547137368520d072f6e87ac3208)
@@ -0,0 +1,174 @@
+# parallelepiped model
+# Note: model title and parameter table are inserted automatically
+r"""
+The form factor is normalized by the particle volume.
+
+For information about polarised and magnetic scattering, click here_.
+
+Definition
+----------
+
+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 the 
+RectangularPrismModel_.
+
+The calculated form factor is:
+
+.. math::
+
+    P(Q) = {\text{scale} \over V} F^2(Q) + \text{background}
+
+where the volume *V* = *A B C* and the averaging < > is applied over all 
+orientations for 1D.
+
+.. image:: img/parallelepiped.jpg
+
+*Figure. 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{\textstyle{scale}}{V}\int_0^1 \phi(\mu \sqrt{1-\sigma^2},a) 
+    [S(\mu c \sigma/2)]^2 d\sigma
+
+with
+
+.. math::
+
+    \phi(\mu,a) = \int_0^1 \{S[\frac{\mu}{2}\cos(\frac{\pi}{2}u)] 
+    S[\frac{\mu a}{2}\sin(\frac{\pi}{2}u)]\}^2 du
+
+    S(x) = \frac{\sin x}{x}
+    
+    \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|; 
+ie, *I(q)* = |phi| *P(q)*\ .
+
+NB: The 2nd virial coefficient of the parallelpiped is calculated based on 
+the averaged effective radius (= sqrt(*short_a* \* *short_b* / |pi|)) and 
+length(= *long_c*) values, and used as the effective radius for
+*S(Q)* when *P(Q)* \* *S(Q)* is applied.
+
+To provide easy access to the orientation of the parallelepiped, we define 
+three angles |theta|, |phi| and |bigpsi|. The definition of |theta| and |phi|
+is the same as for the cylinder model (see also figures below). 
+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.
+
+
+.. _parallelepiped-orientation:
+
+.. figure:: img/orientation.jpg
+
+    Definition of the angles for oriented parallelepipeds.
+
+.. figure:: img/orientation2.jpg
+
+    Examples of the angles for oriented parallelepipeds against the detector plane.
+
+
+Validation
+----------
+
+Validation of the code was done by comparing the output of the 1D calculation 
+to the angular average of the output of a 2D calculation over all possible 
+angles. The Figure below shows the comparison where the solid dot refers to 
+averaged 2D while the line represents the result of the 1D calculation (for 
+the averaging, 76, 180, 76 points are taken for the angles of |theta|, |phi|, 
+and |psi| respectively).
+
+.. _parallelepiped-compare:
+
+.. figure:: img/parallelepiped_compare.jpg
+
+*Figure. Comparison between 1D and averaged 2D.*
+
+This model reimplements the form factor calculations implemented in a c-library 
+provided by the NIST Center for Neutron Research (Kline, 2006).
+
+"""
+
+from numpy import pi, inf
+
+name = "parallelepiped"
+title = "Rectangular parallelepiped with uniform scattering length density."
+description = """
+     P(q)= scale/V*integral from 0 to 1 of ...
+           phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma
+     
+            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(x) = sin(x)/x
+	    mu = q*B
+"""
+
+parameters = [
+#   [ "name", "units", default, [lower, upper], "type",
+#     "description" ],
+    [ "sld", "6e-6/Ang^2", 4, [-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" ],
+    [ "b_side", "Ang",  75, [0, inf], "volume",
+      "Second side of the parallelepiped" ],
+    [ "c_side", "Ang",  400, [0, inf], "volume",
+      "Larger side of the parallelepiped" ],
+    [ "theta", "degrees", 60, [-inf, inf], "orientation",
+      "In plane angle" ],
+    [ "phi", "degrees", 60, [-inf, inf], "orientation",
+      "Out of plane angle" ],
+    [ "psi", "degrees", 60, [-inf, inf], "orientation",
+      "Rotation angle around its own c axis against q plane" ],
+    ]
+
+source = [ "lib/J1.c", "lib/gauss76.c", "parallelepiped.c" ]
+
+def ER(a_side, b_side, c_side):
+
+    # 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
+
+    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, b_side=75, c_side=400,
+    theta=45, phi=30, psi=15,
+    a_side_pd=0.1, a_side_pd_n=10,
+    b_side_pd=0.1, b_side_pd_n=1,
+    c_side_pd=0.1, c_side_pd_n=10,
+    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='ParallelepipedModel'
+oldpars=dict(theta='parallel_theta', phi='parallel_phi', psi='parallel_psi',
+             a_side='short_a', b_side='short_b', c_side='long_c',
+             sld='sldPipe', solvent_sld='sldSolv')
+
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision 250fa25ad675ea9a85cddbc91dd32cd7a1255fb0)
+++ sasmodels/sasview_model.py	(revision 0a82216fe41763d3d9ab2277f58228a0921edcce)
@@ -32,9 +32,12 @@
 
 
-def make_class(kernel_module, dtype='single'):
+def make_class(kernel_module, dtype='single', namestyle='name'):
     """
     Load the sasview model defined in *kernel_module*.
 
     Returns a class that can be used directly as a sasview model.
+
+    Defaults to using the new name for a model. Setting namestyle='name'
+    will produce a class with a name compatible with SasView
     """
     model =  load_model(kernel_module, dtype=dtype)
@@ -42,5 +45,5 @@
         SasviewModel.__init__(self, model)
     attrs = dict(__init__=__init__)
-    ConstructedModel = type(model.info['name'],  (SasviewModel,), attrs)
+    ConstructedModel = type(model.info[namestyle],  (SasviewModel,), attrs)
     return ConstructedModel
 
@@ -54,4 +57,5 @@
 
         self.name = model.info['name']
+        self.oldname = model.info['oldname']
         self.description = model.info['description']
         self.category = None
