Index: sansmodels/src/sans/models/CSParallelepipedModel.py
===================================================================
--- sansmodels/src/sans/models/CSParallelepipedModel.py	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
+++ sansmodels/src/sans/models/CSParallelepipedModel.py	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
@@ -0,0 +1,191 @@
+#!/usr/bin/env python
+
+##############################################################################
+#	This software was developed by the University of Tennessee as part of the
+#	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+#	project funded by the US National Science Foundation.
+#
+#	If you use DANSE applications to do scientific research that leads to
+#	publication, we ask that you acknowledge the use of the software with the
+#	following sentence:
+#
+#	"This work benefited from DANSE software developed under NSF award DMR-0520547."
+#
+#	copyright 2008, University of Tennessee
+##############################################################################
+
+
+""" 
+Provide functionality for a C extension model
+
+:WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+         DO NOT MODIFY THIS FILE, MODIFY ..\c_extensions\csparallelepiped.h
+         AND RE-RUN THE GENERATOR SCRIPT
+
+"""
+
+from sans.models.BaseComponent import BaseComponent
+from sans_extension.c_models import CCSParallelepipedModel
+import copy    
+    
+class CSParallelepipedModel(CCSParallelepipedModel, BaseComponent):
+    """ 
+    Class that evaluates a CSParallelepipedModel model. 
+    This file was auto-generated from ..\c_extensions\csparallelepiped.h.
+    Refer to that file and the structure it contains
+    for details of the model.
+    List of default parameters:
+         scale           = 1.0 
+         shortA          = 35.0 [A]
+         midB            = 75.0 [A]
+         longC           = 400.0 [A]
+         rimA            = 10.0 [A]
+         rimB            = 10.0 [A]
+         rimC            = 10.0 [A]
+         sld_rimA        = 2e-006 [1/A^(2)]
+         sld_rimB        = 4e-006 [1/A^(2)]
+         sld_rimC        = 2e-006 [1/A^(2)]
+         sld_pcore       = 1e-006 [1/A^(2)]
+         sld_solv        = 6e-006 [1/A^(2)]
+         background      = 0.06 [1/cm]
+         parallel_theta  = 0.0 [rad]
+         parallel_phi    = 0.0 [rad]
+         parallel_psi    = 0.0 [rad]
+
+    """
+        
+    def __init__(self):
+        """ Initialization """
+        
+        # Initialize BaseComponent first, then sphere
+        BaseComponent.__init__(self)
+        CCSParallelepipedModel.__init__(self)
+        
+        ## Name of the model
+        self.name = "CSParallelepipedModel"
+        ## Model description
+        self.description =""" Form factor for a rectangular Shell. Below are the Parameters.
+		scale: scale factor
+		shortA: length of short edge  [A]
+		midB: length of another short edge [A]
+		longC: length of long edge  of the parallelepiped [A]
+		rimA: length of short edge  [A]
+		rimB: length of another short edge [A]
+		rimC: length of long edge  of the parallelepiped [A]
+		sld_rimA: sld of rimA [1/A^(2)]
+		sld_rimB: sld of rimB [1/A^(2)]
+		sld_rimC: sld of rimC [1/A^(2)]
+		sld_core: Pipe_sld [1/A^(2)]
+		sld_solv: solvent_sld [1/A^(2)]
+		background: incoherent Background [1/cm]"""
+       
+        ## Parameter details [units, min, max]
+        self.details = {}
+        self.details['scale'] = ['', None, None]
+        self.details['shortA'] = ['[A]', None, None]
+        self.details['midB'] = ['[A]', None, None]
+        self.details['longC'] = ['[A]', None, None]
+        self.details['rimA'] = ['[A]', None, None]
+        self.details['rimB'] = ['[A]', None, None]
+        self.details['rimC'] = ['[A]', None, None]
+        self.details['sld_rimA'] = ['[1/A^(2)]', None, None]
+        self.details['sld_rimB'] = ['[1/A^(2)]', None, None]
+        self.details['sld_rimC'] = ['[1/A^(2)]', None, None]
+        self.details['sld_pcore'] = ['[1/A^(2)]', None, None]
+        self.details['sld_solv'] = ['[1/A^(2)]', None, None]
+        self.details['background'] = ['[1/cm]', None, None]
+        self.details['parallel_theta'] = ['[rad]', None, None]
+        self.details['parallel_phi'] = ['[rad]', None, None]
+        self.details['parallel_psi'] = ['[rad]', None, None]
+
+        ## fittable parameters
+        self.fixed=['shortA.width', 'midB.width', 'longC.width', 'parallel_phi.width', 'parallel_psi.width', 'parallel_theta.width']
+        
+        ## non-fittable parameters
+        self.non_fittable=[]
+        
+        ## parameters with orientation
+        self.orientation_params =['parallel_phi', 'parallel_psi', 'parallel_theta', 'parallel_phi.width', 'parallel_psi.width', 'parallel_theta.width']
+   
+    def clone(self):
+        """ Return a identical copy of self """
+        return self._clone(CSParallelepipedModel())   
+        
+    def __getstate__(self):
+        """
+        return object state for pickling and copying
+        """
+        model_state = {'params': self.params, 'dispersion': self.dispersion, 'log': self.log}
+        
+        return self.__dict__, model_state
+        
+    def __setstate__(self, state):
+        """
+        create object from pickled state
+        
+        :param state: the state of the current model
+        
+        """
+        
+        self.__dict__, model_state = state
+        self.params = model_state['params']
+        self.dispersion = model_state['dispersion']
+        self.log = model_state['log']
+       	
+   
+    def run(self, x=0.0):
+        """ 
+        Evaluate the model
+        
+        :param x: input q, or [q,phi]
+        
+        :return: scattering function P(q)
+        
+        """
+        
+        return CCSParallelepipedModel.run(self, x)
+   
+    def runXY(self, x=0.0):
+        """ 
+        Evaluate the model in cartesian coordinates
+        
+        :param x: input q, or [qx, qy]
+        
+        :return: scattering function P(q)
+        
+        """
+        
+        return CCSParallelepipedModel.runXY(self, x)
+        
+    def evalDistribution(self, x=[]):
+        """ 
+        Evaluate the model in cartesian coordinates
+        
+        :param x: input q[], or [qx[], qy[]]
+        
+        :return: scattering function P(q[])
+        
+        """
+        return CCSParallelepipedModel.evalDistribution(self, x)
+        
+    def calculate_ER(self):
+        """ 
+        Calculate the effective radius for P(q)*S(q)
+        
+        :return: the value of the effective radius
+        
+        """       
+        return CCSParallelepipedModel.calculate_ER(self)
+        
+    def set_dispersion(self, parameter, dispersion):
+        """
+        Set the dispersion object for a model parameter
+        
+        :param parameter: name of the parameter [string]
+        :param dispersion: dispersion object of type DispersionModel
+        
+        """
+        return CCSParallelepipedModel.set_dispersion(self, parameter, dispersion.cdisp)
+        
+   
+# End of file
Index: sansmodels/src/sans/models/CappedCylinderModel.py
===================================================================
--- sansmodels/src/sans/models/CappedCylinderModel.py	(revision 339ce675656e3ab89eff7b84d8292a2e0ffe6109)
+++ sansmodels/src/sans/models/CappedCylinderModel.py	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
@@ -58,6 +58,6 @@
         self.name = "CappedCylinderModel"
         ## Model description
-        self.description ="""Calculates the scattering from a cylinder with spherical
-		section end-caps. That is, a sphereocylinder
+        self.description ="""Calculates the scattering from a cylinder with spherical section end-caps.
+		That is, a sphereocylinder
 		with end caps that have a radius larger than
 		that of the cylinder and the center of the
Index: sansmodels/src/sans/models/c_extensions/bcc.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/bcc.c	(revision 8f20419d5be1d467bd17c7dfade28c96f90656ca)
+++ sansmodels/src/sans/models/c_extensions/bcc.c	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
@@ -110,14 +110,14 @@
     // Compute the angle btw vector q and the a3 axis
     cos_val_a3 = a3_x*q_x + a3_y*q_y + a3_z*q_z;
