Index: sasmodels/models/bcc_paracrystal.c
===================================================================
--- sasmodels/models/bcc_paracrystal.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
+++ sasmodels/models/bcc_paracrystal.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -0,0 +1,170 @@
+double form_volume(double radius);
+double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld);
+double Iqxy(double qx, double qy, double dnn,
+    double d_factor, double radius,double sld, double solvent_sld,
+    double theta, double phi, double psi);
+
+double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
+double _BCCeval(double Theta, double Phi, double temp1, double temp3);
+double _sphereform(double q, double radius, double sld, double solvent_sld);
+
+
+double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) {
+
+	const double Da = d_factor*dnn;
+	const double temp1 = q*q*Da*Da;
+	const double temp3 = q*dnn;
+
+	double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI);
+	return(retVal);
+}
+
+double _BCCeval(double Theta, double Phi, double temp1, double temp3) {
+
+	double result;
+	double sin_theta,cos_theta,sin_phi,cos_phi;
+	SINCOS(Theta, sin_theta, cos_theta);
+	SINCOS(Phi, sin_phi, cos_phi);
+
+	const double temp6 =  sin_theta;
+	const double temp7 =  sin_theta*cos_phi + sin_theta*sin_phi + cos_theta;
+	const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta;
+	const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta;
+	const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9));
+	result = pow(1.0-(temp10*temp10),3)*temp6
+	    / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10)
+	      * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10)
+	      * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10));
+
+	return (result);
+}
+
+double form_volume(double radius){
+    return sphere_volume(radius);
+}
+
+
+double Iq(double q, double dnn,
+  double d_factor, double radius,
+  double sld, double solvent_sld){
+
+	//Volume fraction calculated from lattice symmetry and sphere radius
+	const double s1 = dnn/sqrt(0.75);
+	const double latticescale = 2.0*sphere_volume(radius/s1);
+
+    const double va = 0.0;
+    const double vb = 2.0*M_PI;
+    const double vaj = 0.0;
+    const double vbj = M_PI;
+
+    double summ = 0.0;
+    double answer = 0.0;
+	for(int i=0; i<150; i++) {
+		//setup inner integral over the ellipsoidal cross-section
+		double summj=0.0;
+		const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;		//the outer dummy is phi
+		for(int j=0;j<150;j++) {
+			//20 gauss points for the inner integral
+			double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;		//the inner dummy is theta
+			double yyy = Gauss150Wt[j] * _BCC_Integrand(q,dnn,d_factor,ztheta,zphi);
+			summj += yyy;
+		}
+		//now calculate the value of the inner integral
+		double answer = (vbj-vaj)/2.0*summj;
+
+		//now calculate outer integral
+		summ = summ+(Gauss150Wt[i] * answer);
+	}		//final scaling is done at the end of the function, after the NT_FP64 case
+
+	answer = (vb-va)/2.0*summ;
+	answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale;
+
+    return answer;
+
+
+}
+
+
+double Iqxy(double qx, double qy, double dnn,
+    double d_factor, double radius,double sld, double solvent_sld,
+    double theta, double phi, double psi){
+
+  double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
+  //double q_z;
+  double cos_val_b3, cos_val_b2, cos_val_b1;
+  double a1_dot_q, a2_dot_q,a3_dot_q;
+  double answer;
+  double Zq, Fkq, Fkq_2;
+
+  //convert to q and make scaled values
+  double q = sqrt(qx*qx+qy*qy);
+  double q_x = qx/q;
+  double q_y = qy/q;
+
+  //convert angle degree to radian
+  theta = theta * M_PI_180;
+  phi = phi * M_PI_180;
+  psi = psi * M_PI_180;
+
+  const double Da = d_factor*dnn;
+  const double s1 = dnn/sqrt(0.75);
+
+
+  //the occupied volume of the lattice
+  const double latticescale = 2.0*sphere_volume(radius/s1);
+  // q vector
+  //q_z = 0.0; // for SANS; assuming qz is negligible
+  /// Angles here are respect to detector coordinate
+  ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
+    // b3 axis orientation
+    b3_x = cos(theta) * cos(phi);
+    b3_y = sin(theta);
+    //b3_z = -cos(theta) * sin(phi);
+    cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z;
+
+    //alpha = acos(cos_val_b3);
+    // b1 axis orientation
+    b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
+    b1_y = sin(psi)*cos(theta);
+    cos_val_b1 = b1_x*q_x + b1_y*q_y;
+    // b2 axis orientation
+    b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
+  	b2_y = cos(theta)*cos(psi);
+    cos_val_b2 = b2_x*q_x + b2_y*q_y;
+
+    // The following test should always pass
+    if (fabs(cos_val_b3)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b3 = 1.0;
+    }
+    if (fabs(cos_val_b2)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b2 = 1.0;
+    }
+    if (fabs(cos_val_b1)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b1 = 1.0;
+    }
+    // Compute the angle btw vector q and the a3 axis
+    a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3);
+
+    // a1 axis
+    a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1);
+
+    // a2 axis
+    a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2);
+
+
+    // Get Fkq and Fkq_2
+    Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
+    Fkq_2 = Fkq*Fkq;
+    // Call Zq=Z1*Z2*Z3
+    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2);
+
+  // Use SphereForm directly from libigor
+  answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale;
+
+  return answer;
+ }
Index: sasmodels/models/bcc_paracrystal.py
===================================================================
--- sasmodels/models/bcc_paracrystal.py	(revision 42356c8890c04b71169941b12a36fc0e07e2034b)
+++ sasmodels/models/bcc_paracrystal.py	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -124,5 +124,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc_paracrystal_kernel.c"]
+source = ["lib/sph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc_paracrystal.c"]
 
 # parameters for demo
Index: smodels/models/bcc_paracrystal_kernel.c
===================================================================
--- sasmodels/models/bcc_paracrystal_kernel.c	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
+++ 	(revision )
@@ -1,170 +1,0 @@
-double form_volume(double radius);
-double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld);
-double Iqxy(double qx, double qy, double dnn,
-    double d_factor, double radius,double sld, double solvent_sld,
-    double theta, double phi, double psi);
-
-double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
-double _BCCeval(double Theta, double Phi, double temp1, double temp3);
-double _sphereform(double q, double radius, double sld, double solvent_sld);
-
-
-double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) {
-
-	const double Da = d_factor*dnn;
-	const double temp1 = q*q*Da*Da;
-	const double temp3 = q*dnn;
-
-	double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI);
-	return(retVal);
-}
-
-double _BCCeval(double Theta, double Phi, double temp1, double temp3) {
-
-	double result;
-	double sin_theta,cos_theta,sin_phi,cos_phi;
-	SINCOS(Theta, sin_theta, cos_theta);
-	SINCOS(Phi, sin_phi, cos_phi);
-
-	const double temp6 =  sin_theta;
-	const double temp7 =  sin_theta*cos_phi + sin_theta*sin_phi + cos_theta;
-	const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta;
-	const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta;
-	const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9));
-	result = pow(1.0-(temp10*temp10),3)*temp6
-	    / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10)
-	      * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10)
-	      * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10));
-
-	return (result);
-}
-
-double form_volume(double radius){
-    return sphere_volume(radius);
-}
-
-
-double Iq(double q, double dnn,
-  double d_factor, double radius,
-  double sld, double solvent_sld){
-
-	//Volume fraction calculated from lattice symmetry and sphere radius
-	const double s1 = dnn/sqrt(0.75);
-	const double latticescale = 2.0*sphere_volume(radius/s1);
-
-    const double va = 0.0;
-    const double vb = 2.0*M_PI;
-    const double vaj = 0.0;
-    const double vbj = M_PI;
-
-    double summ = 0.0;
-    double answer = 0.0;
-	for(int i=0; i<150; i++) {
-		//setup inner integral over the ellipsoidal cross-section
-		double summj=0.0;
-		const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;		//the outer dummy is phi
-		for(int j=0;j<150;j++) {
-			//20 gauss points for the inner integral
-			double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;		//the inner dummy is theta
-			double yyy = Gauss150Wt[j] * _BCC_Integrand(q,dnn,d_factor,ztheta,zphi);
-			summj += yyy;
-		}
-		//now calculate the value of the inner integral
-		double answer = (vbj-vaj)/2.0*summj;
-
-		//now calculate outer integral
-		summ = summ+(Gauss150Wt[i] * answer);
-	}		//final scaling is done at the end of the function, after the NT_FP64 case
-
-	answer = (vb-va)/2.0*summ;
-	answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale;
-
-    return answer;
-
-
-}
-
-
-double Iqxy(double qx, double qy, double dnn,
-    double d_factor, double radius,double sld, double solvent_sld,
-    double theta, double phi, double psi){
-
-  double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
-  //double q_z;
-  double cos_val_b3, cos_val_b2, cos_val_b1;
-  double a1_dot_q, a2_dot_q,a3_dot_q;
-  double answer;
-  double Zq, Fkq, Fkq_2;
-
-  //convert to q and make scaled values
-  double q = sqrt(qx*qx+qy*qy);
-  double q_x = qx/q;
-  double q_y = qy/q;
-
-  //convert angle degree to radian
-  theta = theta * M_PI_180;
-  phi = phi * M_PI_180;
-  psi = psi * M_PI_180;
-
-  const double Da = d_factor*dnn;
-  const double s1 = dnn/sqrt(0.75);
-
-
-  //the occupied volume of the lattice
-  const double latticescale = 2.0*sphere_volume(radius/s1);
-  // q vector
-  //q_z = 0.0; // for SANS; assuming qz is negligible
-  /// Angles here are respect to detector coordinate
-  ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
-    // b3 axis orientation
-    b3_x = cos(theta) * cos(phi);
-    b3_y = sin(theta);
-    //b3_z = -cos(theta) * sin(phi);
-    cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z;
-
-    //alpha = acos(cos_val_b3);
-    // b1 axis orientation
-    b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
-    b1_y = sin(psi)*cos(theta);
-    cos_val_b1 = b1_x*q_x + b1_y*q_y;
-    // b2 axis orientation
-    b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
-  	b2_y = cos(theta)*cos(psi);
-    cos_val_b2 = b2_x*q_x + b2_y*q_y;
-
-    // The following test should always pass
-    if (fabs(cos_val_b3)>1.0) {
-      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
-      cos_val_b3 = 1.0;
-    }
-    if (fabs(cos_val_b2)>1.0) {
-      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
-      cos_val_b2 = 1.0;
-    }
-    if (fabs(cos_val_b1)>1.0) {
-      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
-      cos_val_b1 = 1.0;
-    }
-    // Compute the angle btw vector q and the a3 axis
-    a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3);
-
-    // a1 axis
-    a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1);
-
-    // a2 axis
-    a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2);
-
-
-    // Get Fkq and Fkq_2
-    Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
-    Fkq_2 = Fkq*Fkq;
-    // Call Zq=Z1*Z2*Z3
-    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
-    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
-    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2);
-
-  // Use SphereForm directly from libigor
-  answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale;
-
-  return answer;
- }
Index: sasmodels/models/fcc_paracrystal.c
===================================================================
--- sasmodels/models/fcc_paracrystal.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
+++ sasmodels/models/fcc_paracrystal.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -0,0 +1,170 @@
+double form_volume(double radius);
+double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld);
+double Iqxy(double qx, double qy, double dnn,
+    double d_factor, double radius,double sld, double solvent_sld,
+    double theta, double phi, double psi);
+
+double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
+double _FCCeval(double Theta, double Phi, double temp1, double temp3);
+
+
+double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) {
+
+	const double Da = d_factor*dnn;
+	const double temp1 = q*q*Da*Da;
+	const double temp3 = q*dnn;
+
+	double retVal = _FCCeval(theta,phi,temp1,temp3)/(4.0*M_PI);
+	return(retVal);
+}
+
+double _FCCeval(double Theta, double Phi, double temp1, double temp3) {
+
+	double result;
+        double sin_theta,cos_theta,sin_phi,cos_phi;
+	SINCOS(Theta, sin_theta, cos_theta);
+	SINCOS(Phi, sin_phi, cos_phi);
+
+	const double temp6 =  sin_theta;
+	const double temp7 =  sin_theta*sin_phi + cos_theta;
+	const double temp8 = -sin_theta*cos_phi + cos_theta;
+	const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi;
+
+	const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
+	result = pow((1.0-(temp10*temp10)),3)*temp6
+	    / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10)
+	      * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10)
+	      * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10));
+
+	return (result);
+}
+
+double form_volume(double radius){
+    return sphere_volume(radius);
+}
+
+
+double Iq(double q, double dnn,
+  double d_factor, double radius,
+  double sld, double solvent_sld){
+
+	//Volume fraction calculated from lattice symmetry and sphere radius
+	const double s1 = dnn*sqrt(2.0);
+	const double latticescale = 4.0*sphere_volume(radius/s1);
+
+    const double va = 0.0;
+    const double vb = 2.0*M_PI;
+    const double vaj = 0.0;
+    const double vbj = M_PI;
+
+    double summ = 0.0;
+    double answer = 0.0;
+	for(int i=0; i<150; i++) {
+		//setup inner integral over the ellipsoidal cross-section
+		double summj=0.0;
+		const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;		//the outer dummy is phi
+		for(int j=0;j<150;j++) {
+			//20 gauss points for the inner integral
+			double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;		//the inner dummy is theta
+			double yyy = Gauss150Wt[j] * _FCC_Integrand(q,dnn,d_factor,ztheta,zphi);
+			summj += yyy;
+		}
+		//now calculate the value of the inner integral
+		double answer = (vbj-vaj)/2.0*summj;
+
+		//now calculate outer integral
+		summ = summ+(Gauss150Wt[i] * answer);
+	}		//final scaling is done at the end of the function, after the NT_FP64 case
+
+	answer = (vb-va)/2.0*summ;
+	answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale;
+
+    return answer;
+
+
+}
+
+
+double Iqxy(double qx, double qy, double dnn,
+    double d_factor, double radius,double sld, double solvent_sld,
+    double theta, double phi, double psi){
+
+  double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
+  // double q_z;
+  double cos_val_b3, cos_val_b2, cos_val_b1;
+  double a1_dot_q, a2_dot_q,a3_dot_q;
+  double answer;
+  double Zq, Fkq, Fkq_2;
+
+  //convert to q and make scaled values
+  double q = sqrt(qx*qx+qy*qy);
+  double q_x = qx/q;
+  double q_y = qy/q;
+
+  //convert angle degree to radian
+  theta = theta * M_PI_180;
+  phi = phi * M_PI_180;
+  psi = psi * M_PI_180;
+
+  const double Da = d_factor*dnn;
+  const double s1 = dnn/sqrt(0.75);
+
+
+  //the occupied volume of the lattice
+  const double latticescale = 2.0*sphere_volume(radius/s1);
+  // q vector
+  // q_z = 0.0; // for SANS; assuming qz is negligible
+  /// Angles here are respect to detector coordinate
+  ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
+    // b3 axis orientation
+    b3_x = cos(theta) * cos(phi);
+    b3_y = sin(theta);
+    //b3_z = -cos(theta) * sin(phi);
+    cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z;
+
+    //alpha = acos(cos_val_b3);
+    // b1 axis orientation
+    b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
+    b1_y = sin(psi)*cos(theta);
+    cos_val_b1 = b1_x*q_x + b1_y*q_y;
+    // b2 axis orientation
+    b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
+  	b2_y = cos(theta)*cos(psi);
+    cos_val_b2 = b2_x*q_x + b2_y*q_y;
+
+    // The following test should always pass
+    if (fabs(cos_val_b3)>1.0) {
+      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b3 = 1.0;
+    }
+    if (fabs(cos_val_b2)>1.0) {
+      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b2 = 1.0;
+    }
+    if (fabs(cos_val_b1)>1.0) {
+      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b1 = 1.0;
+    }
+    // Compute the angle btw vector q and the a3 axis
+    a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3);
+
+    // a1 axis
+    a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1);
+
+    // a2 axis
+    a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2);
+
+
+    // Get Fkq and Fkq_2
+    Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
+    Fkq_2 = Fkq*Fkq;
+    // Call Zq=Z1*Z2*Z3
+    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2);
+
+  // Use SphereForm directly from libigor
+  answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale;
+
+  return answer;
+ }
Index: sasmodels/models/fcc_paracrystal.py
===================================================================
--- sasmodels/models/fcc_paracrystal.py	(revision 42356c8890c04b71169941b12a36fc0e07e2034b)
+++ sasmodels/models/fcc_paracrystal.py	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -116,5 +116,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc_paracrystal_kernel.c"]
+source = ["lib/sph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc_paracrystal.c"]
 
 # parameters for demo
