Index: sansmodels/src/sans/models/SphereExpShellModel.py
===================================================================
--- sansmodels/src/sans/models/SphereExpShellModel.py	(revision 339ce675656e3ab89eff7b84d8292a2e0ffe6109)
+++ sansmodels/src/sans/models/SphereExpShellModel.py	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -164,5 +164,6 @@
             
             A = self.params['A_shell%s'% str(n)]
-            if A <1.0e-16:
+            from math import fabs
+            if fabs(A) <1.0e-16:
                 # Right side of each shells
                 r0 += self.params['thick_shell%s'% str(n)]
Index: sansmodels/src/sans/models/c_extensions/bcc.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/bcc.c	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
+++ sansmodels/src/sans/models/c_extensions/bcc.c	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -61,5 +61,5 @@
 	double b3_x, b3_y, b3_z, b2_x, b2_y, b1_x, b1_y;
 	double q_z;
-	double alpha, vol, cos_val_a3, cos_val_a2, cos_val_a1;
+	double alpha, vol, cos_val_b3, cos_val_b2, cos_val_b1;
 	double a1_dot_q, a2_dot_q,a3_dot_q;
 	double answer;
@@ -80,5 +80,6 @@
 	//the occupied volume of the lattice
 	latticeScale = 2.0*(4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa/sqrt(3.0/4.0),3.0);
-
+	// 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))
@@ -87,41 +88,29 @@
     b3_y = sin(pars->theta) * sin(pars->phi);
     b3_z = cos(pars->theta);
+    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 = sin(pars->psi);
     b1_y = cos(pars->psi);
+	cos_val_b1 = (b1_x*q_x + b1_y*q_y);
     // b2 axis orientation
-    b2_x = sqrt(1.0-sin(pars->theta)*cos(pars->phi))*cos(pars->psi);
-    b2_y = sqrt(1.0-sin(pars->theta)*cos(pars->phi))*sin(pars->psi);
-
-    // a3 axis orientation
-    a3_x = 0.5*(b2_x + b1_x - b3_x);
-    a3_y = 0.5*(b2_y + b1_y - b3_y);
-    a3_z = 0.0;
-    // a1 axis orientation
-    a1_x = 0.5*(b3_x + b2_x - b1_x);
-    a1_y = 0.5*(b3_y + b2_y - b1_y);
-    // a2 axis orientation
-    a2_x = 0.5*(b3_x + b1_x - b2_x);
-    a2_y = 0.5*(b3_y + b1_y - b2_y);
-
-    // q vector
-    q_z = 0.0; // for SANS; assuming qz is negligible
+	cos_val_b2 = sin(acos(cos_val_b1));
+	// alpha corrections
+	cos_val_b2 *= sin(alpha);
+	cos_val_b1 *= sin(alpha);
 
     // 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);
-    a3_dot_q = aa*q*cos_val_a3;
+    a3_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b1-cos_val_b3);
 
     // a1 axis
-    cos_val_a1 = a1_x*q_x + a1_y*q_y;
-    a1_dot_q = aa*q*cos_val_a1*sin(alpha);
+    a1_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b2-cos_val_b1);
 
     // a2 axis
-    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
+    a2_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b1-cos_val_b2);
 
     // The following test should always pass
