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