-    //alpha = acos(cos_val_a3);
+    alpha = acos(cos_val_a3);
     a3_dot_q = aa*q*cos_val_a3;
 
     // a1 axis
     cos_val_a1 = a1_x*q_x + a1_y*q_y;
-    a1_dot_q = aa*q*cos_val_a1;
+    a1_dot_q = aa*q*cos_val_a1*sin(alpha);
 
     // a2 axis
-    cos_val_a2 = a2_x*q_x + a2_y*q_y; //sin(acos(cos_val_a1)) ;
-    a2_dot_q = aa*q*cos_val_a2;
+    cos_val_a2 = sin(acos(cos_val_a1));//a2_x*q_x + a2_y*q_y;
+    a2_dot_q = aa*q*cos_val_a2*sin(alpha); //aa*q*cos_val_a2
 
     // The following test should always pass
Index: sansmodels/src/sans/models/c_extensions/capcyl.h
===================================================================
--- sansmodels/src/sans/models/c_extensions/capcyl.h	(revision 339ce675656e3ab89eff7b84d8292a2e0ffe6109)
+++ sansmodels/src/sans/models/c_extensions/capcyl.h	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
@@ -7,6 +7,6 @@
  //[PYTHONCLASS] = CappedCylinderModel
  //[DISP_PARAMS] = rad_cyl,len_cyl,rad_cap,phi,  theta
- //[DESCRIPTION] =<text>Calculates the scattering from a cylinder with spherical
- //				section end-caps. That is, a sphereocylinder
+ //[DESCRIPTION] =<text>Calculates the scattering from a cylinder with spherical section end-caps.
+ //				That is, a sphereocylinder
  //				with end caps that have a radius larger than
  //				that of the cylinder and the center of the