-    if (fabs(cos_val_a3)>1.0) {
-    	printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+    if (fabs(cos_val_b3)>1.0) {
+    	printf("bcc_ana_2D: Unexpected error: cos(alpha)>1\n");
      	return 0;
     }
@@ -131,6 +120,6 @@
     // Call Zq=Z1*Z2*Z3
     Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
-    Zq = Zq * (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
-    Zq = Zq * (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_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
Index: sansmodels/src/sans/models/c_extensions/fcc.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/fcc.c	(revision 18f2ca1d86b9efb3c8d2b861df22475a4c46cb07)
+++ sansmodels/src/sans/models/c_extensions/fcc.c	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -61,5 +61,5 @@
 	double b3_x, b3_y, b3_z, b2_x, b2_y, b1_x, b1_y;
 	double q_z;
-	double alpha, vol, cos_val_a3, cos_val_a2, cos_val_a1;
+	double alpha, vol, cos_val_b3, cos_val_b2, cos_val_b1, cos_val_a3, cos_val_a2, cos_val_a1;
 	double a1_dot_q, a2_dot_q,a3_dot_q;
 	double answer;
@@ -80,48 +80,36 @@
 
 	latticeScale = 4.0*(4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa*sqrt(2.0),3.0);
-
+	// 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))
+	///  instead of against q coordinate in PRB 36(46), 3(6), 1754(3854)
     // b3 axis orientation
     b3_x = sin(pars->theta) * cos(pars->phi);//negative sign here???
     b3_y = sin(pars->theta) * sin(pars->phi);
     b3_z = cos(pars->theta);
+    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 = sin(pars->psi);
     b1_y = cos(pars->psi);
+	cos_val_b1 = (b1_x*q_x + b1_y*q_y);
     // b2 axis orientation
-    b2_x = sqrt(1.0-sin(pars->theta)*cos(pars->phi))*cos(pars->psi);
-    b2_y = sqrt(1.0-sin(pars->theta)*cos(pars->phi))*sin(pars->psi);
-
-    // a3 axis orientation
-    a3_x = 0.5*(b2_x + b1_x);
-    a3_y = 0.5*(b2_y + b1_y);
-    a3_z = 0.0;
-    // a1 axis orientation
-    a1_x = 0.5*(b3_x + b2_x);
-    a1_y = 0.5*(b3_y + b2_y);
-    // a2 axis orientation
-    a2_x = 0.5*(b3_x + b1_x);
-    a2_y = 0.5*(b3_y + b1_y);
-
-    // q vector
-    q_z = 0.0; // for SANS; assuming qz is negligible
+	cos_val_b2 = sin(acos(cos_val_b1));
+	// alpha correction
+	cos_val_b2 *= sin(alpha);
+	cos_val_b1 *= sin(alpha);
 
     // 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);
-    a3_dot_q = aa*q*cos_val_a3;
+    a3_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b1);
 
     // a1 axis
-    cos_val_a1 = a1_x*q_x + a1_y*q_y;
-    a1_dot_q = aa*q*cos_val_a1*sin(alpha);
+    a1_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b3);
 
     // a2 axis
-    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);
+    a2_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b1);
 
     // The following test should always pass
-    if (fabs(cos_val_a3)>1.0) {
-    	printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+    if (fabs(cos_val_b3)>1.0) {
+    	printf("fcc_ana_2D: Unexpected error: cos(alpha)>1\n");
      	return 0;
     }
Index: sansmodels/src/sans/models/c_models/barbell.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/barbell.cpp	(revision 339ce675656e3ab89eff7b84d8292a2e0ffe6109)
+++ sansmodels/src/sans/models/c_models/barbell.cpp	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -198,5 +198,5 @@
 						}
 						if (weights_theta.size()>1) {
-							_ptvalue *= sin(weights_theta[l].value);
+							_ptvalue *= fabs(sin(weights_theta[l].value));
 						}
 						sum += _ptvalue;
Index: sansmodels/src/sans/models/c_models/bcc.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/bcc.cpp	(revision d5b6a9d7cd9e26b5f9e28e55e8cf848a40de6ae9)
+++ sansmodels/src/sans/models/c_models/bcc.cpp	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -166,5 +166,5 @@
 					}
 					if (weights_theta.size()>1) {
-						_ptvalue *= sin(weights_theta[j].value);
+						_ptvalue *= fabs(sin(weights_theta[j].value));
 					}
 					sum += _ptvalue;
Index: sansmodels/src/sans/models/c_models/capcyl.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/capcyl.cpp	(revision 339ce675656e3ab89eff7b84d8292a2e0ffe6109)
+++ sansmodels/src/sans/models/c_models/capcyl.cpp	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -198,5 +198,5 @@
 						}
 						if (weights_theta.size()>1) {
-							_ptvalue *= sin(weights_theta[l].value);
+							_ptvalue *= fabs(sin(weights_theta[l].value));
 						}
 						sum += _ptvalue;
Index: sansmodels/src/sans/models/c_models/coreshellcylinder.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/coreshellcylinder.cpp	(revision c451be953dd7ec3d17683e043145313f5c761ac0)
+++ sansmodels/src/sans/models/c_models/coreshellcylinder.cpp	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -199,5 +199,5 @@
 
 					if (weights_theta.size()>1) {
-						_ptvalue *= sin(weights_theta[k].value);
+						_ptvalue *= fabs(sin(weights_theta[k].value));
 					}
 					sum += _ptvalue;
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 85bb87048005735e428588676cad765827cc1f95)
@@ -224,5 +224,5 @@
 
 							if (weights_parallel_theta.size()>1) {
-								_ptvalue *= sin(weights_parallel_theta[l].value);
+								_ptvalue *= fabs(sin(weights_parallel_theta[l].value));
 							}
 							sum += _ptvalue;
