Index: sansmodels/src/c_models/binaryHS.cpp
===================================================================
--- sansmodels/src/c_models/binaryHS.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/binaryHS.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -21,26 +21,25 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "binaryHS.h"
 
 extern "C" {
-	#include "libSphere.h"
-	#include "binaryHS.h"
+#include "libSphere.h"
 }
 
 BinaryHSModel :: BinaryHSModel() {
 
-	l_radius     = Parameter(160.0, true);
-	l_radius.set_min(0.0);
-	s_radius    = Parameter(25.0, true);
-	s_radius.set_min(0.0);
-	vol_frac_ls  = Parameter(0.2);
-	vol_frac_ss  = Parameter(0.1);
-	ls_sld      = Parameter(3.5e-6);
-	ss_sld     = Parameter(5e-7);
-	solvent_sld   = Parameter(6.36e-6);
-	background = Parameter(0.0);
+  l_radius     = Parameter(160.0, true);
+  l_radius.set_min(0.0);
+  s_radius    = Parameter(25.0, true);
+  s_radius.set_min(0.0);
+  vol_frac_ls  = Parameter(0.2);
+  vol_frac_ss  = Parameter(0.1);
+  ls_sld      = Parameter(3.5e-6);
+  ss_sld     = Parameter(5e-7);
+  solvent_sld   = Parameter(6.36e-6);
+  background = Parameter(0.0);
 }
 
@@ -52,44 +51,44 @@
  */
 double BinaryHSModel :: operator()(double q) {
-	double dp[8];
+  double dp[8];
 
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = l_radius();
-	dp[1] = s_radius();
-	dp[2] = vol_frac_ls();
-	dp[3] = vol_frac_ss();
-	dp[4] = ls_sld();
-	dp[5] = ss_sld();
-	dp[6] = solvent_sld();
-	dp[7] = 0.0;
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = l_radius();
+  dp[1] = s_radius();
+  dp[2] = vol_frac_ls();
+  dp[3] = vol_frac_ss();
+  dp[4] = ls_sld();
+  dp[5] = ss_sld();
+  dp[6] = solvent_sld();
+  dp[7] = 0.0;
 
 
-	// Get the dispersion points for the large radius
-	vector<WeightPoint> weights_l_radius;
-	l_radius.get_weights(weights_l_radius);
+  // Get the dispersion points for the large radius
+  vector<WeightPoint> weights_l_radius;
+  l_radius.get_weights(weights_l_radius);
 
-	// Get the dispersion points for the small radius
-	vector<WeightPoint> weights_s_radius;
-	s_radius.get_weights(weights_s_radius);
+  // Get the dispersion points for the small radius
+  vector<WeightPoint> weights_s_radius;
+  s_radius.get_weights(weights_s_radius);
 
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
 
-	// Loop over larger radius weight points
-	for(int i=0; i< (int)weights_l_radius.size(); i++) {
-		dp[0] = weights_l_radius[i].value;
+  // Loop over larger radius weight points
+  for(int i=0; i< (int)weights_l_radius.size(); i++) {
+    dp[0] = weights_l_radius[i].value;
 
-		// Loop over small radius weight points
-		for(int j=0; j< (int)weights_s_radius.size(); j++) {
-			dp[1] = weights_s_radius[j].value;
+    // Loop over small radius weight points
+    for(int j=0; j< (int)weights_s_radius.size(); j++) {
+      dp[1] = weights_s_radius[j].value;
 
 
-			sum += weights_l_radius[i].weight *weights_s_radius[j].weight * BinaryHS(dp, q);
-			norm += weights_l_radius[i].weight *weights_s_radius[j].weight;
-		}
-	}
-	return sum/norm + background();
+      sum += weights_l_radius[i].weight *weights_s_radius[j].weight * BinaryHS(dp, q);
+      norm += weights_l_radius[i].weight *weights_s_radius[j].weight;
+    }
+  }
+  return sum/norm + background();
 }
 