Index: sansmodels/src/sans/models/c_extensions/csparallelepiped.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/csparallelepiped.c	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
+++ sansmodels/src/sans/models/c_extensions/csparallelepiped.c	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
@@ -0,0 +1,239 @@
+/**
+ * Scattering model for a csparallelepiped
+ */
+
+#include "csparallelepiped.h"
+#include <math.h>
+#include "libCylinder.h"
+#include <stdio.h>
+#include <stdlib.h>
+
+
+/**
+ * Function to evaluate 1D scattering function
+ * @param pars: parameters of the CSparallelepiped
+ * @param q: q-value
+ * @return: function value
+ */
+double csparallelepiped_analytical_1D(CSParallelepipedParameters *pars, double q) {
+	double dp[13];
+
+	// Fill paramater array
+	dp[0] = pars->scale;
+	dp[1] = pars->shortA;
+	dp[2] = pars->midB;
+	dp[3] = pars->longC;
+	dp[4] = pars->rimA;
+	dp[5] = pars->rimB;
+	dp[6] = pars->rimC;
+	dp[7] = pars->sld_rimA;
+	dp[8] = pars->sld_rimB;
+	dp[9] = pars->sld_rimC;
+	dp[10] = pars->sld_pcore;
+	dp[11] = pars->sld_solv;
+	dp[12] = pars->background;
+
+	// Call library function to evaluate model
+	//ToDo: Correct this 1d model, CSParallelepiped in libigor (2D corrected).
+	return CSParallelepiped(dp, q);
+}
+
+
+double cspkernel(double dp[],double q, double ala, double alb, double alc){
+    // mu passed in is really mu*sqrt(1-sig^2)
+    double argA,argB,argC,argtA,argtB,argtC,tmp1,tmp2,tmp3,tmpt1,tmpt2,tmpt3;			//local variables
+
+	double aa,bb,cc, ta,tb,tc;
+	double Vin,Vot,V1,V2,V3;
+	double rhoA,rhoB,rhoC, rhoP, rhosolv;
+	double dr0, drA,drB, drC;
+	double Pi,retVal;
+
+	aa = dp[1];
+	bb = dp[2];
+	cc = dp[3];
+	ta = dp[4];
+	tb = dp[5];
+	tc = dp[6];
+	rhoA=dp[7];
+	rhoB=dp[8];
+	rhoC=dp[9];
+	rhoP=dp[10];
+	rhosolv=dp[11];
+	dr0=rhoP-rhosolv;
+	drA=rhoA-rhosolv;
+	drB=rhoB-rhosolv;
+	drC=rhoC-rhosolv;
+	Vin=(aa*bb*cc);
+	Vot=(aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc);
+	V1=(2.0*ta*bb*cc);   //  incorrect V1 (aa*bb*cc+2*ta*bb*cc)
+	V2=(2.0*aa*tb*cc);  // incorrect V2(aa*bb*cc+2*aa*tb*cc)
+	V3=(2.0*aa*bb*tc);
+	//aa = aa/bb;
+	ta=(aa+2.0*ta);///bb;
+	tb=(aa+2.0*tb);///bb;
+	tc=(aa+2.0*tc);
+    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
+    argA = q*aa*ala/2.0;
+    argB = q*bb*alb/2.0;
+    argC = q*cc*alc/2.0;
+    argtA = q*ta*ala/2.0;
+	argtB = q*tb*alb/2.0;
+	argtC = q*tc*alc/2.0;
+
+    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;
+    }
+    // This expression is different from NIST/IGOR package (I strongly believe the IGOR is wrong!!!). 10/15/2010.
+    retVal =( 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
+    //retVal *= (tmp3*tmp3);
+    retVal /= Vot;
+
+    return (retVal);
+
+}//Function cspkernel()
+
+
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the CSparallelepiped
+ * @param q: q-value
+ * @return: function value
+ */
+double csparallelepiped_analytical_2DXY(CSParallelepipedParameters *pars, double qx, double qy) {
+	double q;
+	q = sqrt(qx*qx+qy*qy);
+    return csparallelepiped_analytical_2D_scaled(pars, q, qx/q, qy/q);
+}
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the CSParallelepiped
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double csparallelepiped_analytical_2D(CSParallelepipedParameters *pars, double q, double phi) {
+    return csparallelepiped_analytical_2D_scaled(pars, q, cos(phi), sin(phi));
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the CSparallelepiped
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+double csparallelepiped_analytical_2D_scaled(CSParallelepipedParameters *pars, double q, double q_x, double q_y) {
+	double dp[13];
+
+	// Fill paramater array
+	dp[0] = 1.0;
+	dp[1] = pars->shortA;
+	dp[2] = pars->midB;
+	dp[3] = pars->longC;
+	dp[4] = pars->rimA;
+	dp[5] = pars->rimB;
+	dp[6] = pars->rimC;
+	dp[7] = pars->sld_rimA;
+	dp[8] = pars->sld_rimB;
+	dp[9] = pars->sld_rimC;
+	dp[10] = pars->sld_pcore;
+	dp[11] = pars->sld_solv;
+	dp[12] = 0.0;
+
+	double cparallel_x, cparallel_y, cparallel_z, bparallel_x, bparallel_y, parallel_x, parallel_y, parallel_z;
+	double q_z;
+	double alpha, vol, cos_val_c, cos_val_b, cos_val_a, edgeA, edgeB, edgeC;
+
+	double answer;
+	double pi = 4.0*atan(1.0);
+
+	edgeA = pars->shortA;
+	edgeB = pars->midB;
+	edgeC = pars->longC;
+
+
+    // parallelepiped c axis orientation
+    cparallel_x = sin(pars->parallel_theta) * cos(pars->parallel_phi);
+    cparallel_y = sin(pars->parallel_theta) * sin(pars->parallel_phi);
+    cparallel_z = cos(pars->parallel_theta);
+
+    // q vector
+    q_z = 0.0;
+
+    // Compute the angle btw vector q and the
+    // axis of the parallelepiped
+    cos_val_c = cparallel_x*q_x + cparallel_y*q_y + cparallel_z*q_z;
+    alpha = acos(cos_val_c);
+
+    // parallelepiped a axis orientation
+    parallel_x = sin(pars->parallel_psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)*sin(pars->parallel_psi);
+    parallel_y = cos(pars->parallel_psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)*cos(pars->parallel_psi);
+
+    cos_val_a = parallel_x*q_x + parallel_y*q_y;
+
+
+
+    // parallelepiped b axis orientation
+    bparallel_x = sqrt(1.0-sin(pars->parallel_theta)*cos(pars->parallel_phi))*cos(pars->parallel_psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)* cos(pars->parallel_psi);
+    bparallel_y = sqrt(1.0-sin(pars->parallel_theta)*cos(pars->parallel_phi))*sin(pars->parallel_psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)* sin(pars->parallel_psi);
+    // axis of the parallelepiped
+    cos_val_b = sin(acos(cos_val_a)) ;
+
+
+
+    // The following test should always pass
+    if (fabs(cos_val_c)>1.0) {
+    	printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+     	return 0;
+    }
+
+	// Call the IGOR library function to get the kernel
+	answer = cspkernel( dp,q, sin(alpha)*cos_val_a,sin(alpha)*cos_val_b,cos_val_c);
+
+	//convert to [cm-1]
+	answer *= 1.0e8;
+
+	//Scale
+	answer *= pars->scale;
+
+	// add in the background
+	answer += pars->background;
+
+	return answer;
+}
+
Index: sansmodels/src/sans/models/c_extensions/csparallelepiped.h
===================================================================
--- sansmodels/src/sans/models/c_extensions/csparallelepiped.h	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
+++ sansmodels/src/sans/models/c_extensions/csparallelepiped.h	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
@@ -0,0 +1,92 @@
+/*
+	TODO: Add 2D model
+*/
+
+#if !defined(csparallelepiped_h)
+#define csparallelepiped_h
+/** Structure definition for CSParallelepiped parameters
+ * [PYTHONCLASS] = CSParallelepipedModel
+ * [DISP_PARAMS] = shortA, midB, longC,parallel_phi,parallel_psi, parallel_theta
+   [DESCRIPTION] = <text> Form factor for a rectangular Shell. Below are the Parameters.
+		scale: scale factor
+		shortA: length of short edge  [A]
+		midB: length of another short edge [A]
+		longC: length of long edge  of the parallelepiped [A]
+		rimA: length of short edge  [A]
+		rimB: length of another short edge [A]
+		rimC: length of long edge  of the parallelepiped [A]
+		sld_rimA: sld of rimA [1/A^(2)]
+		sld_rimB: sld of rimB [1/A^(2)]
+		sld_rimC: sld of rimC [1/A^(2)]
+		sld_core: Pipe_sld [1/A^(2)]
+		sld_solv: solvent_sld [1/A^(2)]
+		background: incoherent Background [1/cm]
+		</text>
+	[FIXED]= <text>shortA.width; midB.width; longC.width;parallel_phi.width;parallel_psi.width; parallel_theta.width</text>
+	[ORIENTATION_PARAMS]= <text>parallel_phi;parallel_psi; parallel_theta; parallel_phi.width;parallel_psi.width; parallel_theta.width</text>
+
+
+ **/
+typedef struct {
+    /// Scale factor
+    //  [DEFAULT]=scale=1.0
+    double scale;
+    ///  Length of short edge of the parallelepiped [A]
+    //  [DEFAULT]=shortA=35 [A]
+    double shortA;
+	/// Length of mid edge of the parallelepiped [A]
+    //  [DEFAULT]=midB=75 [A]
+    double midB;
+	/// Length of long edge of the parallelepiped [A]
+    //  [DEFAULT]=longC=400 [A]
+    double longC;
+    ///  Thickness of rimA [A]
+    //  [DEFAULT]=rimA=10 [A]
+    double rimA;
+	/// Thickness of rimB [A] [A]
+    //  [DEFAULT]=rimB=10 [A]
+    double rimB;
+	/// Thickness of rimC [A] [A]
+    //  [DEFAULT]=rimC=10 [A]
+    double rimC;
+    /// SLD_rimA [1/A^(2)]
+    //  [DEFAULT]=sld_rimA=2e-6 [1/A^(2)]
+    double sld_rimA;
+    /// SLD_rimB [1/A^(2)]
+    //  [DEFAULT]=sld_rimB=4e-6 [1/A^(2)]
+    double sld_rimB;
+    /// SLD_rimC [1/A^(2)]
+    //  [DEFAULT]=sld_rimC=2e-6 [1/A^(2)]
+    double sld_rimC;
+    /// SLD_pcore [1/A^(2)]
+    //  [DEFAULT]=sld_pcore=1e-6 [1/A^(2)]
+    double sld_pcore;
+    /// sld_solv [1/A^(2)]
+    //  [DEFAULT]=sld_solv=6e-6 [1/A^(2)]
+    double sld_solv;
+	/// Incoherent Background [1/cm]
+	//  [DEFAULT]=background=0.06 [1/cm]
+	double background;
+    /// Orientation of the parallelepiped axis w/respect incoming beam [rad]
+    //  [DEFAULT]=parallel_theta=0.0 [rad]
+    double parallel_theta;
+    /// Orientation of the longitudinal axis of the parallelepiped in the plane of the detector [rad]
+    //  [DEFAULT]=parallel_phi=0.0 [rad]
+    double parallel_phi;
+    /// Orientation of the cross-sectional minor axis of the parallelepiped in the plane of the detector [rad]
+    //  [DEFAULT]=parallel_psi=0.0 [rad]
+    double parallel_psi;
+
+
+} CSParallelepipedParameters;
+
+
+
+/// 1D scattering function
+double csparallelepiped_analytical_1D(CSParallelepipedParameters *pars, double q);
+
+/// 2D scattering function
+double csparallelepiped_analytical_2D(CSParallelepipedParameters *pars, double q, double phi);
+double csparallelepiped_analytical_2DXY(CSParallelepipedParameters *pars, double qx, double qy);
+double csparallelepiped_analytical_2D_scaled(CSParallelepipedParameters *pars, double q, double q_x, double q_y);
+#endif
Index: sansmodels/src/sans/models/c_extensions/fcc.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/fcc.c	(revision 8f20419d5be1d467bd17c7dfade28c96f90656ca)
+++ sansmodels/src/sans/models/c_extensions/fcc.c	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
@@ -77,5 +77,5 @@
 	Da = pars->d_factor*aa;
 	qDa_2 = pow(q*Da,2.0);
-	contrast = pars->sldSph - pars->sldSolv;
+	//contrast = pars->sldSph - pars->sldSolv;
 
 	latticeScale = 4.0*(4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa*sqrt(2.0),3.0);
@@ -109,14 +109,15 @@
 
     // Compute the angle btw vector q and the a3 axis
-    cos_val_a3 = a3_x*q_x + a3_y*q_y;// + a3_z*q_z;
+    cos_val_a3 = a3_x*q_x + a3_y*q_y + a3_z*q_z;
+    alpha = acos(cos_val_a3);
     a3_dot_q = aa*q*cos_val_a3;
 
     // a1 axis
     cos_val_a1 = a1_x*q_x + a1_y*q_y;
-    a1_dot_q = aa*q*cos_val_a1;
+    a1_dot_q = aa*q*cos_val_a1*sin(alpha);
 
     // a2 axis
-    cos_val_a2 = a2_x*q_x + a2_y*q_y; //sin(acos(cos_val_a1)) ;
-    a2_dot_q = aa*q*cos_val_a2;
+    cos_val_a2 = sin(acos(cos_val_a1));//a2_x*q_x + a2_y*q_y;
+    a2_dot_q = aa*q*cos_val_a2*sin(alpha);
 
     // The following test should always pass
@@ -134,5 +135,5 @@
 
 	// Use SphereForm directly from libigor
-	answer = SphereForm_Paracrystal(pars->radius,contrast,q)*Zq;
+	answer = SphereForm(dp,q)*Zq;
 
 	//consider scales
Index: sansmodels/src/sans/models/c_extensions/sc.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/sc.c	(revision 8f20419d5be1d467bd17c7dfade28c96f90656ca)
+++ sansmodels/src/sans/models/c_extensions/sc.c	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
@@ -89,4 +89,5 @@
     // Compute the angle btw vector q and the a3 axis
     cos_val_a3 = a3_x*q_x + a3_y*q_y + a3_z*q_z;
+    alpha = acos(cos_val_a3);
     //alpha = acos(cos_val_a3);
     a3_dot_q = aa*q*cos_val_a3;
@@ -96,5 +97,5 @@
 
     cos_val_a1 = a1_x*q_x + a1_y*q_y;
-    a1_dot_q = aa*q*cos_val_a1;
+    a1_dot_q = aa*q*cos_val_a1*sin(alpha);
 
     // a2 axis orientation
@@ -102,6 +103,6 @@
     a2_y = sqrt(1.0-sin(pars->theta)*cos(pars->phi))*sin(pars->psi);
     // a2 axis
-    cos_val_a2 =  a2_x*q_x + a2_y*q_y;//sin(acos(cos_val_a1)) ;
-    a2_dot_q = aa*q*cos_val_a2;
+    cos_val_a2 =  sin(acos(cos_val_a1));//a2_x*q_x + a2_y*q_y;
+    a2_dot_q = aa*q*cos_val_a2*sin(alpha);
 
     // The following test should always pass
Index: sansmodels/src/sans/models/c_models/CCSParallelepipedModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CCSParallelepipedModel.cpp	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
+++ sansmodels/src/sans/models/c_models/CCSParallelepipedModel.cpp	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
@@ -0,0 +1,673 @@
+/**
+	This software was developed by the University of Tennessee as part of the
+	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+	project funded by the US National Science Foundation.
+
+	If you use DANSE applications to do scientific research that leads to
+	publication, we ask that you acknowledge the use of the software with the
+	following sentence:
+
+	"This work benefited from DANSE software developed under NSF award DMR-0520547."
+
+	copyright 2008, University of Tennessee
+ */
+
+/** CCSParallelepipedModel
+ *
+ * C extension 
+ *
+ * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ *          DO NOT MODIFY THIS FILE, MODIFY csparallelepiped.h
+ *          AND RE-RUN THE GENERATOR SCRIPT
+ *
+ */
+#define NO_IMPORT_ARRAY
+#define PY_ARRAY_UNIQUE_SYMBOL PyArray_API_sans
+ 
+extern "C" {
+#include <Python.h>
+#include <arrayobject.h>
+#include "structmember.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include "csparallelepiped.h"
+}
+
+#include "models.hh"
+#include "dispersion_visitor.hh"
+
+/// Error object for raised exceptions
+static PyObject * CCSParallelepipedModelError = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    /// Dispersion parameters
+    PyObject * dispersion;
+    /// Underlying model object
+    CSParallelepipedModel * model;
+    /// Log for unit testing
+    PyObject * log;
+} CCSParallelepipedModel;
+
+
+static void
+CCSParallelepipedModel_dealloc(CCSParallelepipedModel* self)
+{
+    Py_DECREF(self->params);
+    Py_DECREF(self->dispersion);
+    Py_DECREF(self->log);
+    delete self->model;
+    self->ob_type->tp_free((PyObject*)self);
+    
+
+}
+
+static PyObject *
+CCSParallelepipedModel_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CCSParallelepipedModel *self;
+    
+    self = (CCSParallelepipedModel *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+CCSParallelepipedModel_init(CCSParallelepipedModel *self, PyObject *args, PyObject *kwds)
+{
+    if (self != NULL) {
+    	
+    	// Create parameters
+        self->params = PyDict_New();
+        self->dispersion = PyDict_New();
+        self->model = new CSParallelepipedModel();
+        
+        // Initialize parameter dictionary
+        PyDict_SetItemString(self->params,"sld_rimB",Py_BuildValue("d",0.000004000000));
+        PyDict_SetItemString(self->params,"sld_rimC",Py_BuildValue("d",0.000002000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.060000000000));
+        PyDict_SetItemString(self->params,"parallel_psi",Py_BuildValue("d",0.000000000000));
+        PyDict_SetItemString(self->params,"midB",Py_BuildValue("d",75.000000000000));
+        PyDict_SetItemString(self->params,"parallel_theta",Py_BuildValue("d",0.000000000000));
+        PyDict_SetItemString(self->params,"parallel_phi",Py_BuildValue("d",0.000000000000));
+        PyDict_SetItemString(self->params,"sld_pcore",Py_BuildValue("d",0.000001000000));
+        PyDict_SetItemString(self->params,"sld_rimA",Py_BuildValue("d",0.000002000000));
+        PyDict_SetItemString(self->params,"shortA",Py_BuildValue("d",35.000000000000));
+        PyDict_SetItemString(self->params,"rimB",Py_BuildValue("d",10.000000000000));
+        PyDict_SetItemString(self->params,"rimC",Py_BuildValue("d",10.000000000000));
+        PyDict_SetItemString(self->params,"longC",Py_BuildValue("d",400.000000000000));
+        PyDict_SetItemString(self->params,"sld_solv",Py_BuildValue("d",0.000006000000));
+        PyDict_SetItemString(self->params,"rimA",Py_BuildValue("d",10.000000000000));
+        // Initialize dispersion / averaging parameter dict
+        DispersionVisitor* visitor = new DispersionVisitor();
+        PyObject * disp_dict;
+        disp_dict = PyDict_New();
+        self->model->shortA.dispersion->accept_as_source(visitor, self->model->shortA.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "shortA", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->midB.dispersion->accept_as_source(visitor, self->model->midB.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "midB", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->longC.dispersion->accept_as_source(visitor, self->model->longC.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "longC", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->parallel_phi.dispersion->accept_as_source(visitor, self->model->parallel_phi.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "parallel_phi", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->parallel_psi.dispersion->accept_as_source(visitor, self->model->parallel_psi.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "parallel_psi", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->parallel_theta.dispersion->accept_as_source(visitor, self->model->parallel_theta.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "parallel_theta", disp_dict);
+
+
+         
+        // Create empty log
+        self->log = PyDict_New();
+        
+        
+
+    }
+    return 0;
+}
+
+static PyMemberDef CCSParallelepipedModel_members[] = {
+    {"params", T_OBJECT, offsetof(CCSParallelepipedModel, params), 0,
+     "Parameters"},
+	{"dispersion", T_OBJECT, offsetof(CCSParallelepipedModel, dispersion), 0,
+	  "Dispersion parameters"},     
+    {"log", T_OBJECT, offsetof(CCSParallelepipedModel, log), 0,
+     "Log"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CCSParallelepipedModel_readDouble(PyObject *p) {
+    if (PyFloat_Check(p)==1) {
+        return (double)(((PyFloatObject *)(p))->ob_fval);
+    } else if (PyInt_Check(p)==1) {
+        return (double)(((PyIntObject *)(p))->ob_ival);
+    } else if (PyLong_Check(p)==1) {
+        return (double)PyLong_AsLong(p);
+    } else {
+        return 0.0;
+    }
+}
+/**
+ * Function to call to evaluate model
+ * @param args: input numpy array q[] 
+ * @return: numpy array object 
+ */
+ 
+static PyObject *evaluateOneDim(CSParallelepipedModel* model, PyArrayObject *q){
+    PyArrayObject *result;
+   
+    // Check validity of array q , q must be of dimension 1, an array of double
+    if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE)
+    {
+        //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num;
+        //PyErr_SetString(PyExc_ValueError , message);
+        return NULL;
+    }
+    result = (PyArrayObject *)PyArray_FromDims(q->nd, (int *)(q->dimensions), 
+										  PyArray_DOUBLE);
+	if (result == NULL) {
+        const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+		return NULL;
+	}
+	 for (int i = 0; i < q->dimensions[0]; i++){
+      double q_value  = *(double *)(q->data + i*q->strides[0]);
+      double *result_value = (double *)(result->data + i*result->strides[0]);
+      *result_value =(*model)(q_value);
+	}
+    return PyArray_Return(result); 
+ }
+
+ /**
+ * Function to call to evaluate model
+ * @param args: input numpy array  [x[],y[]]
+ * @return: numpy array object 
+ */
+ static PyObject * evaluateTwoDimXY( CSParallelepipedModel* model, 
+                              PyArrayObject *x, PyArrayObject *y)
+ {
+    PyArrayObject *result;
+    int i,j, x_len, y_len, dims[1];
+    //check validity of input vectors
+    if (x->nd != 1 || x->descr->type_num != PyArray_DOUBLE
+        || y->nd != 1 || y->descr->type_num != PyArray_DOUBLE
+        || y->dimensions[0] != x->dimensions[0]){
+        const char * message= "evaluateTwoDimXY  expect 2 numpy arrays";
+        PyErr_SetString(PyExc_ValueError , message); 
+        return NULL;
+    }
+   
+	if (PyArray_Check(x) && PyArray_Check(y)) {
+		
+	    x_len = dims[0]= x->dimensions[0];
+        y_len = dims[0]= y->dimensions[0];
+	    
+	    // Make a new double matrix of same dims
+        result=(PyArrayObject *) PyArray_FromDims(1,dims,NPY_DOUBLE);
+        if (result == NULL){
+	    const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+	    return NULL;
+	    }
+       
+        /* Do the calculation. */
+        for ( i=0; i< x_len; i++) {
+            double x_value = *(double *)(x->data + i*x->strides[0]);
+  		    double y_value = *(double *)(y->data + i*y->strides[0]);
+  			double *result_value = (double *)(result->data +
+  			      i*result->strides[0]);
+  			*result_value = (*model)(x_value, y_value);
+        }           
+        return PyArray_Return(result); 
+        
+        }else{
+		    PyErr_SetString(CCSParallelepipedModelError, 
+                   "CCSParallelepipedModel.evaluateTwoDimXY couldn't run.");
+	        return NULL;
+		}      	
+}
+/**
+ *  evalDistribution function evaluate a model function with input vector
+ *  @param args: input q as vector or [qx, qy] where qx, qy are vectors
+ *
+ */ 
+static PyObject * evalDistribution(CCSParallelepipedModel *self, PyObject *args){
+	PyObject *qx, *qy;
+	PyArrayObject * pars;
+	int npars ,mpars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->sld_rimB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_rimB") );
+    self->model->sld_rimC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_rimC") );
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->parallel_psi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_psi") );
+    self->model->midB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "midB") );
+    self->model->parallel_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_theta") );
+    self->model->parallel_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_phi") );
+    self->model->sld_pcore = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_pcore") );
+    self->model->sld_rimA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_rimA") );
+    self->model->shortA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "shortA") );
+    self->model->rimB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rimB") );
+    self->model->rimC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rimC") );
+    self->model->longC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longC") );
+    self->model->sld_solv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_solv") );
+    self->model->rimA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rimA") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "shortA");
+    self->model->shortA.dispersion->accept_as_destination(visitor, self->model->shortA.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "midB");
+    self->model->midB.dispersion->accept_as_destination(visitor, self->model->midB.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "longC");
+    self->model->longC.dispersion->accept_as_destination(visitor, self->model->longC.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_phi");
+    self->model->parallel_phi.dispersion->accept_as_destination(visitor, self->model->parallel_phi.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_psi");
+    self->model->parallel_psi.dispersion->accept_as_destination(visitor, self->model->parallel_psi.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_theta");
+    self->model->parallel_theta.dispersion->accept_as_destination(visitor, self->model->parallel_theta.dispersion, disp_dict);
+
+	
+	// Get input and determine whether we have to supply a 1D or 2D return value.
+	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
+	    PyErr_SetString(CCSParallelepipedModelError, 
+	    	"CCSParallelepipedModel.evalDistribution expects a q value.");
+		return NULL;
+	}
+    // Check params
+	
+    if(PyArray_Check(pars)==1) {
+		
+	    // Length of list should 1 or 2
+	    npars = pars->nd; 
+	    if(npars==1) {
+	        // input is a numpy array
+	        if (PyArray_Check(pars)) {
+		        return evaluateOneDim(self->model, (PyArrayObject*)pars); 
+		    }
+		}else{
+		    PyErr_SetString(CCSParallelepipedModelError, 
+                   "CCSParallelepipedModel.evalDistribution expect numpy array of one dimension.");
+	        return NULL;
+		}
+    }else if( PyList_Check(pars)==1) {
+    	// Length of list should be 2 for I(qx,qy)
+	    mpars = PyList_GET_SIZE(pars); 
+	    if(mpars!=2) {
+	    	PyErr_SetString(CCSParallelepipedModelError, 
+	    		"CCSParallelepipedModel.evalDistribution expects a list of dimension 2.");
+	    	return NULL;
+	    }
+	     qx = PyList_GET_ITEM(pars,0);
+	     qy = PyList_GET_ITEM(pars,1);
+	     if (PyArray_Check(qx) && PyArray_Check(qy)) {
+	         return evaluateTwoDimXY(self->model, (PyArrayObject*)qx,
+		           (PyArrayObject*)qy);
+		 }else{
+		    PyErr_SetString(CCSParallelepipedModelError, 
+                   "CCSParallelepipedModel.evalDistribution expect 2 numpy arrays in list.");
+	        return NULL;
+	     }
+	}
+	PyErr_SetString(CCSParallelepipedModelError, 
+                   "CCSParallelepipedModel.evalDistribution couln't be run.");
+	return NULL;
+	
+}
+
+/**
+ * Function to call to evaluate model
+ * @param args: input q or [q,phi]
+ * @return: function value
+ */
+static PyObject * run(CCSParallelepipedModel *self, PyObject *args) {
+	double q_value, phi_value;
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->sld_rimB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_rimB") );
+    self->model->sld_rimC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_rimC") );
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->parallel_psi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_psi") );
+    self->model->midB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "midB") );
+    self->model->parallel_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_theta") );
+    self->model->parallel_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_phi") );
+    self->model->sld_pcore = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_pcore") );
+    self->model->sld_rimA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_rimA") );
+    self->model->shortA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "shortA") );
+    self->model->rimB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rimB") );
+    self->model->rimC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rimC") );
+    self->model->longC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longC") );
+    self->model->sld_solv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_solv") );
+    self->model->rimA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rimA") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "shortA");
+    self->model->shortA.dispersion->accept_as_destination(visitor, self->model->shortA.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "midB");
+    self->model->midB.dispersion->accept_as_destination(visitor, self->model->midB.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "longC");
+    self->model->longC.dispersion->accept_as_destination(visitor, self->model->longC.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_phi");
+    self->model->parallel_phi.dispersion->accept_as_destination(visitor, self->model->parallel_phi.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_psi");
+    self->model->parallel_psi.dispersion->accept_as_destination(visitor, self->model->parallel_psi.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_theta");
+    self->model->parallel_theta.dispersion->accept_as_destination(visitor, self->model->parallel_theta.dispersion, disp_dict);
+
+	
+	// Get input and determine whether we have to supply a 1D or 2D return value.
+	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
+	    PyErr_SetString(CCSParallelepipedModelError, 
+	    	"CCSParallelepipedModel.run expects a q value.");
+		return NULL;
+	}
+	  
+	// Check params
+	if( PyList_Check(pars)==1) {
+		
+		// Length of list should be 2 for I(q,phi)
+	    npars = PyList_GET_SIZE(pars); 
+	    if(npars!=2) {
+	    	PyErr_SetString(CCSParallelepipedModelError, 
+	    		"CCSParallelepipedModel.run expects a double or a list of dimension 2.");
+	    	return NULL;
+	    }
+	    // We have a vector q, get the q and phi values at which
+	    // to evaluate I(q,phi)
+	    q_value = CCSParallelepipedModel_readDouble(PyList_GET_ITEM(pars,0));
+	    phi_value = CCSParallelepipedModel_readDouble(PyList_GET_ITEM(pars,1));
+	    // Skip zero
+	    if (q_value==0) {
+	    	return Py_BuildValue("d",0.0);
+	    }
+		return Py_BuildValue("d",(*(self->model)).evaluate_rphi(q_value,phi_value));
+
+	} else {
+
+		// We have a scalar q, we will evaluate I(q)
+		q_value = CCSParallelepipedModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",(*(self->model))(q_value));
+	}	
+}
+/**
+ * Function to call to calculate_ER
+ * @return: effective radius value 
+ */
+static PyObject * calculate_ER(CCSParallelepipedModel *self) {
+
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->sld_rimB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_rimB") );
+    self->model->sld_rimC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_rimC") );
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->parallel_psi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_psi") );
+    self->model->midB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "midB") );
+    self->model->parallel_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_theta") );
+    self->model->parallel_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_phi") );
+    self->model->sld_pcore = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_pcore") );
+    self->model->sld_rimA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_rimA") );
+    self->model->shortA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "shortA") );
+    self->model->rimB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rimB") );
+    self->model->rimC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rimC") );
+    self->model->longC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longC") );
+    self->model->sld_solv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_solv") );
+    self->model->rimA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rimA") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "shortA");
+    self->model->shortA.dispersion->accept_as_destination(visitor, self->model->shortA.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "midB");
+    self->model->midB.dispersion->accept_as_destination(visitor, self->model->midB.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "longC");
+    self->model->longC.dispersion->accept_as_destination(visitor, self->model->longC.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_phi");
+    self->model->parallel_phi.dispersion->accept_as_destination(visitor, self->model->parallel_phi.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_psi");
+    self->model->parallel_psi.dispersion->accept_as_destination(visitor, self->model->parallel_psi.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_theta");
+    self->model->parallel_theta.dispersion->accept_as_destination(visitor, self->model->parallel_theta.dispersion, disp_dict);
+
+		
+	return Py_BuildValue("d",(*(self->model)).calculate_ER());
+
+}
+/**
+ * Function to call to evaluate model in cartesian coordinates
+ * @param args: input q or [qx, qy]]
+ * @return: function value
+ */
+static PyObject * runXY(CCSParallelepipedModel *self, PyObject *args) {
+	double qx_value, qy_value;
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->sld_rimB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_rimB") );
+    self->model->sld_rimC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_rimC") );
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->parallel_psi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_psi") );
+    self->model->midB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "midB") );
+    self->model->parallel_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_theta") );
+    self->model->parallel_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_phi") );
+    self->model->sld_pcore = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_pcore") );
+    self->model->sld_rimA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_rimA") );
+    self->model->shortA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "shortA") );
+    self->model->rimB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rimB") );
+    self->model->rimC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rimC") );
+    self->model->longC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longC") );
+    self->model->sld_solv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_solv") );
+    self->model->rimA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rimA") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "shortA");
+    self->model->shortA.dispersion->accept_as_destination(visitor, self->model->shortA.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "midB");
+    self->model->midB.dispersion->accept_as_destination(visitor, self->model->midB.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "longC");
+    self->model->longC.dispersion->accept_as_destination(visitor, self->model->longC.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_phi");
+    self->model->parallel_phi.dispersion->accept_as_destination(visitor, self->model->parallel_phi.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_psi");
+    self->model->parallel_psi.dispersion->accept_as_destination(visitor, self->model->parallel_psi.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_theta");
+    self->model->parallel_theta.dispersion->accept_as_destination(visitor, self->model->parallel_theta.dispersion, disp_dict);
+
+	
+	// Get input and determine whether we have to supply a 1D or 2D return value.
+	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
+	    PyErr_SetString(CCSParallelepipedModelError, 
+	    	"CCSParallelepipedModel.run expects a q value.");
+		return NULL;
+	}
+	  
+	// Check params
+	if( PyList_Check(pars)==1) {
+		
+		// Length of list should be 2 for I(qx, qy))
+	    npars = PyList_GET_SIZE(pars); 
+	    if(npars!=2) {
+	    	PyErr_SetString(CCSParallelepipedModelError, 
+	    		"CCSParallelepipedModel.run expects a double or a list of dimension 2.");
+	    	return NULL;
+	    }
+	    // We have a vector q, get the qx and qy values at which
+	    // to evaluate I(qx,qy)
+	    qx_value = CCSParallelepipedModel_readDouble(PyList_GET_ITEM(pars,0));
+	    qy_value = CCSParallelepipedModel_readDouble(PyList_GET_ITEM(pars,1));
+	    return Py_BuildValue("d",(*(self->model))(qx_value,qy_value));
+
+	} else {
+
+		// We have a scalar q, we will evaluate I(q)
+		qx_value = CCSParallelepipedModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",(*(self->model))(qx_value));
+	}	
+}
+
+static PyObject * reset(CCSParallelepipedModel *self, PyObject *args) {
+    
+
+    return Py_BuildValue("d",0.0);
+}
+
+static PyObject * set_dispersion(CCSParallelepipedModel *self, PyObject *args) {
+	PyObject * disp;
+	const char * par_name;
+
+	if ( !PyArg_ParseTuple(args,"sO", &par_name, &disp) ) {
+	    PyErr_SetString(CCSParallelepipedModelError,
+	    	"CCSParallelepipedModel.set_dispersion expects a DispersionModel object.");
+		return NULL;
+	}
+	void *temp = PyCObject_AsVoidPtr(disp);
+	DispersionModel * dispersion = static_cast<DispersionModel *>(temp);
+
+
+	// Ugliness necessary to go from python to C
+	    // TODO: refactor this
+    if (!strcmp(par_name, "shortA")) {
+        self->model->shortA.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "midB")) {
+        self->model->midB.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "longC")) {
+        self->model->longC.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "parallel_phi")) {
+        self->model->parallel_phi.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "parallel_psi")) {
+        self->model->parallel_psi.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "parallel_theta")) {
+        self->model->parallel_theta.dispersion = dispersion;
+    } else {
+	    PyErr_SetString(CCSParallelepipedModelError,
+	    	"CCSParallelepipedModel.set_dispersion expects a valid parameter name.");
+		return NULL;
+	}
+
+	DispersionVisitor* visitor = new DispersionVisitor();
+	PyObject * disp_dict = PyDict_New();
+	dispersion->accept_as_source(visitor, dispersion, disp_dict);
+	PyDict_SetItemString(self->dispersion, par_name, disp_dict);
+    return Py_BuildValue("i",1);
+}
+
+
+static PyMethodDef CCSParallelepipedModel_methods[] = {
+    {"run",      (PyCFunction)run     , METH_VARARGS,
+      "Evaluate the model at a given Q or Q, phi"},
+    {"runXY",      (PyCFunction)runXY     , METH_VARARGS,
+      "Evaluate the model at a given Q or Qx, Qy"},
+    {"calculate_ER",      (PyCFunction)calculate_ER     , METH_VARARGS,
+      "Evaluate the model at a given Q or Q, phi"},
+      
+    {"evalDistribution",  (PyCFunction)evalDistribution , METH_VARARGS,
+      "Evaluate the model at a given Q or Qx, Qy vector "},
+    {"reset",    (PyCFunction)reset   , METH_VARARGS,
+      "Reset pair correlation"},
+    {"set_dispersion",      (PyCFunction)set_dispersion     , METH_VARARGS,
+      "Set the dispersion model for a given parameter"},
+   {NULL}
+};
+
+static PyTypeObject CCSParallelepipedModelType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CCSParallelepipedModel",             /*tp_name*/
+    sizeof(CCSParallelepipedModel),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CCSParallelepipedModel_dealloc, /*tp_dealloc*/
+    0,                         /*tp_print*/
+    0,                         /*tp_getattr*/
+    0,                         /*tp_setattr*/
+    0,                         /*tp_compare*/
+    0,                         /*tp_repr*/
+    0,                         /*tp_as_number*/
+    0,                         /*tp_as_sequence*/
+    0,                         /*tp_as_mapping*/
+    0,                         /*tp_hash */
+    0,                         /*tp_call*/
+    0,                         /*tp_str*/
+    0,                         /*tp_getattro*/
+    0,                         /*tp_setattro*/
+    0,                         /*tp_as_buffer*/
+    Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
+    "CCSParallelepipedModel objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CCSParallelepipedModel_methods,             /* tp_methods */
+    CCSParallelepipedModel_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CCSParallelepipedModel_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CCSParallelepipedModel_new,                 /* tp_new */
+};
+
+
+//static PyMethodDef module_methods[] = {
+//    {NULL} 
+//};
+
+/**
+ * Function used to add the model class to a module
+ * @param module: module to add the class to
+ */ 
+void addCCSParallelepipedModel(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CCSParallelepipedModelType) < 0)
+        return;
+
+    Py_INCREF(&CCSParallelepipedModelType);
+    PyModule_AddObject(module, "CCSParallelepipedModel", (PyObject *)&CCSParallelepipedModelType);
+    
+    d = PyModule_GetDict(module);
+    CCSParallelepipedModelError = PyErr_NewException("CCSParallelepipedModel.error", NULL, NULL);
+    PyDict_SetItemString(d, "CCSParallelepipedModelError", CCSParallelepipedModelError);
+}
+
Index: sansmodels/src/sans/models/c_models/c_models.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/c_models.cpp	(revision fbd09fa09c7abeb0de3821f31585186902b5f12e)
+++ sansmodels/src/sans/models/c_models/c_models.cpp	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
@@ -16,4 +16,5 @@
 void addCTriaxialEllipsoidModel(PyObject *module);
 void addCParallelepipedModel(PyObject *module);