Index: smodels/models/fcc_paracrystal_kernel.c
===================================================================
--- sasmodels/models/fcc_paracrystal_kernel.c	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
+++ 	(revision )
@@ -1,170 +1,0 @@
-double form_volume(double radius);
-double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld);
-double Iqxy(double qx, double qy, double dnn,
-    double d_factor, double radius,double sld, double solvent_sld,
-    double theta, double phi, double psi);
-
-double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
-double _FCCeval(double Theta, double Phi, double temp1, double temp3);
-
-
-double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) {
-
-	const double Da = d_factor*dnn;
-	const double temp1 = q*q*Da*Da;
-	const double temp3 = q*dnn;
-
-	double retVal = _FCCeval(theta,phi,temp1,temp3)/(4.0*M_PI);
-	return(retVal);
-}
-
-double _FCCeval(double Theta, double Phi, double temp1, double temp3) {
-
-	double result;
-        double sin_theta,cos_theta,sin_phi,cos_phi;
-	SINCOS(Theta, sin_theta, cos_theta);
-	SINCOS(Phi, sin_phi, cos_phi);
-
-	const double temp6 =  sin_theta;
-	const double temp7 =  sin_theta*sin_phi + cos_theta;
-	const double temp8 = -sin_theta*cos_phi + cos_theta;
-	const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi;
-
-	const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
-	result = pow((1.0-(temp10*temp10)),3)*temp6
-	    / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10)
-	      * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10)
-	      * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10));
-
-	return (result);
-}
-
-double form_volume(double radius){
-    return sphere_volume(radius);
-}
-
-
-double Iq(double q, double dnn,
-  double d_factor, double radius,
-  double sld, double solvent_sld){
-
-	//Volume fraction calculated from lattice symmetry and sphere radius
-	const double s1 = dnn*sqrt(2.0);
-	const double latticescale = 4.0*sphere_volume(radius/s1);
-
-    const double va = 0.0;
-    const double vb = 2.0*M_PI;
-    const double vaj = 0.0;
-    const double vbj = M_PI;
-
-    double summ = 0.0;
-    double answer = 0.0;
-	for(int i=0; i<150; i++) {
-		//setup inner integral over the ellipsoidal cross-section
-		double summj=0.0;
-		const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;		//the outer dummy is phi
-		for(int j=0;j<150;j++) {
-			//20 gauss points for the inner integral
-			double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;		//the inner dummy is theta
-			double yyy = Gauss150Wt[j] * _FCC_Integrand(q,dnn,d_factor,ztheta,zphi);
-			summj += yyy;
-		}
-		//now calculate the value of the inner integral
-		double answer = (vbj-vaj)/2.0*summj;
-
-		//now calculate outer integral
-		summ = summ+(Gauss150Wt[i] * answer);
-	}		//final scaling is done at the end of the function, after the NT_FP64 case
-
-	answer = (vb-va)/2.0*summ;
-	answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale;
-
-    return answer;
-
-
-}
-
-
-double Iqxy(double qx, double qy, double dnn,
-    double d_factor, double radius,double sld, double solvent_sld,
-    double theta, double phi, double psi){
-
-  double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
-  // double q_z;
-  double cos_val_b3, cos_val_b2, cos_val_b1;
-  double a1_dot_q, a2_dot_q,a3_dot_q;
-  double answer;
-  double Zq, Fkq, Fkq_2;
-
-  //convert to q and make scaled values
-  double q = sqrt(qx*qx+qy*qy);
-  double q_x = qx/q;
-  double q_y = qy/q;
-
-  //convert angle degree to radian
-  theta = theta * M_PI_180;
-  phi = phi * M_PI_180;
-  psi = psi * M_PI_180;
-
-  const double Da = d_factor*dnn;
-  const double s1 = dnn/sqrt(0.75);
-
-
-  //the occupied volume of the lattice
-  const double latticescale = 2.0*sphere_volume(radius/s1);
-  // q vector
-  // q_z = 0.0; // for SANS; assuming qz is negligible
-  /// Angles here are respect to detector coordinate
-  ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
-    // b3 axis orientation
-    b3_x = cos(theta) * cos(phi);
-    b3_y = sin(theta);
-    //b3_z = -cos(theta) * sin(phi);
-    cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z;
-
-    //alpha = acos(cos_val_b3);
-    // b1 axis orientation
-    b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
-    b1_y = sin(psi)*cos(theta);
-    cos_val_b1 = b1_x*q_x + b1_y*q_y;
-    // b2 axis orientation
-    b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
-  	b2_y = cos(theta)*cos(psi);
-    cos_val_b2 = b2_x*q_x + b2_y*q_y;
-
-    // The following test should always pass
-    if (fabs(cos_val_b3)>1.0) {
-      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
-      cos_val_b3 = 1.0;
-    }
-    if (fabs(cos_val_b2)>1.0) {
-      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
-      cos_val_b2 = 1.0;
-    }
-    if (fabs(cos_val_b1)>1.0) {
-      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
-      cos_val_b1 = 1.0;
-    }
-    // Compute the angle btw vector q and the a3 axis
-    a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3);
-
-    // a1 axis
-    a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1);
-
-    // a2 axis
-    a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2);
-
-
-    // Get Fkq and Fkq_2
-    Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
-    Fkq_2 = Fkq*Fkq;
-    // Call Zq=Z1*Z2*Z3
-    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
-    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
-    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2);
-
-  // Use SphereForm directly from libigor
-  answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale;
-
-  return answer;
- }
Index: sasmodels/models/hayter_msa.c
===================================================================
--- sasmodels/models/hayter_msa.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
+++ sasmodels/models/hayter_msa.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -0,0 +1,544 @@
+// Hayter-Penfold (rescaled) MSA structure factor for screened Coulomb interactions 
+//
+// C99 needs declarations of routines here
+double Iq(double QQ,
+      double radius_effective, double zz, double VolFrac, double Temp, double csalt, double dialec);
+int
+sqcoef(int ir, double gMSAWave[]);
+
+int
+sqfun(int ix, int ir, double gMSAWave[]);
+
+double
+sqhcal(double qq, double gMSAWave[]);
+  
+double Iq(double QQ,
+      double radius_effective, double VolFrac, double zz, double Temp, double csalt, double dialec)
+{
+    double gMSAWave[17]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17};
+	double Elcharge=1.602189e-19;		// electron charge in Coulombs (C)
+	double kB=1.380662e-23;				// Boltzman constant in J/K
+	double FrSpPerm=8.85418782E-12;	//Permittivity of free space in C^2/(N m^2)
+	double SofQ, Qdiam, Vp, ss;
+	double SIdiam, diam, Kappa, cs, IonSt;
+	double  Perm, Beta;
+	double pi, charge;
+	int ierr;
+	
+	pi = M_PI;
+
+	diam=2*radius_effective;		//in A
+
+						////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+						//////////////////////////// convert to USEFUL inputs in SI units                                                //
+						////////////////////////////    NOTE: easiest to do EVERYTHING in SI units                               //
+						////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+	Beta=1.0/(kB*Temp);		// in Joules^-1
+	Perm=dialec*FrSpPerm;	//in C^2/(N  m^2)
+	charge=zz*Elcharge;		//in Coulomb (C)
+	SIdiam = diam*1.0E-10;		//in m
+	Vp=4.0*pi/3.0*(SIdiam/2.0)*(SIdiam/2.0)*(SIdiam/2.0);	//in m^3
+	cs=csalt*6.022E23*1.0E3;	//# salt molecules/m^3
+	
+	//         Compute the derived values of :
+	//			 Ionic strength IonSt (in C^2/m^3)  
+	// 			Kappa (Debye-Huckel screening length in m)
+	//	and		gamma Exp(-k)
+	
+	// the zz*VolFrac/Vp is for the counterions from the micelle, assumed monovalent, the 2.0*cs if for added salt, assumed 1:1 electolyte 
+	IonSt=0.5 * Elcharge*Elcharge*(zz*VolFrac/Vp+2.0*cs);
+	Kappa=sqrt(2*Beta*IonSt/Perm);     //Kappa calc from Ionic strength
+									   //	Kappa=2/SIdiam					// Use to compare with HP paper
+	gMSAWave[5]=Beta*charge*charge/(pi*Perm*SIdiam*pow((2.0+Kappa*SIdiam),2));
+	
+	//         Finally set up dimensionless parameters 
+	Qdiam=QQ*diam;
+	gMSAWave[6] = Kappa*SIdiam;
+	gMSAWave[4] = VolFrac;
+	
+	//Function sqhpa(qq)  {this is where Hayter-Penfold program began}
+	
+	//       FIRST CALCULATE COUPLING
+	
+	ss=pow(gMSAWave[4],(1.0/3.0));
+	gMSAWave[9] = 2.0*ss*gMSAWave[5]*exp(gMSAWave[6]-gMSAWave[6]/ss);
+	
+	//        CALCULATE COEFFICIENTS, CHECK ALL IS WELL
+	//        AND IF SO CALCULATE S(Q*SIG)
+	
+	ierr=0;
+	ierr=sqcoef(ierr, gMSAWave);
+	if (ierr>=0) {
+		SofQ=sqhcal(Qdiam, gMSAWave);
+	}else{
+       	//SofQ=NaN;
+		SofQ=-1.0;
+		//	print "Error Level = ",ierr
+		//      print "Please report HPMSA problem with above error code"
+	}
+	
+	return(SofQ);
+}
+
+
+
+/////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////
+//
+//
+//      CALCULATES RESCALED VOLUME FRACTION AND CORRESPONDING
+//      COEFFICIENTS FOR "SQHPA"
+//
+//      JOHN B. HAYTER   (I.L.L.)    14-SEP-81
+//
+//      ON EXIT:
+//
+//      SETA IS THE RESCALED VOLUME FRACTION
+//      SGEK IS THE RESCALED CONTACT POTENTIAL
+//      SAK IS THE RESCALED SCREENING CONSTANT
+//      A,B,C,F,U,V ARE THE MSA COEFFICIENTS
+//      G1= G(1+) IS THE CONTACT VALUE OF G(R/SIG):
+//      FOR THE GILLAN CONDITION, THE DIFFERENCE FROM
+//      ZERO INDICATES THE COMPUTATIONAL ACCURACY.
+//
+//      IR > 0:    NORMAL EXIT,  IR IS THE NUMBER OF ITERATIONS.
+//         < 0:    FAILED TO CONVERGE
+//
+int
+sqcoef(int ir, double gMSAWave[])
+{	
+	int itm=40,ix,ig,ii;
+	double acc=5.0E-6,del,e1,e2,f1,f2;
+
+	//      WAVE gMSAWave = $"root:HayPenMSA:gMSAWave"
+	f1=0;		//these were never properly initialized...
+	f2=0;
+	
+	ig=1;
+	if (gMSAWave[6]>=(1.0+8.0*gMSAWave[4])) {
+		ig=0;
+		gMSAWave[15]=gMSAWave[14];
+		gMSAWave[16]=gMSAWave[4];
+		ix=1;
+		ir = sqfun(ix,ir,gMSAWave);
+		gMSAWave[14]=gMSAWave[15];
+		gMSAWave[4]=gMSAWave[16];
+		if((ir<0.0) || (gMSAWave[14]>=0.0)) {
+			return ir;
+		}
+	}
+	gMSAWave[10]=fmin(gMSAWave[4],0.20);
+	if ((ig!=1) || ( gMSAWave[9]>=0.15)) {
+		ii=0;                             
+		do {
+			ii=ii+1;
+			if(ii>itm) {
+				ir=-1;
+				return ir;		
+			}
+			if (gMSAWave[10]<=0.0) {
+			    gMSAWave[10]=gMSAWave[4]/ii;
+			}
+			if(gMSAWave[10]>0.6) {
+			    gMSAWave[10] = 0.35/ii;
+			}
+			e1=gMSAWave[10];
+			gMSAWave[15]=f1;
+			gMSAWave[16]=e1;
+			ix=2;
+			ir = sqfun(ix,ir,gMSAWave);
+			f1=gMSAWave[15];
+			e1=gMSAWave[16];
+			e2=gMSAWave[10]*1.01;
+			gMSAWave[15]=f2;
+			gMSAWave[16]=e2;
+			ix=2;
+			ir = sqfun(ix,ir,gMSAWave);
+			f2=gMSAWave[15];
+			e2=gMSAWave[16];
+			e2=e1-(e2-e1)*f1/(f2-f1);
+			gMSAWave[10] = e2;
+			del = fabs((e2-e1)/e1);
+		} while (del>acc);
+		gMSAWave[15]=gMSAWave[14];
+		gMSAWave[16]=e2;
+		ix=4;
+		ir = sqfun(ix,ir,gMSAWave);
+		gMSAWave[14]=gMSAWave[15];
+		e2=gMSAWave[16];
+		ir=ii;
+		if ((ig!=1) || (gMSAWave[10]>=gMSAWave[4])) {
+		    return ir;
+		}
+	}
+	gMSAWave[15]=gMSAWave[14];
+	gMSAWave[16]=gMSAWave[4];
+	ix=3;
+	ir = sqfun(ix,ir,gMSAWave);
+	gMSAWave[14]=gMSAWave[15];
+	gMSAWave[4]=gMSAWave[16];
+	if ((ir>=0) && (gMSAWave[14]<0.0)) {
+		ir=-3;
+	}
+	return ir;
+}
+
+
+int
+sqfun(int ix, int ir, double gMSAWave[])
+{	
+	double acc=1.0e-6;
+	double reta,eta2,eta21,eta22,eta3,eta32,eta2d,eta2d2,eta3d,eta6d,e12,e24,rgek;
+	double rak,ak1,ak2,dak,dak2,dak4,d,d2,dd2,dd4,dd45,ex1,ex2,sk,ck,ckma,skma;
+	double al1,al2,al3,al4,al5,al6,be1,be2,be3,vu1,vu2,vu3,vu4,vu5,ph1,ph2,ta1,ta2,ta3,ta4,ta5;
+	double a1,a2,a3,b1,b2,b3,v1,v2,v3,p1,p2,p3,pp,pp1,pp2,p1p2,t1,t2,t3,um1,um2,um3,um4,um5,um6;
+	double w0,w1,w2,w3,w4,w12,w13,w14,w15,w16,w24,w25,w26,w32,w34,w3425,w35,w3526,w36,w46,w56;
+	double fa,fap,ca,e24g,pwk,qpw,pg,del,fun,fund,g24;
+	int ii,ibig,itm=40;
+	//      WAVE gMSAWave = $"root:HayPenMSA:gMSAWave"
+	a2=0;
+	a3=0;
+	b2=0;
+	b3=0;
+	v2=0;
+	v3=0;
+	p2=0;
+	p3=0;
+	
+	//     CALCULATE CONSTANTS; NOTATION IS HAYTER PENFOLD (1981)
+	
+	reta = gMSAWave[16];                                                
+	eta2 = reta*reta;
+	eta3 = eta2*reta;
+	e12 = 12.0*reta;
+	e24 = e12+e12;
+	gMSAWave[13] = pow( (gMSAWave[4]/gMSAWave[16]),(1.0/3.0));
+	gMSAWave[12]=gMSAWave[6]/gMSAWave[13];
+	ibig=0;
+	if (( gMSAWave[12]>15.0) && (ix==1)) {
+		ibig=1;
+	}
+    
+	gMSAWave[11] = gMSAWave[5]*gMSAWave[13]*exp(gMSAWave[6]- gMSAWave[12]);
+	rgek =  gMSAWave[11];
+	rak =  gMSAWave[12];
+	ak2 = rak*rak;
+	ak1 = 1.0+rak;
+	dak2 = 1.0/ak2;
+	dak4 = dak2*dak2;
+	d = 1.0-reta;
+	d2 = d*d;
+	dak = d/rak;
+	dd2 = 1.0/d2;
+	dd4 = dd2*dd2;
+	dd45 = dd4*2.0e-1;
+	eta3d=3.0*reta;
+	eta6d = eta3d+eta3d;
+	eta32 = eta3+ eta3;
+	eta2d = reta+2.0;
+	eta2d2 = eta2d*eta2d;
+	eta21 = 2.0*reta+1.0;
+	eta22 = eta21*eta21;
+	
+	//     ALPHA(I)
+	
+	al1 = -eta21*dak;
+	al2 = (14.0*eta2-4.0*reta-1.0)*dak2;
+	al3 = 36.0*eta2*dak4;
+	
+	//      BETA(I)
+	
+	be1 = -(eta2+7.0*reta+1.0)*dak;
+	be2 = 9.0*reta*(eta2+4.0*reta-2.0)*dak2;
+	be3 = 12.0*reta*(2.0*eta2+8.0*reta-1.0)*dak4;
+	
+	//      NU(I)
+	
+	vu1 = -(eta3+3.0*eta2+45.0*reta+5.0)*dak;
+	vu2 = (eta32+3.0*eta2+42.0*reta-2.0e1)*dak2;
+	vu3 = (eta32+3.0e1*reta-5.0)*dak4;
+	vu4 = vu1+e24*rak*vu3;
+	vu5 = eta6d*(vu2+4.0*vu3);
+	
+	//      PHI(I)
+	
+	ph1 = eta6d/rak;
+	ph2 = d-e12*dak2;
+	
+	//      TAU(I)
+	
+	ta1 = (reta+5.0)/(5.0*rak);
+	ta2 = eta2d*dak2;
+	ta3 = -e12*rgek*(ta1+ta2);
+	ta4 = eta3d*ak2*(ta1*ta1-ta2*ta2);
+	ta5 = eta3d*(reta+8.0)*1.0e-1-2.0*eta22*dak2;
+	
+	//     double PRECISION SINH(K), COSH(K)
+	
+	ex1 = exp(rak);
+	ex2 = 0.0;
+	if ( gMSAWave[12]<20.0) {
+		ex2=exp(-rak);
+	}
+	sk=0.5*(ex1-ex2);
+	ck = 0.5*(ex1+ex2);
+	ckma = ck-1.0-rak*sk;
+	skma = sk-rak*ck;
+	
+	//      a(I)
+	
+	a1 = (e24*rgek*(al1+al2+ak1*al3)-eta22)*dd4;
+	if (ibig==0) {
+		a2 = e24*(al3*skma+al2*sk-al1*ck)*dd4;
+		a3 = e24*(eta22*dak2-0.5*d2+al3*ckma-al1*sk+al2*ck)*dd4;
+	}
+	
+	//      b(I)
+	
+	b1 = (1.5*reta*eta2d2-e12*rgek*(be1+be2+ak1*be3))*dd4;
+	if (ibig==0) {
+		b2 = e12*(-be3*skma-be2*sk+be1*ck)*dd4;
+		b3 = e12*(0.5*d2*eta2d-eta3d*eta2d2*dak2-be3*ckma+be1*sk-be2*ck)*dd4;
+	}
+	
+	//      V(I)
+	
+	v1 = (eta21*(eta2-2.0*reta+1.0e1)*2.5e-1-rgek*(vu4+vu5))*dd45;
+	if (ibig==0) {
+		v2 = (vu4*ck-vu5*sk)*dd45;
+		v3 = ((eta3-6.0*eta2+5.0)*d-eta6d*(2.0*eta3-3.0*eta2+18.0*reta+1.0e1)*dak2+e24*vu3+vu4*sk-vu5*ck)*dd45;
+	}
+	
+	
+	//       P(I)
+	
+	pp1 = ph1*ph1;
+	pp2 = ph2*ph2;
+	pp = pp1+pp2;
+	p1p2 = ph1*ph2*2.0;
+	p1 = (rgek*(pp1+pp2-p1p2)-0.5*eta2d)*dd2;
+	if (ibig==0) {
+		p2 = (pp*sk+p1p2*ck)*dd2;
+		p3 = (pp*ck+p1p2*sk+pp1-pp2)*dd2;
+	}
+	
+	//       T(I)
+	
+	t1 = ta3+ta4*a1+ta5*b1;
+	if (ibig!=0) {
+		
+		//		VERY LARGE SCREENING:  ASYMPTOTIC SOLUTION
+		
+  		v3 = ((eta3-6.0*eta2+5.0)*d-eta6d*(2.0*eta3-3.0*eta2+18.0*reta+1.0e1)*dak2+e24*vu3)*dd45;
+		t3 = ta4*a3+ta5*b3+e12*ta2 - 4.0e-1*reta*(reta+1.0e1)-1.0;
+		p3 = (pp1-pp2)*dd2;
+		b3 = e12*(0.5*d2*eta2d-eta3d*eta2d2*dak2+be3)*dd4;
+		a3 = e24*(eta22*dak2-0.5*d2-al3)*dd4;
+		um6 = t3*a3-e12*v3*v3;
+		um5 = t1*a3+a1*t3-e24*v1*v3;
+		um4 = t1*a1-e12*v1*v1;
+		al6 = e12*p3*p3;
+		al5 = e24*p1*p3-b3-b3-ak2;
+		al4 = e12*p1*p1-b1-b1;
+		w56 = um5*al6-al5*um6;
+		w46 = um4*al6-al4*um6;
+		fa = -w46/w56;
+		ca = -fa;
+		gMSAWave[3] = fa;
+		gMSAWave[2] = ca;
+		gMSAWave[1] = b1+b3*fa;
+		gMSAWave[0] = a1+a3*fa;
+		gMSAWave[8] = v1+v3*fa;
+		gMSAWave[14] = -(p1+p3*fa);
+		gMSAWave[15] = gMSAWave[14];
+		if (fabs(gMSAWave[15])<1.0e-3) {
+			gMSAWave[15] = 0.0;
+		}
+		gMSAWave[10] = gMSAWave[16];
+		
+	} else {
+        
+		t2 = ta4*a2+ta5*b2+e12*(ta1*ck-ta2*sk);
+		t3 = ta4*a3+ta5*b3+e12*(ta1*sk-ta2*(ck-1.0))-4.0e-1*reta*(reta+1.0e1)-1.0;
+		
+		//		MU(i)
+		
+		um1 = t2*a2-e12*v2*v2;
+		um2 = t1*a2+t2*a1-e24*v1*v2;
+		um3 = t2*a3+t3*a2-e24*v2*v3;
+		um4 = t1*a1-e12*v1*v1;
+		um5 = t1*a3+t3*a1-e24*v1*v3;
+		um6 = t3*a3-e12*v3*v3;
+		
+		//			GILLAN CONDITION ?
+		//
+		//			YES - G(X=1+) = 0
+		//
+		//			COEFFICIENTS AND FUNCTION VALUE
+		//
+		if ((ix==1) || (ix==3)) {
+			
+			//			NO - CALCULATE REMAINING COEFFICIENTS.
+			
+			//			LAMBDA(I)
+			
+			al1 = e12*p2*p2;
+			al2 = e24*p1*p2-b2-b2;
+			al3 = e24*p2*p3;
+			al4 = e12*p1*p1-b1-b1;
+			al5 = e24*p1*p3-b3-b3-ak2;
+			al6 = e12*p3*p3;
+			
+			//			OMEGA(I)
+			
+			w16 = um1*al6-al1*um6;
+			w15 = um1*al5-al1*um5;
+			w14 = um1*al4-al1*um4;
+			w13 = um1*al3-al1*um3;
+			w12 = um1*al2-al1*um2;
+			
+			w26 = um2*al6-al2*um6;
+			w25 = um2*al5-al2*um5;
+			w24 = um2*al4-al2*um4;
+			
+			w36 = um3*al6-al3*um6;
+			w35 = um3*al5-al3*um5;
+			w34 = um3*al4-al3*um4;
+			w32 = um3*al2-al3*um2;
+			//
+			w46 = um4*al6-al4*um6;
+			w56 = um5*al6-al5*um6;
+			w3526 = w35+w26;
+			w3425 = w34+w25;
+			
+			//			QUARTIC COEFFICIENTS
+			
+			w4 = w16*w16-w13*w36;
+			w3 = 2.0*w16*w15-w13*w3526-w12*w36;
+			w2 = w15*w15+2.0*w16*w14-w13*w3425-w12*w3526;
+			w1 = 2.0*w15*w14-w13*w24-w12*w3425;
+			w0 = w14*w14-w12*w24;
+			
+			//			ESTIMATE THE STARTING VALUE OF f
+			
+			if (ix==1) {
+				//				LARGE K
+				fap = (w14-w34-w46)/(w12-w15+w35-w26+w56-w32);
+			} else {
+				//				ASSUME NOT TOO FAR FROM GILLAN CONDITION.
+				//				IF BOTH RGEK AND RAK ARE SMALL, USE P-W ESTIMATE.
+				gMSAWave[14]=0.5*eta2d*dd2*exp(-rgek);
+				if (( gMSAWave[11]<=2.0) && ( gMSAWave[11]>=0.0) && ( gMSAWave[12]<=1.0)) {
+					e24g = e24*rgek*exp(rak);
+					pwk = sqrt(e24g);
+					qpw = (1.0-sqrt(1.0+2.0*d2*d*pwk/eta22))*eta21/d;
+					gMSAWave[14] = -qpw*qpw/e24+0.5*eta2d*dd2;
+				}
+  				pg = p1+gMSAWave[14];
+				ca = ak2*pg+2.0*(b3*pg-b1*p3)+e12*gMSAWave[14]*gMSAWave[14]*p3;
+				ca = -ca/(ak2*p2+2.0*(b3*p2-b2*p3));
+				fap = -(pg+p2*ca)/p3;
+			}
+			
+			//			AND REFINE IT ACCORDING TO NEWTON
+			ii=0;
+			do {
+				ii = ii+1;
+				if (ii>itm) {
+					//					FAILED TO CONVERGE IN ITM ITERATIONS
+					ir=-2;
+					return (ir);
+				}
+				fa = fap;
+				fun = w0+(w1+(w2+(w3+w4*fa)*fa)*fa)*fa;
+				fund = w1+(2.0*w2+(3.0*w3+4.0*w4*fa)*fa)*fa;
+				fap = fa-fun/fund;
+				del=fabs((fap-fa)/fa);
+			} while (del>acc);
+			
+			ir = ir+ii;
+			fa = fap;
+			ca = -(w16*fa*fa+w15*fa+w14)/(w13*fa+w12);
+			gMSAWave[14] = -(p1+p2*ca+p3*fa);
+			gMSAWave[15] = gMSAWave[14];
+			if (fabs(gMSAWave[15])<1.0e-3) {
+				gMSAWave[15] = 0.0;
+			}
+			gMSAWave[10] = gMSAWave[16];
+		} else {
+			ca = ak2*p1+2.0*(b3*p1-b1*p3);
+			ca = -ca/(ak2*p2+2.0*(b3*p2-b2*p3));
+			fa = -(p1+p2*ca)/p3;
+			if (ix==2) {
+				gMSAWave[15] = um1*ca*ca+(um2+um3*fa)*ca+um4+um5*fa+um6*fa*fa;
+			}
+			if (ix==4) {
+				gMSAWave[15] = -(p1+p2*ca+p3*fa);
+			}
+		}
+   		gMSAWave[3] = fa;
+		gMSAWave[2] = ca;
+		gMSAWave[1] = b1+b2*ca+b3*fa;
+		gMSAWave[0] = a1+a2*ca+a3*fa;
+		gMSAWave[8] = (v1+v2*ca+v3*fa)/gMSAWave[0];
+	}
+   	g24 = e24*rgek*ex1;
+	gMSAWave[7] = (rak*ak2*ca-g24)/(ak2*g24);
+	return (ir);
+}
+
+double
+sqhcal(double qq, double gMSAWave[])
+{      	
+    double SofQ,etaz,akz,gekz,e24,x1,x2,ck,sk,ak2,qk,q2k,qk2,qk3,qqk,sink,cosk,asink,qcosk,aqk,inter; 		
+	//	WAVE gMSAWave = $"root:HayPenMSA:gMSAWave"
+
+	etaz = gMSAWave[10];
+	akz =  gMSAWave[12];
+	gekz =  gMSAWave[11];
+	e24 = 24.0*etaz;
+	x1 = exp(akz);
+	x2 = 0.0;
+	if ( gMSAWave[12]<20.0) {
+		x2 = exp(-akz);
+	}
+	ck = 0.5*(x1+x2);
+	sk = 0.5*(x1-x2);
+	ak2 = akz*akz;
+	
+	qk = qq/gMSAWave[13];
+	q2k = qk*qk;
+	if (qk<=1.0e-08) {
+		SofQ = -1.0/gMSAWave[0];
+	} else {
+	// this rescales Q.sigma = 2.Q.Radius, so is hard to predict the value to test the function
+	if (qk<=0.01) {
+		// try Taylor series expansion at small qk (RKH Feb 2016, with help from Mathematica), 
+		// transition point may need to depend on precision of cpu used and ALAS on the values of some of the parameters !
+		// note have adsorbed a factor 24 from SofQ=
+		// needs thorough test over wide range of parameter space!
+		// there seem to be some rounding issues here in single precision, must use double
+		aqk = gMSAWave[0]*(8.0+2.0*etaz) + 6*gMSAWave[1] -12.0*gMSAWave[3] 
+			-24*(gekz*(1.0+akz) -ck*akz*gMSAWave[2] +gMSAWave[3]*(ck-1.0) +(gMSAWave[2]-gMSAWave[3]*akz)*sk )/ak2
+			+q2k*( -(gMSAWave[0]*(48.0+15.0*etaz) +40.0*gMSAWave[1])/60.0 +gMSAWave[3] 
+			+(4.0/ak2)*(gekz*(9.0+7.0*akz) +ck*(9.0*gMSAWave[3] -7.0*gMSAWave[2]*akz) +sk*(9.0*gMSAWave[2] -7.0*gMSAWave[3]*akz)) );
+		SofQ = 1.0/(1.0-gMSAWave[10]*aqk);
+	} else {
+		qk2 = 1.0/q2k;
+		qk3 = qk2/qk;
+		qqk = 1.0/(qk*(q2k+ak2));
+		SINCOS(qk,sink,cosk);
+		asink = akz*sink;
+		qcosk = qk*cosk;
+		aqk = gMSAWave[0]*(sink-qcosk);
+		aqk=aqk+gMSAWave[1]*((2.0*qk2-1.0)*qcosk+2.0*sink-2.0/qk);
+		inter=24.0*qk3+4.0*(1.0-6.0*qk2)*sink;
+		aqk=(aqk+0.5*etaz*gMSAWave[0]*(inter-(1.0-12.0*qk2+24.0*qk2*qk2)*qcosk))*qk3;
+		aqk=aqk +gMSAWave[2]*(ck*asink-sk*qcosk)*qqk;
+		aqk=aqk +gMSAWave[3]*(sk*asink-qk*(ck*cosk-1.0))*qqk;
+		aqk=aqk +gMSAWave[3]*(cosk-1.0)*qk2;
+		aqk=aqk -gekz*(asink+qcosk)*qqk;
+		SofQ = 1.0/(1.0  -e24*aqk);
+	} }
+	return (SofQ);
+}
Index: sasmodels/models/hayter_msa.py
===================================================================
--- sasmodels/models/hayter_msa.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/models/hayter_msa.py	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -81,5 +81,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["hayter_msa_kernel.c"]
+source = ["hayter_msa.c"]
 # No volume normalization despite having a volume parameter
 # This should perhaps be volume normalized?