Index: sansmodels/src/sans/models/c_models/cylinder.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/cylinder.cpp	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/cylinder.cpp	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -167,5 +167,5 @@
 						*pow(weights_rad[i].value,2)*weights_len[j].value;
 					if (weights_theta.size()>1) {
-						_ptvalue *= sin(weights_theta[k].value);
+						_ptvalue *= fabs(sin(weights_theta[k].value));
 					}
 					sum += _ptvalue;
Index: sansmodels/src/sans/models/c_models/ellipsoid.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/ellipsoid.cpp	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/ellipsoid.cpp	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -170,5 +170,5 @@
 						* pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
 					if (weights_theta.size()>1) {
-						_ptvalue *= sin(weights_theta[k].value);
+						_ptvalue *= fabs(sin(weights_theta[k].value));
 					}
 					sum += _ptvalue;
Index: sansmodels/src/sans/models/c_models/ellipticalcylinder.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/ellipticalcylinder.cpp	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/ellipticalcylinder.cpp	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -207,5 +207,5 @@
 						* weights_rat[m].value;
 					if (weights_theta.size()>1) {
-						_ptvalue *= sin(weights_theta[k].value);
+						_ptvalue *= fabs(sin(weights_theta[k].value));
 					}
 					sum += _ptvalue;
Index: sansmodels/src/sans/models/c_models/fcc.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/fcc.cpp	(revision d5b6a9d7cd9e26b5f9e28e55e8cf848a40de6ae9)
+++ sansmodels/src/sans/models/c_models/fcc.cpp	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -165,5 +165,5 @@
 					}
 					if (weights_theta.size()>1) {
-						_ptvalue *= sin(weights_theta[j].value);
+						_ptvalue *= fabs(sin(weights_theta[j].value));
 					}
 					sum += _ptvalue;
Index: sansmodels/src/sans/models/c_models/hollowcylinder.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/hollowcylinder.cpp	(revision 2339249f38c9787cd4816eda4965172695a27647)
+++ sansmodels/src/sans/models/c_models/hollowcylinder.cpp	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -195,5 +195,5 @@
 						* weights_length[j].value;
 					if (weights_theta.size()>1) {
-						_ptvalue *= sin(weights_theta[k].value);
+						_ptvalue *= fabs(sin(weights_theta[k].value));
 					}
 					sum += _ptvalue;
Index: sansmodels/src/sans/models/c_models/parallelepiped.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/parallelepiped.cpp	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/parallelepiped.cpp	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -208,5 +208,5 @@
 
 							if (weights_parallel_theta.size()>1) {
-								_ptvalue *= sin(weights_parallel_theta[l].value);
+								_ptvalue *= fabs(sin(weights_parallel_theta[l].value));
 							}
 							sum += _ptvalue;
Index: sansmodels/src/sans/models/c_models/sc.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/sc.cpp	(revision 94a3f8f49326fa54618105cf9c371ebe39bb3dda)
+++ sansmodels/src/sans/models/c_models/sc.cpp	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -165,5 +165,5 @@
 					}
 					if (weights_theta.size()>1) {
-						_ptvalue *= sin(weights_theta[j].value);
+						_ptvalue *= fabs(sin(weights_theta[j].value));
 					}
 					sum += _ptvalue;
Index: sansmodels/src/sans/models/c_models/spheroid.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/spheroid.cpp	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/spheroid.cpp	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -241,5 +241,5 @@
 								* pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
 							if (weights_theta.size()>1) {
-								_ptvalue *= sin(weights_theta[m].value);
+								_ptvalue *= fabs(sin(weights_theta[m].value));
 							}
 							sum += _ptvalue;
Index: sansmodels/src/sans/models/c_models/stackeddisks.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/stackeddisks.cpp	(revision c451be953dd7ec3d17683e043145313f5c761ac0)
+++ sansmodels/src/sans/models/c_models/stackeddisks.cpp	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -198,5 +198,5 @@
 								*pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value);
 							if (weights_theta.size()>1) {
-								_ptvalue *= sin(weights_theta[l].value);
+								_ptvalue *= fabs(sin(weights_theta[l].value));
 							}
 							sum += _ptvalue;
Index: sansmodels/src/sans/models/c_models/triaxialellipsoid.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/triaxialellipsoid.cpp	(revision 13eb1c4b175eb150bd941382067266672545173d)
+++ sansmodels/src/sans/models/c_models/triaxialellipsoid.cpp	(revision 85bb87048005735e428588676cad765827cc1f95)
@@ -199,5 +199,5 @@
 								* weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
 							if (weights_theta.size()>1) {
-								_ptvalue *= sin(weights_theta[k].value);
+								_ptvalue *= fabs(sin(weights_theta[k].value));
 							}
 							sum += _ptvalue;