+void addCCSParallelepipedModel(PyObject *module);
 void addCSphereModel(PyObject *module);
 void addCOnionModel(PyObject *module);
@@ -227,4 +228,5 @@
 	addCCappedCylinderModel(m);
 	addCParallelepipedModel(m);
+	addCCSParallelepipedModel(m);
 	addCCoreShellCylinderModel(m);
 	addCCoreShellModel(m);
Index: sansmodels/src/sans/models/c_models/csparallelepiped.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/csparallelepiped.cpp	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
+++ sansmodels/src/sans/models/c_models/csparallelepiped.cpp	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
@@ -0,0 +1,343 @@
+/**
+	This software was developed by the University of Tennessee as part of the
+	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+	project funded by the US National Science Foundation.
+
+	If you use DANSE applications to do scientific research that leads to
+	publication, we ask that you acknowledge the use of the software with the
+	following sentence:
+
+	"This work benefited from DANSE software developed under NSF award DMR-0520547."
+
+	copyright 2010, University of Tennessee
+ */
+
+/**
+ * Scattering model classes
+ * The classes use the IGOR library found in
+ *   sansmodels/src/libigor
+ */
+
+#include <math.h>
+#include "models.hh"
+#include "parameters.hh"
+#include <stdio.h>
+using namespace std;
+
+extern "C" {
+	#include "libCylinder.h"
+	#include "libStructureFactor.h"
+	#include "csparallelepiped.h"
+}
+
+CSParallelepipedModel :: CSParallelepipedModel() {
+	scale      = Parameter(1.0);
+	shortA     = Parameter(35.0, true);
+	shortA.set_min(1.0);
+	midB     = Parameter(75.0, true);
+	midB.set_min(1.0);
+	longC    = Parameter(400.0, true);
+	longC.set_min(1.0);
+	rimA     = Parameter(10.0, true);
+	rimB     = Parameter(10.0, true);
+	rimC     = Parameter(10.0, true);
+	sld_rimA     = Parameter(2.0e-6, true);
+	sld_rimB     = Parameter(4.0e-6, true);
+	sld_rimC    = Parameter(2.0e-6, true);
+	sld_pcore   = Parameter(1.0e-6);
+	sld_solv   = Parameter(6.0e-6);
+	background = Parameter(0.06);
+	parallel_theta  = Parameter(0.0, true);
+	parallel_phi    = Parameter(0.0, true);
+	parallel_psi    = Parameter(0.0, true);
+}
+
+/**
+ * Function to evaluate 1D scattering function
+ * The NIST IGOR library is used for the actual calculation.
+ * @param q: q-value
+ * @return: function value
+ */
+double CSParallelepipedModel :: operator()(double q) {
+	double dp[13];
+
+	// Fill parameter array for IGOR library
+	// Add the background after averaging
+	dp[0] = scale();
+	dp[1] = shortA();
+	dp[2] = midB();
+	dp[3] = longC();
+	dp[4] = rimA();
+	dp[5] = rimB();
+	dp[6] = rimC();
+	dp[7] = sld_rimA();
+	dp[8] = sld_rimB();
+	dp[9] = sld_rimC();
+	dp[10] = sld_pcore();
+	dp[11] = sld_solv();
+	dp[12] = 0.0;
+
+	// Get the dispersion points for the short_edgeA
+	vector<WeightPoint> weights_shortA;
+	shortA.get_weights(weights_shortA);
+
+	// Get the dispersion points for the longer_edgeB
+	vector<WeightPoint> weights_midB;
+	midB.get_weights(weights_midB);
+
+	// Get the dispersion points for the longuest_edgeC
+	vector<WeightPoint> weights_longC;
+	longC.get_weights(weights_longC);
+
+
+
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+	double vol = 0.0;
+
+	// Loop over short_edgeA weight points
+	for(int i=0; i< (int)weights_shortA.size(); i++) {
+		dp[1] = weights_shortA[i].value;
+
+		// Loop over longer_edgeB weight points
+		for(int j=0; j< (int)weights_midB.size(); j++) {
+			dp[2] = weights_midB[j].value;
+
+			// Loop over longuest_edgeC weight points
+			for(int k=0; k< (int)weights_longC.size(); k++) {
+				dp[3] = weights_longC[k].value;
+				//Un-normalize  by volume
+				sum += weights_shortA[i].weight * weights_midB[j].weight
+					* weights_longC[k].weight * CSParallelepiped(dp, q)
+					* weights_shortA[i].value*weights_midB[j].value
+					* weights_longC[k].value;
+				//Find average volume
+				vol += weights_shortA[i].weight * weights_midB[j].weight
+					* weights_longC[k].weight
+					* weights_shortA[i].value * weights_midB[j].value
+					* weights_longC[k].value;
+
+				norm += weights_shortA[i].weight
+					 * weights_midB[j].weight * weights_longC[k].weight;
+			}
+		}
+	}
+	if (vol != 0.0 && norm != 0.0) {
+		//Re-normalize by avg volume
+		sum = sum/(vol/norm);}
+
+	return sum/norm + background();
+}
+/**
+ * Function to evaluate 2D scattering function
+ * @param q_x: value of Q along x
+ * @param q_y: value of Q along y
+ * @return: function value
+ */
+double CSParallelepipedModel :: operator()(double qx, double qy) {
+	CSParallelepipedParameters dp;
+	// Fill parameter array
+	dp.scale      = scale();
+	dp.shortA   = shortA();
+	dp.midB   = midB();
+	dp.longC  = longC();
+	dp.rimA   = rimA();
+	dp.rimB   = rimB();
+	dp.rimC  = rimC();
+	dp.sld_rimA   = sld_rimA();
+	dp.sld_rimB   = sld_rimB();
+	dp.sld_rimC  = sld_rimC();
+	dp.sld_pcore   = sld_pcore();
+	dp.sld_solv   = sld_solv();
+	dp.background = 0.0;
+	//dp.background = background();
+	dp.parallel_theta  = parallel_theta();
+	dp.parallel_phi    = parallel_phi();
+	dp.parallel_psi    = parallel_psi();
+
+
+
+	// Get the dispersion points for the short_edgeA
+	vector<WeightPoint> weights_shortA;
+	shortA.get_weights(weights_shortA);
+
+	// Get the dispersion points for the longer_edgeB
+	vector<WeightPoint> weights_midB;
+	midB.get_weights(weights_midB);
+
+	// Get the dispersion points for the longuest_edgeC
+	vector<WeightPoint> weights_longC;
+	longC.get_weights(weights_longC);
+
+	// Get angular averaging for theta
+	vector<WeightPoint> weights_parallel_theta;
+	parallel_theta.get_weights(weights_parallel_theta);
+
+	// Get angular averaging for phi
+	vector<WeightPoint> weights_parallel_phi;
+	parallel_phi.get_weights(weights_parallel_phi);
+
+	// Get angular averaging for psi
+	vector<WeightPoint> weights_parallel_psi;
+	parallel_psi.get_weights(weights_parallel_psi);
+
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+	double norm_vol = 0.0;
+	double vol = 0.0;
+
+	// Loop over radius weight points
+	for(int i=0; i< (int)weights_shortA.size(); i++) {
+		dp.shortA = weights_shortA[i].value;
+
+		// Loop over longer_edgeB weight points
+		for(int j=0; j< (int)weights_midB.size(); j++) {
+			dp.midB = weights_midB[j].value;
+
+			// Average over longuest_edgeC distribution
+			for(int k=0; k< (int)weights_longC.size(); k++) {
+				dp.longC = weights_longC[k].value;
+
+				// Average over theta distribution
+				for(int l=0; l< (int)weights_parallel_theta.size(); l++) {
+				dp.parallel_theta = weights_parallel_theta[l].value;
+
+					// Average over phi distribution
+					for(int m=0; m< (int)weights_parallel_phi.size(); m++) {
+						dp.parallel_phi = weights_parallel_phi[m].value;
+
+						// Average over phi distribution
+						for(int n=0; n< (int)weights_parallel_psi.size(); n++) {
+							dp.parallel_psi = weights_parallel_psi[n].value;
+							//Un-normalize by volume
+							double _ptvalue = weights_shortA[i].weight
+								* weights_midB[j].weight
+								* weights_longC[k].weight
+								* weights_parallel_theta[l].weight
+								* weights_parallel_phi[m].weight
+								* weights_parallel_psi[n].weight
+								* csparallelepiped_analytical_2DXY(&dp, qx, qy)
+								* weights_shortA[i].value*weights_midB[j].value
+								* weights_longC[k].value;
+
+							if (weights_parallel_theta.size()>1) {
+								_ptvalue *= sin(weights_parallel_theta[l].value);
+							}
+							sum += _ptvalue;
+							//Find average volume
+							vol += weights_shortA[i].weight
+								* weights_midB[j].weight
+								* weights_longC[k].weight
+								* weights_shortA[i].value*weights_midB[j].value
+								* weights_longC[k].value;
+							//Find norm for volume
+							norm_vol += weights_shortA[i].weight
+								* weights_midB[j].weight
+								* weights_longC[k].weight;
+
+							norm += weights_shortA[i].weight
+								* weights_midB[j].weight
+								* weights_longC[k].weight
+								* weights_parallel_theta[l].weight
+								* weights_parallel_phi[m].weight
+								* weights_parallel_psi[n].weight;
+						}
+					}
+
+				}
+			}
+		}
+	}
+	// Averaging in theta needs an extra normalization
+	// factor to account for the sin(theta) term in the
+	// integration (see documentation).
+	if (weights_parallel_theta.size()>1) norm = norm / asin(1.0);
+
+	if (vol != 0.0 && norm_vol != 0.0) {
+		//Re-normalize by avg volume
+		sum = sum/(vol/norm_vol);}
+
+	return sum/norm + background();
+}
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double CSParallelepipedModel :: evaluate_rphi(double q, double phi) {
+	double qx = q*cos(phi);
+	double qy = q*sin(phi);
+	return (*this).operator()(qx, qy);
+}
+/**
+ * Function to calculate effective radius
+ * @return: effective radius value
+ */
+double CSParallelepipedModel :: calculate_ER() {
+	CSParallelepipedParameters dp;
+	dp.shortA   = shortA();
+	dp.midB   = midB();
+	dp.longC  = longC();
+	dp.rimA   = rimA();
+	dp.rimB   = rimB();
+	dp.rimC  = rimC();
+
+	double rad_out = 0.0;
+	double pi = 4.0*atan(1.0);
+	double suf_rad = sqrt((dp.shortA*dp.midB+2.0*dp.rimA*dp.midB+2.0*dp.rimA*dp.shortA)/pi);
+	double height =(dp.longC + 2.0*dp.rimC);
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+
+	// Get the dispersion points for the short_edgeA
+	vector<WeightPoint> weights_shortA;
+	shortA.get_weights(weights_shortA);
+
+	// Get the dispersion points for the longer_edgeB
+	vector<WeightPoint> weights_midB;
+	midB.get_weights(weights_midB);
+
+	// Get angular averaging for the longuest_edgeC
+	vector<WeightPoint> weights_longC;
+	longC.get_weights(weights_longC);
+
+	// Loop over radius weight points
+	for(int i=0; i< (int)weights_shortA.size(); i++) {
+		dp.shortA = weights_shortA[i].value;
+
+		// Loop over longer_edgeB weight points
+		for(int j=0; j< (int)weights_midB.size(); j++) {
+			dp.midB = weights_midB[j].value;
+
+			// Average over longuest_edgeC distribution
+			for(int k=0; k< (int)weights_longC.size(); k++) {
+				dp.longC = weights_longC[k].value;
+				//Calculate surface averaged radius
+				//This is rough approximation.
+				suf_rad = sqrt((dp.shortA*dp.midB+2.0*dp.rimA*dp.midB+2.0*dp.rimA*dp.shortA)/pi);
+				height =(dp.longC + 2.0*dp.rimC);
+				//Note: output of "DiamCyl(dp.length,dp.radius)" is a DIAMETER.
+				sum +=weights_shortA[i].weight* weights_midB[j].weight
+					* weights_longC[k].weight*DiamCyl(height, suf_rad)/2.0;
+				norm += weights_shortA[i].weight* weights_midB[j].weight*weights_longC[k].weight;
+			}
+		}
+	}
+
+	if (norm != 0){
+		//return the averaged value
+		rad_out =  sum/norm;}
+	else{
+		//return normal value
+		//Note: output of "DiamCyl(length,radius)" is DIAMETER.
+		rad_out = DiamCyl(dp.longC, suf_rad)/2.0;}
+	return rad_out;
+
+}
Index: sansmodels/src/sans/models/c_models/models.hh
===================================================================
--- sansmodels/src/sans/models/c_models/models.hh	(revision 339ce675656e3ab89eff7b84d8292a2e0ffe6109)
+++ sansmodels/src/sans/models/c_models/models.hh	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
@@ -98,5 +98,4 @@
 class ParallelepipedModel{
 public:
-	// TODO: add 2D
 	// Model parameters
 	Parameter scale;
@@ -121,4 +120,34 @@
 };
 