Index: smodels/models/hayter_msa_kernel.c
===================================================================
--- sasmodels/models/hayter_msa_kernel.c	(revision ce43de0146d69b5f356808e50821fbfa3cf1d0fa)
+++ 	(revision )
@@ -1,544 +1,0 @@
-// Hayter-Penfold (rescaled) MSA structure factor for screened Coulomb interactions 
-//
-// C99 needs declarations of routines here
-double Iq(double QQ,
-      double radius_effective, double zz, double VolFrac, double Temp, double csalt, double dialec);
-int
-sqcoef(int ir, double gMSAWave[]);
-
-int
-sqfun(int ix, int ir, double gMSAWave[]);
-
-double
-sqhcal(double qq, double gMSAWave[]);
-  
-double Iq(double QQ,
-      double radius_effective, double VolFrac, double zz, double Temp, double csalt, double dialec)
-{
-    double gMSAWave[17]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17};
-	double Elcharge=1.602189e-19;		// electron charge in Coulombs (C)
-	double kB=1.380662e-23;				// Boltzman constant in J/K
-	double FrSpPerm=8.85418782E-12;	//Permittivity of free space in C^2/(N m^2)
-	double SofQ, Qdiam, Vp, ss;
-	double SIdiam, diam, Kappa, cs, IonSt;
-	double  Perm, Beta;
-	double pi, charge;
-	int ierr;
-	
-	pi = M_PI;
-
-	diam=2*radius_effective;		//in A
-
-						////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-						//////////////////////////// convert to USEFUL inputs in SI units                                                //
-						////////////////////////////    NOTE: easiest to do EVERYTHING in SI units                               //
-						////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-	Beta=1.0/(kB*Temp);		// in Joules^-1
-	Perm=dialec*FrSpPerm;	//in C^2/(N  m^2)
-	charge=zz*Elcharge;		//in Coulomb (C)
-	SIdiam = diam*1.0E-10;		//in m
-	Vp=4.0*pi/3.0*(SIdiam/2.0)*(SIdiam/2.0)*(SIdiam/2.0);	//in m^3
-	cs=csalt*6.022E23*1.0E3;	//# salt molecules/m^3
-	
-	//         Compute the derived values of :
-	//			 Ionic strength IonSt (in C^2/m^3)  
-	// 			Kappa (Debye-Huckel screening length in m)
-	//	and		gamma Exp(-k)
-	
-	// the zz*VolFrac/Vp is for the counterions from the micelle, assumed monovalent, the 2.0*cs if for added salt, assumed 1:1 electolyte 
-	IonSt=0.5 * Elcharge*Elcharge*(zz*VolFrac/Vp+2.0*cs);
-	Kappa=sqrt(2*Beta*IonSt/Perm);     //Kappa calc from Ionic strength
-									   //	Kappa=2/SIdiam					// Use to compare with HP paper
-	gMSAWave[5]=Beta*charge*charge/(pi*Perm*SIdiam*pow((2.0+Kappa*SIdiam),2));
-	
-	//         Finally set up dimensionless parameters 
-	Qdiam=QQ*diam;
-	gMSAWave[6] = Kappa*SIdiam;
-	gMSAWave[4] = VolFrac;
-	
-	//Function sqhpa(qq)  {this is where Hayter-Penfold program began}
-	
-	//       FIRST CALCULATE COUPLING
-	
-	ss=pow(gMSAWave[4],(1.0/3.0));
-	gMSAWave[9] = 2.0*ss*gMSAWave[5]*exp(gMSAWave[6]-gMSAWave[6]/ss);
-	
-	//        CALCULATE COEFFICIENTS, CHECK ALL IS WELL
-	//        AND IF SO CALCULATE S(Q*SIG)
-	
-	ierr=0;
-	ierr=sqcoef(ierr, gMSAWave);
-	if (ierr>=0) {
-		SofQ=sqhcal(Qdiam, gMSAWave);
-	}else{
-       	//SofQ=NaN;
-		SofQ=-1.0;
-		//	print "Error Level = ",ierr
-		//      print "Please report HPMSA problem with above error code"
-	}
-	
-	return(SofQ);
-}
-
-
-
-/////////////////////////////////////////////////////////////
-/////////////////////////////////////////////////////////////
-//
-//
-//      CALCULATES RESCALED VOLUME FRACTION AND CORRESPONDING
-//      COEFFICIENTS FOR "SQHPA"
-//
-//      JOHN B. HAYTER   (I.L.L.)    14-SEP-81
-//
-//      ON EXIT:
-//
-//      SETA IS THE RESCALED VOLUME FRACTION
-//      SGEK IS THE RESCALED CONTACT POTENTIAL
-//      SAK IS THE RESCALED SCREENING CONSTANT
-//      A,B,C,F,U,V ARE THE MSA COEFFICIENTS
-//      G1= G(1+) IS THE CONTACT VALUE OF G(R/SIG):
-//      FOR THE GILLAN CONDITION, THE DIFFERENCE FROM
-//      ZERO INDICATES THE COMPUTATIONAL ACCURACY.
-//
-//      IR > 0:    NORMAL EXIT,  IR IS THE NUMBER OF ITERATIONS.
-//         < 0:    FAILED TO CONVERGE
-//
-int
-sqcoef(int ir, double gMSAWave[])
-{	
-	int itm=40,ix,ig,ii;
-	double acc=5.0E-6,del,e1,e2,f1,f2;
-
-	//      WAVE gMSAWave = $"root:HayPenMSA:gMSAWave"
-	f1=0;		//these were never properly initialized...
-	f2=0;
-	
-	ig=1;
-	if (gMSAWave[6]>=(1.0+8.0*gMSAWave[4])) {
-		ig=0;
-		gMSAWave[15]=gMSAWave[14];
-		gMSAWave[16]=gMSAWave[4];
-		ix=1;
-		ir = sqfun(ix,ir,gMSAWave);
-		gMSAWave[14]=gMSAWave[15];
-		gMSAWave[4]=gMSAWave[16];
-		if((ir<0.0) || (gMSAWave[14]>=0.0)) {
-			return ir;
-		}
-	}
-	gMSAWave[10]=fmin(gMSAWave[4],0.20);
-	if ((ig!=1) || ( gMSAWave[9]>=0.15)) {
-		ii=0;                             
-		do {
-			ii=ii+1;
-			if(ii>itm) {
-				ir=-1;
-				return ir;		
-			}
-			if (gMSAWave[10]<=0.0) {
-			    gMSAWave[10]=gMSAWave[4]/ii;
-			}
-			if(gMSAWave[10]>0.6) {
-			    gMSAWave[10] = 0.35/ii;
-			}
-			e1=gMSAWave[10];
-			gMSAWave[15]=f1;
-			gMSAWave[16]=e1;
-			ix=2;
-			ir = sqfun(ix,ir,gMSAWave);
-			f1=gMSAWave[15];
-			e1=gMSAWave[16];
-			e2=gMSAWave[10]*1.01;
-			gMSAWave[15]=f2;
-			gMSAWave[16]=e2;
-			ix=2;
-			ir = sqfun(ix,ir,gMSAWave);
-			f2=gMSAWave[15];
-			e2=gMSAWave[16];
-			e2=e1-(e2-e1)*f1/(f2-f1);
-			gMSAWave[10] = e2;
-			del = fabs((e2-e1)/e1);
-		} while (del>acc);
-		gMSAWave[15]=gMSAWave[14];
-		gMSAWave[16]=e2;
-		ix=4;
-		ir = sqfun(ix,ir,gMSAWave);
-		gMSAWave[14]=gMSAWave[15];
-		e2=gMSAWave[16];
-		ir=ii;
-		if ((ig!=1) || (gMSAWave[10]>=gMSAWave[4])) {
-		    return ir;
-		}
-	}
-	gMSAWave[15]=gMSAWave[14];
-	gMSAWave[16]=gMSAWave[4];
-	ix=3;
-	ir = sqfun(ix,ir,gMSAWave);
-	gMSAWave[14]=gMSAWave[15];
-	gMSAWave[4]=gMSAWave[16];
-	if ((ir>=0) && (gMSAWave[14]<0.0)) {
-		ir=-3;
-	}
-	return ir;
-}
-
-
-int
-sqfun(int ix, int ir, double gMSAWave[])
-{	
-	double acc=1.0e-6;
-	double reta,eta2,eta21,eta22,eta3,eta32,eta2d,eta2d2,eta3d,eta6d,e12,e24,rgek;
-	double rak,ak1,ak2,dak,dak2,dak4,d,d2,dd2,dd4,dd45,ex1,ex2,sk,ck,ckma,skma;
-	double al1,al2,al3,al4,al5,al6,be1,be2,be3,vu1,vu2,vu3,vu4,vu5,ph1,ph2,ta1,ta2,ta3,ta4,ta5;
-	double a1,a2,a3,b1,b2,b3,v1,v2,v3,p1,p2,p3,pp,pp1,pp2,p1p2,t1,t2,t3,um1,um2,um3,um4,um5,um6;
-	double w0,w1,w2,w3,w4,w12,w13,w14,w15,w16,w24,w25,w26,w32,w34,w3425,w35,w3526,w36,w46,w56;
-	double fa,fap,ca,e24g,pwk,qpw,pg,del,fun,fund,g24;
-	int ii,ibig,itm=40;
-	//      WAVE gMSAWave = $"root:HayPenMSA:gMSAWave"
-	a2=0;
-	a3=0;
-	b2=0;
-	b3=0;
-	v2=0;
-	v3=0;
-	p2=0;
-	p3=0;
-	
-	//     CALCULATE CONSTANTS; NOTATION IS HAYTER PENFOLD (1981)
-	
-	reta = gMSAWave[16];                                                
-	eta2 = reta*reta;
-	eta3 = eta2*reta;
-	e12 = 12.0*reta;
-	e24 = e12+e12;
-	gMSAWave[13] = pow( (gMSAWave[4]/gMSAWave[16]),(1.0/3.0));
-	gMSAWave[12]=gMSAWave[6]/gMSAWave[13];
-	ibig=0;
-	if (( gMSAWave[12]>15.0) && (ix==1)) {
-		ibig=1;
-	}
-    
-	gMSAWave[11] = gMSAWave[5]*gMSAWave[13]*exp(gMSAWave[6]- gMSAWave[12]);
-	rgek =  gMSAWave[11];
-	rak =  gMSAWave[12];
-	ak2 = rak*rak;
-	ak1 = 1.0+rak;
-	dak2 = 1.0/ak2;
-	dak4 = dak2*dak2;
-	d = 1.0-reta;
-	d2 = d*d;
-	dak = d/rak;
-	dd2 = 1.0/d2;
-	dd4 = dd2*dd2;
-	dd45 = dd4*2.0e-1;
-	eta3d=3.0*reta;
-	eta6d = eta3d+eta3d;
-	eta32 = eta3+ eta3;
-	eta2d = reta+2.0;
-	eta2d2 = eta2d*eta2d;
-	eta21 = 2.0*reta+1.0;
-	eta22 = eta21*eta21;
-	
-	//     ALPHA(I)
-	
-	al1 = -eta21*dak;
-	al2 = (14.0*eta2-4.0*reta-1.0)*dak2;
-	al3 = 36.0*eta2*dak4;
-	
-	//      BETA(I)
-	
-	be1 = -(eta2+7.0*reta+1.0)*dak;
-	be2 = 9.0*reta*(eta2+4.0*reta-2.0)*dak2;
-	be3 = 12.0*reta*(2.0*eta2+8.0*reta-1.0)*dak4;
-	
-	//      NU(I)
-	
-	vu1 = -(eta3+3.0*eta2+45.0*reta+5.0)*dak;
-	vu2 = (eta32+3.0*eta2+42.0*reta-2.0e1)*dak2;
-	vu3 = (eta32+3.0e1*reta-5.0)*dak4;
-	vu4 = vu1+e24*rak*vu3;
-	vu5 = eta6d*(vu2+4.0*vu3);
-	
-	//      PHI(I)
-	
-	ph1 = eta6d/rak;
-	ph2 = d-e12*dak2;
-	
-	//      TAU(I)
-	
-	ta1 = (reta+5.0)/(5.0*rak);
-	ta2 = eta2d*dak2;
-	ta3 = -e12*rgek*(ta1+ta2);
-	ta4 = eta3d*ak2*(ta1*ta1-ta2*ta2);
-	ta5 = eta3d*(reta+8.0)*1.0e-1-2.0*eta22*dak2;
-	
-	//     double PRECISION SINH(K), COSH(K)
-	
-	ex1 = exp(rak);
-	ex2 = 0.0;
-	if ( gMSAWave[12]<20.0) {
-		ex2=exp(-rak);
-	}
-	sk=0.5*(ex1-ex2);
-	ck = 0.5*(ex1+ex2);
-	ckma = ck-1.0-rak*sk;
-	skma = sk-rak*ck;
-	
-	//      a(I)
-	
-	a1 = (e24*rgek*(al1+al2+ak1*al3)-eta22)*dd4;
-	if (ibig==0) {
-		a2 = e24*(al3*skma+al2*sk-al1*ck)*dd4;
-		a3 = e24*(eta22*dak2-0.5*d2+al3*ckma-al1*sk+al2*ck)*dd4;
-	}
-	
-	//      b(I)
-	
-	b1 = (1.5*reta*eta2d2-e12*rgek*(be1+be2+ak1*be3))*dd4;
-	if (ibig==0) {
-		b2 = e12*(-be3*skma-be2*sk+be1*ck)*dd4;
-		b3 = e12*(0.5*d2*eta2d-eta3d*eta2d2*dak2-be3*ckma+be1*sk-be2*ck)*dd4;
-	}
-	
-	//      V(I)
-	
-	v1 = (eta21*(eta2-2.0*reta+1.0e1)*2.5e-1-rgek*(vu4+vu5))*dd45;
-	if (ibig==0) {
-		v2 = (vu4*ck-vu5*sk)*dd45;
-		v3 = ((eta3-6.0*eta2+5.0)*d-eta6d*(2.0*eta3-3.0*eta2+18.0*reta+1.0e1)*dak2+e24*vu3+vu4*sk-vu5*ck)*dd45;
-	}
-	
-	
-	//       P(I)
-	
-	pp1 = ph1*ph1;
-	pp2 = ph2*ph2;
-	pp = pp1+pp2;
-	p1p2 = ph1*ph2*2.0;
-	p1 = (rgek*(pp1+pp2-p1p2)-0.5*eta2d)*dd2;
-	if (ibig==0) {
-		p2 = (pp*sk+p1p2*ck)*dd2;
-		p3 = (pp*ck+p1p2*sk+pp1-pp2)*dd2;
-	}
-	
-	//       T(I)
-	
-	t1 = ta3+ta4*a1+ta5*b1;
-	if (ibig!=0) {
-		
-		//		VERY LARGE SCREENING:  ASYMPTOTIC SOLUTION
-		
-  		v3 = ((eta3-6.0*eta2+5.0)*d-eta6d*(2.0*eta3-3.0*eta2+18.0*reta+1.0e1)*dak2+e24*vu3)*dd45;
-		t3 = ta4*a3+ta5*b3+e12*ta2 - 4.0e-1*reta*(reta+1.0e1)-1.0;
-		p3 = (pp1-pp2)*dd2;
-		b3 = e12*(0.5*d2*eta2d-eta3d*eta2d2*dak2+be3)*dd4;
-		a3 = e24*(eta22*dak2-0.5*d2-al3)*dd4;
-		um6 = t3*a3-e12*v3*v3;
-		um5 = t1*a3+a1*t3-e24*v1*v3;
-		um4 = t1*a1-e12*v1*v1;
-		al6 = e12*p3*p3;
-		al5 = e24*p1*p3-b3-b3-ak2;
-		al4 = e12*p1*p1-b1-b1;
-		w56 = um5*al6-al5*um6;
-		w46 = um4*al6-al4*um6;
-		fa = -w46/w56;
-		ca = -fa;
-		gMSAWave[3] = fa;
-		gMSAWave[2] = ca;
-		gMSAWave[1] = b1+b3*fa;
-		gMSAWave[0] = a1+a3*fa;
-		gMSAWave[8] = v1+v3*fa;
-		gMSAWave[14] = -(p1+p3*fa);
-		gMSAWave[15] = gMSAWave[14];
-		if (fabs(gMSAWave[15])<1.0e-3) {
-			gMSAWave[15] = 0.0;
-		}
-		gMSAWave[10] = gMSAWave[16];
-		
-	} else {
-        
-		t2 = ta4*a2+ta5*b2+e12*(ta1*ck-ta2*sk);
-		t3 = ta4*a3+ta5*b3+e12*(ta1*sk-ta2*(ck-1.0))-4.0e-1*reta*(reta+1.0e1)-1.0;
-		
-		//		MU(i)
-		
-		um1 = t2*a2-e12*v2*v2;
-		um2 = t1*a2+t2*a1-e24*v1*v2;
-		um3 = t2*a3+t3*a2-e24*v2*v3;
-		um4 = t1*a1-e12*v1*v1;
-		um5 = t1*a3+t3*a1-e24*v1*v3;
-		um6 = t3*a3-e12*v3*v3;
-		
-		//			GILLAN CONDITION ?
-		//
-		//			YES - G(X=1+) = 0
-		//
-		//			COEFFICIENTS AND FUNCTION VALUE
-		//
-		if ((ix==1) || (ix==3)) {
-			
-			//			NO - CALCULATE REMAINING COEFFICIENTS.
-			
-			//			LAMBDA(I)
-			
-			al1 = e12*p2*p2;
-			al2 = e24*p1*p2-b2-b2;
-			al3 = e24*p2*p3;
-			al4 = e12*p1*p1-b1-b1;
-			al5 = e24*p1*p3-b3-b3-ak2;
-			al6 = e12*p3*p3;
-			
-			//			OMEGA(I)
-			
-			w16 = um1*al6-al1*um6;
-			w15 = um1*al5-al1*um5;
-			w14 = um1*al4-al1*um4;
-			w13 = um1*al3-al1*um3;
-			w12 = um1*al2-al1*um2;
-			
-			w26 = um2*al6-al2*um6;
-			w25 = um2*al5-al2*um5;
-			w24 = um2*al4-al2*um4;
-			
-			w36 = um3*al6-al3*um6;
-			w35 = um3*al5-al3*um5;
-			w34 = um3*al4-al3*um4;
-			w32 = um3*al2-al3*um2;
-			//
-			w46 = um4*al6-al4*um6;
-			w56 = um5*al6-al5*um6;
-			w3526 = w35+w26;
-			w3425 = w34+w25;
-			
-			//			QUARTIC COEFFICIENTS
-			
-			w4 = w16*w16-w13*w36;
-			w3 = 2.0*w16*w15-w13*w3526-w12*w36;
-			w2 = w15*w15+2.0*w16*w14-w13*w3425-w12*w3526;
-			w1 = 2.0*w15*w14-w13*w24-w12*w3425;
-			w0 = w14*w14-w12*w24;
-			
-			//			ESTIMATE THE STARTING VALUE OF f
-			
-			if (ix==1) {
-				//				LARGE K
-				fap = (w14-w34-w46)/(w12-w15+w35-w26+w56-w32);
-			} else {
-				//				ASSUME NOT TOO FAR FROM GILLAN CONDITION.
-				//				IF BOTH RGEK AND RAK ARE SMALL, USE P-W ESTIMATE.
-				gMSAWave[14]=0.5*eta2d*dd2*exp(-rgek);
-				if (( gMSAWave[11]<=2.0) && ( gMSAWave[11]>=0.0) && ( gMSAWave[12]<=1.0)) {
-					e24g = e24*rgek*exp(rak);
-					pwk = sqrt(e24g);
-					qpw = (1.0-sqrt(1.0+2.0*d2*d*pwk/eta22))*eta21/d;
-					gMSAWave[14] = -qpw*qpw/e24+0.5*eta2d*dd2;
-				}
-  				pg = p1+gMSAWave[14];
-				ca = ak2*pg+2.0*(b3*pg-b1*p3)+e12*gMSAWave[14]*gMSAWave[14]*p3;
-				ca = -ca/(ak2*p2+2.0*(b3*p2-b2*p3));
-				fap = -(pg+p2*ca)/p3;
-			}
-			
-			//			AND REFINE IT ACCORDING TO NEWTON
-			ii=0;
-			do {
-				ii = ii+1;
-				if (ii>itm) {
-					//					FAILED TO CONVERGE IN ITM ITERATIONS
-					ir=-2;
-					return (ir);
-				}
-				fa = fap;
-				fun = w0+(w1+(w2+(w3+w4*fa)*fa)*fa)*fa;
-				fund = w1+(2.0*w2+(3.0*w3+4.0*w4*fa)*fa)*fa;
-				fap = fa-fun/fund;
-				del=fabs((fap-fa)/fa);
-			} while (del>acc);
-			
-			ir = ir+ii;
-			fa = fap;
-			ca = -(w16*fa*fa+w15*fa+w14)/(w13*fa+w12);
-			gMSAWave[14] = -(p1+p2*ca+p3*fa);
-			gMSAWave[15] = gMSAWave[14];
-			if (fabs(gMSAWave[15])<1.0e-3) {
-				gMSAWave[15] = 0.0;
-			}
-			gMSAWave[10] = gMSAWave[16];
-		} else {
-			ca = ak2*p1+2.0*(b3*p1-b1*p3);
-			ca = -ca/(ak2*p2+2.0*(b3*p2-b2*p3));
-			fa = -(p1+p2*ca)/p3;
-			if (ix==2) {
-				gMSAWave[15] = um1*ca*ca+(um2+um3*fa)*ca+um4+um5*fa+um6*fa*fa;
-			}
-			if (ix==4) {
-				gMSAWave[15] = -(p1+p2*ca+p3*fa);
-			}
-		}
-   		gMSAWave[3] = fa;
-		gMSAWave[2] = ca;
-		gMSAWave[1] = b1+b2*ca+b3*fa;
-		gMSAWave[0] = a1+a2*ca+a3*fa;
-		gMSAWave[8] = (v1+v2*ca+v3*fa)/gMSAWave[0];
-	}
-   	g24 = e24*rgek*ex1;
-	gMSAWave[7] = (rak*ak2*ca-g24)/(ak2*g24);
-	return (ir);
-}
-
-double
-sqhcal(double qq, double gMSAWave[])
-{      	
-    double SofQ,etaz,akz,gekz,e24,x1,x2,ck,sk,ak2,qk,q2k,qk2,qk3,qqk,sink,cosk,asink,qcosk,aqk,inter; 		
-	//	WAVE gMSAWave = $"root:HayPenMSA:gMSAWave"
-
-	etaz = gMSAWave[10];
-	akz =  gMSAWave[12];
-	gekz =  gMSAWave[11];
-	e24 = 24.0*etaz;
-	x1 = exp(akz);
-	x2 = 0.0;
-	if ( gMSAWave[12]<20.0) {
-		x2 = exp(-akz);
-	}
-	ck = 0.5*(x1+x2);
-	sk = 0.5*(x1-x2);
-	ak2 = akz*akz;
-	
-	qk = qq/gMSAWave[13];
-	q2k = qk*qk;
-	if (qk<=1.0e-08) {
-		SofQ = -1.0/gMSAWave[0];
-	} else {
-	// this rescales Q.sigma = 2.Q.Radius, so is hard to predict the value to test the function
-	if (qk<=0.01) {
-		// try Taylor series expansion at small qk (RKH Feb 2016, with help from Mathematica), 
-		// transition point may need to depend on precision of cpu used and ALAS on the values of some of the parameters !
-		// note have adsorbed a factor 24 from SofQ=
-		// needs thorough test over wide range of parameter space!
-		// there seem to be some rounding issues here in single precision, must use double
-		aqk = gMSAWave[0]*(8.0+2.0*etaz) + 6*gMSAWave[1] -12.0*gMSAWave[3] 
-			-24*(gekz*(1.0+akz) -ck*akz*gMSAWave[2] +gMSAWave[3]*(ck-1.0) +(gMSAWave[2]-gMSAWave[3]*akz)*sk )/ak2
-			+q2k*( -(gMSAWave[0]*(48.0+15.0*etaz) +40.0*gMSAWave[1])/60.0 +gMSAWave[3] 
-			+(4.0/ak2)*(gekz*(9.0+7.0*akz) +ck*(9.0*gMSAWave[3] -7.0*gMSAWave[2]*akz) +sk*(9.0*gMSAWave[2] -7.0*gMSAWave[3]*akz)) );
-		SofQ = 1.0/(1.0-gMSAWave[10]*aqk);
-	} else {
-		qk2 = 1.0/q2k;
-		qk3 = qk2/qk;
-		qqk = 1.0/(qk*(q2k+ak2));
-		SINCOS(qk,sink,cosk);
-		asink = akz*sink;
-		qcosk = qk*cosk;
-		aqk = gMSAWave[0]*(sink-qcosk);
-		aqk=aqk+gMSAWave[1]*((2.0*qk2-1.0)*qcosk+2.0*sink-2.0/qk);
-		inter=24.0*qk3+4.0*(1.0-6.0*qk2)*sink;
-		aqk=(aqk+0.5*etaz*gMSAWave[0]*(inter-(1.0-12.0*qk2+24.0*qk2*qk2)*qcosk))*qk3;
-		aqk=aqk +gMSAWave[2]*(ck*asink-sk*qcosk)*qqk;
-		aqk=aqk +gMSAWave[3]*(sk*asink-qk*(ck*cosk-1.0))*qqk;
-		aqk=aqk +gMSAWave[3]*(cosk-1.0)*qk2;
-		aqk=aqk -gekz*(asink+qcosk)*qqk;
-		SofQ = 1.0/(1.0  -e24*aqk);
-	} }
-	return (SofQ);
-}
Index: sasmodels/models/lamellar_hg_stack_caille.c
===================================================================
--- sasmodels/models/lamellar_hg_stack_caille.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
+++ sasmodels/models/lamellar_hg_stack_caille.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -0,0 +1,75 @@
+/*	LamellarCailleHG kernel - allows for name changes of passed parameters ...
+    Maths identical to LamellarCaille apart from the line for P(Q)
+*/
+
+double Iq(double qval,
+      double tail_length,
+      double head_length,
+      double Nlayers, 
+      double dd,
+      double Cp,
+      double tail_sld,
+      double head_sld,
+      double solvent_sld);
+
+double Iq(double qval,
+      double tail_length,
+      double head_length,
+      double Nlayers, 
+      double dd,
+      double Cp,
+      double tail_sld,
+      double head_sld,
+      double solvent_sld)
+{
+  double NN;   //local variables of coefficient wave
+  double inten,Pq,Sq,alpha,temp,t2;
+  //double dQ, dQDefault, t1, t3;
+  int ii,NNint;
+  // from wikipedia 0.577215664901532860606512090082402431042159335
+  const double Euler = 0.577215664901533;   // Euler's constant, increased sig figs for new models Feb 2015
+  //dQDefault = 0.0;    //[=] 1/A, q-resolution, default value
+  //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015
+
+  NN = trunc(Nlayers);    //be sure that NN is an integer
+  
+  Pq = (head_sld-solvent_sld)*(sin(qval*(head_length+tail_length))-sin(qval*tail_length)) +
+              (tail_sld-solvent_sld)*sin(qval*tail_length);
+  Pq *= Pq;
+  Pq *= 4.0/(qval*qval);
+
+  NNint = (int)NN;    //cast to an integer for the loop
+  ii=0;
+  Sq = 0.0;
+  for(ii=1;ii<=(NNint-1);ii+=1) {
+
+    //fii = (double)ii;   //do I really need to do this? - unused variable, removed 18Feb2015
+
+    temp = 0.0;
+    alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler);
+    //t1 = 2.0*dQ*dQ*dd*dd*alpha;
+    t2 = 2.0*qval*qval*dd*dd*alpha;
+    //t3 = dQ*dQ*dd*dd*ii*ii;
+
+    temp = 1.0-ii/NN;
+    //temp *= cos(dd*qval*ii/(1.0+t1));
+    temp *= cos(dd*qval*ii);
+    //if (temp < 0) printf("q=%g: ii=%d, cos(dd*q*ii)=cos(%g) < 0\n",qval,ii,dd*qval*ii);
+    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
+    temp *= exp(-t2/2.0 );
+    //temp /= sqrt(1.0+t1);
+
+    Sq += temp;
+  }
+
+  Sq *= 2.0;
+  Sq += 1.0;
+
+  //if (Sq < 0) printf("q=%g: S(q) =%g\n", qval, Sq);
+
+  inten = 2.0*M_PI*Pq*Sq/(dd*qval*qval);
+
+  inten *= 1.0e-04;   // 1/A to 1/cm
+  return(inten);
+}
+
Index: sasmodels/models/lamellar_hg_stack_caille.py
===================================================================
--- sasmodels/models/lamellar_hg_stack_caille.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/models/lamellar_hg_stack_caille.py	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -112,5 +112,5 @@
     ]
 