@@ -101,6 +100,6 @@
  */
 double BinaryHSModel :: operator()(double qx, double qy) {
-	double q = sqrt(qx*qx + qy*qy);
-	return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -113,5 +112,5 @@
  */
 double BinaryHSModel :: evaluate_rphi(double q, double phi) {
-	return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 
@@ -121,6 +120,6 @@
  */
 double BinaryHSModel :: calculate_ER() {
-	//NOT implemented yet!!!
-	return 0.0;
+  //NOT implemented yet!!!
+  return 0.0;
 }
 
Index: sansmodels/src/c_models/binaryHS_PSF11.cpp
===================================================================
--- sansmodels/src/c_models/binaryHS_PSF11.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/binaryHS_PSF11.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -21,5 +21,4 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
@@ -29,4 +28,26 @@
 	#include "libSphere.h"
 }
+
+class BinaryHSPSF11Model{
+public:
+  // Model parameters
+  Parameter l_radius;
+  Parameter s_radius;
+  Parameter vol_frac_ls;
+  Parameter vol_frac_ss;
+  Parameter ls_sld;
+  Parameter ss_sld;
+  Parameter solvent_sld;
+  Parameter background;
+
+  //Constructor
+  BinaryHSPSF11Model();
+
+  //Operators to get I(Q)
+  double operator()(double q);
+  double operator()(double qx , double qy);
+  double calculate_ER();
+  double evaluate_rphi(double q, double phi);
+};
 
 BinaryHSPSF11Model :: BinaryHSPSF11Model() {
Index: sansmodels/src/c_models/c_models.cpp
===================================================================
--- sansmodels/src/c_models/c_models.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/c_models.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -68,8 +68,8 @@
 	void addDisperser(PyObject *module);
 	void addCGaussian(PyObject *module);
-	void addCLorentzian(PyObject *module);
 	void addCLogNormal(PyObject *module);
 	void addCSchulz(PyObject *module);
 }
+void addCLorentzian(PyObject *module);
 
 /**
Index: sansmodels/src/c_models/coreshellcylinder.cpp
===================================================================
--- sansmodels/src/c_models/coreshellcylinder.cpp	(revision 011e0e4569d018a5f39b70032285911f67750383)
+++ sansmodels/src/c_models/coreshellcylinder.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -27,6 +27,6 @@
 
 extern "C" {
-	#include "libCylinder.h"
-	#include "libStructureFactor.h"
+#include "libCylinder.h"
+#include "libStructureFactor.h"
 }
 
@@ -63,21 +63,21 @@
   double phi = pars->axis_phi * pi/180.0;
 
-    // Cylinder orientation
-    cyl_x = sin(theta) * cos(phi);
-    cyl_y = sin(theta) * sin(phi);
-    cyl_z = cos(theta);
-
-    // q vector
-    q_z = 0;
-
-    // Compute the angle btw vector q and the
-    // axis of the cylinder
-    cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
-
-    // The following test should always pass
-    if (fabs(cos_val)>1.0) {
-      printf("core_shell_cylinder_analytical_2D: Unexpected error: cos(alpha)=%g\n", cos_val);
-      return 0;
-    }
+  // Cylinder orientation
+  cyl_x = sin(theta) * cos(phi);
+  cyl_y = sin(theta) * sin(phi);
+  cyl_z = cos(theta);
+
+  // q vector
+  q_z = 0;
+
+  // Compute the angle btw vector q and the
+  // axis of the cylinder
+  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
+
+  // The following test should always pass
+  if (fabs(cos_val)>1.0) {
+    printf("core_shell_cylinder_analytical_2D: Unexpected error: cos(alpha)=%g\n", cos_val);
+    return 0;
+  }
 
   alpha = acos( cos_val );
@@ -85,11 +85,11 @@
   // Call the IGOR library function to get the kernel
   answer = CoreShellCylKernel(q, pars->radius, pars->thickness,
-                pars->core_sld,pars->shell_sld,
-                pars->solvent_sld, pars->length/2.0, alpha) / fabs(sin(alpha));
+      pars->core_sld,pars->shell_sld,
+      pars->solvent_sld, pars->length/2.0, alpha) / fabs(sin(alpha));
 
   //normalize by cylinder volume
   vol=pi*(pars->radius+pars->thickness)
-      *(pars->radius+pars->thickness)
-      *(pars->length+2.0*pars->thickness);
+          *(pars->radius+pars->thickness)
+          *(pars->length+2.0*pars->thickness);
   answer /= vol;
 
@@ -115,22 +115,22 @@
   double q;
   q = sqrt(qx*qx+qy*qy);
-    return core_shell_cylinder_analytical_2D_scaled(pars, q, qx/q, qy/q);
+  return core_shell_cylinder_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
 
 CoreShellCylinderModel :: CoreShellCylinderModel() {
-	scale      = Parameter(1.0);
-	radius     = Parameter(20.0, true);
-	radius.set_min(0.0);
-	thickness  = Parameter(10.0, true);
-	thickness.set_min(0.0);
-	length     = Parameter(400.0, true);
-	length.set_min(0.0);
-	core_sld   = Parameter(1.e-6);
-	shell_sld  = Parameter(4.e-6);
-	solvent_sld= Parameter(1.e-6);
-	background = Parameter(0.0);
-	axis_theta = Parameter(90.0, true);
-	axis_phi   = Parameter(0.0, true);
+  scale      = Parameter(1.0);
+  radius     = Parameter(20.0, true);
+  radius.set_min(0.0);
+  thickness  = Parameter(10.0, true);
+  thickness.set_min(0.0);
+  length     = Parameter(400.0, true);
+  length.set_min(0.0);
+  core_sld   = Parameter(1.e-6);
+  shell_sld  = Parameter(4.e-6);
+  solvent_sld= Parameter(1.e-6);
+  background = Parameter(0.0);
+  axis_theta = Parameter(90.0, true);
+  axis_phi   = Parameter(0.0, true);
 }
 
@@ -142,68 +142,68 @@
  */
 double CoreShellCylinderModel :: operator()(double q) {
-	double dp[8];
-
-	dp[0] = scale();
-	dp[1] = radius();
-	dp[2] = thickness();
-	dp[3] = length();
-	dp[4] = core_sld();
-	dp[5] = shell_sld();
-	dp[6] = solvent_sld();
-	dp[7] = 0.0;
-
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_rad;
-	radius.get_weights(weights_rad);
-
-	// Get the dispersion points for the thickness
-	vector<WeightPoint> weights_thick;
-	thickness.get_weights(weights_thick);
-
-	// Get the dispersion points for the length
-	vector<WeightPoint> weights_len;
-	length.get_weights(weights_len);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
-
-	// Loop over radius weight points
-	for(size_t i=0; i<weights_rad.size(); i++) {
-		dp[1] = weights_rad[i].value;
-
-		// Loop over length weight points
-		for(size_t j=0; j<weights_len.size(); j++) {
-			dp[3] = weights_len[j].value;
-
-			// Loop over thickness weight points
-			for(size_t k=0; k<weights_thick.size(); k++) {
-				dp[2] = weights_thick[k].value;
-				//Un-normalize by volume
-				sum += weights_rad[i].weight
-					* weights_len[j].weight
-					* weights_thick[k].weight
-					* CoreShellCylinder(dp, q)
-					* pow(weights_rad[i].value+weights_thick[k].value,2)
-					*(weights_len[j].value+2.0*weights_thick[k].value);
-				//Find average volume
-				vol += weights_rad[i].weight
-					* weights_len[j].weight
-					* weights_thick[k].weight
-					* pow(weights_rad[i].value+weights_thick[k].value,2)
-					*(weights_len[j].value+2.0*weights_thick[k].value);
-				norm += weights_rad[i].weight
-				* weights_len[j].weight
-				* weights_thick[k].weight;
-			}
-		}
-	}
-
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
-
-	return sum/norm + background();
+  double dp[8];
+
+  dp[0] = scale();
+  dp[1] = radius();
+  dp[2] = thickness();
+  dp[3] = length();
+  dp[4] = core_sld();
+  dp[5] = shell_sld();
+  dp[6] = solvent_sld();
+  dp[7] = 0.0;
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad;
+  radius.get_weights(weights_rad);
+
+  // Get the dispersion points for the thickness
+  vector<WeightPoint> weights_thick;
+  thickness.get_weights(weights_thick);
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_len;
+  length.get_weights(weights_len);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp[1] = weights_rad[i].value;
+
+    // Loop over length weight points
+    for(size_t j=0; j<weights_len.size(); j++) {
+      dp[3] = weights_len[j].value;
+
+      // Loop over thickness weight points
+      for(size_t k=0; k<weights_thick.size(); k++) {
+        dp[2] = weights_thick[k].value;
+        //Un-normalize by volume
+        sum += weights_rad[i].weight
+            * weights_len[j].weight
+            * weights_thick[k].weight
+            * CoreShellCylinder(dp, q)
+        * pow(weights_rad[i].value+weights_thick[k].value,2)
+        *(weights_len[j].value+2.0*weights_thick[k].value);
+        //Find average volume
+        vol += weights_rad[i].weight
+            * weights_len[j].weight
+            * weights_thick[k].weight
+            * pow(weights_rad[i].value+weights_thick[k].value,2)
+        *(weights_len[j].value+2.0*weights_thick[k].value);
+        norm += weights_rad[i].weight
+            * weights_len[j].weight
+            * weights_thick[k].weight;
+      }
+    }
+  }
+
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -215,110 +215,110 @@
  */
 double CoreShellCylinderModel :: operator()(double qx, double qy) {
-	CoreShellCylinderParameters dp;
-	// Fill parameter array
-	dp.scale      = scale();
-	dp.radius     = radius();
-	dp.thickness  = thickness();
-	dp.length     = length();
-	dp.core_sld   = core_sld();
-	dp.shell_sld  = shell_sld();
-	dp.solvent_sld= solvent_sld();
-	dp.background = 0.0;
-	dp.axis_theta = axis_theta();
-	dp.axis_phi   = axis_phi();
-
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_rad;
-	radius.get_weights(weights_rad);
-
-	// Get the dispersion points for the thickness
-	vector<WeightPoint> weights_thick;
-	thickness.get_weights(weights_thick);
-
-	// Get the dispersion points for the length
-	vector<WeightPoint> weights_len;
-	length.get_weights(weights_len);
-
-	// Get angular averaging for theta
-	vector<WeightPoint> weights_theta;
-	axis_theta.get_weights(weights_theta);
-
-	// Get angular averaging for phi
-	vector<WeightPoint> weights_phi;
-	axis_phi.get_weights(weights_phi);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double norm_vol = 0.0;
-	double vol = 0.0;
-	double pi = 4.0*atan(1.0);
-	// Loop over radius weight points
-	for(size_t i=0; i<weights_rad.size(); i++) {
-		dp.radius = weights_rad[i].value;
-
-
-		// Loop over length weight points
-		for(size_t j=0; j<weights_len.size(); j++) {
-			dp.length = weights_len[j].value;
-
-			// Loop over thickness weight points
-			for(size_t m=0; m<weights_thick.size(); m++) {
-				dp.thickness = weights_thick[m].value;
-
-			// Average over theta distribution
-			for(size_t k=0; k<weights_theta.size(); k++) {
-				dp.axis_theta = weights_theta[k].value;
-
-				// Average over phi distribution
-				for(size_t l=0; l<weights_phi.size(); l++) {
-					dp.axis_phi = weights_phi[l].value;
-					//Un-normalize by volume
-					double _ptvalue = weights_rad[i].weight
-						* weights_len[j].weight
-						* weights_thick[m].weight
-						* weights_theta[k].weight
-						* weights_phi[l].weight
-						* core_shell_cylinder_analytical_2DXY(&dp, qx, qy)
-						* pow(weights_rad[i].value+weights_thick[m].value,2)
-						*(weights_len[j].value+2.0*weights_thick[m].value);
-
-					if (weights_theta.size()>1) {
-						_ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
-					}
-					sum += _ptvalue;
-
-					//Find average volume
-					vol += weights_rad[i].weight
-							* weights_len[j].weight
-							* weights_thick[m].weight
-							* pow(weights_rad[i].value+weights_thick[m].value,2)
-							*(weights_len[j].value+2.0*weights_thick[m].value);
-					//Find norm for volume
-					norm_vol += weights_rad[i].weight
-								* weights_len[j].weight
-								* weights_thick[m].weight;
-
-					norm += weights_rad[i].weight
-						* weights_len[j].weight
-						* weights_thick[m].weight
-						* weights_theta[k].weight
-						* weights_phi[l].weight;
-
-				}
-			}
-			}
-		}
-	}
-	// Averaging in theta needs an extra normalization
-	// factor to account for the sin(theta) term in the
-	// integration (see documentation).
-	if (weights_theta.size()>1) norm = norm / asin(1.0);
-
-	if (vol != 0.0 && norm_vol != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm_vol);}
-
-	return sum/norm + background();
+  CoreShellCylinderParameters dp;
+  // Fill parameter array
+  dp.scale      = scale();
+  dp.radius     = radius();
+  dp.thickness  = thickness();
+  dp.length     = length();
+  dp.core_sld   = core_sld();
+  dp.shell_sld  = shell_sld();
+  dp.solvent_sld= solvent_sld();
+  dp.background = 0.0;
+  dp.axis_theta = axis_theta();
+  dp.axis_phi   = axis_phi();
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad;
+  radius.get_weights(weights_rad);
+
+  // Get the dispersion points for the thickness
+  vector<WeightPoint> weights_thick;
+  thickness.get_weights(weights_thick);
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_len;
+  length.get_weights(weights_len);
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_theta;
+  axis_theta.get_weights(weights_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_phi;
+  axis_phi.get_weights(weights_phi);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double norm_vol = 0.0;
+  double vol = 0.0;
+  double pi = 4.0*atan(1.0);
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp.radius = weights_rad[i].value;
+
+
+    // Loop over length weight points
+    for(size_t j=0; j<weights_len.size(); j++) {
+      dp.length = weights_len[j].value;
+
+      // Loop over thickness weight points
+      for(size_t m=0; m<weights_thick.size(); m++) {
+        dp.thickness = weights_thick[m].value;
+
+        // Average over theta distribution
+        for(size_t k=0; k<weights_theta.size(); k++) {
+          dp.axis_theta = weights_theta[k].value;
+
+          // Average over phi distribution
+          for(size_t l=0; l<weights_phi.size(); l++) {
+            dp.axis_phi = weights_phi[l].value;
+            //Un-normalize by volume
+            double _ptvalue = weights_rad[i].weight
+                * weights_len[j].weight
+                * weights_thick[m].weight
+                * weights_theta[k].weight
+                * weights_phi[l].weight
+                * core_shell_cylinder_analytical_2DXY(&dp, qx, qy)
+            * pow(weights_rad[i].value+weights_thick[m].value,2)
+            *(weights_len[j].value+2.0*weights_thick[m].value);
+
+            if (weights_theta.size()>1) {
+              _ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
+            }
+            sum += _ptvalue;
+
+            //Find average volume
+            vol += weights_rad[i].weight
+                * weights_len[j].weight
+                * weights_thick[m].weight
+                * pow(weights_rad[i].value+weights_thick[m].value,2)
+            *(weights_len[j].value+2.0*weights_thick[m].value);
+            //Find norm for volume
+            norm_vol += weights_rad[i].weight
+                * weights_len[j].weight
+                * weights_thick[m].weight;
+
+            norm += weights_rad[i].weight
+                * weights_len[j].weight
+                * weights_thick[m].weight
+                * weights_theta[k].weight
+                * weights_phi[l].weight;
+
+          }
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_theta.size()>1) norm = norm / asin(1.0);
+
+  if (vol != 0.0 && norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm_vol);}
+
+  return sum/norm + background();
 }
 
@@ -331,7 +331,7 @@
  */
 double CoreShellCylinderModel :: evaluate_rphi(double q, double phi) {
-	double qx = q*cos(phi);
-	double qy = q*sin(phi);
-	return (*this).operator()(qx, qy);
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
 }
 /**
@@ -340,49 +340,49 @@
  */
 double CoreShellCylinderModel :: calculate_ER() {
-	CoreShellCylinderParameters dp;
-
-	dp.radius     = radius();
-	dp.thickness  = thickness();
-	dp.length     = length();
-	double rad_out = 0.0;
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-
-	// Get the dispersion points for the length
-	vector<WeightPoint> weights_length;
-	length.get_weights(weights_length);
-
-	// Get the dispersion points for the thickness
-	vector<WeightPoint> weights_thickness;
-	thickness.get_weights(weights_thickness);
-
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_radius ;
-	radius.get_weights(weights_radius);
-
-	// Loop over major shell weight points
-	for(int i=0; i< (int)weights_length.size(); i++) {
-		dp.length = weights_length[i].value;
-		for(int j=0; j< (int)weights_thickness.size(); j++) {
-			dp.thickness = weights_thickness[j].value;
-			for(int k=0; k< (int)weights_radius.size(); k++) {
-				dp.radius = weights_radius[k].value;
-				//Note: output of "DiamCyl( )" is DIAMETER.
-				sum +=weights_length[i].weight * weights_thickness[j].weight
-					* weights_radius[k].weight*DiamCyl(dp.length+2.0*dp.thickness,dp.radius+dp.thickness)/2.0;
-				norm += weights_length[i].weight* weights_thickness[j].weight* weights_radius[k].weight;
-			}
-		}
-	}
-	if (norm != 0){
-		//return the averaged value
-		rad_out =  sum/norm;}
-	else{
-		//return normal value
-		//Note: output of "DiamCyl()" is DIAMETER.
-		rad_out = DiamCyl(dp.length+2.0*dp.thickness,dp.radius+dp.thickness)/2.0;}
-
-	return rad_out;
-}
+  CoreShellCylinderParameters dp;
+
+  dp.radius     = radius();
+  dp.thickness  = thickness();
+  dp.length     = length();
+  double rad_out = 0.0;
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_length;
+  length.get_weights(weights_length);
+
+  // Get the dispersion points for the thickness
+  vector<WeightPoint> weights_thickness;
+  thickness.get_weights(weights_thickness);
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_radius ;
+  radius.get_weights(weights_radius);
+
+  // Loop over major shell weight points
+  for(int i=0; i< (int)weights_length.size(); i++) {
+    dp.length = weights_length[i].value;
+    for(int j=0; j< (int)weights_thickness.size(); j++) {
+      dp.thickness = weights_thickness[j].value;
+      for(int k=0; k< (int)weights_radius.size(); k++) {
+        dp.radius = weights_radius[k].value;
+        //Note: output of "DiamCyl( )" is DIAMETER.
+        sum +=weights_length[i].weight * weights_thickness[j].weight
+            * weights_radius[k].weight*DiamCyl(dp.length+2.0*dp.thickness,dp.radius+dp.thickness)/2.0;
+        norm += weights_length[i].weight* weights_thickness[j].weight* weights_radius[k].weight;
+      }
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    //Note: output of "DiamCyl()" is DIAMETER.
+    rad_out = DiamCyl(dp.length+2.0*dp.thickness,dp.radius+dp.thickness)/2.0;}
+
+  return rad_out;
+}
Index: sansmodels/src/c_models/ellipsoid.cpp
===================================================================
--- sansmodels/src/c_models/ellipsoid.cpp	(revision 011e0e4569d018a5f39b70032285911f67750383)
+++ sansmodels/src/c_models/ellipsoid.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -27,6 +27,6 @@
 
 extern "C" {
-	#include "libCylinder.h"
-	#include "libStructureFactor.h"
+#include "libCylinder.h"
+#include "libStructureFactor.h"
 }
 
@@ -60,23 +60,23 @@
   double phi = pars->axis_phi * pi/180.0;
 
-    // Ellipsoid orientation
-    cyl_x = sin(theta) * cos(phi);
-    cyl_y = sin(theta) * sin(phi);
-    cyl_z = cos(theta);
-
-    // q vector
-    q_z = 0.0;
-
-    // Compute the angle btw vector q and the
-    // axis of the cylinder
-    cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
-
-    // The following test should always pass
-    if (fabs(cos_val)>1.0) {
-      printf("ellipsoid_ana_2D: Unexpected error: cos(alpha)>1\n");
-      return 0;
-    }
-
-    // Angle between rotation axis and q vector
+  // Ellipsoid orientation
+  cyl_x = sin(theta) * cos(phi);
+  cyl_y = sin(theta) * sin(phi);
+  cyl_z = cos(theta);
+
+  // q vector
+  q_z = 0.0;
+
+  // Compute the angle btw vector q and the
+  // axis of the cylinder
+  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
+
+  // The following test should always pass
+  if (fabs(cos_val)>1.0) {
+    printf("ellipsoid_ana_2D: Unexpected error: cos(alpha)>1\n");
+    return 0;
+  }
+
+  // Angle between rotation axis and q vector
   alpha = acos( cos_val );
 
@@ -88,5 +88,5 @@
 
   //normalize by cylinder volume
-    vol = 4.0/3.0 * acos(-1.0) * pars->radius_b * pars->radius_b * pars->radius_a;
+  vol = 4.0/3.0 * acos(-1.0) * pars->radius_b * pars->radius_b * pars->radius_a;
   answer *= vol;
 
@@ -112,18 +112,18 @@
   double q;
   q = sqrt(qx*qx+qy*qy);
-    return ellipsoid_analytical_2D_scaled(pars, q, qx/q, qy/q);
+  return ellipsoid_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
 EllipsoidModel :: EllipsoidModel() {
-	scale      = Parameter(1.0);
-	radius_a   = Parameter(20.0, true);
-	radius_a.set_min(0.0);
-	radius_b   = Parameter(400.0, true);
-	radius_b.set_min(0.0);
-	sldEll   = Parameter(4.e-6);
-	sldSolv   = Parameter(1.e-6);
-	background = Parameter(0.0);
-	axis_theta  = Parameter(57.325, true);
-	axis_phi    = Parameter(0.0, true);
+  scale      = Parameter(1.0);
+  radius_a   = Parameter(20.0, true);
+  radius_a.set_min(0.0);
+  radius_b   = Parameter(400.0, true);
+  radius_b.set_min(0.0);
+  sldEll   = Parameter(4.e-6);
+  sldSolv   = Parameter(1.e-6);
+  background = Parameter(0.0);
+  axis_theta  = Parameter(57.325, true);
+  axis_phi    = Parameter(0.0, true);
 }
 
@@ -135,55 +135,55 @@
  */
 double EllipsoidModel :: operator()(double q) {
-	double dp[6];
-
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = scale();
-	dp[1] = radius_a();
-	dp[2] = radius_b();
-	dp[3] = sldEll();
-	dp[4] = sldSolv();
-	dp[5] = 0.0;
-
-	// Get the dispersion points for the radius_a
-	vector<WeightPoint> weights_rad_a;
-	radius_a.get_weights(weights_rad_a);
-
-	// Get the dispersion points for the radius_b
-	vector<WeightPoint> weights_rad_b;
-	radius_b.get_weights(weights_rad_b);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
-
-	// Loop over radius_a weight points
-	for(size_t i=0; i<weights_rad_a.size(); i++) {
-		dp[1] = weights_rad_a[i].value;
-
-		// Loop over radius_b weight points
-		for(size_t j=0; j<weights_rad_b.size(); j++) {
-			dp[2] = weights_rad_b[j].value;
-			//Un-normalize  by volume
-			sum += weights_rad_a[i].weight
-				* weights_rad_b[j].weight * EllipsoidForm(dp, q)
-				* pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
-
-			//Find average volume
-			vol += weights_rad_a[i].weight
-				* weights_rad_b[j].weight
-				* pow(weights_rad_b[j].value,2)
-				* weights_rad_a[i].value;
-			norm += weights_rad_a[i].weight
-				* weights_rad_b[j].weight;
-		}
-	}
-
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
-
-	return sum/norm + background();
+  double dp[6];
+
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = radius_a();
+  dp[2] = radius_b();
+  dp[3] = sldEll();
+  dp[4] = sldSolv();
+  dp[5] = 0.0;
+
+  // Get the dispersion points for the radius_a
+  vector<WeightPoint> weights_rad_a;
+  radius_a.get_weights(weights_rad_a);
+
+  // Get the dispersion points for the radius_b
+  vector<WeightPoint> weights_rad_b;
+  radius_b.get_weights(weights_rad_b);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over radius_a weight points
+  for(size_t i=0; i<weights_rad_a.size(); i++) {
+    dp[1] = weights_rad_a[i].value;
+
+    // Loop over radius_b weight points
+    for(size_t j=0; j<weights_rad_b.size(); j++) {
+      dp[2] = weights_rad_b[j].value;
+      //Un-normalize  by volume
+      sum += weights_rad_a[i].weight
+          * weights_rad_b[j].weight * EllipsoidForm(dp, q)
+      * pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
+
+      //Find average volume
+      vol += weights_rad_a[i].weight
+          * weights_rad_b[j].weight
+          * pow(weights_rad_b[j].value,2)
+      * weights_rad_a[i].value;
+      norm += weights_rad_a[i].weight
+          * weights_rad_b[j].weight;
+    }
+  }
+
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -195,91 +195,91 @@
  */
 double EllipsoidModel :: operator()(double qx, double qy) {
-	EllipsoidParameters dp;
-	// Fill parameter array
-	dp.scale      = scale();
-	dp.radius_a   = radius_a();
-	dp.radius_b   = radius_b();
-	dp.sldEll   = sldEll();
-	dp.sldSolv   = sldSolv();
-	dp.background = 0.0;
-	dp.axis_theta = axis_theta();
-	dp.axis_phi   = axis_phi();
-
-	// Get the dispersion points for the radius_a
-	vector<WeightPoint> weights_rad_a;
-	radius_a.get_weights(weights_rad_a);
-
-	// Get the dispersion points for the radius_b
-	vector<WeightPoint> weights_rad_b;
-	radius_b.get_weights(weights_rad_b);
-
-	// Get angular averaging for theta
-	vector<WeightPoint> weights_theta;
-	axis_theta.get_weights(weights_theta);
-
-	// Get angular averaging for phi
-	vector<WeightPoint> weights_phi;
-	axis_phi.get_weights(weights_phi);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double norm_vol = 0.0;
-	double vol = 0.0;
-	double pi = 4.0*atan(1.0);
-	// Loop over radius weight points
-	for(size_t i=0; i<weights_rad_a.size(); i++) {
-		dp.radius_a = weights_rad_a[i].value;
-
-
-		// Loop over length weight points
-		for(size_t j=0; j<weights_rad_b.size(); j++) {
-			dp.radius_b = weights_rad_b[j].value;
-
-			// Average over theta distribution
-			for(size_t k=0; k<weights_theta.size(); k++) {
-				dp.axis_theta = weights_theta[k].value;
-
-				// Average over phi distribution
-				for(size_t l=0; l<weights_phi.size(); l++) {
-					dp.axis_phi = weights_phi[l].value;
-					//Un-normalize by volume
-					double _ptvalue = weights_rad_a[i].weight
-						* weights_rad_b[j].weight
-						* weights_theta[k].weight
-						* weights_phi[l].weight
-						* ellipsoid_analytical_2DXY(&dp, qx, qy)
-						* pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
-					if (weights_theta.size()>1) {
-						_ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
-					}
-					sum += _ptvalue;
-					//Find average volume
-					vol += weights_rad_a[i].weight
-						* weights_rad_b[j].weight
-						* pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
-					//Find norm for volume
-					norm_vol += weights_rad_a[i].weight
-						* weights_rad_b[j].weight;
-
-					norm += weights_rad_a[i].weight
-						* weights_rad_b[j].weight
-						* weights_theta[k].weight
-						* weights_phi[l].weight;
-
-				}
-			}
-		}
-	}
-	// Averaging in theta needs an extra normalization
-	// factor to account for the sin(theta) term in the
-	// integration (see documentation).
-	if (weights_theta.size()>1) norm = norm / asin(1.0);
-
-	if (vol != 0.0 && norm_vol != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm_vol);}
-
-	return sum/norm + background();
+  EllipsoidParameters dp;
+  // Fill parameter array
+  dp.scale      = scale();
+  dp.radius_a   = radius_a();
+  dp.radius_b   = radius_b();
+  dp.sldEll   = sldEll();
+  dp.sldSolv   = sldSolv();
+  dp.background = 0.0;
+  dp.axis_theta = axis_theta();
+  dp.axis_phi   = axis_phi();
+
+  // Get the dispersion points for the radius_a
+  vector<WeightPoint> weights_rad_a;
+  radius_a.get_weights(weights_rad_a);
+
+  // Get the dispersion points for the radius_b
+  vector<WeightPoint> weights_rad_b;
+  radius_b.get_weights(weights_rad_b);
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_theta;
+  axis_theta.get_weights(weights_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_phi;
+  axis_phi.get_weights(weights_phi);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double norm_vol = 0.0;
+  double vol = 0.0;
+  double pi = 4.0*atan(1.0);
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad_a.size(); i++) {
+    dp.radius_a = weights_rad_a[i].value;
+
+
+    // Loop over length weight points
+    for(size_t j=0; j<weights_rad_b.size(); j++) {
+      dp.radius_b = weights_rad_b[j].value;
+
+      // Average over theta distribution
+      for(size_t k=0; k<weights_theta.size(); k++) {
+        dp.axis_theta = weights_theta[k].value;
+
+        // Average over phi distribution
+        for(size_t l=0; l<weights_phi.size(); l++) {
+          dp.axis_phi = weights_phi[l].value;
+          //Un-normalize by volume
+          double _ptvalue = weights_rad_a[i].weight
+              * weights_rad_b[j].weight
+              * weights_theta[k].weight
+              * weights_phi[l].weight
+              * ellipsoid_analytical_2DXY(&dp, qx, qy)
+          * pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
+          if (weights_theta.size()>1) {
+            _ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
+          }
+          sum += _ptvalue;
+          //Find average volume
+          vol += weights_rad_a[i].weight
+              * weights_rad_b[j].weight
+              * pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
+          //Find norm for volume
+          norm_vol += weights_rad_a[i].weight
+              * weights_rad_b[j].weight;
+
+          norm += weights_rad_a[i].weight
+              * weights_rad_b[j].weight
+              * weights_theta[k].weight
+              * weights_phi[l].weight;
+
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_theta.size()>1) norm = norm / asin(1.0);
+
+  if (vol != 0.0 && norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm_vol);}
+
+  return sum/norm + background();
 }
 
@@ -292,7 +292,7 @@
  */
 double EllipsoidModel :: evaluate_rphi(double q, double phi) {
-	double qx = q*cos(phi);
-	double qy = q*sin(phi);
-	return (*this).operator()(qx, qy);
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
 }
 
@@ -302,40 +302,40 @@
  */
 double EllipsoidModel :: calculate_ER() {
-	EllipsoidParameters dp;
-
-	dp.radius_a = radius_a();
-	dp.radius_b = radius_b();
-
-	double rad_out = 0.0;
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-
-	// Get the dispersion points for the major shell
-	vector<WeightPoint> weights_radius_a;
-	radius_a.get_weights(weights_radius_a);
-
-	// Get the dispersion points for the minor shell
-	vector<WeightPoint> weights_radius_b;
-	radius_b.get_weights(weights_radius_b);
-
-	// Loop over major shell weight points
-	for(int i=0; i< (int)weights_radius_b.size(); i++) {
-		dp.radius_b = weights_radius_b[i].value;
-		for(int k=0; k< (int)weights_radius_a.size(); k++) {
-			dp.radius_a = weights_radius_a[k].value;
-			sum +=weights_radius_b[i].weight
-				* weights_radius_a[k].weight*DiamEllip(dp.radius_a,dp.radius_b)/2.0;
-			norm += weights_radius_b[i].weight* weights_radius_a[k].weight;
-		}
-	}
-	if (norm != 0){
-		//return the averaged value
-		rad_out =  sum/norm;}
-	else{
-		//return normal value
-		rad_out = DiamEllip(dp.radius_a,dp.radius_b)/2.0;}
-
-	return rad_out;
-}
+  EllipsoidParameters dp;
+
+  dp.radius_a = radius_a();
+  dp.radius_b = radius_b();
+
+  double rad_out = 0.0;
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the major shell
+  vector<WeightPoint> weights_radius_a;
+  radius_a.get_weights(weights_radius_a);
+
+  // Get the dispersion points for the minor shell
+  vector<WeightPoint> weights_radius_b;
+  radius_b.get_weights(weights_radius_b);
+
+  // Loop over major shell weight points
+  for(int i=0; i< (int)weights_radius_b.size(); i++) {
+    dp.radius_b = weights_radius_b[i].value;
+    for(int k=0; k< (int)weights_radius_a.size(); k++) {
+      dp.radius_a = weights_radius_a[k].value;
+      sum +=weights_radius_b[i].weight
+          * weights_radius_a[k].weight*DiamEllip(dp.radius_a,dp.radius_b)/2.0;
+      norm += weights_radius_b[i].weight* weights_radius_a[k].weight;
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    rad_out = DiamEllip(dp.radius_a,dp.radius_b)/2.0;}
+
+  return rad_out;
+}
Index: sansmodels/src/c_models/ellipticalcylinder.cpp
===================================================================
--- sansmodels/src/c_models/ellipticalcylinder.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/ellipticalcylinder.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -18,33 +18,173 @@
  *   sansmodels/src/libigor
  *
- *	TODO: refactor so that we pull in the old sansmodels.c_extensions
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
+#include <stdlib.h>
 using namespace std;
+#include "elliptical_cylinder.h"
 
 extern "C" {
-	#include "libCylinder.h"
-	#include "libStructureFactor.h"
-	#include "elliptical_cylinder.h"
+#include "libCylinder.h"
+#include "libStructureFactor.h"
+}
+
+typedef struct {
+  double scale;
+  double r_minor;
+  double r_ratio;
+  double length;
+  double sldCyl;
+  double sldSolv;
+  double background;
+  double cyl_theta;
+  double cyl_phi;
+  double cyl_psi;
+} EllipticalCylinderParameters;
+
+
+static double elliptical_cylinder_kernel(EllipticalCylinderParameters *pars, double q, double alpha, double nu) {
+  double qr;
+  double qL;
+  double Be,Si;
+  double r_major;
+  double kernel;
+
+  r_major = pars->r_ratio * pars->r_minor;
+
+  qr = q*sin(alpha)*sqrt( r_major*r_major*sin(nu)*sin(nu) + pars->r_minor*pars->r_minor*cos(nu)*cos(nu) );
+  qL = q*pars->length*cos(alpha)/2.0;
+
+  if (qr==0){
+    Be = 0.5;
+  }else{
+    Be = NR_BessJ1(qr)/qr;
+  }
+  if (qL==0){
+    Si = 1.0;
+  }else{
+    Si = sin(qL)/qL;
+  }
+
+
+  kernel = 2.0*Be * Si;
+  return kernel*kernel;
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double elliptical_cylinder_analytical_2D_scaled(EllipticalCylinderParameters *pars, double q, double q_x, double q_y) {
+  double cyl_x, cyl_y, cyl_z;
+  double ell_x, ell_y;
+  double q_z;
+  double alpha, vol, cos_val;
+  double nu, cos_nu;
+  double answer;
+  //convert angle degree to radian
+  double pi = 4.0*atan(1.0);
+  double theta = pars->cyl_theta * pi/180.0;
+  double phi = pars->cyl_phi * pi/180.0;
+  double psi = pars->cyl_psi * pi/180.0;
+
+  //Cylinder orientation
+  cyl_x = sin(theta) * cos(phi);
+  cyl_y = sin(theta) * sin(phi);
+  cyl_z = cos(theta);
+
+  // q vector
+  q_z = 0;
+
+  // Compute the angle btw vector q and the
+  // axis of the cylinder
+  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
+
+  // The following test should always pass
+  if (fabs(cos_val)>1.0) {
+    printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
+    return 0;
+  }
+
+  // Note: cos(alpha) = 0 and 1 will get an
+  // undefined value from CylKernel
+  alpha = acos( cos_val );
+
+  //ellipse orientation:
+  // the elliptical corss section was transformed and projected
+  // into the detector plane already through sin(alpha)and furthermore psi remains as same
+  // on the detector plane.
+  // So, all we need is to calculate the angle (nu) of the minor axis of the ellipse wrt
+  // the wave vector q.
+
+  //x- y- component on the detector plane.
+  ell_x =  cos(psi);
+  ell_y =  sin(psi);
+
+  // calculate the axis of the ellipse wrt q-coord.
+  cos_nu = ell_x*q_x + ell_y*q_y;
+  nu = acos(cos_nu);
+
+  // The following test should always pass
+  if (fabs(cos_nu)>1.0) {
+    printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n");
+    return 0;
+  }
+
+  answer = elliptical_cylinder_kernel(pars, q, alpha,nu);
+
+  // Multiply by contrast^2
+  answer *= (pars->sldCyl - pars->sldSolv) * (pars->sldCyl - pars->sldSolv);
+
+  //normalize by cylinder volume
+  //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
+  vol = acos(-1.0) * pars->r_minor * pars->r_minor * pars->r_ratio * pars->length;
+  answer *= vol;
+
+  //convert to [cm-1]
+  answer *= 1.0e8;
+
+  //Scale
+  answer *= pars->scale;
+
+  // add in the background
+  answer += pars->background;
+
+  return answer;
+}
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @return: function value
+ */
+static double elliptical_cylinder_analytical_2DXY(EllipticalCylinderParameters *pars, double qx, double qy) {
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return elliptical_cylinder_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
 EllipticalCylinderModel :: EllipticalCylinderModel() {
-	scale      = Parameter(1.0);
-	r_minor    = Parameter(20.0, true);
-	r_minor.set_min(0.0);
-	r_ratio    = Parameter(1.5, true);
-	r_ratio.set_min(0.0);
-	length     = Parameter(400.0, true);
-	length.set_min(0.0);
-	sldCyl   = Parameter(4.e-6);
-	sldSolv   = Parameter(1.e-6);
-	background = Parameter(0.0);
-	cyl_theta  = Parameter(57.325, true);
-	cyl_phi    = Parameter(0.0, true);
-	cyl_psi    = Parameter(0.0, true);
+  scale      = Parameter(1.0);
+  r_minor    = Parameter(20.0, true);
+  r_minor.set_min(0.0);
+  r_ratio    = Parameter(1.5, true);
+  r_ratio.set_min(0.0);
+  length     = Parameter(400.0, true);
+  length.set_min(0.0);
+  sldCyl   = Parameter(4.e-6);
+  sldSolv   = Parameter(1.e-6);
+  background = Parameter(0.0);
+  cyl_theta  = Parameter(57.325, true);
+  cyl_phi    = Parameter(0.0, true);
+  cyl_psi    = Parameter(0.0, true);
 }
 
@@ -56,67 +196,67 @@
  */
 double EllipticalCylinderModel :: operator()(double q) {
-	double dp[7];
-
-	dp[0] = scale();
-	dp[1] = r_minor();
-	dp[2] = r_ratio();
-	dp[3] = length();
-	dp[4] = sldCyl();
-	dp[5] = sldSolv();
-	dp[6] = 0.0;
-
-	// Get the dispersion points for the r_minor
-	vector<WeightPoint> weights_rad;
-	r_minor.get_weights(weights_rad);
-
-	// Get the dispersion points for the r_ratio
-	vector<WeightPoint> weights_rat;
-	r_ratio.get_weights(weights_rat);
-
-	// Get the dispersion points for the length
-	vector<WeightPoint> weights_len;
-	length.get_weights(weights_len);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
-
-	// Loop over r_minor weight points
-	for(size_t i=0; i<weights_rad.size(); i++) {
-		dp[1] = weights_rad[i].value;
-
-		// Loop over r_ratio weight points
-		for(size_t j=0; j<weights_rat.size(); j++) {
-			dp[2] = weights_rat[j].value;
-
-			// Loop over length weight points
-			for(size_t k=0; k<weights_len.size(); k++) {
-				dp[3] = weights_len[k].value;
-				//Un-normalize  by volume
-				sum += weights_rad[i].weight
-					* weights_len[k].weight
-					* weights_rat[j].weight
-					* EllipCyl20(dp, q)
-					* pow(weights_rad[i].value,2) * weights_rat[j].value
-					* weights_len[k].value;
-				//Find average volume
-				vol += weights_rad[i].weight
-					* weights_len[k].weight
-					* weights_rat[j].weight
-					* pow(weights_rad[i].value,2) * weights_rat[j].value
-					* weights_len[k].value;
-				norm += weights_rad[i].weight
-					* weights_len[k].weight
-					* weights_rat[j].weight;
-			}
-		}
-	}
-
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
-
-	return sum/norm + background();
+  double dp[7];
+
+  dp[0] = scale();
+  dp[1] = r_minor();
+  dp[2] = r_ratio();
+  dp[3] = length();
+  dp[4] = sldCyl();
+  dp[5] = sldSolv();
+  dp[6] = 0.0;
+
+  // Get the dispersion points for the r_minor
+  vector<WeightPoint> weights_rad;
+  r_minor.get_weights(weights_rad);
+
+  // Get the dispersion points for the r_ratio
+  vector<WeightPoint> weights_rat;
+  r_ratio.get_weights(weights_rat);
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_len;
+  length.get_weights(weights_len);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over r_minor weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp[1] = weights_rad[i].value;
+
+    // Loop over r_ratio weight points
+    for(size_t j=0; j<weights_rat.size(); j++) {
+      dp[2] = weights_rat[j].value;
+
+      // Loop over length weight points
+      for(size_t k=0; k<weights_len.size(); k++) {
+        dp[3] = weights_len[k].value;
+        //Un-normalize  by volume
+        sum += weights_rad[i].weight
+            * weights_len[k].weight
+            * weights_rat[j].weight
+            * EllipCyl20(dp, q)
+        * pow(weights_rad[i].value,2) * weights_rat[j].value
+        * weights_len[k].value;
+        //Find average volume
+        vol += weights_rad[i].weight
+            * weights_len[k].weight
+            * weights_rat[j].weight
+            * pow(weights_rad[i].value,2) * weights_rat[j].value
+            * weights_len[k].value;
+        norm += weights_rad[i].weight
+            * weights_len[k].weight
+            * weights_rat[j].weight;
+      }
+    }
+  }
+
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -128,121 +268,121 @@
  */
 double EllipticalCylinderModel :: operator()(double qx, double qy) {
-	EllipticalCylinderParameters dp;
-	// Fill parameter array
-	dp.scale      = scale();
-	dp.r_minor    = r_minor();
-	dp.r_ratio    = r_ratio();
-	dp.length     = length();
-	dp.sldCyl   = sldCyl();
-	dp.sldSolv   = sldSolv();
-	dp.background = 0.0;
-	dp.cyl_theta  = cyl_theta();
-	dp.cyl_phi    = cyl_phi();
-	dp.cyl_psi    = cyl_psi();
-
-	// Get the dispersion points for the r_minor
-	vector<WeightPoint> weights_rad;
-	r_minor.get_weights(weights_rad);
-
-	// Get the dispersion points for the r_ratio
-	vector<WeightPoint> weights_rat;
-	r_ratio.get_weights(weights_rat);
-
-	// Get the dispersion points for the length
-	vector<WeightPoint> weights_len;
-	length.get_weights(weights_len);
-
-	// Get angular averaging for theta
-	vector<WeightPoint> weights_theta;
-	cyl_theta.get_weights(weights_theta);
-
-	// Get angular averaging for phi
-	vector<WeightPoint> weights_phi;
-	cyl_phi.get_weights(weights_phi);
-
-	// Get angular averaging for psi
-	vector<WeightPoint> weights_psi;
-	cyl_psi.get_weights(weights_psi);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double norm_vol = 0.0;
-	double vol = 0.0;
-	double pi = 4.0*atan(1.0);
-	// Loop over minor radius weight points
-	for(size_t i=0; i<weights_rad.size(); i++) {
-		dp.r_minor = weights_rad[i].value;
-
-
-		// Loop over length weight points
-		for(size_t j=0; j<weights_len.size(); j++) {
-			dp.length = weights_len[j].value;
-
-			// Loop over r_ration weight points
-			for(size_t m=0; m<weights_rat.size(); m++) {
-				dp.r_ratio = weights_rat[m].value;
-
-			// Average over theta distribution
-			for(size_t k=0; k<weights_theta.size(); k++) {
-				dp.cyl_theta = weights_theta[k].value;
-
-				// Average over phi distribution
-				for(size_t l=0; l<weights_phi.size(); l++) {
-					dp.cyl_phi = weights_phi[l].value;
-
-				// Average over phi distribution
-				for(size_t o=0; o<weights_psi.size(); o++) {
-					dp.cyl_psi = weights_psi[o].value;
-					//Un-normalize by volume
-					double _ptvalue = weights_rad[i].weight
-						* weights_len[j].weight
-						* weights_rat[m].weight
-						* weights_theta[k].weight
-						* weights_phi[l].weight
-						* weights_psi[o].weight
-						* elliptical_cylinder_analytical_2DXY(&dp, qx, qy)
-						* pow(weights_rad[i].value,2)
-						* weights_len[j].value
-						* weights_rat[m].value;
-					if (weights_theta.size()>1) {
-						_ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
-					}
-					sum += _ptvalue;
-					//Find average volume
-					vol += weights_rad[i].weight
-						* weights_len[j].weight
-						* weights_rat[m].weight
-						* pow(weights_rad[i].value,2)
-						* weights_len[j].value
-						* weights_rat[m].value;
-					//Find norm for volume
-					norm_vol += weights_rad[i].weight
-						* weights_len[j].weight
-						* weights_rat[m].weight;
-
-					norm += weights_rad[i].weight
-						* weights_len[j].weight
-						* weights_rat[m].weight
-						* weights_theta[k].weight
-						* weights_phi[l].weight
-						* weights_psi[o].weight;
-
-				}
-				}
-			}
-			}
-		}
-	}
-	// Averaging in theta needs an extra normalization
-	// factor to account for the sin(theta) term in the
-	// integration (see documentation).
-	if (weights_theta.size()>1) norm = norm / asin(1.0);
-
-	if (vol != 0.0 && norm_vol != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm_vol);}
-
-	return sum/norm + background();
+  EllipticalCylinderParameters dp;
+  // Fill parameter array
+  dp.scale      = scale();
+  dp.r_minor    = r_minor();
+  dp.r_ratio    = r_ratio();
+  dp.length     = length();
+  dp.sldCyl   = sldCyl();
+  dp.sldSolv   = sldSolv();
+  dp.background = 0.0;
+  dp.cyl_theta  = cyl_theta();
+  dp.cyl_phi    = cyl_phi();
+  dp.cyl_psi    = cyl_psi();
+
+  // Get the dispersion points for the r_minor
+  vector<WeightPoint> weights_rad;
+  r_minor.get_weights(weights_rad);
+
+  // Get the dispersion points for the r_ratio
+  vector<WeightPoint> weights_rat;
+  r_ratio.get_weights(weights_rat);
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_len;
+  length.get_weights(weights_len);
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_theta;
+  cyl_theta.get_weights(weights_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_phi;
+  cyl_phi.get_weights(weights_phi);
+
+  // Get angular averaging for psi
+  vector<WeightPoint> weights_psi;
+  cyl_psi.get_weights(weights_psi);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double norm_vol = 0.0;
+  double vol = 0.0;
+  double pi = 4.0*atan(1.0);
+  // Loop over minor radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp.r_minor = weights_rad[i].value;
+
+
+    // Loop over length weight points
+    for(size_t j=0; j<weights_len.size(); j++) {
+      dp.length = weights_len[j].value;
+
+      // Loop over r_ration weight points
+      for(size_t m=0; m<weights_rat.size(); m++) {
+        dp.r_ratio = weights_rat[m].value;
+
+        // Average over theta distribution
+        for(size_t k=0; k<weights_theta.size(); k++) {
+          dp.cyl_theta = weights_theta[k].value;
+
+          // Average over phi distribution
+          for(size_t l=0; l<weights_phi.size(); l++) {
+            dp.cyl_phi = weights_phi[l].value;
+
+            // Average over phi distribution
+            for(size_t o=0; o<weights_psi.size(); o++) {
+              dp.cyl_psi = weights_psi[o].value;
+              //Un-normalize by volume
+              double _ptvalue = weights_rad[i].weight
+                  * weights_len[j].weight
+                  * weights_rat[m].weight
+                  * weights_theta[k].weight
+                  * weights_phi[l].weight
+                  * weights_psi[o].weight
+                  * elliptical_cylinder_analytical_2DXY(&dp, qx, qy)
+              * pow(weights_rad[i].value,2)
+              * weights_len[j].value
+              * weights_rat[m].value;
+              if (weights_theta.size()>1) {
+                _ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
+              }
+              sum += _ptvalue;
+              //Find average volume
+              vol += weights_rad[i].weight
+                  * weights_len[j].weight
+                  * weights_rat[m].weight
+                  * pow(weights_rad[i].value,2)
+              * weights_len[j].value
+              * weights_rat[m].value;
+              //Find norm for volume
+              norm_vol += weights_rad[i].weight
+                  * weights_len[j].weight
+                  * weights_rat[m].weight;
+
+              norm += weights_rad[i].weight
+                  * weights_len[j].weight
+                  * weights_rat[m].weight
+                  * weights_theta[k].weight
+                  * weights_phi[l].weight
+                  * weights_psi[o].weight;
+
+            }
+          }
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_theta.size()>1) norm = norm / asin(1.0);
+
+  if (vol != 0.0 && norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm_vol);}
+
+  return sum/norm + background();
 
 }
@@ -256,7 +396,7 @@
  */
 double EllipticalCylinderModel :: evaluate_rphi(double q, double phi) {
-	double qx = q*cos(phi);
-	double qy = q*sin(phi);
-	return (*this).operator()(qx, qy);
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
 }
 /**
@@ -265,56 +405,56 @@
  */
 double EllipticalCylinderModel :: calculate_ER() {
-	EllipticalCylinderParameters dp;
-	dp.r_minor    = r_minor();
-	dp.r_ratio    = r_ratio();
-	dp.length     = length();
-	double rad_out = 0.0;
-	double suf_rad = sqrt(dp.r_minor*dp.r_minor*dp.r_ratio);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-
-	// Get the dispersion points for the r_minor
-	vector<WeightPoint> weights_rad;
-	r_minor.get_weights(weights_rad);
-
-	// Get the dispersion points for the r_ratio
-	vector<WeightPoint> weights_rat;
-	r_ratio.get_weights(weights_rat);
-
-	// Get the dispersion points for the length
-	vector<WeightPoint> weights_len;
-	length.get_weights(weights_len);
-
-	// Loop over minor radius weight points
-	for(size_t i=0; i<weights_rad.size(); i++) {
-		dp.r_minor = weights_rad[i].value;
-
-		// Loop over length weight points
-		for(size_t j=0; j<weights_len.size(); j++) {
-			dp.length = weights_len[j].value;
-
-			// Loop over r_ration weight points
-			for(size_t m=0; m<weights_rat.size(); m++) {
-				dp.r_ratio = weights_rat[m].value;
-				//Calculate surface averaged radius
-				suf_rad = sqrt(dp.r_minor * dp.r_minor * dp.r_ratio);
-
-				//Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-				sum +=weights_rad[i].weight *weights_len[j].weight
-					* weights_rat[m].weight*DiamCyl(dp.length, suf_rad)/2.0;
-				norm += weights_rad[i].weight *weights_len[j].weight* weights_rat[m].weight;
-			}
-		}
-	}
-	if (norm != 0){
-		//return the averaged value
-		rad_out =  sum/norm;}
-	else{
-		//return normal value
-		//Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-		rad_out = DiamCyl(dp.length, suf_rad)/2.0;}
-
-	return rad_out;
-}
+  EllipticalCylinderParameters dp;
+  dp.r_minor    = r_minor();
+  dp.r_ratio    = r_ratio();
+  dp.length     = length();
+  double rad_out = 0.0;
+  double suf_rad = sqrt(dp.r_minor*dp.r_minor*dp.r_ratio);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the r_minor
+  vector<WeightPoint> weights_rad;
+  r_minor.get_weights(weights_rad);
+
+  // Get the dispersion points for the r_ratio
+  vector<WeightPoint> weights_rat;
+  r_ratio.get_weights(weights_rat);
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_len;
+  length.get_weights(weights_len);
+
+  // Loop over minor radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp.r_minor = weights_rad[i].value;
+
+    // Loop over length weight points
+    for(size_t j=0; j<weights_len.size(); j++) {
+      dp.length = weights_len[j].value;
+
+      // Loop over r_ration weight points
+      for(size_t m=0; m<weights_rat.size(); m++) {
+        dp.r_ratio = weights_rat[m].value;
+        //Calculate surface averaged radius
+        suf_rad = sqrt(dp.r_minor * dp.r_minor * dp.r_ratio);
+
+        //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+        sum +=weights_rad[i].weight *weights_len[j].weight
+            * weights_rat[m].weight*DiamCyl(dp.length, suf_rad)/2.0;
+        norm += weights_rad[i].weight *weights_len[j].weight* weights_rat[m].weight;
+      }
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+    rad_out = DiamCyl(dp.length, suf_rad)/2.0;}
+
+  return rad_out;
+}
Index: sansmodels/src/c_models/flexcyl_ellipX.cpp
===================================================================
--- sansmodels/src/c_models/flexcyl_ellipX.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/flexcyl_ellipX.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -21,28 +21,38 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "flexcyl_ellipX.h"
 
 extern "C" {
-	#include "libCylinder.h"
-	#include "libStructureFactor.h"
-	#include "flexcyl_ellipX.h"
-}
+#include "libCylinder.h"
+#include "libStructureFactor.h"
+}
+
+typedef struct {
+  double scale;
+  double length;
+  double kuhn_length;
+  double radius;
+  double axis_ratio;
+  double sldCyl;
+  double sldSolv;
+  double background;
+} FlexCylEXParameters;
 
 FlexCylEllipXModel :: FlexCylEllipXModel() {
-	scale      = Parameter(1.0);
-	length     = Parameter(1000.0, true);
-	length.set_min(0.0);
-	kuhn_length = Parameter(100.0, true);
-	kuhn_length.set_min(0.0);
-	radius  = Parameter(20.0, true);
-	radius.set_min(0.0);
-	axis_ratio = Parameter(1.5);
-	axis_ratio.set_min(0.0);
-	sldCyl   = Parameter(1.0e-6);
-	sldSolv   = Parameter(6.3e-6);
-	background = Parameter(0.0001);
+  scale      = Parameter(1.0);
+  length     = Parameter(1000.0, true);
+  length.set_min(0.0);
+  kuhn_length = Parameter(100.0, true);
+  kuhn_length.set_min(0.0);
+  radius  = Parameter(20.0, true);
+  radius.set_min(0.0);
+  axis_ratio = Parameter(1.5);
+  axis_ratio.set_min(0.0);
+  sldCyl   = Parameter(1.0e-6);
+  sldSolv   = Parameter(6.3e-6);
+  background = Parameter(0.0001);
 }
 
@@ -54,72 +64,72 @@
  */
 double FlexCylEllipXModel :: operator()(double q) {
-	double dp[8];
-
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = scale();
-	dp[1] = length();
-	dp[2] = kuhn_length();
-	dp[3] = radius();
-	dp[4] = axis_ratio();
-	dp[5] = sldCyl();
-	dp[6] = sldSolv();
-	dp[7] = 0.0;
-
-	// Get the dispersion points for the length
-	vector<WeightPoint> weights_len;
-	length.get_weights(weights_len);
-
-	// Get the dispersion points for the kuhn_length
-	vector<WeightPoint> weights_kuhn;
-	kuhn_length.get_weights(weights_kuhn);
-
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_rad;
-	radius.get_weights(weights_rad);
-
-	// Get the dispersion points for the axis_ratio
-	vector<WeightPoint> weights_ratio;
-	axis_ratio.get_weights(weights_ratio);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
-
-	// Loop over length weight points
-	for(int i=0; i< (int)weights_len.size(); i++) {
-		dp[1] = weights_len[i].value;
-
-		// Loop over kuhn_length weight points
-		for(int j=0; j< (int)weights_kuhn.size(); j++) {
-			dp[2] = weights_kuhn[j].value;
-
-			// Loop over radius weight points
-			for(int k=0; k< (int)weights_rad.size(); k++) {
-				dp[3] = weights_rad[k].value;
-				// Loop over axis_ratio weight points
-				for(int l=0; l< (int)weights_ratio.size(); l++) {
-					dp[4] = weights_ratio[l].value;
-
-					//Un-normalize by volume
-					sum += weights_len[i].weight * weights_kuhn[j].weight*weights_rad[k].weight
-						* weights_ratio[l].weight * FlexCyl_Ellip(dp, q)
-						* (pow(weights_rad[k].value,2.0) * weights_ratio[l].value * weights_len[i].value);
-					//Find weighted volume
-					vol += weights_rad[k].weight * weights_kuhn[j].weight
-						* weights_len[i].weight * weights_ratio[l].weight
-						*pow(weights_rad[k].value,2.0)* weights_ratio[l].weight*weights_len[i].value;
-					norm += weights_len[i].weight * weights_kuhn[j].weight
-						*weights_rad[k].weight* weights_ratio[l].weight;
-				}
-			}
-		}
-	}
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
-
-	return sum/norm + background();
+  double dp[8];
+
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = length();
+  dp[2] = kuhn_length();
+  dp[3] = radius();
+  dp[4] = axis_ratio();
+  dp[5] = sldCyl();
+  dp[6] = sldSolv();
+  dp[7] = 0.0;
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_len;
+  length.get_weights(weights_len);
+
+  // Get the dispersion points for the kuhn_length
+  vector<WeightPoint> weights_kuhn;
+  kuhn_length.get_weights(weights_kuhn);
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad;
+  radius.get_weights(weights_rad);
+
+  // Get the dispersion points for the axis_ratio
+  vector<WeightPoint> weights_ratio;
+  axis_ratio.get_weights(weights_ratio);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over length weight points
+  for(int i=0; i< (int)weights_len.size(); i++) {
+    dp[1] = weights_len[i].value;
+
+    // Loop over kuhn_length weight points
+    for(int j=0; j< (int)weights_kuhn.size(); j++) {
+      dp[2] = weights_kuhn[j].value;
+
+      // Loop over radius weight points
+      for(int k=0; k< (int)weights_rad.size(); k++) {
+        dp[3] = weights_rad[k].value;
+        // Loop over axis_ratio weight points
+        for(int l=0; l< (int)weights_ratio.size(); l++) {
+          dp[4] = weights_ratio[l].value;
+
+          //Un-normalize by volume
+          sum += weights_len[i].weight * weights_kuhn[j].weight*weights_rad[k].weight
+              * weights_ratio[l].weight * FlexCyl_Ellip(dp, q)
+          * (pow(weights_rad[k].value,2.0) * weights_ratio[l].value * weights_len[i].value);
+          //Find weighted volume
+          vol += weights_rad[k].weight * weights_kuhn[j].weight
+              * weights_len[i].weight * weights_ratio[l].weight
+              *pow(weights_rad[k].value,2.0)* weights_ratio[l].weight*weights_len[i].value;
+          norm += weights_len[i].weight * weights_kuhn[j].weight
+              *weights_rad[k].weight* weights_ratio[l].weight;
+        }
+      }
+    }
+  }
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -131,6 +141,6 @@
  */
 double FlexCylEllipXModel :: operator()(double qx, double qy) {
-	double q = sqrt(qx*qx + qy*qy);
-	return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -143,7 +153,7 @@
  */
 double FlexCylEllipXModel :: evaluate_rphi(double q, double phi) {
-	//double qx = q*cos(phi);
-	//double qy = q*sin(phi);
-	return (*this).operator()(q);
+  //double qx = q*cos(phi);
+  //double qy = q*sin(phi);
+  return (*this).operator()(q);
 }
 /**
@@ -152,52 +162,52 @@
  */
 double FlexCylEllipXModel :: calculate_ER() {
-	FlexCylEXParameters dp;
-
-	dp.radius  = radius();
-	dp.length     = length();
-	dp.axis_ratio  = axis_ratio();
-
-	double rad_out = 0.0;
-	double suf_rad = sqrt(dp.radius*dp.radius*dp.axis_ratio );
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-
-	// Get the dispersion points for the total length
-	vector<WeightPoint> weights_length;
-	length.get_weights(weights_length);
-
-	// Get the dispersion points for minor radius
-	vector<WeightPoint> weights_radius ;
-	radius.get_weights(weights_radius);
-
-	// Get the dispersion points for axis ratio = major_radius/minor_radius
-	vector<WeightPoint> weights_ratio ;
-	axis_ratio.get_weights(weights_ratio);
-
-	// Loop over major shell weight points
-	for(int i=0; i< (int)weights_length.size(); i++) {
-		dp.length = weights_length[i].value;
-		for(int k=0; k< (int)weights_radius.size(); k++) {
-			dp.radius = weights_radius[k].value;
-			// Loop over axis_ratio weight points
-			for(int l=0; l< (int)weights_ratio.size(); l++) {
-				dp.axis_ratio = weights_ratio[l].value;
-				suf_rad = sqrt(dp.radius  * dp.radius  * dp.axis_ratio);
-				//Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-				sum +=weights_length[i].weight * weights_radius[k].weight
-						* weights_ratio[l].weight *DiamCyl(dp.length,suf_rad)/2.0;
-				norm += weights_length[i].weight* weights_radius[k].weight* weights_ratio[l].weight;
-			}
-		}
-	}
-	if (norm != 0){
-		//return the averaged value
-		rad_out =  sum/norm;}
-	else{
-		//return normal value
-		//Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-		rad_out = DiamCyl(dp.length,suf_rad)/2.0;}
-
-	return rad_out;
-}
+  FlexCylEXParameters dp;
+
+  dp.radius  = radius();
+  dp.length     = length();
+  dp.axis_ratio  = axis_ratio();
+
+  double rad_out = 0.0;
+  double suf_rad = sqrt(dp.radius*dp.radius*dp.axis_ratio );
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the total length
+  vector<WeightPoint> weights_length;
+  length.get_weights(weights_length);
+
+  // Get the dispersion points for minor radius
+  vector<WeightPoint> weights_radius ;
+  radius.get_weights(weights_radius);
+
+  // Get the dispersion points for axis ratio = major_radius/minor_radius
+  vector<WeightPoint> weights_ratio ;
+  axis_ratio.get_weights(weights_ratio);
+
+  // Loop over major shell weight points
+  for(int i=0; i< (int)weights_length.size(); i++) {
+    dp.length = weights_length[i].value;
+    for(int k=0; k< (int)weights_radius.size(); k++) {
+      dp.radius = weights_radius[k].value;
+      // Loop over axis_ratio weight points
+      for(int l=0; l< (int)weights_ratio.size(); l++) {
+        dp.axis_ratio = weights_ratio[l].value;
+        suf_rad = sqrt(dp.radius  * dp.radius  * dp.axis_ratio);
+        //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+        sum +=weights_length[i].weight * weights_radius[k].weight
+            * weights_ratio[l].weight *DiamCyl(dp.length,suf_rad)/2.0;
+        norm += weights_length[i].weight* weights_radius[k].weight* weights_ratio[l].weight;
+      }
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+    rad_out = DiamCyl(dp.length,suf_rad)/2.0;}
+
+  return rad_out;
+}
Index: sansmodels/src/c_models/flexiblecylinder.cpp
===================================================================
--- sansmodels/src/c_models/flexiblecylinder.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/flexiblecylinder.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -18,19 +18,27 @@
  *   sansmodels/src/libigor
  *
- *	TODO: refactor so that we pull in the old sansmodels.c_extensions
  *	TODO: add 2d
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "flexible_cylinder.h"
 
 extern "C" {
 	#include "libCylinder.h"
 	#include "libStructureFactor.h"
-	#include "flexible_cylinder.h"
 }
+
+typedef struct {
+  double scale;
+  double length;
+  double kuhn_length;
+  double radius;
+  double sldCyl;
+  double sldSolv;
+  double background;
+} FlexibleCylinderParameters;
 
 FlexibleCylinderModel :: FlexibleCylinderModel() {
Index: sansmodels/src/c_models/fractal.cpp
===================================================================
--- sansmodels/src/c_models/fractal.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/fractal.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -21,23 +21,22 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "fractal.h"
 
 extern "C" {
-	#include "libTwoPhase.h"
-	#include "fractal.h"
+#include "libTwoPhase.h"
 }
 
 FractalModel :: FractalModel() {
-	scale      = Parameter(1.0);
-    radius = Parameter(5.0, true);
-    radius.set_min(0.0);
-    fractal_dim = Parameter(2.0);
-    cor_length = Parameter(100.0);
-    sldBlock = Parameter(2.0e-6);
-    sldSolv= Parameter(6.35e-6);
-	background = Parameter(0.0);
+  scale      = Parameter(1.0);
+  radius = Parameter(5.0, true);
+  radius.set_min(0.0);
+  fractal_dim = Parameter(2.0);
+  cor_length = Parameter(100.0);
+  sldBlock = Parameter(2.0e-6);
+  sldSolv= Parameter(6.35e-6);
+  background = Parameter(0.0);
 }
 
@@ -49,43 +48,43 @@
  */
 double FractalModel :: operator()(double q) {
-	double dp[7];
+  double dp[7];
 
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = scale();
-	dp[1] = radius();
-	dp[2] = fractal_dim();
-	dp[3] = cor_length();
-	dp[4] = sldBlock();
-	dp[5] = sldSolv();
-	dp[6] = 0.0;
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = radius();
+  dp[2] = fractal_dim();
+  dp[3] = cor_length();
+  dp[4] = sldBlock();
+  dp[5] = sldSolv();
+  dp[6] = 0.0;
 
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_rad;
-	radius.get_weights(weights_rad);
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad;
+  radius.get_weights(weights_rad);
 
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
 
-	// Loop over radius weight points
-	for(size_t i=0; i<weights_rad.size(); i++) {
-		dp[1] = weights_rad[i].value;
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp[1] = weights_rad[i].value;
 
-		//Un-normalize Fractal by volume
-		sum += weights_rad[i].weight
-			* Fractal(dp, fabs(q)) * pow(weights_rad[i].value,3);
-		//Find average volume
-		vol += weights_rad[i].weight
-			* pow(weights_rad[i].value,3);
+    //Un-normalize Fractal by volume
+    sum += weights_rad[i].weight
+        * Fractal(dp, fabs(q)) * pow(weights_rad[i].value,3);
+    //Find average volume
+    vol += weights_rad[i].weight
+        * pow(weights_rad[i].value,3);
 
-		norm += weights_rad[i].weight;
-	}
+    norm += weights_rad[i].weight;
+  }
 
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
-	return sum/norm + background();
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+  return sum/norm + background();
 }
 
@@ -97,6 +96,6 @@
  */
 double FractalModel :: operator()(double qx, double qy) {
-	double q = sqrt(qx*qx + qy*qy);
-	return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -109,5 +108,5 @@
  */
 double FractalModel :: evaluate_rphi(double q, double phi) {
-	return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 
@@ -117,5 +116,5 @@
  */
 double FractalModel :: calculate_ER() {
-	//NOT implemented yet!!! 'cause None shape Model
-	return 0.0;
+  //NOT implemented yet!!! 'cause None shape Model
+  return 0.0;
 }
Index: sansmodels/src/c_models/hollowcylinder.cpp
===================================================================
--- sansmodels/src/c_models/hollowcylinder.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/hollowcylinder.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -18,32 +18,117 @@
  *   sansmodels/src/libigor
  *
- *	TODO: refactor so that we pull in the old sansmodels.c_extensions
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "hollow_cylinder.h"
 
 extern "C" {
-	#include "libCylinder.h"
-	#include "libStructureFactor.h"
-	#include "hollow_cylinder.h"
+#include "libCylinder.h"
+#include "libStructureFactor.h"
+}
+
+typedef struct {
+  double scale;
+  double core_radius;
+  double radius;
+  double length;
+  double sldCyl;
+  double sldSolv;
+  double background;
+  double axis_theta;
+  double axis_phi;
+} HollowCylinderParameters;
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the hollow cylinder
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double hollow_cylinder_analytical_2D_scaled(HollowCylinderParameters *pars, double q, double q_x, double q_y) {
+  double cyl_x, cyl_y, cyl_z;
+  double q_z;
+  double  alpha,vol, cos_val;
+  double answer;
+  //convert angle degree to radian
+  double pi = 4.0*atan(1.0);
+  double theta = pars->axis_theta * pi/180.0;
+  double phi = pars->axis_phi * pi/180.0;
+
+  // Cylinder orientation
+  cyl_x = sin(theta) * cos(phi);
+  cyl_y = sin(theta) * sin(phi);
+  cyl_z = cos(theta);
+
+  // q vector
+  q_z = 0;
+
+  // Compute the angle btw vector q and the
+  // axis of the cylinder
+  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
+
+  // The following test should always pass
+  if (fabs(cos_val)>1.0) {
+    printf("core_shell_cylinder_analytical_2D: Unexpected error: cos(alpha)=%g\n", cos_val);
+    return 0;
+  }
+
+  alpha = acos( cos_val );
+
+  // Call the IGOR library function to get the kernel
+  answer = HolCylKernel(q, pars->core_radius, pars->radius, pars->length, cos_val);
+
+  // Multiply by contrast^2
+  answer *= (pars->sldCyl - pars->sldSolv)*(pars->sldCyl - pars->sldSolv);
+
+  //normalize by cylinder volume
+  vol=pi*((pars->radius*pars->radius)-(pars->core_radius *pars->core_radius))
+          *(pars->length);
+  answer *= vol;
+
+  //convert to [cm-1]
+  answer *= 1.0e8;
+
+  //Scale
+  answer *= pars->scale;
+
+  // add in the background
+  answer += pars->background;
+
+  return answer;
+}
+
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the Hollow cylinder
+ * @param q: q-value
+ * @return: function value
+ */
+static double hollow_cylinder_analytical_2DXY(HollowCylinderParameters *pars, double qx, double qy) {
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return hollow_cylinder_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
 HollowCylinderModel :: HollowCylinderModel() {
-	scale      = Parameter(1.0);
-	core_radius = Parameter(20.0, true);
-	core_radius.set_min(0.0);
-	radius  = Parameter(30.0, true);
-	radius.set_min(0.0);
-	length     = Parameter(400.0, true);
-	length.set_min(0.0);
-	sldCyl  = Parameter(6.3e-6);
-	sldSolv  = Parameter(1.0e-6);
-	background = Parameter(0.0);
-	axis_theta = Parameter(0.0, true);
-	axis_phi   = Parameter(0.0, true);
+  scale      = Parameter(1.0);
+  core_radius = Parameter(20.0, true);
+  core_radius.set_min(0.0);
+  radius  = Parameter(30.0, true);
+  radius.set_min(0.0);
+  length     = Parameter(400.0, true);
+  length.set_min(0.0);
+  sldCyl  = Parameter(6.3e-6);
+  sldSolv  = Parameter(1.0e-6);
+  background = Parameter(0.0);
+  axis_theta = Parameter(0.0, true);
+  axis_phi   = Parameter(0.0, true);
 }
 
@@ -55,67 +140,67 @@
  */
 double HollowCylinderModel :: operator()(double q) {
-	double dp[7];
-
-	dp[0] = scale();
-	dp[1] = core_radius();
-	dp[2] = radius();
-	dp[3] = length();
-	dp[4] = sldCyl();
-	dp[5] = sldSolv();
-	dp[6] = 0.0;
-
-	// Get the dispersion points for the core radius
-	vector<WeightPoint> weights_core_radius;
-	core_radius.get_weights(weights_core_radius);
-
-	// Get the dispersion points for the shell radius
-	vector<WeightPoint> weights_radius;
-	radius.get_weights(weights_radius);
-
-	// Get the dispersion points for the length
-	vector<WeightPoint> weights_length;
-	length.get_weights(weights_length);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
-
-	// Loop over core radius weight points
-	for(int i=0; i< (int)weights_core_radius.size(); i++) {
-		dp[1] = weights_core_radius[i].value;
-
-		// Loop over length weight points
-		for(int j=0; j< (int)weights_length.size(); j++) {
-			dp[3] = weights_length[j].value;
-
-			// Loop over shell radius weight points
-			for(int k=0; k< (int)weights_radius.size(); k++) {
-				dp[2] = weights_radius[k].value;
-				//Un-normalize  by volume
-				sum += weights_core_radius[i].weight
-					* weights_length[j].weight
-					* weights_radius[k].weight
-					* HollowCylinder(dp, q)
-					* (pow(weights_radius[k].value,2)-pow(weights_core_radius[i].value,2))
-					* weights_length[j].value;
-				//Find average volume
-				vol += weights_core_radius[i].weight
-					* weights_length[j].weight
-					* weights_radius[k].weight
-					* (pow(weights_radius[k].value,2)-pow(weights_core_radius[i].value,2))
-					* weights_length[j].value;
-
-				norm += weights_core_radius[i].weight
-					* weights_length[j].weight
-					* weights_radius[k].weight;
-			}
-		}
-	}
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
-
-	return sum/norm + background();
+  double dp[7];
+
+  dp[0] = scale();
+  dp[1] = core_radius();
+  dp[2] = radius();
+  dp[3] = length();
+  dp[4] = sldCyl();
+  dp[5] = sldSolv();
+  dp[6] = 0.0;
+
+  // Get the dispersion points for the core radius
+  vector<WeightPoint> weights_core_radius;
+  core_radius.get_weights(weights_core_radius);
+
+  // Get the dispersion points for the shell radius
+  vector<WeightPoint> weights_radius;
+  radius.get_weights(weights_radius);
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_length;
+  length.get_weights(weights_length);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over core radius weight points
+  for(int i=0; i< (int)weights_core_radius.size(); i++) {
+    dp[1] = weights_core_radius[i].value;
+
+    // Loop over length weight points
+    for(int j=0; j< (int)weights_length.size(); j++) {
+      dp[3] = weights_length[j].value;
+
+      // Loop over shell radius weight points
+      for(int k=0; k< (int)weights_radius.size(); k++) {
+        dp[2] = weights_radius[k].value;
+        //Un-normalize  by volume
+        sum += weights_core_radius[i].weight
+            * weights_length[j].weight
+            * weights_radius[k].weight
+            * HollowCylinder(dp, q)
+        * (pow(weights_radius[k].value,2)-pow(weights_core_radius[i].value,2))
+        * weights_length[j].value;
+        //Find average volume
+        vol += weights_core_radius[i].weight
+            * weights_length[j].weight
+            * weights_radius[k].weight
+            * (pow(weights_radius[k].value,2)-pow(weights_core_radius[i].value,2))
+            * weights_length[j].value;
+
+        norm += weights_core_radius[i].weight
+            * weights_length[j].weight
+            * weights_radius[k].weight;
+      }
+    }
+  }
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -127,105 +212,105 @@
  */
 double HollowCylinderModel :: operator()(double qx, double qy) {
-	HollowCylinderParameters dp;
-	// Fill parameter array
-	dp.scale      = scale();
-	dp.core_radius     = core_radius();
-	dp.radius  = radius();
-	dp.length     = length();
-	dp.sldCyl   = sldCyl();
-	dp.sldSolv  = sldSolv();
-	dp.background = 0.0;
-	dp.axis_theta = axis_theta();
-	dp.axis_phi   = axis_phi();
-
-	// Get the dispersion points for the core radius
-	vector<WeightPoint> weights_core_radius;
-	core_radius.get_weights(weights_core_radius);
-
-	// Get the dispersion points for the shell radius
-	vector<WeightPoint> weights_radius;
-	radius.get_weights(weights_radius);
-
-	// Get the dispersion points for the length
-	vector<WeightPoint> weights_length;
-	length.get_weights(weights_length);
-
-	// Get angular averaging for theta
-	vector<WeightPoint> weights_theta;
-	axis_theta.get_weights(weights_theta);
-
-	// Get angular averaging for phi
-	vector<WeightPoint> weights_phi;
-	axis_phi.get_weights(weights_phi);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double norm_vol = 0.0;
-	double vol = 0.0;
-	double pi = 4.0*atan(1.0);
-	// Loop over core radius weight points
-	for(int i=0; i<(int)weights_core_radius.size(); i++) {
-		dp.core_radius = weights_core_radius[i].value;
-
-
-		// Loop over length weight points
-		for(int j=0; j<(int)weights_length.size(); j++) {
-			dp.length = weights_length[j].value;
-
-			// Loop over shell radius weight points
-			for(int m=0; m< (int)weights_radius.size(); m++) {
-				dp.radius = weights_radius[m].value;
-
-			// Average over theta distribution
-			for(int k=0; k< (int)weights_theta.size(); k++) {
-				dp.axis_theta = weights_theta[k].value;
-
-				// Average over phi distribution
-				for(int l=0; l< (int)weights_phi.size(); l++) {
-					dp.axis_phi = weights_phi[l].value;
-					//Un-normalize by volume
-					double _ptvalue = weights_core_radius[i].weight
-						* weights_length[j].weight
-						* weights_radius[m].weight
-						* weights_theta[k].weight
-						* weights_phi[l].weight
-						* hollow_cylinder_analytical_2DXY(&dp, qx, qy)
-						/ ((pow(weights_radius[m].value,2)-pow(weights_core_radius[i].value,2))
-						* weights_length[j].value);
-					if (weights_theta.size()>1) {
-						_ptvalue *= fabs(sin(weights_theta[k].value * pi/180.0));
-					}
-					sum += _ptvalue;
-					//Find average volume
-					vol += weights_core_radius[i].weight
-						* weights_length[j].weight
-						* weights_radius[m].weight
-						* (pow(weights_radius[m].value,2)-pow(weights_core_radius[i].value,2))
-						* weights_length[j].value;
-					//Find norm for volume
-					norm_vol += weights_core_radius[i].weight
-						* weights_length[j].weight
-						* weights_radius[m].weight;
-
-					norm += weights_core_radius[i].weight
-						* weights_length[j].weight
-						* weights_radius[m].weight
-						* weights_theta[k].weight
-						* weights_phi[l].weight;
-
-				}
-			}
-			}
-		}
-	}
-	// Averaging in theta needs an extra normalization
-	// factor to account for the sin(theta) term in the
-	// integration (see documentation).
-	if (weights_theta.size()>1) norm = norm/asin(1.0);
-	if (vol != 0.0 || norm_vol != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum*(vol/norm_vol);}
-	return sum/norm + background();
+  HollowCylinderParameters dp;
+  // Fill parameter array
+  dp.scale      = scale();
+  dp.core_radius     = core_radius();
+  dp.radius  = radius();
+  dp.length     = length();
+  dp.sldCyl   = sldCyl();
+  dp.sldSolv  = sldSolv();
+  dp.background = 0.0;
+  dp.axis_theta = axis_theta();
+  dp.axis_phi   = axis_phi();
+
+  // Get the dispersion points for the core radius
+  vector<WeightPoint> weights_core_radius;
+  core_radius.get_weights(weights_core_radius);
+
+  // Get the dispersion points for the shell radius
+  vector<WeightPoint> weights_radius;
+  radius.get_weights(weights_radius);
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_length;
+  length.get_weights(weights_length);
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_theta;
+  axis_theta.get_weights(weights_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_phi;
+  axis_phi.get_weights(weights_phi);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double norm_vol = 0.0;
+  double vol = 0.0;
+  double pi = 4.0*atan(1.0);
+  // Loop over core radius weight points
+  for(int i=0; i<(int)weights_core_radius.size(); i++) {
+    dp.core_radius = weights_core_radius[i].value;
+
+
+    // Loop over length weight points
+    for(int j=0; j<(int)weights_length.size(); j++) {
+      dp.length = weights_length[j].value;
+
+      // Loop over shell radius weight points
+      for(int m=0; m< (int)weights_radius.size(); m++) {
+        dp.radius = weights_radius[m].value;
+
+        // Average over theta distribution
+        for(int k=0; k< (int)weights_theta.size(); k++) {
+          dp.axis_theta = weights_theta[k].value;
+
+          // Average over phi distribution
+          for(int l=0; l< (int)weights_phi.size(); l++) {
+            dp.axis_phi = weights_phi[l].value;
+            //Un-normalize by volume
+            double _ptvalue = weights_core_radius[i].weight
+                * weights_length[j].weight
+                * weights_radius[m].weight
+                * weights_theta[k].weight
+                * weights_phi[l].weight
+                * hollow_cylinder_analytical_2DXY(&dp, qx, qy)
+            / ((pow(weights_radius[m].value,2)-pow(weights_core_radius[i].value,2))
+                * weights_length[j].value);
+            if (weights_theta.size()>1) {
+              _ptvalue *= fabs(sin(weights_theta[k].value * pi/180.0));
+            }
+            sum += _ptvalue;
+            //Find average volume
+            vol += weights_core_radius[i].weight
+                * weights_length[j].weight
+                * weights_radius[m].weight
+                * (pow(weights_radius[m].value,2)-pow(weights_core_radius[i].value,2))
+                * weights_length[j].value;
+            //Find norm for volume
+            norm_vol += weights_core_radius[i].weight
+                * weights_length[j].weight
+                * weights_radius[m].weight;
+
+            norm += weights_core_radius[i].weight
+                * weights_length[j].weight
+                * weights_radius[m].weight
+                * weights_theta[k].weight
+                * weights_phi[l].weight;
+
+          }
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_theta.size()>1) norm = norm/asin(1.0);
+  if (vol != 0.0 || norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum*(vol/norm_vol);}
+  return sum/norm + background();
 }
 
@@ -238,7 +323,7 @@
  */
 double HollowCylinderModel :: evaluate_rphi(double q, double phi) {
-	double qx = q*cos(phi);
-	double qy = q*sin(phi);
-	return (*this).operator()(qx, qy);
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
 }
 /**
@@ -247,42 +332,42 @@
  */
 double HollowCylinderModel :: calculate_ER() {
-	HollowCylinderParameters dp;
-
-	dp.radius  = radius();
-	dp.length     = length();
-
-	double rad_out = 0.0;
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-
-	// Get the dispersion points for the major shell
-	vector<WeightPoint> weights_length;
-	length.get_weights(weights_length);
-
-	// Get the dispersion points for the minor shell
-	vector<WeightPoint> weights_radius ;
-	radius.get_weights(weights_radius);
-
-	// Loop over major shell weight points
-	for(int i=0; i< (int)weights_length.size(); i++) {
-		dp.length = weights_length[i].value;
-		for(int k=0; k< (int)weights_radius.size(); k++) {
-			dp.radius = weights_radius[k].value;
-			//Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-			sum +=weights_length[i].weight
-				* weights_radius[k].weight*DiamCyl(dp.length,dp.radius)/2.0;
-			norm += weights_length[i].weight* weights_radius[k].weight;
-		}
-	}
-	if (norm != 0){
-		//return the averaged value
-		rad_out =  sum/norm;}
-	else{
-		//return normal value
-		//Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-		rad_out = DiamCyl(dp.length,dp.radius)/2.0;}
-
-	return rad_out;
-}
+  HollowCylinderParameters dp;
+
+  dp.radius  = radius();
+  dp.length     = length();
+
+  double rad_out = 0.0;
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the major shell
+  vector<WeightPoint> weights_length;
+  length.get_weights(weights_length);
+
+  // Get the dispersion points for the minor shell
+  vector<WeightPoint> weights_radius ;
+  radius.get_weights(weights_radius);
+
+  // Loop over major shell weight points
+  for(int i=0; i< (int)weights_length.size(); i++) {
+    dp.length = weights_length[i].value;
+    for(int k=0; k< (int)weights_radius.size(); k++) {
+      dp.radius = weights_radius[k].value;
+      //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+      sum +=weights_length[i].weight
+          * weights_radius[k].weight*DiamCyl(dp.length,dp.radius)/2.0;
+      norm += weights_length[i].weight* weights_radius[k].weight;
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+    rad_out = DiamCyl(dp.length,dp.radius)/2.0;}
+
+  return rad_out;
+}
Index: sansmodels/src/c_models/lamellar.cpp
===================================================================
--- sansmodels/src/c_models/lamellar.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/lamellar.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -18,16 +18,39 @@
  *   sansmodels/src/libigor
  *
- *	TODO: refactor so that we pull in the old sansmodels.c_extensions
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "lamellar.h"
 
-extern "C" {
-//	#include "libCylinder.h"
-	#include "lamellar.h"
+/*  LamellarFFX  :  calculates the form factor of a lamellar structure - no S(q) effects included
+            -NO polydispersion included
+*/
+static double lamellar_kernel(double dp[], double q){
+  double scale,del,sld_bi,sld_sol,contr,bkg;    //local variables of coefficient wave
+  double inten, qval,Pq;
+  double Pi;
+
+
+  Pi = 4.0*atan(1.0);
+  scale = dp[0];
+  del = dp[1];
+  sld_bi = dp[2];
+  sld_sol = dp[3];
+  bkg = dp[4];
+  qval = q;
+  contr = sld_bi -sld_sol;
+
+  Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del));
+
+  inten = 2.0*Pi*scale*Pq/(qval*qval);    //this is now dimensionless...
+
+  inten /= del;     //normalize by the thickness (in A)
+
+  inten *= 1.0e8;   // 1/A to 1/cm
+
+  return(inten+bkg);
 }
 
@@ -39,5 +62,4 @@
 	sld_sol    = Parameter(6.3e-6);
 	background = Parameter(0.0);
-
 }
 
Index: sansmodels/src/c_models/lamellarFF_HG.cpp
===================================================================
--- sansmodels/src/c_models/lamellarFF_HG.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/lamellarFF_HG.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -18,28 +18,26 @@
  *   sansmodels/src/libigor
  *
- *	TODO: refactor so that we pull in the old sansmodels.c_extensions
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "lamellarFF_HG.h"
 
 extern "C" {
-	#include "libCylinder.h"
-	#include "lamellarFF_HG.h"
+#include "libCylinder.h"
 }
 
 LamellarFFHGModel :: LamellarFFHGModel() {
-	scale      = Parameter(1.0);
-	t_length     = Parameter(15.0, true);
-	t_length.set_min(0.0);
-	h_thickness    = Parameter(10.0, true);
-	h_thickness.set_min(0.0);
-	sld_tail   = Parameter(4e-7);
-	sld_head  = Parameter(3e-6);
-	sld_solvent    = Parameter(6e-6);
-	background = Parameter(0.0);
+  scale      = Parameter(1.0);
+  t_length     = Parameter(15.0, true);
+  t_length.set_min(0.0);
+  h_thickness    = Parameter(10.0, true);
+  h_thickness.set_min(0.0);
+  sld_tail   = Parameter(4e-7);
+  sld_head  = Parameter(3e-6);
+  sld_solvent    = Parameter(6e-6);
+  background = Parameter(0.0);
 
 }
@@ -52,41 +50,41 @@
  */
 double LamellarFFHGModel :: operator()(double q) {
-	double dp[7];
+  double dp[7];
 
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = scale();
-	dp[1] = t_length();
-	dp[2] = h_thickness();
-	dp[3] = sld_tail();
-	dp[4] = sld_head();
-	dp[5] = sld_solvent();
-	dp[6] = 0.0;
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = t_length();
+  dp[2] = h_thickness();
+  dp[3] = sld_tail();
+  dp[4] = sld_head();
+  dp[5] = sld_solvent();
+  dp[6] = 0.0;
 
-	// Get the dispersion points for the tail length
-	vector<WeightPoint> weights_t_length;
-	t_length.get_weights(weights_t_length);
+  // Get the dispersion points for the tail length
+  vector<WeightPoint> weights_t_length;
+  t_length.get_weights(weights_t_length);
 
-	// Get the dispersion points for the head thickness
-	vector<WeightPoint> weights_h_thickness;
-	h_thickness.get_weights(weights_h_thickness);
+  // Get the dispersion points for the head thickness
+  vector<WeightPoint> weights_h_thickness;
+  h_thickness.get_weights(weights_h_thickness);
 
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
 
-	// Loop over semi axis A weight points
-	for(int i=0; i< (int)weights_t_length.size(); i++) {
-		dp[1] = weights_t_length[i].value;
+  // Loop over semi axis A weight points
+  for(int i=0; i< (int)weights_t_length.size(); i++) {
+    dp[1] = weights_t_length[i].value;
 
-		for (int j=0; j< (int)weights_h_thickness.size(); j++){
-			dp[2] = weights_h_thickness[j].value;
+    for (int j=0; j< (int)weights_h_thickness.size(); j++){
+      dp[2] = weights_h_thickness[j].value;
 
-			sum += weights_t_length[i].weight* weights_h_thickness[j].weight*LamellarFF_HG(dp, q);
-			norm += weights_t_length[i].weight* weights_h_thickness[j].weight;
-		}
+      sum += weights_t_length[i].weight* weights_h_thickness[j].weight*LamellarFF_HG(dp, q);
+      norm += weights_t_length[i].weight* weights_h_thickness[j].weight;
+    }
 
-	}
-	return sum/norm + background();
+  }
+  return sum/norm + background();
 }
 
@@ -99,6 +97,6 @@
 
 double LamellarFFHGModel :: operator()(double qx, double qy) {
-	double q = sqrt(qx*qx + qy*qy);
-	return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -111,5 +109,5 @@
  */
 double LamellarFFHGModel :: evaluate_rphi(double q, double phi) {
-	return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 /**
@@ -118,5 +116,5 @@
  */
 double LamellarFFHGModel :: calculate_ER() {
-//NOT implemented yet!!!
-	return 0.0;
+  //NOT implemented yet!!!
+  return 0.0;
 }
Index: sansmodels/src/c_models/lamellarPC.cpp
===================================================================
--- sansmodels/src/c_models/lamellarPC.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/lamellarPC.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -20,25 +20,24 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "lamellarPC.h"
 
 extern "C" {
-	#include "libCylinder.h"
-	#include "lamellarPC.h"
+#include "libCylinder.h"
 }
 
 LamellarPCrystalModel :: LamellarPCrystalModel() {
-	scale      = Parameter(1.0);
-	thickness     = Parameter(33.0, true);
-	thickness.set_min(0.0);
-	Nlayers    = Parameter(20.0, true);
-	Nlayers.set_min(0.0);
-	spacing   = Parameter(250);
-	pd_spacing   = Parameter(0.0);
-	sld_layer  = Parameter(1.0e-6);
-	sld_solvent    = Parameter(6.34e-6);
-	background = Parameter(0.0);
+  scale      = Parameter(1.0);
+  thickness     = Parameter(33.0, true);
+  thickness.set_min(0.0);
+  Nlayers    = Parameter(20.0, true);
+  Nlayers.set_min(0.0);
+  spacing   = Parameter(250);
+  pd_spacing   = Parameter(0.0);
+  sld_layer  = Parameter(1.0e-6);
+  sld_solvent    = Parameter(6.34e-6);
+  background = Parameter(0.0);
 
 }
@@ -51,41 +50,41 @@
  */
 double LamellarPCrystalModel :: operator()(double q) {
-	double dp[8];
+  double dp[8];
 
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = scale();
-	dp[1] = thickness();
-	dp[2] = Nlayers();
-	dp[3] = spacing();
-	dp[4] = pd_spacing();
-	dp[5] = sld_layer();
-	dp[6] = sld_solvent();
-	dp[7] = 0.0; // Do not apply background here.
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = thickness();
+  dp[2] = Nlayers();
+  dp[3] = spacing();
+  dp[4] = pd_spacing();
+  dp[5] = sld_layer();
+  dp[6] = sld_solvent();
+  dp[7] = 0.0; // Do not apply background here.
 
-	// Get the dispersion points for the head thickness
-	vector<WeightPoint> weights_thickness;
-	thickness.get_weights(weights_thickness);
+  // Get the dispersion points for the head thickness
+  vector<WeightPoint> weights_thickness;
+  thickness.get_weights(weights_thickness);
 
-	// Let's provide from the func which is more accurate especially for small q region.
-	// Get the dispersion points for the tail length
-	//vector<WeightPoint> weights_spacing;
-	//spacing.get_weights(weights_spacing);
+  // Let's provide from the func which is more accurate especially for small q region.
+  // Get the dispersion points for the tail length
+  //vector<WeightPoint> weights_spacing;
+  //spacing.get_weights(weights_spacing);
 
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
 
-	// Loop over thickness and spacing weight points
-	for(int i=0; i< (int)weights_thickness.size(); i++) {
-		dp[1] = weights_thickness[i].value;
-		//for (int j=0; j< (int)weights_spacing.size(); j++){
-			//dp[3] = weights_spacing[j].value;
-			sum += weights_thickness[i].weight*Lamellar_ParaCrystal(dp, q);
-			norm += weights_thickness[i].weight;
-		//}
-	}
-	//apply norm and background
-	return sum/norm + background();
+  // Loop over thickness and spacing weight points
+  for(int i=0; i< (int)weights_thickness.size(); i++) {
+    dp[1] = weights_thickness[i].value;
+    //for (int j=0; j< (int)weights_spacing.size(); j++){
+    //dp[3] = weights_spacing[j].value;
+    sum += weights_thickness[i].weight*Lamellar_ParaCrystal(dp, q);
+    norm += weights_thickness[i].weight;
+    //}
+  }
+  //apply norm and background
+  return sum/norm + background();
 }
 
@@ -98,6 +97,6 @@
 
 double LamellarPCrystalModel :: operator()(double qx, double qy) {
-	double q = sqrt(qx*qx + qy*qy);
-	return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -110,5 +109,5 @@
  */
 double LamellarPCrystalModel :: evaluate_rphi(double q, double phi) {
-	return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 /**
@@ -117,5 +116,5 @@
  */
 double LamellarPCrystalModel :: calculate_ER() {
-//NOT implemented yet!!!
-	return 0.0;
+  //NOT implemented yet!!!
+  return 0.0;
 }
Index: sansmodels/src/c_models/lamellarPS.cpp
===================================================================
--- sansmodels/src/c_models/lamellarPS.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/lamellarPS.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -23,25 +23,88 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
+#include <stdlib.h>
 using namespace std;
+#include "lamellarPS.h"
 
 extern "C" {
-	#include "libCylinder.h"
-	#include "lamellarPS.h"
+#include "libCylinder.h"
+}
+
+/*LamellarPS_kernel() was moved from libigor to get rid of polydipersity in del(thickness) that we provide from control panel.
+  LamellarPSX  :  calculates the form factor of a lamellar structure - with S(q) effects included
+-------
+------- resolution effects ARE NOT included, but only a CONSTANT default value, not the real q-dependent resolution!!
+
+ */
+static double LamellarPS_kernel(double dp[], double q)
+{
+  double scale,dd,del,sld_bi,sld_sol,contr,NN,Cp,bkg;   //local variables of coefficient wave
+  double inten, qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ;
+  double Pi,Euler,dQDefault,fii;
+  int ii,NNint;
+  Euler = 0.5772156649;   // Euler's constant
+  dQDefault = 0.0;    //[=] 1/A, q-resolution, default value
+  dQ = dQDefault;
+
+  Pi = 4.0*atan(1.0);
+  qval = q;
+
+  scale = dp[0];
+  dd = dp[1];
+  del = dp[2];
+  sld_bi = dp[3];
+  sld_sol = dp[4];
+  NN = trunc(dp[5]);    //be sure that NN is an integer
+  Cp = dp[6];
+  bkg = dp[7];
+
+  contr = sld_bi - sld_sol;
+
+  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;
+  for(ii=1;ii<(NNint-1);ii+=1) {
+
+    fii = (double)ii;   //do I really need to do this?
+
+    temp = 0.0;
+    alpha = Cp/4.0/Pi/Pi*(log(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 *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
+    temp /= sqrt(1.0+t1);
+
+    Sq += temp;
+  }
+
+  Sq *= 2.0;
+  Sq += 1.0;
+
+  inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval);
+
+  inten *= 1.0e8;   // 1/A to 1/cm
+
+  return(inten+bkg);
 }
 
 LamellarPSModel :: LamellarPSModel() {
-	scale      = Parameter(1.0);
-	spacing    = Parameter(400.0, true);
-	spacing.set_min(0.0);
-	delta     = Parameter(30.0);
-	delta.set_min(0.0);
-	sld_bi   = Parameter(6.3e-6);
-	sld_sol   = Parameter(1.0e-6);
-	n_plates     = Parameter(20.0);
-	caille = Parameter(0.1);
-	background = Parameter(0.0);
+  scale      = Parameter(1.0);
+  spacing    = Parameter(400.0, true);
+  spacing.set_min(0.0);
+  delta     = Parameter(30.0);
+  delta.set_min(0.0);
+  sld_bi   = Parameter(6.3e-6);
+  sld_sol   = Parameter(1.0e-6);
+  n_plates     = Parameter(20.0);
+  caille = Parameter(0.1);
+  background = Parameter(0.0);
 
 }
@@ -54,39 +117,39 @@
  */
 double LamellarPSModel :: operator()(double q) {
-	double dp[8];
+  double dp[8];
 
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = scale();
-	dp[1] = spacing();
-	dp[2] = delta();
-	dp[3] = sld_bi();
-	dp[4] = sld_sol();
-	dp[5] = n_plates();
-	dp[6] = caille();
-	dp[7] = 0.0;
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = spacing();
+  dp[2] = delta();
+  dp[3] = sld_bi();
+  dp[4] = sld_sol();
+  dp[5] = n_plates();
+  dp[6] = caille();
+  dp[7] = 0.0;
 
 
-	// Get the dispersion points for spacing and delta (thickness)
-	vector<WeightPoint> weights_spacing;
-	spacing.get_weights(weights_spacing);
-	vector<WeightPoint> weights_delta;
-	delta.get_weights(weights_delta);
+  // Get the dispersion points for spacing and delta (thickness)
+  vector<WeightPoint> weights_spacing;
+  spacing.get_weights(weights_spacing);
+  vector<WeightPoint> weights_delta;
+  delta.get_weights(weights_delta);
 
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
 
-	// Loop over short_edgeA weight points
-	for(int i=0; i< (int)weights_spacing.size(); i++) {
-		dp[1] = weights_spacing[i].value;
-		for(int j=0; j< (int)weights_spacing.size(); j++) {
-			dp[2] = weights_delta[i].value;
+  // Loop over short_edgeA weight points
+  for(int i=0; i< (int)weights_spacing.size(); i++) {
+    dp[1] = weights_spacing[i].value;
+    for(int j=0; j< (int)weights_spacing.size(); j++) {
+      dp[2] = weights_delta[i].value;
 
-			sum += weights_spacing[i].weight * weights_delta[j].weight * LamellarPS_kernel(dp, q);
-			norm += weights_spacing[i].weight * weights_delta[j].weight;
-		}
-	}
-	return sum/norm + background();
+      sum += weights_spacing[i].weight * weights_delta[j].weight * LamellarPS_kernel(dp, q);
+      norm += weights_spacing[i].weight * weights_delta[j].weight;
+    }
+  }
+  return sum/norm + background();
 }
 /**
@@ -97,6 +160,6 @@
  */
 double LamellarPSModel :: operator()(double qx, double qy) {
-	double q = sqrt(qx*qx + qy*qy);
-	return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -109,5 +172,5 @@
  */
 double LamellarPSModel :: evaluate_rphi(double q, double phi) {
-	return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 /**
@@ -116,6 +179,6 @@
  */
 double LamellarPSModel :: calculate_ER() {
-//NOT implemented yet!!!
-	return 0.0;
+  //NOT implemented yet!!!
+  return 0.0;
 }
 
Index: sansmodels/src/c_models/lamellarPS_HG.cpp
===================================================================
--- sansmodels/src/c_models/lamellarPS_HG.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/lamellarPS_HG.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -18,17 +18,15 @@
  *   sansmodels/src/libigor
  *
- *	TODO: refactor so that we pull in the old sansmodels.c_extensions
  *	TODO: add 2D function
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "lamellarPS_HG.h"
 
 extern "C" {
 	#include "libCylinder.h"
-	#include "lamellarPS_HG.h"
 }
 
@@ -137,47 +135,2 @@
 	return 0.0;
 }
-
-/*
-double LamellarPSHGModel :: operator()(double qx, double qy) {
-	LamellarPSHGParameters dp;
-	// Fill parameter array
-	dp.scale      = scale();
-	dp.spacing   = spacing();
-	dp.deltaT  = deltaT();
-	dp.deltaH = deltaH();
-	dp.sld_tail   = sld_tail();
-	dp.sld_head = sld_head();
-	dp.sld_solvent   = sld_solvent();
-	dp.n_plates = n_plates();
-	dp.caille = caille();
-	dp.background    = background();
-
-	// Get the dispersion points for the deltaT
-	vector<WeightPoint> weights_deltaT;
-	deltaT.get_weights(weights_deltaT);
-
-	// Get the dispersion points for the deltaH
-	vector<WeightPoint> weights_deltaH;
-	deltaH.get_weights(weights_deltaH);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-
-	// Loop over deltaT weight points
-	for(int i=0; i< (int)weights_deltaT.size(); i++) {
-		dp.deltaT = weights_deltaT[i].value;
-
-		// Loop over deltaH weight points
-		for(int j=0; j< (int)weights_deltaH.size(); j++) {
-			dp.deltaH = weights_deltaH[j].value;
-
-			sum += weights_deltaT[i].weight *weights_deltaH[j].weight *lamellarPS_HG_analytical_2DXY(&dp, qx, qy);
-			norm += weights_deltaT[i].weight * weights_deltaH[j].weight;
-		}
-	}
-	return sum/norm + background();
-}
-*/
-
-
Index: sansmodels/src/c_models/lorentzian.cpp
===================================================================
--- sansmodels/src/c_models/lorentzian.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
+++ sansmodels/src/c_models/lorentzian.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -0,0 +1,67 @@
+/**
+	This software was developed by the University of Tennessee as part of the
+	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+	project funded by the US National Science Foundation.
+
+	If you use DANSE applications to do scientific research that leads to
+	publication, we ask that you acknowledge the use of the software with the
+	following sentence:
+
+	"This work benefited from DANSE software developed under NSF award DMR-0520547."
+
+	copyright 2008, University of Tennessee
+ */
+
+#include <math.h>
+#include "parameters.hh"
+#include <stdio.h>
+using namespace std;
+#include "lorentzian.h"
+
+Lorentzian :: Lorentzian() {
+  scale  = Parameter(1.0, true);
+  gamma  = Parameter(1.0, true);
+  center = Parameter(0.0, true);
+}
+
+/**
+ * Function to evaluate 1D scattering function
+  * @param q: q-value
+ * @return: function value
+ */
+double Lorentzian :: operator()(double q) {
+  double delta2 = pow( (q-center()), 2);
+  return scale() / acos(-1.0) * 0.5*gamma() /
+    (delta2 + 0.25*gamma()*gamma());
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param q_x: value of Q along x
+ * @param q_y: value of Q along y
+ * @return: function value
+ */
+double Lorentzian :: operator()(double qx, double qy) {
+  return (*this).operator()(qx) * (*this).operator()(qy);
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double Lorentzian :: evaluate_rphi(double q, double phi) {
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
+}
+/**
+ * Function to calculate effective radius
+ * @return: effective radius value
+ */
+double Lorentzian :: calculate_ER() {
+  //NOT implemented yet!!!
+  return 0.0;
+}
Index: sansmodels/src/c_models/models.hh
===================================================================
--- sansmodels/src/c_models/models.hh	(revision 011e0e4569d018a5f39b70032285911f67750383)
+++ 	(revision )
@@ -1,431 +1,0 @@
-/**
-	This software was developed by the University of Tennessee as part of the
-	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
-	project funded by the US National Science Foundation.
-
-	If you use DANSE applications to do scientific research that leads to
-	publication, we ask that you acknowledge the use of the software with the
-	following sentence:
-
-	"This work benefited from DANSE software developed under NSF award DMR-0520547."
-
-	copyright 2008, University of Tennessee
- */
-#ifndef MODEL_CLASS_H
-#define MODEL_CLASS_H
-
-#include <vector>
-#include "parameters.hh"
-extern "C" {
-	#include "cylinder.h"
-	#include "parallelepiped.h"
-	#include "lamellarPS.h"
-	#include "lamellar.h"
-	#include "fuzzysphere.h"
-}
-
-using namespace std;
-
-class EllipticalCylinderModel{
-public:
-	// Model parameters
-	Parameter r_minor;
-	Parameter scale;
-	Parameter r_ratio;
-	Parameter length;
-	Parameter sldCyl;
-	Parameter sldSolv;
-	Parameter background;
-	Parameter cyl_theta;
-	Parameter cyl_phi;
-	Parameter cyl_psi;
-
-	// Constructor
-	EllipticalCylinderModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-class TriaxialEllipsoidModel{
-public:
-	// Model parameters
-	Parameter scale;
-	Parameter semi_axisA;
-	Parameter semi_axisB;
-	Parameter semi_axisC;
-	Parameter sldEll;
-	Parameter sldSolv;
-	Parameter background;
-	Parameter axis_theta;
-	Parameter axis_phi;
-	Parameter axis_psi;
-
-	// Constructor
-	TriaxialEllipsoidModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class FlexibleCylinderModel{
-public:
-	// Model parameters
-	Parameter scale;
-	Parameter length;
-	Parameter kuhn_length;
-	Parameter radius;
-	Parameter sldCyl;
-	Parameter sldSolv;
-	Parameter background;
-
-	// Constructor
-	FlexibleCylinderModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class FlexCylEllipXModel{
-public:
-	// Model parameters
-	Parameter scale;
-	Parameter length;
-	Parameter kuhn_length;
-	Parameter radius;
-	Parameter axis_ratio;
-	Parameter sldCyl;
-	Parameter sldSolv;
-	Parameter background;
-
-	// Constructor
-	FlexCylEllipXModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class StackedDisksModel{
-public:
-	// Model parameters
-	Parameter scale;
-	Parameter core_thick;
-	Parameter radius;
-	Parameter layer_thick;
-	Parameter core_sld;
-	Parameter layer_sld;
-	Parameter solvent_sld;
-	Parameter n_stacking;
-	Parameter sigma_d;
-	Parameter background;
-	Parameter axis_theta;
-	Parameter axis_phi;
-
-	// Constructor
-	StackedDisksModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class LamellarModel{
-public:
-	// Model parameters
-	Parameter scale;
-	Parameter bi_thick;
-	Parameter sld_bi;
-	Parameter sld_sol;
-	Parameter background;
-
-	// Constructor
-	LamellarModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-
-};
-
-class LamellarFFHGModel{
-public:
-	// Model parameters
-	Parameter scale;
-	Parameter t_length;
-	Parameter h_thickness;
-	Parameter sld_tail;
-	Parameter sld_head;
-	Parameter sld_solvent;
-	Parameter background;
-
-	// Constructor
-	LamellarFFHGModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-
-};
-
-
-
-class LamellarPSModel{
-public:
-	// Model parameters
-	Parameter scale;
-	Parameter spacing;
-	Parameter delta;
-	Parameter sld_bi;
-	Parameter sld_sol;
-	Parameter n_plates;
-	Parameter caille;
-	Parameter background;
-
-	// Constructor
-	LamellarPSModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class LamellarPSHGModel{
-public:
-	// Model parameters
-	Parameter scale;
-	Parameter spacing;
-	Parameter deltaT;
-	Parameter deltaH;
-	Parameter sld_tail;
-	Parameter sld_head;
-	Parameter sld_solvent;
-	Parameter n_plates;
-	Parameter caille;
-	Parameter background;
-
-	// Constructor
-	LamellarPSHGModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-
-class LamellarPCrystalModel{
-public:
-	// Model parameters
-	Parameter scale;
-	Parameter thickness;
-	Parameter Nlayers;
-	Parameter spacing;
-	Parameter pd_spacing;
-	Parameter sld_layer;
-	Parameter sld_solvent;
-	Parameter background;
-
-	// Constructor
-	LamellarPCrystalModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class CoreShellEllipsoidModel{
-public:
-	// Model parameters
-	Parameter scale;
-	Parameter equat_core;
-	Parameter polar_core;
-	Parameter equat_shell;
-	Parameter polar_shell;
-	Parameter sld_core;
-	Parameter sld_shell;
-	Parameter sld_solvent;
-	Parameter background;
-	Parameter axis_theta;
-	Parameter axis_phi;
-
-	// Constructor
-	CoreShellEllipsoidModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class HollowCylinderModel{
-public:
-	// Model parameters
-	Parameter scale;
-	Parameter core_radius;
-	Parameter radius;
-	Parameter length;
-	Parameter sldCyl;
-	Parameter sldSolv;
-	Parameter background;
-	Parameter axis_theta;
-	Parameter axis_phi;
-
-	//Constructor
-	HollowCylinderModel();
-
-	//Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx , double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class MultiShellModel{
-public:
-	// Model parameters
-	Parameter scale;
-	Parameter core_radius;
-	Parameter s_thickness;
-	Parameter w_thickness;
-	Parameter core_sld;
-	Parameter shell_sld;
-	Parameter n_pairs;
-	Parameter background;
-
-	//Constructor
-	MultiShellModel();
-
-	//Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx , double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class VesicleModel{
-public:
-	// Model parameters
-	Parameter scale;
-	Parameter radius;
-	Parameter thickness;
-	Parameter core_sld;
-	Parameter shell_sld;
-	Parameter background;
-
-	//Constructor
-	VesicleModel();
-
-	//Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx , double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class BinaryHSModel{
-public:
-	// Model parameters
-	Parameter l_radius;
-	Parameter s_radius;
-	Parameter vol_frac_ls;
-	Parameter vol_frac_ss;
-	Parameter ls_sld;
-	Parameter ss_sld;
-	Parameter solvent_sld;
-	Parameter background;
-
-	//Constructor
-	BinaryHSModel();
-
-	//Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx , double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class BinaryHSPSF11Model{
-public:
-	// Model parameters
-	Parameter l_radius;
-	Parameter s_radius;
-	Parameter vol_frac_ls;
-	Parameter vol_frac_ss;
-	Parameter ls_sld;
-	Parameter ss_sld;
-	Parameter solvent_sld;
-	Parameter background;
-
-	//Constructor
-	BinaryHSPSF11Model();
-
-	//Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx , double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class Poly_GaussCoil{
-public:
-	// Model parameters
-	Parameter rg;
-	Parameter scale;
-	Parameter poly_m;
-	Parameter background;
-
-	// Constructor
-	Poly_GaussCoil();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class FractalModel{
-public:
-	// Model parameters
-	Parameter radius;
-	Parameter scale;
-	Parameter fractal_dim;
-	Parameter cor_length;
-	Parameter sldBlock;
-	Parameter sldSolv;
-	Parameter background;
-
-	// Constructor
-	FractalModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-#endif
Index: sansmodels/src/c_models/multishell.cpp
===================================================================
--- sansmodels/src/c_models/multishell.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/multishell.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -21,26 +21,37 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "multishell.h"
 
 extern "C" {
-	#include "libSphere.h"
-	#include "multishell.h"
+#include "libSphere.h"
 }
 
+typedef struct {
+  double scale;
+  double core_radius;
+  double s_thickness;
+  double w_thickness;
+  double core_sld;
+  double shell_sld;
+  double n_pairs;
+  double background;
+
+} MultiShellParameters;
+
 MultiShellModel :: MultiShellModel() {
-	scale      = Parameter(1.0);
-	core_radius     = Parameter(60.0, true);
-	core_radius.set_min(0.0);
-	s_thickness  = Parameter(10.0, true);
-	s_thickness.set_min(0.0);
-	w_thickness   = Parameter(10.0, true);
-	w_thickness.set_min(0.0);
-	core_sld   = Parameter(6.4e-6);
-	shell_sld   = Parameter(4.0e-7);
-	n_pairs   = Parameter(2);
-	background = Parameter(0.0);
+  scale      = Parameter(1.0);
+  core_radius     = Parameter(60.0, true);
+  core_radius.set_min(0.0);
+  s_thickness  = Parameter(10.0, true);
+  s_thickness.set_min(0.0);
+  w_thickness   = Parameter(10.0, true);
+  w_thickness.set_min(0.0);
+  core_sld   = Parameter(6.4e-6);
+  shell_sld   = Parameter(4.0e-7);
+  n_pairs   = Parameter(2);
+  background = Parameter(0.0);
 }
 
@@ -52,58 +63,58 @@
  */
 double MultiShellModel :: operator()(double q) {
-	double dp[8];
+  double dp[8];
 
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = scale();
-	dp[1] = core_radius();
-	dp[2] = s_thickness();
-	dp[3] = w_thickness();
-	dp[4] = core_sld();
-	dp[5] = shell_sld();
-	dp[6] = n_pairs();
-	dp[7] = 0.0;
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = core_radius();
+  dp[2] = s_thickness();
+  dp[3] = w_thickness();
+  dp[4] = core_sld();
+  dp[5] = shell_sld();
+  dp[6] = n_pairs();
+  dp[7] = 0.0;
 
-	// Get the dispersion points for the core radius
-	vector<WeightPoint> weights_core_radius;
-	core_radius.get_weights(weights_core_radius);
+  // Get the dispersion points for the core radius
+  vector<WeightPoint> weights_core_radius;
+  core_radius.get_weights(weights_core_radius);
 
-	// Get the dispersion points for the s_thickness
-	vector<WeightPoint> weights_s_thickness;
-	s_thickness.get_weights(weights_s_thickness);
+  // Get the dispersion points for the s_thickness
+  vector<WeightPoint> weights_s_thickness;
+  s_thickness.get_weights(weights_s_thickness);
 
-	// Get the dispersion points for the w_thickness
-	vector<WeightPoint> weights_w_thickness;
-	w_thickness.get_weights(weights_w_thickness);
+  // Get the dispersion points for the w_thickness
+  vector<WeightPoint> weights_w_thickness;
+  w_thickness.get_weights(weights_w_thickness);
 
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
 
-	// Loop over radius weight points
-	for(int i=0; i< (int)weights_core_radius.size(); i++) {
-		dp[1] = weights_core_radius[i].value;
-		for(int j=0; j< (int)weights_s_thickness.size(); j++){
-			dp[2] = weights_s_thickness[j].value;
-			for(int k=0; k< (int)weights_w_thickness.size(); k++){
-				dp[3] = weights_w_thickness[k].value;
-				//Un-normalize SphereForm by volume
-				sum += weights_core_radius[i].weight*weights_s_thickness[j].weight
-					*weights_w_thickness[k].weight* MultiShell(dp, q)
-					*pow(weights_core_radius[i].value+dp[6]*weights_s_thickness[j].value+(dp[6]-1)*weights_w_thickness[k].value,3);
-				//Find average volume
-				vol += weights_core_radius[i].weight*weights_s_thickness[j].weight
-					*weights_w_thickness[k].weight
-					*pow(weights_core_radius[i].value+dp[6]*weights_s_thickness[j].value+(dp[6]-1)*weights_w_thickness[k].value,3);
-				norm += weights_core_radius[i].weight*weights_s_thickness[j].weight
-					*weights_w_thickness[k].weight;
-			}
-		}
-	}
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
-	return sum/norm + background();
+  // Loop over radius weight points
+  for(int i=0; i< (int)weights_core_radius.size(); i++) {
+    dp[1] = weights_core_radius[i].value;
+    for(int j=0; j< (int)weights_s_thickness.size(); j++){
+      dp[2] = weights_s_thickness[j].value;
+      for(int k=0; k< (int)weights_w_thickness.size(); k++){
+        dp[3] = weights_w_thickness[k].value;
+        //Un-normalize SphereForm by volume
+        sum += weights_core_radius[i].weight*weights_s_thickness[j].weight
+            *weights_w_thickness[k].weight* MultiShell(dp, q)
+        *pow(weights_core_radius[i].value+dp[6]*weights_s_thickness[j].value+(dp[6]-1)*weights_w_thickness[k].value,3);
+        //Find average volume
+        vol += weights_core_radius[i].weight*weights_s_thickness[j].weight
+            *weights_w_thickness[k].weight
+            *pow(weights_core_radius[i].value+dp[6]*weights_s_thickness[j].value+(dp[6]-1)*weights_w_thickness[k].value,3);
+        norm += weights_core_radius[i].weight*weights_s_thickness[j].weight
+            *weights_w_thickness[k].weight;
+      }
+    }
+  }
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+  return sum/norm + background();
 }
 
@@ -115,6 +126,6 @@
  */
 double MultiShellModel :: operator()(double qx, double qy) {
-	double q = sqrt(qx*qx + qy*qy);
-	return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -127,5 +138,5 @@
  */
 double MultiShellModel :: evaluate_rphi(double q, double phi) {
-	return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 /**
@@ -134,52 +145,52 @@
  */
 double MultiShellModel :: calculate_ER() {
-	MultiShellParameters dp;
+  MultiShellParameters dp;
 
-	dp.core_radius     = core_radius();
-	dp.s_thickness  = s_thickness();
-	dp.w_thickness  = w_thickness();
-	dp.n_pairs = n_pairs();
+  dp.core_radius     = core_radius();
+  dp.s_thickness  = s_thickness();
+  dp.w_thickness  = w_thickness();
+  dp.n_pairs = n_pairs();
 
-	double rad_out = 0.0;
+  double rad_out = 0.0;
 
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	if (dp.n_pairs <= 0.0 ){
-		dp.n_pairs = 0.0;
-	}
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  if (dp.n_pairs <= 0.0 ){
+    dp.n_pairs = 0.0;
+  }
 
-	// Get the dispersion points for the core radius
-	vector<WeightPoint> weights_core_radius;
-	core_radius.get_weights(weights_core_radius);
+  // Get the dispersion points for the core radius
+  vector<WeightPoint> weights_core_radius;
+  core_radius.get_weights(weights_core_radius);
 
-	// Get the dispersion points for the s_thickness
-	vector<WeightPoint> weights_s_thickness;
-	s_thickness.get_weights(weights_s_thickness);
+  // Get the dispersion points for the s_thickness
+  vector<WeightPoint> weights_s_thickness;
+  s_thickness.get_weights(weights_s_thickness);
 
-	// Get the dispersion points for the w_thickness
-	vector<WeightPoint> weights_w_thickness;
-	w_thickness.get_weights(weights_w_thickness);
+  // Get the dispersion points for the w_thickness
+  vector<WeightPoint> weights_w_thickness;
+  w_thickness.get_weights(weights_w_thickness);
 
-	// Loop over major shell weight points
-	for(int i=0; i< (int)weights_s_thickness.size(); i++) {
-		dp.s_thickness = weights_s_thickness[i].value;
-		for(int j=0; j< (int)weights_w_thickness.size(); j++) {
-			dp.w_thickness = weights_w_thickness[j].value;
-			for(int k=0; k< (int)weights_core_radius.size(); k++) {
-				dp.core_radius = weights_core_radius[k].value;
-				sum += weights_s_thickness[i].weight*weights_w_thickness[j].weight
-					* weights_core_radius[k].weight*(dp.core_radius+dp.n_pairs*dp.s_thickness+(dp.n_pairs-1.0)*dp.w_thickness);
-				norm += weights_s_thickness[i].weight*weights_w_thickness[j].weight* weights_core_radius[k].weight;
-			}
-		}
-	}
-	if (norm != 0){
-		//return the averaged value
-		rad_out =  sum/norm;}
-	else{
-		//return normal value
-		rad_out = (dp.core_radius+dp.n_pairs*dp.s_thickness+(dp.n_pairs-1.0)*dp.w_thickness);}
+  // Loop over major shell weight points
+  for(int i=0; i< (int)weights_s_thickness.size(); i++) {
+    dp.s_thickness = weights_s_thickness[i].value;
+    for(int j=0; j< (int)weights_w_thickness.size(); j++) {
+      dp.w_thickness = weights_w_thickness[j].value;
+      for(int k=0; k< (int)weights_core_radius.size(); k++) {
+        dp.core_radius = weights_core_radius[k].value;
+        sum += weights_s_thickness[i].weight*weights_w_thickness[j].weight
+            * weights_core_radius[k].weight*(dp.core_radius+dp.n_pairs*dp.s_thickness+(dp.n_pairs-1.0)*dp.w_thickness);
+        norm += weights_s_thickness[i].weight*weights_w_thickness[j].weight* weights_core_radius[k].weight;
+      }
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    rad_out = (dp.core_radius+dp.n_pairs*dp.s_thickness+(dp.n_pairs-1.0)*dp.w_thickness);}
 
-	return rad_out;
+  return rad_out;
 }
Index: sansmodels/src/c_models/parallelepiped.cpp
===================================================================
--- sansmodels/src/c_models/parallelepiped.cpp	(revision 8343e18b86d18130057f1d8c582e4042eec09721)
+++ sansmodels/src/c_models/parallelepiped.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -22,5 +22,4 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
Index: sansmodels/src/c_models/polygausscoil.cpp
===================================================================
--- sansmodels/src/c_models/polygausscoil.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/polygausscoil.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -21,12 +21,11 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "polygausscoil.h"
 
 extern "C" {
 	#include "libTwoPhase.h"
-	#include "polygausscoil.h"
 }
 
Index: sansmodels/src/c_models/sc.cpp
===================================================================
--- sansmodels/src/c_models/sc.cpp	(revision 6e10cffb542a4045bc3c41ac60c4abf3ed220bb2)
+++ sansmodels/src/c_models/sc.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -21,5 +21,4 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
Index: sansmodels/src/c_models/sphere.cpp
===================================================================
--- sansmodels/src/c_models/sphere.cpp	(revision a8eab1c358d21162446eae5171596814b4f89422)
+++ sansmodels/src/c_models/sphere.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -21,5 +21,4 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
Index: sansmodels/src/c_models/spheroid.cpp
===================================================================
--- sansmodels/src/c_models/spheroid.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/spheroid.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -18,36 +18,131 @@
  *   sansmodels/src/libigor
  *
- *	TODO: refactor so that we pull in the old sansmodels.c_extensions
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
+#include <stdlib.h>
 using namespace std;
+#include "spheroid.h"
 
 extern "C" {
-	#include "libCylinder.h"
-	#include "libStructureFactor.h"
-	#include "spheroid.h"
+#include "libCylinder.h"
+#include "libStructureFactor.h"
+}
+
+typedef struct {
+  double scale;
+  double equat_core;
+  double polar_core;
+  double equat_shell;
+  double polar_shell;
+  double sld_core;
+  double sld_shell;
+  double sld_solvent;
+  double background;
+  double axis_theta;
+  double axis_phi;
+
+} SpheroidParameters;
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the prolate
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double spheroid_analytical_2D_scaled(SpheroidParameters *pars, double q, double q_x, double q_y) {
+
+  double cyl_x, cyl_y, cyl_z;
+  double q_z;
+  double alpha, vol, cos_val;
+  double answer;
+  double Pi = 4.0*atan(1.0);
+  double sldcs,sldss;
+
+  //convert angle degree to radian
+  double theta = pars->axis_theta * Pi/180.0;
+  double phi = pars->axis_phi * Pi/180.0;
+
+
+  // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis.
+  cyl_x = sin(theta) * cos(phi);
+  cyl_y = sin(theta) * sin(phi);
+  cyl_z = cos(theta);
+  //del sld
+  sldcs = pars->sld_core - pars->sld_shell;
+  sldss = pars->sld_shell- pars->sld_solvent;
+
+  // q vector
+  q_z = 0;
+
+  // Compute the angle btw vector q and the
+  // axis of the cylinder
+  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
+
+  // The following test should always pass
+  if (fabs(cos_val)>1.0) {
+    printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
+    return 0;
+  }
+
+  // Note: cos(alpha) = 0 and 1 will get an
+  // undefined value from CylKernel
+  alpha = acos( cos_val );
+
+  // Call the IGOR library function to get the kernel: MUST use gfn4 not gf2 because of the def of params.
+  answer = gfn4(cos_val,pars->equat_core,pars->polar_core,pars->equat_shell,pars->polar_shell,sldcs,sldss,q);
+  //It seems that it should be normalized somehow. How???
+
+  //normalize by cylinder volume
+  //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
+  vol = 4.0*Pi/3.0*pars->equat_shell*pars->equat_shell*pars->polar_shell;
+  answer /= vol;
+
+  //convert to [cm-1]
+  answer *= 1.0e8;
+
+  //Scale
+  answer *= pars->scale;
+
+  // add in the background
+  answer += pars->background;
+
+  return answer;
 }
 
 CoreShellEllipsoidModel :: CoreShellEllipsoidModel() {
-	scale      = Parameter(1.0);
-	equat_core     = Parameter(200.0, true);
-	equat_core.set_min(0.0);
-	polar_core     = Parameter(20.0, true);
-	polar_core.set_min(0.0);
-	equat_shell   = Parameter(250.0, true);
-	equat_shell.set_min(0.0);
-	polar_shell    = Parameter(30.0, true);
-	polar_shell.set_min(0.0);
-	sld_core   = Parameter(2e-6);
-	sld_shell  = Parameter(1e-6);
-	sld_solvent = Parameter(6.3e-6);
-	background = Parameter(0.0);
-	axis_theta  = Parameter(0.0, true);
-	axis_phi    = Parameter(0.0, true);
-
+  scale      = Parameter(1.0);
+  equat_core     = Parameter(200.0, true);
+  equat_core.set_min(0.0);
+  polar_core     = Parameter(20.0, true);
+  polar_core.set_min(0.0);
+  equat_shell   = Parameter(250.0, true);
+  equat_shell.set_min(0.0);
+  polar_shell    = Parameter(30.0, true);
+  polar_shell.set_min(0.0);
+  sld_core   = Parameter(2e-6);
+  sld_shell  = Parameter(1e-6);
+  sld_solvent = Parameter(6.3e-6);
+  background = Parameter(0.0);
+  axis_theta  = Parameter(0.0, true);
+  axis_phi    = Parameter(0.0, true);
+
+}
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the prolate
+ * @param q: q-value
+ * @return: function value
+ */
+static double spheroid_analytical_2DXY(SpheroidParameters *pars, double qx, double qy) {
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return spheroid_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
@@ -59,74 +154,74 @@
  */
 double CoreShellEllipsoidModel :: operator()(double q) {
-	double dp[9];
-
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = scale();
-	dp[1] = equat_core();
-	dp[2] = polar_core();
-	dp[3] = equat_shell();
-	dp[4] = polar_shell();
-	dp[5] = sld_core();
-	dp[6] = sld_shell();
-	dp[7] = sld_solvent();
-	dp[8] = 0.0;
-
-	// Get the dispersion points for the major core
-	vector<WeightPoint> weights_equat_core;
-	equat_core.get_weights(weights_equat_core);
-
-	// Get the dispersion points for the minor core
-	vector<WeightPoint> weights_polar_core;
-	polar_core.get_weights(weights_polar_core);
-
-	// Get the dispersion points for the major shell
-	vector<WeightPoint> weights_equat_shell;
-	equat_shell.get_weights(weights_equat_shell);
-
-	// Get the dispersion points for the minor_shell
-	vector<WeightPoint> weights_polar_shell;
-	polar_shell.get_weights(weights_polar_shell);
-
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
-
-	// Loop over major core weight points
-	for(int i=0; i<(int)weights_equat_core.size(); i++) {
-		dp[1] = weights_equat_core[i].value;
-
-		// Loop over minor core weight points
-		for(int j=0; j<(int)weights_polar_core.size(); j++) {
-			dp[2] = weights_polar_core[j].value;
-
-			// Loop over major shell weight points
-			for(int k=0; k<(int)weights_equat_shell.size(); k++) {
-				dp[3] = weights_equat_shell[k].value;
-
-				// Loop over minor shell weight points
-				for(int l=0; l<(int)weights_polar_shell.size(); l++) {
-					dp[4] = weights_polar_shell[l].value;
-					//Un-normalize  by volume
-					sum += weights_equat_core[i].weight* weights_polar_core[j].weight * weights_equat_shell[k].weight
-						* weights_polar_shell[l].weight * OblateForm(dp, q)
-						* pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
-					//Find average volume
-					vol += weights_equat_core[i].weight* weights_polar_core[j].weight
-						* weights_equat_shell[k].weight
-						* weights_polar_shell[l].weight
-						* pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
-					norm += weights_equat_core[i].weight* weights_polar_core[j].weight * weights_equat_shell[k].weight
-							* weights_polar_shell[l].weight;
-				}
-			}
-		}
-	}
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
-	return sum/norm + background();
+  double dp[9];
+
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = equat_core();
+  dp[2] = polar_core();
+  dp[3] = equat_shell();
+  dp[4] = polar_shell();
+  dp[5] = sld_core();
+  dp[6] = sld_shell();
+  dp[7] = sld_solvent();
+  dp[8] = 0.0;
+
+  // Get the dispersion points for the major core
+  vector<WeightPoint> weights_equat_core;
+  equat_core.get_weights(weights_equat_core);
+
+  // Get the dispersion points for the minor core
+  vector<WeightPoint> weights_polar_core;
+  polar_core.get_weights(weights_polar_core);
+
+  // Get the dispersion points for the major shell
+  vector<WeightPoint> weights_equat_shell;
+  equat_shell.get_weights(weights_equat_shell);
+
+  // Get the dispersion points for the minor_shell
+  vector<WeightPoint> weights_polar_shell;
+  polar_shell.get_weights(weights_polar_shell);
+
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over major core weight points
+  for(int i=0; i<(int)weights_equat_core.size(); i++) {
+    dp[1] = weights_equat_core[i].value;
+
+    // Loop over minor core weight points
+    for(int j=0; j<(int)weights_polar_core.size(); j++) {
+      dp[2] = weights_polar_core[j].value;
+
+      // Loop over major shell weight points
+      for(int k=0; k<(int)weights_equat_shell.size(); k++) {
+        dp[3] = weights_equat_shell[k].value;
+
+        // Loop over minor shell weight points
+        for(int l=0; l<(int)weights_polar_shell.size(); l++) {
+          dp[4] = weights_polar_shell[l].value;
+          //Un-normalize  by volume
+          sum += weights_equat_core[i].weight* weights_polar_core[j].weight * weights_equat_shell[k].weight
+              * weights_polar_shell[l].weight * OblateForm(dp, q)
+          * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
+          //Find average volume
+          vol += weights_equat_core[i].weight* weights_polar_core[j].weight
+              * weights_equat_shell[k].weight
+              * weights_polar_shell[l].weight
+              * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
+          norm += weights_equat_core[i].weight* weights_polar_core[j].weight * weights_equat_shell[k].weight
+              * weights_polar_shell[l].weight;
+        }
+      }
+    }
+  }
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+  return sum/norm + background();
 }
 
@@ -143,5 +238,5 @@
 	return (*this).operator()(q);
 }
-*/
+ */
 
 /**
@@ -153,7 +248,7 @@
  */
 double CoreShellEllipsoidModel :: evaluate_rphi(double q, double phi) {
-	double qx = q*cos(phi);
-	double qy = q*sin(phi);
-	return (*this).operator()(qx, qy);
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
 }
 
@@ -165,110 +260,110 @@
  */
 double CoreShellEllipsoidModel :: operator()(double qx, double qy) {
-	SpheroidParameters dp;
-	// Fill parameter array
-	dp.scale      = scale();
-	dp.equat_core = equat_core();
-	dp.polar_core = polar_core();
-	dp.equat_shell = equat_shell();
-	dp.polar_shell = polar_shell();
-	dp.sld_core = sld_core();
-	dp.sld_shell = sld_shell();
-	dp.sld_solvent = sld_solvent();
-	dp.background = 0.0;
-	dp.axis_theta = axis_theta();
-	dp.axis_phi = axis_phi();
-
-	// Get the dispersion points for the major core
-	vector<WeightPoint> weights_equat_core;
-	equat_core.get_weights(weights_equat_core);
-
-	// Get the dispersion points for the minor core
-	vector<WeightPoint> weights_polar_core;
-	polar_core.get_weights(weights_polar_core);
-
-	// Get the dispersion points for the major shell
-	vector<WeightPoint> weights_equat_shell;
-	equat_shell.get_weights(weights_equat_shell);
-
-	// Get the dispersion points for the minor shell
-	vector<WeightPoint> weights_polar_shell;
-	polar_shell.get_weights(weights_polar_shell);
-
-
-	// Get angular averaging for theta
-	vector<WeightPoint> weights_theta;
-	axis_theta.get_weights(weights_theta);
-
-	// Get angular averaging for phi
-	vector<WeightPoint> weights_phi;
-	axis_phi.get_weights(weights_phi);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double norm_vol = 0.0;
-	double vol = 0.0;
-	double pi = 4.0*atan(1.0);
-	// Loop over major core weight points
-	for(int i=0; i< (int)weights_equat_core.size(); i++) {
-		dp.equat_core = weights_equat_core[i].value;
-
-		// Loop over minor core weight points
-		for(int j=0; j< (int)weights_polar_core.size(); j++) {
-			dp.polar_core = weights_polar_core[j].value;
-
-			// Loop over major shell weight points
-			for(int k=0; k< (int)weights_equat_shell.size(); k++) {
-				dp.equat_shell = weights_equat_shell[i].value;
-
-				// Loop over minor shell weight points
-				for(int l=0; l< (int)weights_polar_shell.size(); l++) {
-					dp.polar_shell = weights_polar_shell[l].value;
-
-					// Average over theta distribution
-					for(int m=0; m< (int)weights_theta.size(); m++) {
-						dp.axis_theta = weights_theta[m].value;
-
-						// Average over phi distribution
-						for(int n=0; n< (int)weights_phi.size(); n++) {
-							dp.axis_phi = weights_phi[n].value;
-							//Un-normalize by volume
-							double _ptvalue = weights_equat_core[i].weight *weights_polar_core[j].weight
-								* weights_equat_shell[k].weight * weights_polar_shell[l].weight
-								* weights_theta[m].weight
-								* weights_phi[n].weight
-								* spheroid_analytical_2DXY(&dp, qx, qy)
-								* pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
-							if (weights_theta.size()>1) {
-								_ptvalue *= fabs(sin(weights_theta[m].value*pi/180.0));
-							}
-							sum += _ptvalue;
-							//Find average volume
-							vol += weights_equat_shell[k].weight
-								* weights_polar_shell[l].weight
-								* pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
-							//Find norm for volume
-							norm_vol += weights_equat_shell[k].weight
-								* weights_polar_shell[l].weight;
-
-							norm += weights_equat_core[i].weight *weights_polar_core[j].weight
-								* weights_equat_shell[k].weight * weights_polar_shell[l].weight
-								* weights_theta[m].weight * weights_phi[n].weight;
-						}
-					}
-				}
-			}
-		}
-	}
-	// Averaging in theta needs an extra normalization
-	// factor to account for the sin(theta) term in the
-	// integration (see documentation).
-	if (weights_theta.size()>1) norm = norm / asin(1.0);
-
-	if (vol != 0.0 && norm_vol != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm_vol);}
-
-	return sum/norm + background();
+  SpheroidParameters dp;
+  // Fill parameter array
+  dp.scale      = scale();
+  dp.equat_core = equat_core();
+  dp.polar_core = polar_core();
+  dp.equat_shell = equat_shell();
+  dp.polar_shell = polar_shell();
+  dp.sld_core = sld_core();
+  dp.sld_shell = sld_shell();
+  dp.sld_solvent = sld_solvent();
+  dp.background = 0.0;
+  dp.axis_theta = axis_theta();
+  dp.axis_phi = axis_phi();
+
+  // Get the dispersion points for the major core
+  vector<WeightPoint> weights_equat_core;
+  equat_core.get_weights(weights_equat_core);
+
+  // Get the dispersion points for the minor core
+  vector<WeightPoint> weights_polar_core;
+  polar_core.get_weights(weights_polar_core);
+
+  // Get the dispersion points for the major shell
+  vector<WeightPoint> weights_equat_shell;
+  equat_shell.get_weights(weights_equat_shell);
+
+  // Get the dispersion points for the minor shell
+  vector<WeightPoint> weights_polar_shell;
+  polar_shell.get_weights(weights_polar_shell);
+
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_theta;
+  axis_theta.get_weights(weights_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_phi;
+  axis_phi.get_weights(weights_phi);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double norm_vol = 0.0;
+  double vol = 0.0;
+  double pi = 4.0*atan(1.0);
+  // Loop over major core weight points
+  for(int i=0; i< (int)weights_equat_core.size(); i++) {
+    dp.equat_core = weights_equat_core[i].value;
+
+    // Loop over minor core weight points
+    for(int j=0; j< (int)weights_polar_core.size(); j++) {
+      dp.polar_core = weights_polar_core[j].value;
+
+      // Loop over major shell weight points
+      for(int k=0; k< (int)weights_equat_shell.size(); k++) {
+        dp.equat_shell = weights_equat_shell[i].value;
+
+        // Loop over minor shell weight points
+        for(int l=0; l< (int)weights_polar_shell.size(); l++) {
+          dp.polar_shell = weights_polar_shell[l].value;
+
+          // Average over theta distribution
+          for(int m=0; m< (int)weights_theta.size(); m++) {
+            dp.axis_theta = weights_theta[m].value;
+
+            // Average over phi distribution
+            for(int n=0; n< (int)weights_phi.size(); n++) {
+              dp.axis_phi = weights_phi[n].value;
+              //Un-normalize by volume
+              double _ptvalue = weights_equat_core[i].weight *weights_polar_core[j].weight
+                  * weights_equat_shell[k].weight * weights_polar_shell[l].weight
+                  * weights_theta[m].weight
+                  * weights_phi[n].weight
+                  * spheroid_analytical_2DXY(&dp, qx, qy)
+              * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
+              if (weights_theta.size()>1) {
+                _ptvalue *= fabs(sin(weights_theta[m].value*pi/180.0));
+              }
+              sum += _ptvalue;
+              //Find average volume
+              vol += weights_equat_shell[k].weight
+                  * weights_polar_shell[l].weight
+                  * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
+              //Find norm for volume
+              norm_vol += weights_equat_shell[k].weight
+                  * weights_polar_shell[l].weight;
+
+              norm += weights_equat_core[i].weight *weights_polar_core[j].weight
+                  * weights_equat_shell[k].weight * weights_polar_shell[l].weight
+                  * weights_theta[m].weight * weights_phi[n].weight;
+            }
+          }
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_theta.size()>1) norm = norm / asin(1.0);
+
+  if (vol != 0.0 && norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm_vol);}
+
+  return sum/norm + background();
 }
 
@@ -278,42 +373,42 @@
  */
 double CoreShellEllipsoidModel :: calculate_ER() {
-	SpheroidParameters dp;
-
-	dp.equat_shell = equat_shell();
-	dp.polar_shell = polar_shell();
-
-	double rad_out = 0.0;
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-
-	// Get the dispersion points for the major shell
-	vector<WeightPoint> weights_equat_shell;
-	equat_shell.get_weights(weights_equat_shell);
-
-	// Get the dispersion points for the minor shell
-	vector<WeightPoint> weights_polar_shell;
-	polar_shell.get_weights(weights_polar_shell);
-
-	// Loop over major shell weight points
-	for(int i=0; i< (int)weights_equat_shell.size(); i++) {
-		dp.equat_shell = weights_equat_shell[i].value;
-		for(int k=0; k< (int)weights_polar_shell.size(); k++) {
-			dp.polar_shell = weights_polar_shell[k].value;
-			//Note: output of "DiamEllip(dp.polar_shell,dp.equat_shell)" is DIAMETER.
-			sum +=weights_equat_shell[i].weight
-				* weights_polar_shell[k].weight*DiamEllip(dp.polar_shell,dp.equat_shell)/2.0;
-			norm += weights_equat_shell[i].weight* weights_polar_shell[k].weight;
-		}
-	}
-	if (norm != 0){
-		//return the averaged value
-		rad_out =  sum/norm;}
-	else{
-		//return normal value
-		//Note: output of "DiamEllip(dp.polar_shell,dp.equat_shell)" is DIAMETER.
-		rad_out = DiamEllip(dp.polar_shell,dp.equat_shell)/2.0;}
-
-	return rad_out;
-}
+  SpheroidParameters dp;
+
+  dp.equat_shell = equat_shell();
+  dp.polar_shell = polar_shell();
+
+  double rad_out = 0.0;
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the major shell
+  vector<WeightPoint> weights_equat_shell;
+  equat_shell.get_weights(weights_equat_shell);
+
+  // Get the dispersion points for the minor shell
+  vector<WeightPoint> weights_polar_shell;
+  polar_shell.get_weights(weights_polar_shell);
+
+  // Loop over major shell weight points
+  for(int i=0; i< (int)weights_equat_shell.size(); i++) {
+    dp.equat_shell = weights_equat_shell[i].value;
+    for(int k=0; k< (int)weights_polar_shell.size(); k++) {
+      dp.polar_shell = weights_polar_shell[k].value;
+      //Note: output of "DiamEllip(dp.polar_shell,dp.equat_shell)" is DIAMETER.
+      sum +=weights_equat_shell[i].weight
+          * weights_polar_shell[k].weight*DiamEllip(dp.polar_shell,dp.equat_shell)/2.0;
+      norm += weights_equat_shell[i].weight* weights_polar_shell[k].weight;
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    //Note: output of "DiamEllip(dp.polar_shell,dp.equat_shell)" is DIAMETER.
+    rad_out = DiamEllip(dp.polar_shell,dp.equat_shell)/2.0;}
+
+  return rad_out;
+}
Index: sansmodels/src/c_models/stackeddisks.cpp
===================================================================
--- sansmodels/src/c_models/stackeddisks.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/stackeddisks.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -23,31 +23,123 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "stacked_disks.h"
 
 extern "C" {
-	#include "libCylinder.h"
-	#include "libStructureFactor.h"
-	#include "stacked_disks.h"
+#include "libCylinder.h"
+#include "libStructureFactor.h"
+}
+
+typedef struct {
+  double scale;
+  double radius;
+  double core_thick;
+  double layer_thick;
+  double core_sld;
+  double layer_sld;
+  double solvent_sld;
+  double n_stacking;
+  double sigma_d;
+  double background;
+  double axis_theta;
+  double axis_phi;
+} StackedDisksParameters;
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the staked disks
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double stacked_disks_analytical_2D_scaled(StackedDisksParameters *pars, double q, double q_x, double q_y) {
+  double cyl_x, cyl_y, cyl_z;
+  double q_z;
+  double alpha, vol, cos_val;
+  double d, dum, halfheight;
+  double answer;
+
+
+
+  // parallelepiped orientation
+  cyl_x = sin(pars->axis_theta) * cos(pars->axis_phi);
+  cyl_y = sin(pars->axis_theta) * sin(pars->axis_phi);
+  cyl_z = cos(pars->axis_theta);
+
+  // q vector
+  q_z = 0;
+
+  // Compute the angle btw vector q and the
+  // axis of the parallelepiped
+  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
+
+  // The following test should always pass
+  if (fabs(cos_val)>1.0) {
+    printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+    return 0;
+  }
+
+  // Note: cos(alpha) = 0 and 1 will get an
+  // undefined value from Stackdisc_kern
+  alpha = acos( cos_val );
+
+  // Call the IGOR library function to get the kernel
+  d = 2 * pars->layer_thick + pars->core_thick;
+  halfheight = pars->core_thick/2.0;
+  dum =alpha ;
+  answer = Stackdisc_kern(q, pars->radius, pars->core_sld,pars->layer_sld,
+      pars->solvent_sld, halfheight, pars->layer_thick, dum, pars->sigma_d, d, pars->n_stacking)/sin(alpha);
+
+  // Multiply by contrast^2
+  //answer *= pars->contrast*pars->contrast;
+
+  //normalize by staked disks volume
+  vol = acos(-1.0) * pars->radius * pars->radius * d * pars->n_stacking;
+  answer /= vol;
+
+  //convert to [cm-1]
+  answer *= 1.0e8;
+
+  //Scale
+  answer *= pars->scale;
+
+  // add in the background
+  answer += pars->background;
+
+  return answer;
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the staked disks
+ * @param q: q-value
+ * @return: function value
+ */
+static double stacked_disks_analytical_2DXY(StackedDisksParameters *pars, double qx, double qy) {
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return stacked_disks_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
 StackedDisksModel :: StackedDisksModel() {
-	scale      = Parameter(1.0);
-	radius     = Parameter(3000.0, true);
-	radius.set_min(0.0);
-	core_thick  = Parameter(10.0, true);
-	core_thick.set_min(0.0);
-	layer_thick     = Parameter(15.0);
-	layer_thick.set_min(0.0);
-	core_sld = Parameter(4.0e-6);
-	layer_sld  = Parameter(-4.0e-7);
-	solvent_sld  = Parameter(5.0e-6);
-	n_stacking   = Parameter(1);
-	sigma_d   = Parameter(0);
-	background = Parameter(0.001);
-	axis_theta  = Parameter(0.0, true);
-	axis_phi    = Parameter(0.0, true);
+  scale      = Parameter(1.0);
+  radius     = Parameter(3000.0, true);
+  radius.set_min(0.0);
+  core_thick  = Parameter(10.0, true);
+  core_thick.set_min(0.0);
+  layer_thick     = Parameter(15.0);
+  layer_thick.set_min(0.0);
+  core_sld = Parameter(4.0e-6);
+  layer_sld  = Parameter(-4.0e-7);
+  solvent_sld  = Parameter(5.0e-6);
+  n_stacking   = Parameter(1);
+  sigma_d   = Parameter(0);
+  background = Parameter(0.001);
+  axis_theta  = Parameter(0.0, true);
+  axis_phi    = Parameter(0.0, true);
 }
 
@@ -59,65 +151,65 @@
  */
 double StackedDisksModel :: operator()(double q) {
-	double dp[10];
-
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = scale();
-	dp[1] = radius();
-	dp[2] = core_thick();
-	dp[3] = layer_thick();
-	dp[4] = core_sld();
-	dp[5] = layer_sld();
-	dp[6] = solvent_sld();
-	dp[7] = n_stacking();
-	dp[8] = sigma_d();
-	dp[9] = 0.0;
-
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_radius;
-	radius.get_weights(weights_radius);
-
-	// Get the dispersion points for the core_thick
-	vector<WeightPoint> weights_core_thick;
-	core_thick.get_weights(weights_core_thick);
-
-	// Get the dispersion points for the layer_thick
-	vector<WeightPoint> weights_layer_thick;
-	layer_thick.get_weights(weights_layer_thick);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
-
-	// Loop over length weight points
-	for(int i=0; i< (int)weights_radius.size(); i++) {
-		dp[1] = weights_radius[i].value;
-
-		// Loop over radius weight points
-		for(int j=0; j< (int)weights_core_thick.size(); j++) {
-			dp[2] = weights_core_thick[j].value;
-
-			// Loop over thickness weight points
-			for(int k=0; k< (int)weights_layer_thick.size(); k++) {
-				dp[3] = weights_layer_thick[k].value;
-				//Un-normalize by volume
-				sum += weights_radius[i].weight
-					* weights_core_thick[j].weight * weights_layer_thick[k].weight* StackedDiscs(dp, q)
-					*pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value);
-				//Find average volume
-				vol += weights_radius[i].weight
-					* weights_core_thick[j].weight * weights_layer_thick[k].weight
-					*pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value);
-				norm += weights_radius[i].weight
-					* weights_core_thick[j].weight* weights_layer_thick[k].weight;
-			}
-		}
-	}
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
-
-	return sum/norm + background();
+  double dp[10];
+
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = radius();
+  dp[2] = core_thick();
+  dp[3] = layer_thick();
+  dp[4] = core_sld();
+  dp[5] = layer_sld();
+  dp[6] = solvent_sld();
+  dp[7] = n_stacking();
+  dp[8] = sigma_d();
+  dp[9] = 0.0;
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_radius;
+  radius.get_weights(weights_radius);
+
+  // Get the dispersion points for the core_thick
+  vector<WeightPoint> weights_core_thick;
+  core_thick.get_weights(weights_core_thick);
+
+  // Get the dispersion points for the layer_thick
+  vector<WeightPoint> weights_layer_thick;
+  layer_thick.get_weights(weights_layer_thick);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over length weight points
+  for(int i=0; i< (int)weights_radius.size(); i++) {
+    dp[1] = weights_radius[i].value;
+
+    // Loop over radius weight points
+    for(int j=0; j< (int)weights_core_thick.size(); j++) {
+      dp[2] = weights_core_thick[j].value;
+
+      // Loop over thickness weight points
+      for(int k=0; k< (int)weights_layer_thick.size(); k++) {
+        dp[3] = weights_layer_thick[k].value;
+        //Un-normalize by volume
+        sum += weights_radius[i].weight
+            * weights_core_thick[j].weight * weights_layer_thick[k].weight* StackedDiscs(dp, q)
+        *pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value);
+        //Find average volume
+        vol += weights_radius[i].weight
+            * weights_core_thick[j].weight * weights_layer_thick[k].weight
+            *pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value);
+        norm += weights_radius[i].weight
+            * weights_core_thick[j].weight* weights_layer_thick[k].weight;
+      }
+    }
+  }
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -129,102 +221,102 @@
  */
 double StackedDisksModel :: operator()(double qx, double qy) {
-	StackedDisksParameters dp;
-	// Fill parameter array
-	dp.scale      = scale();
-	dp.core_thick    = core_thick();
-	dp.radius	  = radius();
-	dp.layer_thick  = layer_thick();
-	dp.core_sld   = core_sld();
-	dp.layer_sld  = layer_sld();
-	dp.solvent_sld= solvent_sld();
-	dp.n_stacking	  = n_stacking();
-	dp.sigma_d   = sigma_d();
-	dp.background = 0.0;
-	dp.axis_theta = axis_theta();
-	dp.axis_phi   = axis_phi();
-
-	// Get the dispersion points for the length
-	vector<WeightPoint> weights_core_thick;
-	core_thick.get_weights(weights_core_thick);
-
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_radius;
-	radius.get_weights(weights_radius);
-
-	// Get the dispersion points for the thickness
-	vector<WeightPoint> weights_layer_thick;
-	layer_thick.get_weights(weights_layer_thick);
-
-	// Get angular averaging for theta
-	vector<WeightPoint> weights_theta;
-	axis_theta.get_weights(weights_theta);
-
-	// Get angular averaging for phi
-	vector<WeightPoint> weights_phi;
-	axis_phi.get_weights(weights_phi);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double norm_vol = 0.0;
-	double vol = 0.0;
-
-	// Loop over length weight points
-	for(int i=0; i< (int)weights_core_thick.size(); i++) {
-		dp.core_thick = weights_core_thick[i].value;
-
-		// Loop over radius weight points
-		for(int j=0; j< (int)weights_radius.size(); j++) {
-			dp.radius = weights_radius[j].value;
-
-				// Loop over thickness weight points
-				for(int k=0; k< (int)weights_layer_thick.size(); k++) {
-				dp.layer_thick = weights_layer_thick[k].value;
-
-					for(int l=0; l< (int)weights_theta.size(); l++) {
-					dp.axis_theta = weights_theta[l].value;
-
-						// Average over phi distribution
-						for(int m=0; m <(int)weights_phi.size(); m++) {
-							dp.axis_phi = weights_phi[m].value;
-
-							//Un-normalize by volume
-							double _ptvalue = weights_core_thick[i].weight
-								* weights_radius[j].weight
-								* weights_layer_thick[k].weight
-								* weights_theta[l].weight
-								* weights_phi[m].weight
-								* stacked_disks_analytical_2DXY(&dp, qx, qy)
-								*pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value);
-							if (weights_theta.size()>1) {
-								_ptvalue *= fabs(sin(weights_theta[l].value));
-							}
-							sum += _ptvalue;
-							//Find average volume
-							vol += weights_radius[j].weight
-								* weights_core_thick[i].weight * weights_layer_thick[k].weight
-								*pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value);
-							//Find norm for volume
-							norm_vol += weights_radius[j].weight
-								* weights_core_thick[i].weight * weights_layer_thick[k].weight;
-
-							norm += weights_core_thick[i].weight
-								* weights_radius[j].weight
-								* weights_layer_thick[k].weight
-								* weights_theta[l].weight
-								* weights_phi[m].weight;
-						}
-				}
-			}
-		}
-	}
-	// Averaging in theta needs an extra normalization
-	// factor to account for the sin(theta) term in the
-	// integration (see documentation).
-	if (weights_theta.size()>1) norm = norm / asin(1.0);
-	if (vol != 0.0 && norm_vol != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm_vol);}
-	return sum/norm + background();
+  StackedDisksParameters dp;
+  // Fill parameter array
+  dp.scale      = scale();
+  dp.core_thick    = core_thick();
+  dp.radius	  = radius();
+  dp.layer_thick  = layer_thick();
+  dp.core_sld   = core_sld();
+  dp.layer_sld  = layer_sld();
+  dp.solvent_sld= solvent_sld();
+  dp.n_stacking	  = n_stacking();
+  dp.sigma_d   = sigma_d();
+  dp.background = 0.0;
+  dp.axis_theta = axis_theta();
+  dp.axis_phi   = axis_phi();
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_core_thick;
+  core_thick.get_weights(weights_core_thick);
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_radius;
+  radius.get_weights(weights_radius);
+
+  // Get the dispersion points for the thickness
+  vector<WeightPoint> weights_layer_thick;
+  layer_thick.get_weights(weights_layer_thick);
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_theta;
+  axis_theta.get_weights(weights_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_phi;
+  axis_phi.get_weights(weights_phi);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double norm_vol = 0.0;
+  double vol = 0.0;
+
+  // Loop over length weight points
+  for(int i=0; i< (int)weights_core_thick.size(); i++) {
+    dp.core_thick = weights_core_thick[i].value;
+
+    // Loop over radius weight points
+    for(int j=0; j< (int)weights_radius.size(); j++) {
+      dp.radius = weights_radius[j].value;
+
+      // Loop over thickness weight points
+      for(int k=0; k< (int)weights_layer_thick.size(); k++) {
+        dp.layer_thick = weights_layer_thick[k].value;
+
+        for(int l=0; l< (int)weights_theta.size(); l++) {
+          dp.axis_theta = weights_theta[l].value;
+
+          // Average over phi distribution
+          for(int m=0; m <(int)weights_phi.size(); m++) {
+            dp.axis_phi = weights_phi[m].value;
+
+            //Un-normalize by volume
+            double _ptvalue = weights_core_thick[i].weight
+                * weights_radius[j].weight
+                * weights_layer_thick[k].weight
+                * weights_theta[l].weight
+                * weights_phi[m].weight
+                * stacked_disks_analytical_2DXY(&dp, qx, qy)
+            *pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value);
+            if (weights_theta.size()>1) {
+              _ptvalue *= fabs(sin(weights_theta[l].value));
+            }
+            sum += _ptvalue;
+            //Find average volume
+            vol += weights_radius[j].weight
+                * weights_core_thick[i].weight * weights_layer_thick[k].weight
+                *pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value);
+            //Find norm for volume
+            norm_vol += weights_radius[j].weight
+                * weights_core_thick[i].weight * weights_layer_thick[k].weight;
+
+            norm += weights_core_thick[i].weight
+                * weights_radius[j].weight
+                * weights_layer_thick[k].weight
+                * weights_theta[l].weight
+                * weights_phi[m].weight;
+          }
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_theta.size()>1) norm = norm / asin(1.0);
+  if (vol != 0.0 && norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm_vol);}
+  return sum/norm + background();
 }
 
@@ -237,7 +329,7 @@
  */
 double StackedDisksModel :: evaluate_rphi(double q, double phi) {
-	double qx = q*cos(phi);
-	double qy = q*sin(phi);
-	return (*this).operator()(qx, qy);
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
 }
 /**
@@ -246,54 +338,54 @@
  */
 double StackedDisksModel :: calculate_ER() {
-	StackedDisksParameters dp;
-
-	dp.core_thick    = core_thick();
-	dp.radius	  = radius();
-	dp.layer_thick  = layer_thick();
-	dp.n_stacking	  = n_stacking();
-
-	double rad_out = 0.0;
-	if (dp.n_stacking <= 0.0){
-		return rad_out;
-	}
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-
-	// Get the dispersion points for the length
-	vector<WeightPoint> weights_core_thick;
-	core_thick.get_weights(weights_core_thick);
-
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_radius;
-	radius.get_weights(weights_radius);
-
-	// Get the dispersion points for the thickness
-	vector<WeightPoint> weights_layer_thick;
-	layer_thick.get_weights(weights_layer_thick);
-
-	// Loop over major shell weight points
-	for(int i=0; i< (int)weights_core_thick.size(); i++) {
-		dp.core_thick = weights_core_thick[i].value;
-		for(int j=0; j< (int)weights_layer_thick.size(); j++) {
-			dp.layer_thick = weights_layer_thick[j].value;
-			for(int k=0; k< (int)weights_radius.size(); k++) {
-				dp.radius = weights_radius[k].value;
-				//Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-				sum +=weights_core_thick[i].weight*weights_layer_thick[j].weight
-					* weights_radius[k].weight*DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0;
-				norm += weights_core_thick[i].weight*weights_layer_thick[j].weight* weights_radius[k].weight;
-			}
-		}
-	}
-	if (norm != 0){
-		//return the averaged value
-		rad_out =  sum/norm;}
-	else{
-		//return normal value
-		//Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-		rad_out = DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0;}
-
-	return rad_out;
-}
+  StackedDisksParameters dp;
+
+  dp.core_thick    = core_thick();
+  dp.radius	  = radius();
+  dp.layer_thick  = layer_thick();
+  dp.n_stacking	  = n_stacking();
+
+  double rad_out = 0.0;
+  if (dp.n_stacking <= 0.0){
+    return rad_out;
+  }
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_core_thick;
+  core_thick.get_weights(weights_core_thick);
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_radius;
+  radius.get_weights(weights_radius);
+
+  // Get the dispersion points for the thickness
+  vector<WeightPoint> weights_layer_thick;
+  layer_thick.get_weights(weights_layer_thick);
+
+  // Loop over major shell weight points
+  for(int i=0; i< (int)weights_core_thick.size(); i++) {
+    dp.core_thick = weights_core_thick[i].value;
+    for(int j=0; j< (int)weights_layer_thick.size(); j++) {
+      dp.layer_thick = weights_layer_thick[j].value;
+      for(int k=0; k< (int)weights_radius.size(); k++) {
+        dp.radius = weights_radius[k].value;
+        //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+        sum +=weights_core_thick[i].weight*weights_layer_thick[j].weight
+            * weights_radius[k].weight*DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0;
+        norm += weights_core_thick[i].weight*weights_layer_thick[j].weight* weights_radius[k].weight;
+      }
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+    rad_out = DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0;}
+
+  return rad_out;
+}
Index: sansmodels/src/c_models/triaxialellipsoid.cpp
===================================================================
--- sansmodels/src/c_models/triaxialellipsoid.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/triaxialellipsoid.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -18,33 +18,159 @@
  *   sansmodels/src/libigor
  *
- *	TODO: refactor so that we pull in the old sansmodels.c_extensions
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
+#include <stdlib.h>
 using namespace std;
+#include "triaxial_ellipsoid.h"
 
 extern "C" {
-	#include "libCylinder.h"
-	#include "libStructureFactor.h"
-	#include "triaxial_ellipsoid.h"
-}
+#include "libCylinder.h"
+#include "libStructureFactor.h"
+}
+
+typedef struct {
+  double scale;
+  double semi_axisA;
+  double semi_axisB;
+  double semi_axisC;
+  double sldEll;
+  double sldSolv;
+  double background;
+  double axis_theta;
+  double axis_phi;
+  double axis_psi;
+
+} TriaxialEllipsoidParameters;
+
+static double triaxial_ellipsoid_kernel(TriaxialEllipsoidParameters *pars, double q, double alpha, double nu) {
+  double t,a,b,c;
+  double kernel;
+
+  a = pars->semi_axisA ;
+  b = pars->semi_axisB ;
+  c = pars->semi_axisC ;
+
+  t = q * sqrt(a*a*cos(nu)*cos(nu)+b*b*sin(nu)*sin(nu)*sin(alpha)*sin(alpha)+c*c*cos(alpha)*cos(alpha));
+  if (t==0.0){
+    kernel  = 1.0;
+  }else{
+    kernel  = 3.0*(sin(t)-t*cos(t))/(t*t*t);
+  }
+  return kernel*kernel;
+}
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the triaxial ellipsoid
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double triaxial_ellipsoid_analytical_2D_scaled(TriaxialEllipsoidParameters *pars, double q, double q_x, double q_y) {
+  double cyl_x, cyl_y, cyl_z, ell_x, ell_y;
+  double q_z;
+  double cos_nu,nu;
+  double alpha, vol, cos_val;
+  double answer;
+  double pi = 4.0*atan(1.0);
+
+  //convert angle degree to radian
+  double theta = pars->axis_theta * pi/180.0;
+  double phi = pars->axis_phi * pi/180.0;
+  double psi = pars->axis_psi * pi/180.0;
+
+  // Cylinder orientation
+  cyl_x = sin(theta) * cos(phi);
+  cyl_y = sin(theta) * sin(phi);
+  cyl_z = cos(theta);
+
+  // q vector
+  q_z = 0.0;
+
+  //dx = 1.0;
+  //dy = 1.0;
+  // Compute the angle btw vector q and the
+  // axis of the cylinder
+  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
+
+  // The following test should always pass
+  if (fabs(cos_val)>1.0) {
+    printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
+    return 0;
+  }
+
+  // Note: cos(alpha) = 0 and 1 will get an
+  // undefined value from CylKernel
+  alpha = acos( cos_val );
+
+  //ellipse orientation:
+  // the elliptical corss section was transformed and projected
+  // into the detector plane already through sin(alpha)and furthermore psi remains as same
+  // on the detector plane.
+  // So, all we need is to calculate the angle (nu) of the minor axis of the ellipse wrt
+  // the wave vector q.
+
+  //x- y- component on the detector plane.
+  ell_x =  cos(psi);
+  ell_y =  sin(psi);
+
+  // calculate the axis of the ellipse wrt q-coord.
+  cos_nu = ell_x*q_x + ell_y*q_y;
+  nu = acos(cos_nu);
+
+  // Call the IGOR library function to get the kernel
+  answer = triaxial_ellipsoid_kernel(pars, q, alpha, nu);
+
+  // Multiply by contrast^2
+  answer *= (pars->sldEll- pars->sldSolv)*(pars->sldEll- pars->sldSolv);
+
+  //normalize by cylinder volume
+  //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
+  vol = 4.0* pi/3.0  * pars->semi_axisA * pars->semi_axisB * pars->semi_axisC;
+  answer *= vol;
+  //convert to [cm-1]
+  answer *= 1.0e8;
+  //Scale
+  answer *= pars->scale;
+
+  // add in the background
+  answer += pars->background;
+
+  return answer;
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the triaxial ellipsoid
+ * @param q: q-value
+ * @return: function value
+ */
+static double triaxial_ellipsoid_analytical_2DXY(TriaxialEllipsoidParameters *pars, double qx, double qy) {
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return triaxial_ellipsoid_analytical_2D_scaled(pars, q, qx/q, qy/q);
+}
+
+
 
 TriaxialEllipsoidModel :: TriaxialEllipsoidModel() {
-	scale      = Parameter(1.0);
-	semi_axisA     = Parameter(35.0, true);
-	semi_axisA.set_min(0.0);
-	semi_axisB     = Parameter(100.0, true);
-	semi_axisB.set_min(0.0);
-	semi_axisC  = Parameter(400.0, true);
-	semi_axisC.set_min(0.0);
-	sldEll   = Parameter(1.0e-6);
-	sldSolv   = Parameter(6.3e-6);
-	background = Parameter(0.0);
-	axis_theta  = Parameter(57.325, true);
-	axis_phi    = Parameter(57.325, true);
-	axis_psi    = Parameter(0.0, true);
+  scale      = Parameter(1.0);
+  semi_axisA     = Parameter(35.0, true);
+  semi_axisA.set_min(0.0);
+  semi_axisB     = Parameter(100.0, true);
+  semi_axisB.set_min(0.0);
+  semi_axisC  = Parameter(400.0, true);
+  semi_axisC.set_min(0.0);
+  sldEll   = Parameter(1.0e-6);
+  sldSolv   = Parameter(6.3e-6);
+  background = Parameter(0.0);
+  axis_theta  = Parameter(57.325, true);
+  axis_phi    = Parameter(57.325, true);
+  axis_psi    = Parameter(0.0, true);
 }
 
@@ -56,63 +182,63 @@
  */
 double TriaxialEllipsoidModel :: operator()(double q) {
-	double dp[7];
-
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = scale();
-	dp[1] = semi_axisA();
-	dp[2] = semi_axisB();
-	dp[3] = semi_axisC();
-	dp[4] = sldEll();
-	dp[5] = sldSolv();
-	dp[6] = 0.0;
-
-	// Get the dispersion points for the semi axis A
-	vector<WeightPoint> weights_semi_axisA;
-	semi_axisA.get_weights(weights_semi_axisA);
-
-	// Get the dispersion points for the semi axis B
-	vector<WeightPoint> weights_semi_axisB;
-	semi_axisB.get_weights(weights_semi_axisB);
-
-	// Get the dispersion points for the semi axis C
-	vector<WeightPoint> weights_semi_axisC;
-	semi_axisC.get_weights(weights_semi_axisC);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
-
-	// Loop over semi axis A weight points
-	for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
-		dp[1] = weights_semi_axisA[i].value;
-
-		// Loop over semi axis B weight points
-		for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
-			dp[2] = weights_semi_axisB[j].value;
-
-			// Loop over semi axis C weight points
-			for(int k=0; k< (int)weights_semi_axisC.size(); k++) {
-				dp[3] = weights_semi_axisC[k].value;
-				//Un-normalize  by volume
-				sum += weights_semi_axisA[i].weight
-					* weights_semi_axisB[j].weight * weights_semi_axisC[k].weight* TriaxialEllipsoid(dp, q)
-					* weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
-				//Find average volume
-				vol += weights_semi_axisA[i].weight
-					* weights_semi_axisB[j].weight * weights_semi_axisC[k].weight
-					* weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
-
-				norm += weights_semi_axisA[i].weight
-					* weights_semi_axisB[j].weight * weights_semi_axisC[k].weight;
-			}
-		}
-	}
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
-
-	return sum/norm + background();
+  double dp[7];
+
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = semi_axisA();
+  dp[2] = semi_axisB();
+  dp[3] = semi_axisC();
+  dp[4] = sldEll();
+  dp[5] = sldSolv();
+  dp[6] = 0.0;
+
+  // Get the dispersion points for the semi axis A
+  vector<WeightPoint> weights_semi_axisA;
+  semi_axisA.get_weights(weights_semi_axisA);
+
+  // Get the dispersion points for the semi axis B
+  vector<WeightPoint> weights_semi_axisB;
+  semi_axisB.get_weights(weights_semi_axisB);
+
+  // Get the dispersion points for the semi axis C
+  vector<WeightPoint> weights_semi_axisC;
+  semi_axisC.get_weights(weights_semi_axisC);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over semi axis A weight points
+  for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
+    dp[1] = weights_semi_axisA[i].value;
+
+    // Loop over semi axis B weight points
+    for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
+      dp[2] = weights_semi_axisB[j].value;
+
+      // Loop over semi axis C weight points
+      for(int k=0; k< (int)weights_semi_axisC.size(); k++) {
+        dp[3] = weights_semi_axisC[k].value;
+        //Un-normalize  by volume
+        sum += weights_semi_axisA[i].weight
+            * weights_semi_axisB[j].weight * weights_semi_axisC[k].weight* TriaxialEllipsoid(dp, q)
+        * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
+        //Find average volume
+        vol += weights_semi_axisA[i].weight
+            * weights_semi_axisB[j].weight * weights_semi_axisC[k].weight
+            * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
+
+        norm += weights_semi_axisA[i].weight
+            * weights_semi_axisB[j].weight * weights_semi_axisC[k].weight;
+      }
+    }
+  }
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -124,115 +250,115 @@
  */
 double TriaxialEllipsoidModel :: operator()(double qx, double qy) {
-	TriaxialEllipsoidParameters dp;
-	// Fill parameter array
-	dp.scale      = scale();
-	dp.semi_axisA   = semi_axisA();
-	dp.semi_axisB     = semi_axisB();
-	dp.semi_axisC     = semi_axisC();
-	dp.sldEll   = sldEll();
-	dp.sldSolv   = sldSolv();
-	dp.background = 0.0;
-	dp.axis_theta  = axis_theta();
-	dp.axis_phi    = axis_phi();
-	dp.axis_psi    = axis_psi();
-
-	// Get the dispersion points for the semi_axis A
-	vector<WeightPoint> weights_semi_axisA;
-	semi_axisA.get_weights(weights_semi_axisA);
-
-	// Get the dispersion points for the semi_axis B
-	vector<WeightPoint> weights_semi_axisB;
-	semi_axisB.get_weights(weights_semi_axisB);
-
-	// Get the dispersion points for the semi_axis C
-	vector<WeightPoint> weights_semi_axisC;
-	semi_axisC.get_weights(weights_semi_axisC);
-
-	// Get angular averaging for theta
-	vector<WeightPoint> weights_theta;
-	axis_theta.get_weights(weights_theta);
-
-	// Get angular averaging for phi
-	vector<WeightPoint> weights_phi;
-	axis_phi.get_weights(weights_phi);
-
-	// Get angular averaging for psi
-	vector<WeightPoint> weights_psi;
-	axis_psi.get_weights(weights_psi);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double norm_vol = 0.0;
-	double vol = 0.0;
-	double pi = 4.0*atan(1.0);
-	// Loop over semi axis A weight points
-	for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
-		dp.semi_axisA = weights_semi_axisA[i].value;
-
-		// Loop over semi axis B weight points
-		for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
-			dp.semi_axisB = weights_semi_axisB[j].value;
-
-			// Loop over semi axis C weight points
-			for(int k=0; k < (int)weights_semi_axisC.size(); k++) {
-			dp.semi_axisC = weights_semi_axisC[k].value;
-
-				// Average over theta distribution
-				for(int l=0; l< (int)weights_theta.size(); l++) {
-					dp.axis_theta = weights_theta[l].value;
-
-					// Average over phi distribution
-					for(int m=0; m <(int)weights_phi.size(); m++) {
-						dp.axis_phi = weights_phi[m].value;
-						// Average over psi distribution
-						for(int n=0; n <(int)weights_psi.size(); n++) {
-							dp.axis_psi = weights_psi[n].value;
-							//Un-normalize  by volume
-							double _ptvalue = weights_semi_axisA[i].weight
-								* weights_semi_axisB[j].weight
-								* weights_semi_axisC[k].weight
-								* weights_theta[l].weight
-								* weights_phi[m].weight
-								* weights_psi[n].weight
-								* triaxial_ellipsoid_analytical_2DXY(&dp, qx, qy)
-								* weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
-							if (weights_theta.size()>1) {
-								_ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
-							}
-							sum += _ptvalue;
-							//Find average volume
-							vol += weights_semi_axisA[i].weight
-								* weights_semi_axisB[j].weight
-								* weights_semi_axisC[k].weight
-								* weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
-							//Find norm for volume
-							norm_vol += weights_semi_axisA[i].weight
-								* weights_semi_axisB[j].weight
-								* weights_semi_axisC[k].weight;
-
-							norm += weights_semi_axisA[i].weight
-								* weights_semi_axisB[j].weight
-								* weights_semi_axisC[k].weight
-								* weights_theta[l].weight
-								* weights_phi[m].weight
-								* weights_psi[n].weight;
-						}
-					}
-
-				}
-			}
-		}
-	}
-	// Averaging in theta needs an extra normalization
-	// factor to account for the sin(theta) term in the
-	// integration (see documentation).
-	if (weights_theta.size()>1) norm = norm / asin(1.0);
-
-	if (vol != 0.0 && norm_vol != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm_vol);}
-
-	return sum/norm + background();
+  TriaxialEllipsoidParameters dp;
+  // Fill parameter array
+  dp.scale      = scale();
+  dp.semi_axisA   = semi_axisA();
+  dp.semi_axisB     = semi_axisB();
+  dp.semi_axisC     = semi_axisC();
+  dp.sldEll   = sldEll();
+  dp.sldSolv   = sldSolv();
+  dp.background = 0.0;
+  dp.axis_theta  = axis_theta();
+  dp.axis_phi    = axis_phi();
+  dp.axis_psi    = axis_psi();
+
+  // Get the dispersion points for the semi_axis A
+  vector<WeightPoint> weights_semi_axisA;
+  semi_axisA.get_weights(weights_semi_axisA);
+
+  // Get the dispersion points for the semi_axis B
+  vector<WeightPoint> weights_semi_axisB;
+  semi_axisB.get_weights(weights_semi_axisB);
+
+  // Get the dispersion points for the semi_axis C
+  vector<WeightPoint> weights_semi_axisC;
+  semi_axisC.get_weights(weights_semi_axisC);
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_theta;
+  axis_theta.get_weights(weights_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_phi;
+  axis_phi.get_weights(weights_phi);
+
+  // Get angular averaging for psi
+  vector<WeightPoint> weights_psi;
+  axis_psi.get_weights(weights_psi);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double norm_vol = 0.0;
+  double vol = 0.0;
+  double pi = 4.0*atan(1.0);
+  // Loop over semi axis A weight points
+  for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
+    dp.semi_axisA = weights_semi_axisA[i].value;
+
+    // Loop over semi axis B weight points
+    for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
+      dp.semi_axisB = weights_semi_axisB[j].value;
+
+      // Loop over semi axis C weight points
+      for(int k=0; k < (int)weights_semi_axisC.size(); k++) {
+        dp.semi_axisC = weights_semi_axisC[k].value;
+
+        // Average over theta distribution
+        for(int l=0; l< (int)weights_theta.size(); l++) {
+          dp.axis_theta = weights_theta[l].value;
+
+          // Average over phi distribution
+          for(int m=0; m <(int)weights_phi.size(); m++) {
+            dp.axis_phi = weights_phi[m].value;
+            // Average over psi distribution
+            for(int n=0; n <(int)weights_psi.size(); n++) {
+              dp.axis_psi = weights_psi[n].value;
+              //Un-normalize  by volume
+              double _ptvalue = weights_semi_axisA[i].weight
+                  * weights_semi_axisB[j].weight
+                  * weights_semi_axisC[k].weight
+                  * weights_theta[l].weight
+                  * weights_phi[m].weight
+                  * weights_psi[n].weight
+                  * triaxial_ellipsoid_analytical_2DXY(&dp, qx, qy)
+              * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
+              if (weights_theta.size()>1) {
+                _ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
+              }
+              sum += _ptvalue;
+              //Find average volume
+              vol += weights_semi_axisA[i].weight
+                  * weights_semi_axisB[j].weight
+                  * weights_semi_axisC[k].weight
+                  * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
+              //Find norm for volume
+              norm_vol += weights_semi_axisA[i].weight
+                  * weights_semi_axisB[j].weight
+                  * weights_semi_axisC[k].weight;
+
+              norm += weights_semi_axisA[i].weight
+                  * weights_semi_axisB[j].weight
+                  * weights_semi_axisC[k].weight
+                  * weights_theta[l].weight
+                  * weights_phi[m].weight
+                  * weights_psi[n].weight;
+            }
+          }
+
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_theta.size()>1) norm = norm / asin(1.0);
+
+  if (vol != 0.0 && norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm_vol);}
+
+  return sum/norm + background();
 }
 
@@ -245,7 +371,7 @@
  */
 double TriaxialEllipsoidModel :: evaluate_rphi(double q, double phi) {
-	double qx = q*cos(phi);
-	double qy = q*sin(phi);
-	return (*this).operator()(qx, qy);
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
 }
 /**
@@ -254,63 +380,63 @@
  */
 double TriaxialEllipsoidModel :: calculate_ER() {
-	TriaxialEllipsoidParameters dp;
-
-	dp.semi_axisA   = semi_axisA();
-	dp.semi_axisB     = semi_axisB();
-	//polar axis C
-	dp.semi_axisC     = semi_axisC();
-
-	double rad_out = 0.0;
-	//Surface average radius at the equat. cross section.
-	double suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-
-	// Get the dispersion points for the semi_axis A
-	vector<WeightPoint> weights_semi_axisA;
-	semi_axisA.get_weights(weights_semi_axisA);
-
-	// Get the dispersion points for the semi_axis B
-	vector<WeightPoint> weights_semi_axisB;
-	semi_axisB.get_weights(weights_semi_axisB);
-
-	// Get the dispersion points for the semi_axis C
-	vector<WeightPoint> weights_semi_axisC;
-	semi_axisC.get_weights(weights_semi_axisC);
-
-	// Loop over semi axis A weight points
-	for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
-		dp.semi_axisA = weights_semi_axisA[i].value;
-
-		// Loop over semi axis B weight points
-		for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
-			dp.semi_axisB = weights_semi_axisB[j].value;
-
-			// Loop over semi axis C weight points
-			for(int k=0; k < (int)weights_semi_axisC.size(); k++) {
-				dp.semi_axisC = weights_semi_axisC[k].value;
-
-				//Calculate surface averaged radius
-				suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB);
-
-				//Sum
-				sum += weights_semi_axisA[i].weight
-					* weights_semi_axisB[j].weight
-					* weights_semi_axisC[k].weight * DiamEllip(dp.semi_axisC, suf_rad)/2.0;
-				//Norm
-				norm += weights_semi_axisA[i].weight* weights_semi_axisB[j].weight
-					* weights_semi_axisC[k].weight;
-			}
-		}
-	}
-	if (norm != 0){
-		//return the averaged value
-		rad_out =  sum/norm;}
-	else{
-		//return normal value
-		rad_out = DiamEllip(dp.semi_axisC, suf_rad)/2.0;}
-
-	return rad_out;
-}
+  TriaxialEllipsoidParameters dp;
+
+  dp.semi_axisA   = semi_axisA();
+  dp.semi_axisB     = semi_axisB();
+  //polar axis C
+  dp.semi_axisC     = semi_axisC();
+
+  double rad_out = 0.0;
+  //Surface average radius at the equat. cross section.
+  double suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the semi_axis A
+  vector<WeightPoint> weights_semi_axisA;
+  semi_axisA.get_weights(weights_semi_axisA);
+
+  // Get the dispersion points for the semi_axis B
+  vector<WeightPoint> weights_semi_axisB;
+  semi_axisB.get_weights(weights_semi_axisB);
+
+  // Get the dispersion points for the semi_axis C
+  vector<WeightPoint> weights_semi_axisC;
+  semi_axisC.get_weights(weights_semi_axisC);
+
+  // Loop over semi axis A weight points
+  for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
+    dp.semi_axisA = weights_semi_axisA[i].value;
+
+    // Loop over semi axis B weight points
+    for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
+      dp.semi_axisB = weights_semi_axisB[j].value;
+
+      // Loop over semi axis C weight points
+      for(int k=0; k < (int)weights_semi_axisC.size(); k++) {
+        dp.semi_axisC = weights_semi_axisC[k].value;
+
+        //Calculate surface averaged radius
+        suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB);
+
+        //Sum
+        sum += weights_semi_axisA[i].weight
+            * weights_semi_axisB[j].weight
+            * weights_semi_axisC[k].weight * DiamEllip(dp.semi_axisC, suf_rad)/2.0;
+        //Norm
+        norm += weights_semi_axisA[i].weight* weights_semi_axisB[j].weight
+            * weights_semi_axisC[k].weight;
+      }
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    rad_out = DiamEllip(dp.semi_axisC, suf_rad)/2.0;}
+
+  return rad_out;
+}
Index: sansmodels/src/c_models/vesicle.cpp
===================================================================
--- sansmodels/src/c_models/vesicle.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/vesicle.cpp	(revision 82c11d354f0cffa7741bd2e74e0641e042eeb034)
@@ -21,23 +21,31 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "vesicle.h"
 
 extern "C" {
-	#include "libSphere.h"
-	#include "vesicle.h"
+#include "libSphere.h"
 }
 
+typedef struct {
+  double scale;
+  double radius;
+  double thickness;
+  double core_sld;
+  double shell_sld;
+  double background;
+} VesicleParameters;
+
 VesicleModel :: VesicleModel() {
-	scale      = Parameter(1.0);
-	radius     = Parameter(100.0, true);
-	radius.set_min(0.0);
-	thickness  = Parameter(30.0, true);
-	thickness.set_min(0.0);
-	core_sld   = Parameter(6.36e-6);
-	shell_sld   = Parameter(5.0e-7);
-	background = Parameter(0.0);
+  scale      = Parameter(1.0);
+  radius     = Parameter(100.0, true);
+  radius.set_min(0.0);
+  thickness  = Parameter(30.0, true);
+  thickness.set_min(0.0);
+  core_sld   = Parameter(6.36e-6);
+  shell_sld   = Parameter(5.0e-7);
+  background = Parameter(0.0);
 }
 
@@ -49,47 +57,47 @@
  */
 double VesicleModel :: operator()(double q) {
-	double dp[6];
+  double dp[6];
 
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = scale();
-	dp[1] = radius();
-	dp[2] = thickness();
-	dp[3] = core_sld();
-	dp[4] = shell_sld();
-	dp[5] = 0.0;
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = radius();
+  dp[2] = thickness();
+  dp[3] = core_sld();
+  dp[4] = shell_sld();
+  dp[5] = 0.0;
 
 
-	// Get the dispersion points for the core radius
-	vector<WeightPoint> weights_radius;
-	radius.get_weights(weights_radius);
-	// Get the dispersion points for the thickness
-	vector<WeightPoint> weights_thickness;
-	thickness.get_weights(weights_thickness);
+  // Get the dispersion points for the core radius
+  vector<WeightPoint> weights_radius;
+  radius.get_weights(weights_radius);
+  // Get the dispersion points for the thickness
+  vector<WeightPoint> weights_thickness;
+  thickness.get_weights(weights_thickness);
 
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
 
-	// Loop over radius weight points
-	for(int i=0; i< (int)weights_radius.size(); i++) {
-		dp[1] = weights_radius[i].value;
-		for(int j=0; j< (int)weights_thickness.size(); j++) {
-			dp[2] = weights_thickness[j].value;
-			sum += weights_radius[i].weight
-				* weights_thickness[j].weight * VesicleForm(dp, q)
-				*(pow(weights_radius[i].value+weights_thickness[j].value,3)-pow(weights_radius[i].value,3));
-			//Find average volume
-			vol += weights_radius[i].weight * weights_thickness[j].weight
-				*(pow(weights_radius[i].value+weights_thickness[j].value,3)-pow(weights_radius[i].value,3));
-			norm += weights_radius[i].weight * weights_thickness[j].weight;
-		}
-	}
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
+  // Loop over radius weight points
+  for(int i=0; i< (int)weights_radius.size(); i++) {
+    dp[1] = weights_radius[i].value;
+    for(int j=0; j< (int)weights_thickness.size(); j++) {
+      dp[2] = weights_thickness[j].value;
+      sum += weights_radius[i].weight
+          * weights_thickness[j].weight * VesicleForm(dp, q)
+      *(pow(weights_radius[i].value+weights_thickness[j].value,3)-pow(weights_radius[i].value,3));
+      //Find average volume
+      vol += weights_radius[i].weight * weights_thickness[j].weight
+          *(pow(weights_radius[i].value+weights_thickness[j].value,3)-pow(weights_radius[i].value,3));
+      norm += weights_radius[i].weight * weights_thickness[j].weight;
+    }
+  }
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
 
-	return sum/norm + background();
+  return sum/norm + background();
 }
 
@@ -101,6 +109,6 @@
  */
 double VesicleModel :: operator()(double qx, double qy) {
-	double q = sqrt(qx*qx + qy*qy);
-	return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -113,5 +121,5 @@
  */
 double VesicleModel :: evaluate_rphi(double q, double phi) {
-	return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 /**
@@ -120,41 +128,41 @@
  */
 double VesicleModel :: calculate_ER() {
-	VesicleParameters dp;
+  VesicleParameters dp;
 
-	dp.radius     = radius();
-	dp.thickness  = thickness();
+  dp.radius     = radius();
+  dp.thickness  = thickness();
 
-	double rad_out = 0.0;
+  double rad_out = 0.0;
 
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
 
 
-	// Get the dispersion points for the major shell
-	vector<WeightPoint> weights_thickness;
-	thickness.get_weights(weights_thickness);
+  // Get the dispersion points for the major shell
+  vector<WeightPoint> weights_thickness;
+  thickness.get_weights(weights_thickness);
 
-	// Get the dispersion points for the minor shell
-	vector<WeightPoint> weights_radius ;
-	radius.get_weights(weights_radius);
+  // Get the dispersion points for the minor shell
+  vector<WeightPoint> weights_radius ;
+  radius.get_weights(weights_radius);
 
-	// Loop over major shell weight points
-	for(int j=0; j< (int)weights_thickness.size(); j++) {
-		dp.thickness = weights_thickness[j].value;
-		for(int k=0; k< (int)weights_radius.size(); k++) {
-			dp.radius = weights_radius[k].value;
-			sum += weights_thickness[j].weight
-				* weights_radius[k].weight*(dp.radius+dp.thickness);
-			norm += weights_thickness[j].weight* weights_radius[k].weight;
-		}
-	}
-	if (norm != 0){
-		//return the averaged value
-		rad_out =  sum/norm;}
-	else{
-		//return normal value
-		rad_out = (dp.radius+dp.thickness);}
+  // Loop over major shell weight points
+  for(int j=0; j< (int)weights_thickness.size(); j++) {
+    dp.thickness = weights_thickness[j].value;
+    for(int k=0; k< (int)weights_radius.size(); k++) {
+      dp.radius = weights_radius[k].value;
+      sum += weights_thickness[j].weight
+          * weights_radius[k].weight*(dp.radius+dp.thickness);
+      norm += weights_thickness[j].weight* weights_radius[k].weight;
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    rad_out = (dp.radius+dp.thickness);}
 
-	return rad_out;
+  return rad_out;
 }