+
+class CSParallelepipedModel{
+public:
+	// Model parameters
+	Parameter scale;
+	Parameter shortA;
+	Parameter midB;
+	Parameter longC;
+	Parameter rimA;
+	Parameter rimB;
+	Parameter rimC;
+	Parameter sld_rimA;
+	Parameter sld_rimB;
+	Parameter sld_rimC;
+	Parameter sld_pcore;
+	Parameter sld_solv;
+	Parameter background;
+	Parameter parallel_theta;
+	Parameter parallel_phi;
+	Parameter parallel_psi;
+
+	// Constructor
+	CSParallelepipedModel();
+
+	// Operators to get I(Q)
+	double operator()(double q);
+	double operator()(double qx, double qy);
+	double calculate_ER();
+	double evaluate_rphi(double q, double phi);
+};
 
 class OnionModel{
Index: sansmodels/src/sans/models/test/utest_extra_models.py
===================================================================
--- sansmodels/src/sans/models/test/utest_extra_models.py	(revision 495c377a84ff414342be48069ee14551a6e609aa)
+++ sansmodels/src/sans/models/test/utest_extra_models.py	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
@@ -346,5 +346,20 @@
         self.assertAlmostEqual(self.model4.runXY(0.4425284), 0.00306386, 8)  
                                                                               
-                                                                                                                                                                                 
+
+class TestCSPP(unittest.TestCase):
+    """
+    Unit tests for CSParallelepiped Model function
+    """
+    def setUp(self):
+        from sans.models.CSParallelepipedModel import CSParallelepipedModel
+        self.model= CSParallelepipedModel()
+        
+    def test1D(self):        
+        # the values are from Igor pro calculation  
+        # the different digits are due to precision of q values  
+        self.assertAlmostEqual(self.model.run(0.001), 1383.96, 2)
+        self.assertAlmostEqual(self.model.run(0.1501412), 2.51932, 4)
+        self.assertAlmostEqual(self.model.runXY(0.4425284), 0.0736735, 6)
+                                                                                                                                                          
 if __name__ == '__main__':
     unittest.main()