-source = ["lamellar_hg_stack_caille_kernel.c"]
+source = ["lamellar_hg_stack_caille.c"]
 
 # No volume normalization despite having a volume parameter
Index: smodels/models/lamellar_hg_stack_caille_kernel.c
===================================================================
--- sasmodels/models/lamellar_hg_stack_caille_kernel.c	(revision 6ab4ed8677894e77dfd90dc3fdd5bdac2df2eff8)
+++ 	(revision )
@@ -1,75 +1,0 @@
-/*	LamellarCailleHG kernel - allows for name changes of passed parameters ...
-    Maths identical to LamellarCaille apart from the line for P(Q)
-*/
-
-double Iq(double qval,
-      double tail_length,
-      double head_length,
-      double Nlayers, 
-      double dd,
-      double Cp,
-      double tail_sld,
-      double head_sld,
-      double solvent_sld);
-
-double Iq(double qval,
-      double tail_length,
-      double head_length,
-      double Nlayers, 
-      double dd,
-      double Cp,
-      double tail_sld,
-      double head_sld,
-      double solvent_sld)
-{
-  double NN;   //local variables of coefficient wave
-  double inten,Pq,Sq,alpha,temp,t2;
-  //double dQ, dQDefault, t1, t3;
-  int ii,NNint;
-  // from wikipedia 0.577215664901532860606512090082402431042159335
-  const double Euler = 0.577215664901533;   // Euler's constant, increased sig figs for new models Feb 2015
-  //dQDefault = 0.0;    //[=] 1/A, q-resolution, default value
-  //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015
-
-  NN = trunc(Nlayers);    //be sure that NN is an integer
-  
-  Pq = (head_sld-solvent_sld)*(sin(qval*(head_length+tail_length))-sin(qval*tail_length)) +
-              (tail_sld-solvent_sld)*sin(qval*tail_length);
-  Pq *= Pq;
-  Pq *= 4.0/(qval*qval);
-
-  NNint = (int)NN;    //cast to an integer for the loop
-  ii=0;
-  Sq = 0.0;
-  for(ii=1;ii<=(NNint-1);ii+=1) {
-
-    //fii = (double)ii;   //do I really need to do this? - unused variable, removed 18Feb2015
-
-    temp = 0.0;
-    alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler);
-    //t1 = 2.0*dQ*dQ*dd*dd*alpha;
-    t2 = 2.0*qval*qval*dd*dd*alpha;
-    //t3 = dQ*dQ*dd*dd*ii*ii;
-
-    temp = 1.0-ii/NN;
-    //temp *= cos(dd*qval*ii/(1.0+t1));
-    temp *= cos(dd*qval*ii);
-    //if (temp < 0) printf("q=%g: ii=%d, cos(dd*q*ii)=cos(%g) < 0\n",qval,ii,dd*qval*ii);
-    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
-    temp *= exp(-t2/2.0 );
-    //temp /= sqrt(1.0+t1);
-
-    Sq += temp;
-  }
-
-  Sq *= 2.0;
-  Sq += 1.0;
-
-  //if (Sq < 0) printf("q=%g: S(q) =%g\n", qval, Sq);
-
-  inten = 2.0*M_PI*Pq*Sq/(dd*qval*qval);
-
-  inten *= 1.0e-04;   // 1/A to 1/cm
-  return(inten);
-}
-
Index: sasmodels/models/lamellar_stack_caille.c
===================================================================
--- sasmodels/models/lamellar_stack_caille.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
+++ sasmodels/models/lamellar_stack_caille.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -0,0 +1,69 @@
+/*	LamellarCaille kernel - allows for name changes of passed parameters ...
+
+*/
+
+double Iq(double qval,
+      double del,
+      double Nlayers, 
+      double dd,
+	  double Cp, 
+      double sld,
+      double solvent_sld);
+
+double Iq(double qval,
+      double del,
+      double Nlayers, 
+      double dd,
+	  double Cp, 
+      double sld,
+      double solvent_sld)
+{
+  double contr,NN;   //local variables of coefficient wave
+  double inten,Pq,Sq,alpha,temp,t2;
+  //double dQ, dQDefault, t1, t3;
+  int ii,NNint;
+  // from wikipedia 0.577215664901532860606512090082402431042159335
+  const double Euler = 0.577215664901533;   // Euler's constant, increased sig figs for new models Feb 2015
+  //dQDefault = 0.0;    //[=] 1/A, q-resolution, default value
+  //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015
+
+  NN = trunc(Nlayers);    //be sure that NN is an integer
+  
+  contr = sld - solvent_sld;
+
+  Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del));
+
+  NNint = (int)NN;    //cast to an integer for the loop
+  ii=0;
+  Sq = 0.0;
+  // the vital "=" in ii<=  added March 2015
+  for(ii=1;ii<=(NNint-1);ii+=1) {
+
+    //fii = (double)ii;   //do I really need to do this? - unused variable, removed 18Feb2015
+
+    temp = 0.0;
+    alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler);
+    //t1 = 2.0*dQ*dQ*dd*dd*alpha;
+    t2 = 2.0*qval*qval*dd*dd*alpha;
+    //t3 = dQ*dQ*dd*dd*ii*ii;
+
+    temp = 1.0-ii/NN;
+    //temp *= cos(dd*qval*ii/(1.0+t1));
+    temp *= cos(dd*qval*ii);
+    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
+    temp *= exp(-t2/2.0 );
+    //temp /= sqrt(1.0+t1);
+
+    Sq += temp;
+  }
+
+  Sq *= 2.0;
+  Sq += 1.0;
+
+  inten = 2.0*M_PI*Pq*Sq/(dd*qval*qval);
+
+  inten *= 1.0e-04;   // 1/A to 1/cm
+
+  return(inten);
+}
+
Index: sasmodels/models/lamellar_stack_caille.py
===================================================================
--- sasmodels/models/lamellar_stack_caille.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/models/lamellar_stack_caille.py	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -96,5 +96,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lamellar_stack_caille_kernel.c"]
+source = ["lamellar_stack_caille.c"]
 
 # No volume normalization despite having a volume parameter
Index: smodels/models/lamellar_stack_caille_kernel.c
===================================================================
--- sasmodels/models/lamellar_stack_caille_kernel.c	(revision 6ab4ed8677894e77dfd90dc3fdd5bdac2df2eff8)
+++ 	(revision )
@@ -1,69 +1,0 @@
-/*	LamellarCaille kernel - allows for name changes of passed parameters ...
-
-*/
-
-double Iq(double qval,
-      double del,
-      double Nlayers, 
-      double dd,
-	  double Cp, 
-      double sld,
-      double solvent_sld);
-
-double Iq(double qval,
-      double del,
-      double Nlayers, 
-      double dd,
-	  double Cp, 
-      double sld,
-      double solvent_sld)
-{
-  double contr,NN;   //local variables of coefficient wave
-  double inten,Pq,Sq,alpha,temp,t2;
-  //double dQ, dQDefault, t1, t3;
-  int ii,NNint;
-  // from wikipedia 0.577215664901532860606512090082402431042159335
-  const double Euler = 0.577215664901533;   // Euler's constant, increased sig figs for new models Feb 2015
-  //dQDefault = 0.0;    //[=] 1/A, q-resolution, default value
-  //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015
-
-  NN = trunc(Nlayers);    //be sure that NN is an integer
-  
-  contr = sld - solvent_sld;
-
-  Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del));
-
-  NNint = (int)NN;    //cast to an integer for the loop
-  ii=0;
-  Sq = 0.0;
-  // the vital "=" in ii<=  added March 2015
-  for(ii=1;ii<=(NNint-1);ii+=1) {
-
-    //fii = (double)ii;   //do I really need to do this? - unused variable, removed 18Feb2015
-
-    temp = 0.0;
-    alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler);
-    //t1 = 2.0*dQ*dQ*dd*dd*alpha;
-    t2 = 2.0*qval*qval*dd*dd*alpha;
-    //t3 = dQ*dQ*dd*dd*ii*ii;
-
-    temp = 1.0-ii/NN;
-    //temp *= cos(dd*qval*ii/(1.0+t1));
-    temp *= cos(dd*qval*ii);
-    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
-    temp *= exp(-t2/2.0 );
-    //temp /= sqrt(1.0+t1);
-
-    Sq += temp;
-  }
-
-  Sq *= 2.0;
-  Sq += 1.0;
-
-  inten = 2.0*M_PI*Pq*Sq/(dd*qval*qval);
-
-  inten *= 1.0e-04;   // 1/A to 1/cm
-
-  return(inten);
-}
-
Index: sasmodels/models/lamellar_stack_paracrystal.c
===================================================================
--- sasmodels/models/lamellar_stack_paracrystal.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
+++ sasmodels/models/lamellar_stack_paracrystal.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -0,0 +1,86 @@
+/*	Lamellar_ParaCrystal - Pedersen's model
+
+*/
+double Iq(double qval,
+      double th,
+      double Nlayers, 
+	  double davg, 
+	  double pd,
+      double sld,
+      double solvent_sld);
+double paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an);
+double paraCryst_an(double ww, double qval, double davg, long Nlayers);
+
+double Iq(double qval,
+      double th,
+      double Nlayers, 
+	  double davg, 
+	  double pd,
+      double sld,
+      double solvent_sld)
+{
+    
+	double inten,contr,xn;
+	double xi,ww,Pbil,Znq,Snq,an;
+	long n1,n2;
+	
+	contr = sld - solvent_sld;
+	//get the fractional part of Nlayers, to determine the "mixing" of N's
+	
+	n1 = (long)trunc(Nlayers);		//rounds towards zero
+	n2 = n1 + 1;
+	xn = (double)n2 - Nlayers;			//fractional contribution of n1
+	
+	ww = exp(-qval*qval*pd*pd*davg*davg/2.0);
+
+	//calculate the n1 contribution
+	an = paraCryst_an(ww,qval,davg,n1);
+	Snq = paraCryst_sn(ww,qval,davg,n1,an);
+	
+	Znq = xn*Snq;
+	
+	//calculate the n2 contribution
+	an = paraCryst_an(ww,qval,davg,n2);
+	Snq = paraCryst_sn(ww,qval,davg,n2,an);
+
+	Znq += (1.0-xn)*Snq;
+	
+	//and the independent contribution
+	Znq += (1.0-ww*ww)/(1.0+ww*ww-2.0*ww*cos(qval*davg));
+	
+	//the limit when Nlayers approaches infinity
+//	Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg))
+	
+	xi = th/2.0;		//use 1/2 the bilayer thickness
+	Pbil = (sin(qval*xi)/(qval*xi))*(sin(qval*xi)/(qval*xi));
+	
+	inten = 2.0*M_PI*contr*contr*Pbil*Znq/(qval*qval);
+	inten *= 1.0e-04;
+//printf("q=%.7e wwm1=%g ww=%.5e an=% 12.5e Snq=% 12.5e Znq=% 12.5e Pbil=% 12.5e\n",qval,wwm1,ww,an,Snq,Znq,Pbil);
+	return(inten);
+}
+
+// functions for the lamellar paracrystal model
+double
+paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an) {
+	
+	double Snq;
+
+	Snq = an/( (double)Nlayers*pow((1.0+ww*ww-2.0*ww*cos(qval*davg)),2) );
+	
+	return(Snq);
+}
+
+double
+paraCryst_an(double ww, double qval, double davg, long Nlayers) {
+	
+	double an;
+	
+	an = 4.0*ww*ww - 2.0*(ww*ww*ww+ww)*cos(qval*davg);
+	an -= 4.0*pow(ww,(Nlayers+2))*cos((double)Nlayers*qval*davg);
+	an += 2.0*pow(ww,(Nlayers+3))*cos((double)(Nlayers-1)*qval*davg);
+	an += 2.0*pow(ww,(Nlayers+1))*cos((double)(Nlayers+1)*qval*davg);
+	
+	return(an);
+}
+
Index: sasmodels/models/lamellar_stack_paracrystal.py
===================================================================
--- sasmodels/models/lamellar_stack_paracrystal.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/models/lamellar_stack_paracrystal.py	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -126,5 +126,5 @@
 
 
-source = ["lamellar_stack_paracrystal_kernel.c"]
+source = ["lamellar_stack_paracrystal.c"]
 
 form_volume = """
Index: smodels/models/lamellar_stack_paracrystal_kernel.c
===================================================================
--- sasmodels/models/lamellar_stack_paracrystal_kernel.c	(revision 6ab4ed8677894e77dfd90dc3fdd5bdac2df2eff8)
+++ 	(revision )
@@ -1,86 +1,0 @@
-/*	Lamellar_ParaCrystal - Pedersen's model
-
-*/
-double Iq(double qval,
-      double th,
-      double Nlayers, 
-	  double davg, 
-	  double pd,
-      double sld,
-      double solvent_sld);
-double paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an);
-double paraCryst_an(double ww, double qval, double davg, long Nlayers);
-
-double Iq(double qval,
-      double th,
-      double Nlayers, 
-	  double davg, 
-	  double pd,
-      double sld,
-      double solvent_sld)
-{
-    
-	double inten,contr,xn;
-	double xi,ww,Pbil,Znq,Snq,an;
-	long n1,n2;
-	
-	contr = sld - solvent_sld;
-	//get the fractional part of Nlayers, to determine the "mixing" of N's
-	
-	n1 = (long)trunc(Nlayers);		//rounds towards zero
-	n2 = n1 + 1;
-	xn = (double)n2 - Nlayers;			//fractional contribution of n1
-	
-	ww = exp(-qval*qval*pd*pd*davg*davg/2.0);
-
-	//calculate the n1 contribution
-	an = paraCryst_an(ww,qval,davg,n1);
-	Snq = paraCryst_sn(ww,qval,davg,n1,an);
-	
-	Znq = xn*Snq;
-	
-	//calculate the n2 contribution
-	an = paraCryst_an(ww,qval,davg,n2);
-	Snq = paraCryst_sn(ww,qval,davg,n2,an);
-
-	Znq += (1.0-xn)*Snq;
-	
-	//and the independent contribution
-	Znq += (1.0-ww*ww)/(1.0+ww*ww-2.0*ww*cos(qval*davg));
-	
-	//the limit when Nlayers approaches infinity
-//	Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg))
-	
-	xi = th/2.0;		//use 1/2 the bilayer thickness
-	Pbil = (sin(qval*xi)/(qval*xi))*(sin(qval*xi)/(qval*xi));
-	
-	inten = 2.0*M_PI*contr*contr*Pbil*Znq/(qval*qval);
-	inten *= 1.0e-04;
-//printf("q=%.7e wwm1=%g ww=%.5e an=% 12.5e Snq=% 12.5e Znq=% 12.5e Pbil=% 12.5e\n",qval,wwm1,ww,an,Snq,Znq,Pbil);
-	return(inten);
-}
-
-// functions for the lamellar paracrystal model
-double
-paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an) {
-	
-	double Snq;
-
-	Snq = an/( (double)Nlayers*pow((1.0+ww*ww-2.0*ww*cos(qval*davg)),2) );
-	
-	return(Snq);
-}
-
-double
-paraCryst_an(double ww, double qval, double davg, long Nlayers) {
-	
-	double an;
-	
-	an = 4.0*ww*ww - 2.0*(ww*ww*ww+ww)*cos(qval*davg);
-	an -= 4.0*pow(ww,(Nlayers+2))*cos((double)Nlayers*qval*davg);
-	an += 2.0*pow(ww,(Nlayers+3))*cos((double)(Nlayers-1)*qval*davg);
-	an += 2.0*pow(ww,(Nlayers+1))*cos((double)(Nlayers+1)*qval*davg);
-	
-	return(an);
-}
-
Index: sasmodels/models/polymer_micelle.c
===================================================================
--- sasmodels/models/polymer_micelle.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
+++ sasmodels/models/polymer_micelle.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -0,0 +1,87 @@
+double Iq(double q,
+        double ndensity,
+        double v_core,
+        double v_corona,
+        double solvent_sld,
+        double core_sld,
+        double corona_sld,
+        double radius_core,
+        double radius_gyr,
+        double d_penetration,
+        double n_aggreg);
+
+static double micelle_spherical_kernel(double q,
+        double ndensity,
+        double v_core,
+        double v_corona,
+        double solvent_sld,
+        double core_sld,
+        double corona_sld,
+        double radius_core,
+        double radius_gyr,
+        double d_penetration,
+        double n_aggreg)
+{
+    const double rho_solv = solvent_sld;     // sld of solvent [1/A^2]
+    const double rho_core = core_sld;        // sld of core [1/A^2]
+    const double rho_corona = corona_sld;    // sld of corona [1/A^2]
+
+    const double beta_core = v_core * (rho_core - rho_solv);
+    const double beta_corona = v_corona * (rho_corona - rho_solv);
+
+    // Self-correlation term of the core
+    const double bes_core = sph_j1c(q*radius_core);
+    const double term1 = n_aggreg*n_aggreg*beta_core*beta_core*bes_core*bes_core;
+
+    // Self-correlation term of the chains
+    const double qrg2 = q*radius_gyr*q*radius_gyr;
+    const double debye_chain = (qrg2 == 0.0) ? 1.0 : 2.0*(expm1(-qrg2)+qrg2)/(qrg2*qrg2);
+    const double term2 = n_aggreg * beta_corona * beta_corona * debye_chain;
+
+    // Interference cross-term between core and chains
+    const double chain_ampl = (qrg2 == 0.0) ? 1.0 : -expm1(-qrg2)/qrg2;
+    const double bes_corona = sinc(q*(radius_core + d_penetration * radius_gyr));
+    const double term3 = 2 * n_aggreg * n_aggreg * beta_core * beta_corona *
+                 bes_core * chain_ampl * bes_corona;
+
+    // Interference cross-term between chains
+    const double term4 = n_aggreg * (n_aggreg - 1.0) * beta_corona * beta_corona *
+                 chain_ampl * chain_ampl * bes_corona * bes_corona;
+
+    // I(q)_micelle : Sum of 4 terms computed above
+    double i_micelle = term1 + term2 + term3 + term4;
+
+    // rescale from [A^2] to [cm^2]
+    i_micelle *= 1.0e-13;
+
+    // "normalize" by number density --> intensity in [cm-1]
+    i_micelle *= ndensity;
+
+    return(i_micelle);
+
+}
+
+double Iq(double q,
+        double ndensity,
+        double v_core,
+        double v_corona,
+        double solvent_sld,
+        double core_sld,
+        double corona_sld,
+        double radius_core,
+        double radius_gyr,
+        double d_penetration,
+        double n_aggreg)
+{
+    return micelle_spherical_kernel(q,
+            ndensity,
+            v_core,
+            v_corona,
+            solvent_sld,
+            core_sld,
+            corona_sld,
+            radius_core,
+            radius_gyr,
+            d_penetration,
+            n_aggreg);
+}
Index: sasmodels/models/polymer_micelle.py
===================================================================
--- sasmodels/models/polymer_micelle.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/models/polymer_micelle.py	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -54,5 +54,5 @@
 single = False
 
-source = ["lib/sph_j1c.c", "polymer_micelle_kernel.c"]
+source = ["lib/sph_j1c.c", "polymer_micelle.c"]
 
 demo = dict(scale=1, background=0,
Index: smodels/models/polymer_micelle_kernel.c
===================================================================
--- sasmodels/models/polymer_micelle_kernel.c	(revision c91537379690f2cf6ec329b94f79d3dabf73fb8c)
+++ 	(revision )
@@ -1,87 +1,0 @@
-double Iq(double q,
-        double ndensity,
-        double v_core,
-        double v_corona,
-        double solvent_sld,
-        double core_sld,
-        double corona_sld,
-        double radius_core,
-        double radius_gyr,
-        double d_penetration,
-        double n_aggreg);
-
-static double micelle_spherical_kernel(double q,
-        double ndensity,
-        double v_core,
-        double v_corona,
-        double solvent_sld,
-        double core_sld,
-        double corona_sld,
-        double radius_core,
-        double radius_gyr,
-        double d_penetration,
-        double n_aggreg)
-{
-    const double rho_solv = solvent_sld;     // sld of solvent [1/A^2]
-    const double rho_core = core_sld;        // sld of core [1/A^2]
-    const double rho_corona = corona_sld;    // sld of corona [1/A^2]
-
-    const double beta_core = v_core * (rho_core - rho_solv);
-    const double beta_corona = v_corona * (rho_corona - rho_solv);
-
-    // Self-correlation term of the core
-    const double bes_core = sph_j1c(q*radius_core);
-    const double term1 = n_aggreg*n_aggreg*beta_core*beta_core*bes_core*bes_core;
-
-    // Self-correlation term of the chains
-    const double qrg2 = q*radius_gyr*q*radius_gyr;
-    const double debye_chain = (qrg2 == 0.0) ? 1.0 : 2.0*(expm1(-qrg2)+qrg2)/(qrg2*qrg2);
-    const double term2 = n_aggreg * beta_corona * beta_corona * debye_chain;
-
-    // Interference cross-term between core and chains
-    const double chain_ampl = (qrg2 == 0.0) ? 1.0 : -expm1(-qrg2)/qrg2;
-    const double bes_corona = sinc(q*(radius_core + d_penetration * radius_gyr));
-    const double term3 = 2 * n_aggreg * n_aggreg * beta_core * beta_corona *
-                 bes_core * chain_ampl * bes_corona;
-
-    // Interference cross-term between chains
-    const double term4 = n_aggreg * (n_aggreg - 1.0) * beta_corona * beta_corona *
-                 chain_ampl * chain_ampl * bes_corona * bes_corona;
-
-    // I(q)_micelle : Sum of 4 terms computed above
-    double i_micelle = term1 + term2 + term3 + term4;
-
-    // rescale from [A^2] to [cm^2]
-    i_micelle *= 1.0e-13;
-
-    // "normalize" by number density --> intensity in [cm-1]
-    i_micelle *= ndensity;
-
-    return(i_micelle);
-
-}
-
-double Iq(double q,
-        double ndensity,
-        double v_core,
-        double v_corona,
-        double solvent_sld,
-        double core_sld,
-        double corona_sld,
-        double radius_core,
-        double radius_gyr,
-        double d_penetration,
-        double n_aggreg)
-{
-    return micelle_spherical_kernel(q,
-            ndensity,
-            v_core,
-            v_corona,
-            solvent_sld,
-            core_sld,
-            corona_sld,
-            radius_core,
-            radius_gyr,
-            d_penetration,
-            n_aggreg);
-}
Index: sasmodels/models/sc_paracrystal.c
===================================================================
--- sasmodels/models/sc_paracrystal.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
+++ sasmodels/models/sc_paracrystal.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -0,0 +1,228 @@
+double form_volume(double radius);
+
+double Iq(double q,
+          double dnn,
+          double d_factor,
+          double radius,
+          double sphere_sld,
+          double solvent_sld);
+
+double Iqxy(double qx, double qy,
+            double dnn,
+            double d_factor,
+            double radius,
+            double sphere_sld,
+            double solvent_sld,
+            double theta,
+            double phi,
+            double psi);
+
+double form_volume(double radius)
+{
+    return sphere_volume(radius);
+}
+
+static double
+sc_eval(double theta, double phi, double temp3, double temp4, double temp5)
+{
+    double cnt, snt;
+    SINCOS(theta, cnt, snt);
+
+    double cnp, snp;
+    SINCOS(phi, cnp, snp);
+
+	double temp6 = snt;
+	double temp7 = -1.0*temp3*snt*cnp;
+	double temp8 = temp3*snt*snp;
+	double temp9 = temp3*cnt;
+	double result = temp6/((1.0-temp4*cos((temp7))+temp5)*
+	                       (1.0-temp4*cos((temp8))+temp5)*
+	                       (1.0-temp4*cos((temp9))+temp5));
+	return (result);
+}
+
+static double
+sc_integrand(double dnn, double d_factor, double qq, double xx, double yy)
+{
+    //Function to calculate integrand values for simple cubic structure
+
+	double da = d_factor*dnn;
+	double temp1 = qq*qq*da*da;
+	double temp2 = pow( 1.0-exp(-1.0*temp1) ,3);
+	double temp3 = qq*dnn;
+	double temp4 = 2.0*exp(-0.5*temp1);
+	double temp5 = exp(-1.0*temp1);
+
+	double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5);
+	integrand *= 2.0/M_PI;
+
+	return(integrand);
+}
+
+static
+double sc_crystal_kernel(double q,
+          double dnn,
+          double d_factor,
+          double radius,
+          double sphere_sld,
+          double solvent_sld)
+{
+	const double va = 0.0;
+	const double vb = M_PI/2.0; //orientation average, outer integral
+
+    double summ=0.0;
+    double answer=0.0;
+
+	for(int i=0;i<150;i++) {
+		//setup inner integral over the ellipsoidal cross-section
+		double summj=0.0;
+		double zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;
+		for(int j=0;j<150;j++) {
+			//150 gauss points for the inner integral
+			double zij = ( Gauss150Z[j]*(vb-va) + va + vb )/2.0;
+			double tmp = Gauss150Wt[j] * sc_integrand(dnn,d_factor,q,zi,zij);
+			summj += tmp;
+		}
+		//now calculate the value of the inner integral
+		answer = (vb-va)/2.0*summj;
+
+		//now calculate outer integral
+		double tmp = Gauss150Wt[i] * answer;
+		summ += tmp;
+	}		//final scaling is done at the end of the function, after the NT_FP64 case
+
+	answer = (vb-va)/2.0*summ;
+
+	//Volume fraction calculated from lattice symmetry and sphere radius
+	// NB: 4/3 pi r^3 / dnn^3 = 4/3 pi(r/dnn)^3
+	const double latticeScale = sphere_volume(radius/dnn);
+
+	answer *= sphere_form(q, radius, sphere_sld, solvent_sld)*latticeScale;
+
+	return answer;
+}
+
+static
+double sc_crystal_kernel_2d(double q, double q_x, double q_y,
+          double dnn,
+          double d_factor,
+          double radius,
+          double sphere_sld,
+          double solvent_sld,
+          double theta,
+          double phi,
+          double psi)
+{
+    //convert angle degree to radian
+    theta = theta * M_PI_180;
+    phi = phi * M_PI_180;
+    psi = psi * M_PI_180;
+
+    const double qda_2 = pow(q*d_factor*dnn,2.0);
+
+    double snt, cnt;
+    SINCOS(theta, snt, cnt);
+
+    double snp, cnp;
+    SINCOS(phi, snp, cnp);
+
+    double sns, cns;
+    SINCOS(psi, sns, cns);
+
+    /// Angles here are respect to detector coordinate instead of against
+    //  q coordinate(PRB 36, 3, 1754)
+    // a3 axis orientation
+
+    const double a3_x = cnt * cnp;
+    const double a3_y = snt;
+
+    // Compute the angle btw vector q and the a3 axis
+    double cos_val_a3 = a3_x*q_x + a3_y*q_y;
+
+    // a1 axis orientation
+    const double a1_x = -cnp*sns * snt+snp*cns;
+    const double a1_y = sns*cnt;
+
+    double cos_val_a1 = a1_x*q_x + a1_y*q_y;
+
+    // a2 axis orientation
+    const double a2_x = -snt*cns*cnp-sns*snp;
+    const double a2_y = cnt*cns;
+
+    // a2 axis
+    const double cos_val_a2 =  a2_x*q_x + a2_y*q_y;
+
+    // The following test should always pass
+    if (fabs(cos_val_a3)>1.0) {
+        //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+        cos_val_a3 = 1.0;
+    }
+    if (fabs(cos_val_a1)>1.0) {
+        //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+        cos_val_a1 = 1.0;
+    }
+    if (fabs(cos_val_a2)>1.0) {
+        //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+        cos_val_a3 = 1.0;
+    }
+
+    const double a3_dot_q = dnn*q*cos_val_a3;
+    const double a1_dot_q = dnn*q*cos_val_a1;
+    const double a2_dot_q = dnn*q*cos_val_a2;
+
+    // Call Zq=Z1*Z2*Z3
+    double Zq = (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a1_dot_q)+exp(-qda_2));
+    Zq *= (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a2_dot_q)+exp(-qda_2));
+    Zq *= (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a3_dot_q)+exp(-qda_2));
+
+    // Use SphereForm directly from libigor
+    double answer = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq;
+
+    //consider scales
+    const double latticeScale = sphere_volume(radius/dnn);
+    answer *= latticeScale;
+
+    return answer;
+}
+
+double Iq(double q,
+          double dnn,
+          double d_factor,
+          double radius,
+          double sphere_sld,
+          double solvent_sld)
+{
+    return sc_crystal_kernel(q,
+              dnn,
+              d_factor,
+              radius,
+              sphere_sld,
+              solvent_sld);
+}
+
+// Iqxy is never called since no orientation or magnetic parameters.
+double Iqxy(double qx, double qy,
+            double dnn,
+            double d_factor,
+            double radius,
+            double sphere_sld,
+            double solvent_sld,
+            double theta,
+            double phi,
+            double psi)
+{
+    double q = sqrt(qx*qx + qy*qy);
+
+
+    return sc_crystal_kernel_2d(q, qx/q, qy/q,
+                  dnn,
+                  d_factor,
+                  radius,
+                  sphere_sld,
+                  solvent_sld,
+                  theta,
+                  phi,
+                  psi);
+
+}
+
Index: sasmodels/models/sc_paracrystal.py
===================================================================
--- sasmodels/models/sc_paracrystal.py	(revision 42356c8890c04b71169941b12a36fc0e07e2034b)
+++ sasmodels/models/sc_paracrystal.py	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
@@ -133,5 +133,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/sphere_form.c", "lib/gauss150.c", "sc_paracrystal_kernel.c"]
+source = ["lib/sph_j1c.c", "lib/sphere_form.c", "lib/gauss150.c", "sc_paracrystal.c"]
 
 demo = dict(scale=1, background=0,
Index: smodels/models/sc_paracrystal_kernel.c
===================================================================
--- sasmodels/models/sc_paracrystal_kernel.c	(revision 02a09200b0a227653f44e340545f09531f7e69f2)
+++ 	(revision )
@@ -1,228 +1,0 @@
-double form_volume(double radius);
-
-double Iq(double q,
-          double dnn,
-          double d_factor,
-          double radius,
-          double sphere_sld,
-          double solvent_sld);
-
-double Iqxy(double qx, double qy,
-            double dnn,
-            double d_factor,
-            double radius,
-            double sphere_sld,
-            double solvent_sld,
-            double theta,
-            double phi,
-            double psi);
-
-double form_volume(double radius)
-{
-    return sphere_volume(radius);
-}
-
-static double
-sc_eval(double theta, double phi, double temp3, double temp4, double temp5)
-{
-    double cnt, snt;
-    SINCOS(theta, cnt, snt);
-
-    double cnp, snp;
-    SINCOS(phi, cnp, snp);
-
-	double temp6 = snt;
-	double temp7 = -1.0*temp3*snt*cnp;
-	double temp8 = temp3*snt*snp;
-	double temp9 = temp3*cnt;
-	double result = temp6/((1.0-temp4*cos((temp7))+temp5)*
-	                       (1.0-temp4*cos((temp8))+temp5)*
-	                       (1.0-temp4*cos((temp9))+temp5));
-	return (result);
-}
-
-static double
-sc_integrand(double dnn, double d_factor, double qq, double xx, double yy)
-{
-    //Function to calculate integrand values for simple cubic structure
-
-	double da = d_factor*dnn;
-	double temp1 = qq*qq*da*da;
-	double temp2 = pow( 1.0-exp(-1.0*temp1) ,3);
-	double temp3 = qq*dnn;
-	double temp4 = 2.0*exp(-0.5*temp1);
-	double temp5 = exp(-1.0*temp1);
-
-	double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5);
-	integrand *= 2.0/M_PI;
-
-	return(integrand);
-}
-
-static
-double sc_crystal_kernel(double q,
-          double dnn,
-          double d_factor,
-          double radius,
-          double sphere_sld,
-          double solvent_sld)
-{
-	const double va = 0.0;
-	const double vb = M_PI/2.0; //orientation average, outer integral
-
-    double summ=0.0;
-    double answer=0.0;
-
-	for(int i=0;i<150;i++) {
-		//setup inner integral over the ellipsoidal cross-section
-		double summj=0.0;
-		double zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;
-		for(int j=0;j<150;j++) {
-			//150 gauss points for the inner integral
-			double zij = ( Gauss150Z[j]*(vb-va) + va + vb )/2.0;
-			double tmp = Gauss150Wt[j] * sc_integrand(dnn,d_factor,q,zi,zij);
-			summj += tmp;
-		}
-		//now calculate the value of the inner integral
-		answer = (vb-va)/2.0*summj;
-
-		//now calculate outer integral
-		double tmp = Gauss150Wt[i] * answer;
-		summ += tmp;
-	}		//final scaling is done at the end of the function, after the NT_FP64 case
-
-	answer = (vb-va)/2.0*summ;
-
-	//Volume fraction calculated from lattice symmetry and sphere radius
-	// NB: 4/3 pi r^3 / dnn^3 = 4/3 pi(r/dnn)^3
-	const double latticeScale = sphere_volume(radius/dnn);
-
-	answer *= sphere_form(q, radius, sphere_sld, solvent_sld)*latticeScale;
-
-	return answer;
-}
-
-static
-double sc_crystal_kernel_2d(double q, double q_x, double q_y,
-          double dnn,
-          double d_factor,
-          double radius,
-          double sphere_sld,
-          double solvent_sld,
-          double theta,
-          double phi,
-          double psi)
-{
-    //convert angle degree to radian
-    theta = theta * M_PI_180;
-    phi = phi * M_PI_180;
-    psi = psi * M_PI_180;
-
-    const double qda_2 = pow(q*d_factor*dnn,2.0);
-
-    double snt, cnt;
-    SINCOS(theta, snt, cnt);
-
-    double snp, cnp;
-    SINCOS(phi, snp, cnp);
-
-    double sns, cns;
-    SINCOS(psi, sns, cns);
-
-    /// Angles here are respect to detector coordinate instead of against
-    //  q coordinate(PRB 36, 3, 1754)
-    // a3 axis orientation
-
-    const double a3_x = cnt * cnp;
-    const double a3_y = snt;
-
-    // Compute the angle btw vector q and the a3 axis
-    double cos_val_a3 = a3_x*q_x + a3_y*q_y;
-
-    // a1 axis orientation
-    const double a1_x = -cnp*sns * snt+snp*cns;
-    const double a1_y = sns*cnt;
-
-    double cos_val_a1 = a1_x*q_x + a1_y*q_y;
-
-    // a2 axis orientation
-    const double a2_x = -snt*cns*cnp-sns*snp;
-    const double a2_y = cnt*cns;
-
-    // a2 axis
-    const double cos_val_a2 =  a2_x*q_x + a2_y*q_y;
-
-    // The following test should always pass
-    if (fabs(cos_val_a3)>1.0) {
-        //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
-        cos_val_a3 = 1.0;
-    }
-    if (fabs(cos_val_a1)>1.0) {
-        //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
-        cos_val_a1 = 1.0;
-    }
-    if (fabs(cos_val_a2)>1.0) {
-        //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
-        cos_val_a3 = 1.0;
-    }
-
-    const double a3_dot_q = dnn*q*cos_val_a3;
-    const double a1_dot_q = dnn*q*cos_val_a1;
-    const double a2_dot_q = dnn*q*cos_val_a2;
-
-    // Call Zq=Z1*Z2*Z3
-    double Zq = (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a1_dot_q)+exp(-qda_2));
-    Zq *= (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a2_dot_q)+exp(-qda_2));
-    Zq *= (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a3_dot_q)+exp(-qda_2));
-
-    // Use SphereForm directly from libigor
-    double answer = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq;
-
-    //consider scales
-    const double latticeScale = sphere_volume(radius/dnn);
-    answer *= latticeScale;
-
-    return answer;
-}
-
-double Iq(double q,
-          double dnn,
-          double d_factor,
-          double radius,
-          double sphere_sld,
-          double solvent_sld)
-{
-    return sc_crystal_kernel(q,
-              dnn,
-              d_factor,
-              radius,
-              sphere_sld,
-              solvent_sld);
-}
-
-// Iqxy is never called since no orientation or magnetic parameters.
-double Iqxy(double qx, double qy,
-            double dnn,
-            double d_factor,
-            double radius,
-            double sphere_sld,
-            double solvent_sld,
-            double theta,
-            double phi,
-            double psi)
-{
-    double q = sqrt(qx*qx + qy*qy);
-
-
-    return sc_crystal_kernel_2d(q, qx/q, qy/q,
-                  dnn,
-                  d_factor,
-                  radius,
-                  sphere_sld,
-                  solvent_sld,
-                  theta,
-                  phi,
-                  psi);
-
-}
-
