Index: nsmodels/src/c_extensions/spheresld.c
===================================================================
--- sansmodels/src/c_extensions/spheresld.c	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ 	(revision )
@@ -1,258 +1,0 @@
-/**
- * spheresld model
- */
-#include <math.h>
-#include "spheresld.h"
-#include "libmultifunc/librefl.h"
-#include <stdio.h>
-#include <stdlib.h>
-
-#define lamda 4.62
-
-
-double sphere_sld_kernel(double dp[], double q) {
-	int n = dp[0];
-	int i,j,k;
-
-	double scale = dp[1];
-	double thick_inter_core = dp[2];
-	double sld_core = dp[4];
-	double sld_solv = dp[5];
-	double background = dp[6];
-	double npts = dp[57]; //number of sub_layers in each interface
-  double nsl=npts;//21.0; //nsl = Num_sub_layer:  MUST ODD number in double //no other number works now
-  int n_s;
-
-  double sld_i,sld_f,dz,bes,fun,f,vol,vol_pre,vol_sub,qr,r,contr,f2;
-  double sign,slope=0.0;
-  double pi;
-
-  int* fun_type;
-  double* sld;
-  double* thick_inter;
-  double* thick;
-  double* fun_coef;
-
-	double total_thick=0.0;
-
-	fun_type = (int*)malloc((n+2)*sizeof(int));
-  sld = (double*)malloc((n+2)*sizeof(double));
-  thick_inter = (double*)malloc((n+2)*sizeof(double));
-  thick = (double*)malloc((n+2)*sizeof(double));
-  fun_coef = (double*)malloc((n+2)*sizeof(double));
-
-	fun_type[0] = dp[3];
-	fun_coef[0] = fabs(dp[58]);
-	for (i =1; i<=n; i++){
-		sld[i] = dp[i+6];
-		thick_inter[i]= dp[i+16];
-		thick[i] = dp[i+26];
-		fun_type[i] = dp[i+36];
-		fun_coef[i] = fabs(dp[i+46]);
-		total_thick += thick[i] + thick_inter[i];
-	}
-	sld[0] = sld_core;
-	sld[n+1] = sld_solv;
-	thick[0] = dp[59];
-	thick[n+1] = total_thick/5.0;
-	thick_inter[0] = thick_inter_core;
-	thick_inter[n+1] = 0.0;
-	fun_coef[n+1] = 0.0;
-
-    pi = 4.0*atan(1.0);
-    f = 0.0;
-    r = 0.0;
-    vol = 0.0;
-    vol_pre = 0.0;
-    vol_sub = 0.0;
-    sld_f = sld_core;
-
-	//floor_nsl = floor(nsl/2.0);
-
-    dz = 0.0;
-    // iteration for # of shells + core + solvent
-    for (i=0;i<=n+1; i++){
-		//iteration for N sub-layers
-    	//if (fabs(thick[i]) <= 1e-24){
-    	//   continue;
-    	//}
-    	// iteration for flat and interface
-		for (j=0;j<2;j++){
-			// iteration for sub_shells in the interface
-			// starts from #1 sub-layer
-			for (n_s=1;n_s<=nsl; n_s++){
-				// for solvent, it doesn't have an interface
-				if (i==n+1 && j==1)
-					break;
-				// for flat layers
-				if (j==0){
-					dz = thick[i];
-					sld_i = sld[i];
-					slope = 0.0;
-				}
-				// for interfacial sub_shells
-				else{
-					dz = thick_inter[i]/nsl;
-					// find sld_i at the outer boundary of sub-layer #n_s
-					sld_i = intersldfunc(fun_type[i],nsl, n_s, fun_coef[i], sld[i], sld[i+1]);
-					// calculate slope
-					slope= (sld_i -sld_f)/dz;
-				}
-				contr = sld_f-slope*r;
-				// iteration for the left and right boundary of the shells(or sub_shells)
-				for (k=0; k<2; k++){
-					// At r=0, the contribution to I is zero, so skip it.
-					if ( i == 0 && j == 0 && k == 0){
-						continue;
-						}
-					// On the top of slovent there is no interface; skip it.
-					if (i == n+1 && k == 1){
-						continue;
-						}
-					// At the right side (outer) boundary
-					if ( k == 1){
-						sign = 1.0;
-						r += dz;
-						}
-					// At the left side (inner) boundary
-					else{
-						sign = -1.0;
-						}
-					qr = q * r;
-					fun = 0.0;
-					if(qr == 0.0){
-						// sigular point
-						bes = sign * 1.0;
-						}
-					else{
-						// for flat sub-layer
-						bes = sign *  3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
-						// with linear slope
-						if (fabs(slope) > 0.0 ){
-							fun = sign * 3.0 * r * (2.0*qr*sin(qr)-((qr*qr)-2.0)*cos(qr))/(qr * qr * qr * qr);
-							}
-						}
-					// update total volume
-					vol = 4.0 * pi / 3.0 * r * r * r;
-					// we won't do the following volume correction for now.
-					// substrate empty area of volume
-					//if (k == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && fun_type[i]==0){
-					//	vol_sub += (vol_pre - vol);
-					//}
-					f += vol * (bes * contr + fun * slope);
-				}
-				// remember this sld as sld_f
-				sld_f =sld_i;
-				// no sub-layer iteration (n_s loop) for the flat layer
-				if (j==0)
-					break;
-			}
-		}
-    }
-    //vol += vol_sub;
-    f2 = f * f / vol * 1.0e8;
-	f2 *= scale;
-	f2 += background;
-
-  free(fun_type);
-  free(sld);
-  free(thick_inter);
-  free(thick);
-  free(fun_coef);
-
-    return (f2);
-}
-
-/**
- * Function to evaluate spheresld function
- * @param pars: parameters of spheresld
- * @param q: q-value
- * @return: function value
- */
-
-double sphere_sld_analytical_1D(SphereSLDParameters *pars, double q) {
-	double dp[60];
-
-	dp[0] = pars->n_shells;
-	dp[1] = pars->scale;
-	dp[2] = pars->thick_inter0;
-	dp[3] = pars->func_inter0;
-	dp[4] = pars->sld_core0;
-	dp[5] = pars->sld_solv;
-	dp[6] = pars->background;
-
-	dp[7] = pars->sld_flat1;
-	dp[8] = pars->sld_flat2;
-	dp[9] = pars->sld_flat3;
-	dp[10] = pars->sld_flat4;
-	dp[11] = pars->sld_flat5;
-	dp[12] = pars->sld_flat6;
-	dp[13] = pars->sld_flat7;
-	dp[14] = pars->sld_flat8;
-	dp[15] = pars->sld_flat9;
-	dp[16] = pars->sld_flat10;
-
-	dp[17] = pars->thick_inter1;
-	dp[18] = pars->thick_inter2;
-	dp[19] = pars->thick_inter3;
-	dp[20] = pars->thick_inter4;
-	dp[21] = pars->thick_inter5;
-	dp[22] = pars->thick_inter6;
-	dp[23] = pars->thick_inter7;
-	dp[24] = pars->thick_inter8;
-	dp[25] = pars->thick_inter9;
-	dp[26] = pars->thick_inter10;
-
-	dp[27] = pars->thick_flat1;
-	dp[28] = pars->thick_flat2;
-	dp[29] = pars->thick_flat3;
-	dp[30] = pars->thick_flat4;
-	dp[31] = pars->thick_flat5;
-	dp[32] = pars->thick_flat6;
-	dp[33] = pars->thick_flat7;
-	dp[34] = pars->thick_flat8;
-	dp[35] = pars->thick_flat9;
-	dp[36] = pars->thick_flat10;
-
-	dp[37] = pars->func_inter1;
-	dp[38] = pars->func_inter2;
-	dp[39] = pars->func_inter3;
-	dp[40] = pars->func_inter4;
-	dp[41] = pars->func_inter5;
-	dp[42] = pars->func_inter6;
-	dp[43] = pars->func_inter7;
-	dp[44] = pars->func_inter8;
-	dp[45] = pars->func_inter9;
-	dp[46] = pars->func_inter10;
-
-	dp[47] = pars->nu_inter1;
-	dp[48] = pars->nu_inter2;
-	dp[49] = pars->nu_inter3;
-	dp[50] = pars->nu_inter4;
-	dp[51] = pars->nu_inter5;
-	dp[52] = pars->nu_inter6;
-	dp[53] = pars->nu_inter7;
-	dp[54] = pars->nu_inter8;
-	dp[55] = pars->nu_inter9;
-	dp[56] = pars->nu_inter10;
-
-	dp[57] = pars->npts_inter;
-	dp[58] = pars->nu_inter0;
-	dp[59] = pars->rad_core0;
-
-	return sphere_sld_kernel(dp, q);
-}
-
-/**
- * Function to evaluate spheresld function
- * @param pars: parameters of spheresld
- * @param q: q-value
- * @return: function value
- */
-double sphere_sld_analytical_2D(SphereSLDParameters *pars, double q, double phi) {
-	return sphere_sld_analytical_1D(pars,q);
-}
-
-double sphere_sld_analytical_2DXY(SphereSLDParameters *pars, double qx, double qy){
-	return sphere_sld_analytical_1D(pars,sqrt(qx*qx+qy*qy));
-}
Index: sansmodels/src/c_extensions/spheresld.h
===================================================================
--- sansmodels/src/c_extensions/spheresld.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/spheresld.h	(revision 37805e993ee77aefcc3df740a88439050390d7ce)
@@ -1,4 +1,5 @@
 #if !defined(o_h)
 #define sphere_sld_h
+#include "parameters.hh"
 
 /**
@@ -30,148 +31,151 @@
  //[ORIENTATION_PARAMS]= <text> </text>
 
-typedef struct {
-	/// number of shells
-	//  [DEFAULT]=n_shells=1
-	int n_shells;
+class SphereSLDModel{
+public:
+  // Model parameters
+  /// number of shells
+  //  [DEFAULT]=n_shells=1
+  Parameter n_shells;
     /// Scale factor
     //  [DEFAULT]=scale= 1.0
-	double scale;
-    ///	thick_inter0 [A]
+  Parameter scale;
+    /// thick_inter0 [A]
     //  [DEFAULT]=thick_inter0=50.0 [A]
-	double thick_inter0;
-	///	func_inter0
-	//  [DEFAULT]=func_inter0= 0
-	double func_inter0;
-	///	sld_core0 [1/A^(2)]
-	//  [DEFAULT]=sld_core0= 2.07e-6 [1/A^(2)]
-	double sld_core0;
-	///	sld_solv [1/A^(2)]
-	//  [DEFAULT]=sld_solv= 1.0e-6 [1/A^(2)]
-	double sld_solv;
-	/// Background
-	//  [DEFAULT]=background=0
-	double background;
+  Parameter thick_inter0;
+  /// func_inter0
+  //  [DEFAULT]=func_inter0= 0
+  Parameter func_inter0;
+  /// sld_core0 [1/A^(2)]
+  //  [DEFAULT]=sld_core0= 2.07e-6 [1/A^(2)]
+  Parameter sld_core0;
+  /// sld_solv [1/A^(2)]
+  //  [DEFAULT]=sld_solv= 1.0e-6 [1/A^(2)]
+  Parameter sld_solv;
+  /// Background
+  //  [DEFAULT]=background=0
+  Parameter background;
 
     //  [DEFAULT]=sld_flat1=4.0e-06 [1/A^(2)]
-    double sld_flat1;
+    Parameter sld_flat1;
     //  [DEFAULT]=sld_flat2=3.5e-06 [1/A^(2)]
-    double sld_flat2;
+    Parameter sld_flat2;
     //  [DEFAULT]=sld_flat3=4.0e-06 [1/A^(2)]
-    double sld_flat3;
+    Parameter sld_flat3;
     //  [DEFAULT]=sld_flat4=3.5e-06 [1/A^(2)]
-    double sld_flat4;
+    Parameter sld_flat4;
     //  [DEFAULT]=sld_flat5=4.0e-06 [1/A^(2)]
-    double sld_flat5;
+    Parameter sld_flat5;
     //  [DEFAULT]=sld_flat6=3.5e-06 [1/A^(2)]
-    double sld_flat6;
+    Parameter sld_flat6;
     //  [DEFAULT]=sld_flat7=4.0e-06 [1/A^(2)]
-    double sld_flat7;
+    Parameter sld_flat7;
     //  [DEFAULT]=sld_flat8=3.5e-06 [1/A^(2)]
-    double sld_flat8;
+    Parameter sld_flat8;
     //  [DEFAULT]=sld_flat9=4.0e-06 [1/A^(2)]
-    double sld_flat9;
+    Parameter sld_flat9;
     //  [DEFAULT]=sld_flat10=3.5e-06 [1/A^(2)]
-    double sld_flat10;
+    Parameter sld_flat10;
 
     //  [DEFAULT]=thick_inter1=50.0 [A]
-    double thick_inter1;
+    Parameter thick_inter1;
     //  [DEFAULT]=thick_inter2=50.0 [A]
-    double thick_inter2;
+    Parameter thick_inter2;
     //  [DEFAULT]=thick_inter3=50.0 [A]
-    double thick_inter3;
+    Parameter thick_inter3;
     //  [DEFAULT]=thick_inter4=50.0 [A]
-    double thick_inter4;
+    Parameter thick_inter4;
     //  [DEFAULT]=thick_inter5=50.0 [A]
-    double thick_inter5;
+    Parameter thick_inter5;
     //  [DEFAULT]=thick_inter6=50.0 [A]
-    double thick_inter6;
+    Parameter thick_inter6;
     //  [DEFAULT]=thick_inter7=50.0 [A]
-    double thick_inter7;
+    Parameter thick_inter7;
     //  [DEFAULT]=thick_inter8=50.0 [A]
-    double thick_inter8;
+    Parameter thick_inter8;
     //  [DEFAULT]=thick_inter9=50.0 [A]
-    double thick_inter9;
+    Parameter thick_inter9;
     //  [DEFAULT]=thick_inter10=50.0 [A]
-    double thick_inter10;
+    Parameter thick_inter10;
 
     //  [DEFAULT]=thick_flat1=100 [A]
-    double thick_flat1;
+    Parameter thick_flat1;
     //  [DEFAULT]=thick_flat2=100 [A]
-    double thick_flat2;
+    Parameter thick_flat2;
     //  [DEFAULT]=thick_flat3=100 [A]
-    double thick_flat3;
+    Parameter thick_flat3;
     //  [DEFAULT]=thick_flat4=100 [A]
-    double thick_flat4;
+    Parameter thick_flat4;
     //  [DEFAULT]=thick_flat5=100 [A]
-    double thick_flat5;
+    Parameter thick_flat5;
     //  [DEFAULT]=thick_flat6=100 [A]
-    double thick_flat6;
+    Parameter thick_flat6;
     //  [DEFAULT]=thick_flat7=100 [A]
-    double thick_flat7;
+    Parameter thick_flat7;
     //  [DEFAULT]=thick_flat8=100 [A]
-    double thick_flat8;
+    Parameter thick_flat8;
     //  [DEFAULT]=thick_flat9=100 [A]
-    double thick_flat9;
+    Parameter thick_flat9;
     //  [DEFAULT]=thick_flat10=100 [A]
-    double thick_flat10;
+    Parameter thick_flat10;
 
     //  [DEFAULT]=func_inter1=0
-    double func_inter1;
+    Parameter func_inter1;
     //  [DEFAULT]=func_inter2=0
-    double func_inter2;
+    Parameter func_inter2;
     //  [DEFAULT]=func_inter3=0
-    double func_inter3;
+    Parameter func_inter3;
     //  [DEFAULT]=func_inter4=0
-    double func_inter4;
+    Parameter func_inter4;
     //  [DEFAULT]=func_inter5=0
-    double func_inter5;
+    Parameter func_inter5;
     //  [DEFAULT]=func_inter6=0
-    double func_inter6;
+    Parameter func_inter6;
     //  [DEFAULT]=func_inter7=0
-    double func_inter7;
+    Parameter func_inter7;
     //  [DEFAULT]=func_inter8=0
-    double func_inter8;
+    Parameter func_inter8;
     //  [DEFAULT]=func_inter9=0
-    double func_inter9;
+    Parameter func_inter9;
     //  [DEFAULT]=func_inter10=0
-    double func_inter10;
+    Parameter func_inter10;
 
     //  [DEFAULT]=nu_inter1=2.5
-    double nu_inter1;
+    Parameter nu_inter1;
     //  [DEFAULT]=nu_inter2=2.5
-    double nu_inter2;
+    Parameter nu_inter2;
     //  [DEFAULT]=nu_inter3=2.5
-    double nu_inter3;
+    Parameter nu_inter3;
     //  [DEFAULT]=nu_inter4=2.5
-    double nu_inter4;
+    Parameter nu_inter4;
     //  [DEFAULT]=nu_inter5=2.5
-    double nu_inter5;
+    Parameter nu_inter5;
     //  [DEFAULT]=nu_inter6=2.5
-    double nu_inter6;
+    Parameter nu_inter6;
     //  [DEFAULT]=nu_inter7=2.5
-    double nu_inter7;
+    Parameter nu_inter7;
     //  [DEFAULT]=nu_inter8=2.5
-    double nu_inter8;
+    Parameter nu_inter8;
     //  [DEFAULT]=nu_inter9=2.5
-    double nu_inter9;
+    Parameter nu_inter9;
     //  [DEFAULT]=nu_inter10=2.5
-    double nu_inter10;
+    Parameter nu_inter10;
 
     //  [DEFAULT]=npts_inter=35.0
-    double npts_inter;
+    Parameter npts_inter;
     //  [DEFAULT]=nu_inter0=2.5
-    double nu_inter0;
+    Parameter nu_inter0;
     //  [DEFAULT]=rad_core0=50.0 [A]
-    double rad_core0;
-} SphereSLDParameters;
+    Parameter rad_core0;
 
-double sphere_sld_kernel(double dq[], double q);
+  // Constructor
+  SphereSLDModel();
 
-/// 1D scattering function
-double sphere_sld_analytical_1D(SphereSLDParameters *pars, double q);
+  // 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);
+};
 
-/// 2D scattering function
-double sphere_sld_analytical_2D(SphereSLDParameters *pars, double q, double phi);
-double sphere_sld_analytical_2DXY(SphereSLDParameters *pars, double qx, double qy);
 
 #endif
Index: sansmodels/src/c_models/models.hh
===================================================================
--- sansmodels/src/c_models/models.hh	(revision a8eab1c358d21162446eae5171596814b4f89422)
+++ sansmodels/src/c_models/models.hh	(revision 37805e993ee77aefcc3df740a88439050390d7ce)
@@ -30,84 +30,4 @@
 
 
-
-class SphereSLDModel{
-public:
-	// Model parameters
-	Parameter n_shells;
-	Parameter scale;
-	Parameter thick_inter0;
-	Parameter func_inter0;
-	Parameter sld_core0;
-	Parameter sld_solv;
-	Parameter background;
-
-	Parameter sld_flat1;
-	Parameter sld_flat2;
-	Parameter sld_flat3;
-	Parameter sld_flat4;
-	Parameter sld_flat5;
-	Parameter sld_flat6;
-	Parameter sld_flat7;
-	Parameter sld_flat8;
-	Parameter sld_flat9;
-	Parameter sld_flat10;
-
-	Parameter thick_inter1;
-	Parameter thick_inter2;
-	Parameter thick_inter3;
-	Parameter thick_inter4;
-	Parameter thick_inter5;
-	Parameter thick_inter6;
-	Parameter thick_inter7;
-	Parameter thick_inter8;
-	Parameter thick_inter9;
-	Parameter thick_inter10;
-
-	Parameter thick_flat1;
-	Parameter thick_flat2;
-	Parameter thick_flat3;
-	Parameter thick_flat4;
-	Parameter thick_flat5;
-	Parameter thick_flat6;
-	Parameter thick_flat7;
-	Parameter thick_flat8;
-	Parameter thick_flat9;
-	Parameter thick_flat10;
-
-	Parameter func_inter1;
-	Parameter func_inter2;
-	Parameter func_inter3;
-	Parameter func_inter4;
-	Parameter func_inter5;
-	Parameter func_inter6;
-	Parameter func_inter7;
-	Parameter func_inter8;
-	Parameter func_inter9;
-	Parameter func_inter10;
-
-	Parameter nu_inter1;
-	Parameter nu_inter2;
-	Parameter nu_inter3;
-	Parameter nu_inter4;
-	Parameter nu_inter5;
-	Parameter nu_inter6;
-	Parameter nu_inter7;
-	Parameter nu_inter8;
-	Parameter nu_inter9;
-	Parameter nu_inter10;
-
-	Parameter npts_inter;
-	Parameter nu_inter0;
-	Parameter rad_core0;
-
-	// Constructor
-	SphereSLDModel();
-
-	// 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);
-};
 
 
Index: sansmodels/src/c_models/refl.cpp
===================================================================
--- sansmodels/src/c_models/refl.cpp	(revision 2d1b7008c4fcc2791722e1581f47a3559d64f820)
+++ sansmodels/src/c_models/refl.cpp	(revision 37805e993ee77aefcc3df740a88439050390d7ce)
@@ -13,5 +13,5 @@
 #define lamda 4.62
 
-double re_kernel(double dp[], double q) {
+static double re_kernel(double dp[], double q) {
   int n = dp[0];
   int i,j;
Index: sansmodels/src/c_models/spheresld.cpp
===================================================================
--- sansmodels/src/c_models/spheresld.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/spheresld.cpp	(revision 37805e993ee77aefcc3df740a88439050390d7ce)
@@ -1,86 +1,381 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
+#include <stdlib.h>
+
+extern "C" {
+#include "libmultifunc/librefl.h"
+}
+
+#include "spheresld.h"
+
 using namespace std;
 
-extern "C" {
-	#include "spheresld.h"
+// Convenience structure
+typedef struct {
+  /// number of shells
+  //  [DEFAULT]=n_shells=1
+  int n_shells;
+  /// Scale factor
+  //  [DEFAULT]=scale= 1.0
+  double scale;
+  /// thick_inter0 [A]
+  //  [DEFAULT]=thick_inter0=50.0 [A]
+  double thick_inter0;
+  /// func_inter0
+  //  [DEFAULT]=func_inter0= 0
+  double func_inter0;
+  /// sld_core0 [1/A^(2)]
+  //  [DEFAULT]=sld_core0= 2.07e-6 [1/A^(2)]
+  double sld_core0;
+  /// sld_solv [1/A^(2)]
+  //  [DEFAULT]=sld_solv= 1.0e-6 [1/A^(2)]
+  double sld_solv;
+  /// Background
+  //  [DEFAULT]=background=0
+  double background;
+
+  //  [DEFAULT]=sld_flat1=4.0e-06 [1/A^(2)]
+  double sld_flat1;
+  //  [DEFAULT]=sld_flat2=3.5e-06 [1/A^(2)]
+  double sld_flat2;
+  //  [DEFAULT]=sld_flat3=4.0e-06 [1/A^(2)]
+  double sld_flat3;
+  //  [DEFAULT]=sld_flat4=3.5e-06 [1/A^(2)]
+  double sld_flat4;
+  //  [DEFAULT]=sld_flat5=4.0e-06 [1/A^(2)]
+  double sld_flat5;
+  //  [DEFAULT]=sld_flat6=3.5e-06 [1/A^(2)]
+  double sld_flat6;
+  //  [DEFAULT]=sld_flat7=4.0e-06 [1/A^(2)]
+  double sld_flat7;
+  //  [DEFAULT]=sld_flat8=3.5e-06 [1/A^(2)]
+  double sld_flat8;
+  //  [DEFAULT]=sld_flat9=4.0e-06 [1/A^(2)]
+  double sld_flat9;
+  //  [DEFAULT]=sld_flat10=3.5e-06 [1/A^(2)]
+  double sld_flat10;
+
+  //  [DEFAULT]=thick_inter1=50.0 [A]
+  double thick_inter1;
+  //  [DEFAULT]=thick_inter2=50.0 [A]
+  double thick_inter2;
+  //  [DEFAULT]=thick_inter3=50.0 [A]
+  double thick_inter3;
+  //  [DEFAULT]=thick_inter4=50.0 [A]
+  double thick_inter4;
+  //  [DEFAULT]=thick_inter5=50.0 [A]
+  double thick_inter5;
+  //  [DEFAULT]=thick_inter6=50.0 [A]
+  double thick_inter6;
+  //  [DEFAULT]=thick_inter7=50.0 [A]
+  double thick_inter7;
+  //  [DEFAULT]=thick_inter8=50.0 [A]
+  double thick_inter8;
+  //  [DEFAULT]=thick_inter9=50.0 [A]
+  double thick_inter9;
+  //  [DEFAULT]=thick_inter10=50.0 [A]
+  double thick_inter10;
+
+  //  [DEFAULT]=thick_flat1=100 [A]
+  double thick_flat1;
+  //  [DEFAULT]=thick_flat2=100 [A]
+  double thick_flat2;
+  //  [DEFAULT]=thick_flat3=100 [A]
+  double thick_flat3;
+  //  [DEFAULT]=thick_flat4=100 [A]
+  double thick_flat4;
+  //  [DEFAULT]=thick_flat5=100 [A]
+  double thick_flat5;
+  //  [DEFAULT]=thick_flat6=100 [A]
+  double thick_flat6;
+  //  [DEFAULT]=thick_flat7=100 [A]
+  double thick_flat7;
+  //  [DEFAULT]=thick_flat8=100 [A]
+  double thick_flat8;
+  //  [DEFAULT]=thick_flat9=100 [A]
+  double thick_flat9;
+  //  [DEFAULT]=thick_flat10=100 [A]
+  double thick_flat10;
+
+  //  [DEFAULT]=func_inter1=0
+  double func_inter1;
+  //  [DEFAULT]=func_inter2=0
+  double func_inter2;
+  //  [DEFAULT]=func_inter3=0
+  double func_inter3;
+  //  [DEFAULT]=func_inter4=0
+  double func_inter4;
+  //  [DEFAULT]=func_inter5=0
+  double func_inter5;
+  //  [DEFAULT]=func_inter6=0
+  double func_inter6;
+  //  [DEFAULT]=func_inter7=0
+  double func_inter7;
+  //  [DEFAULT]=func_inter8=0
+  double func_inter8;
+  //  [DEFAULT]=func_inter9=0
+  double func_inter9;
+  //  [DEFAULT]=func_inter10=0
+  double func_inter10;
+
+  //  [DEFAULT]=nu_inter1=2.5
+  double nu_inter1;
+  //  [DEFAULT]=nu_inter2=2.5
+  double nu_inter2;
+  //  [DEFAULT]=nu_inter3=2.5
+  double nu_inter3;
+  //  [DEFAULT]=nu_inter4=2.5
+  double nu_inter4;
+  //  [DEFAULT]=nu_inter5=2.5
+  double nu_inter5;
+  //  [DEFAULT]=nu_inter6=2.5
+  double nu_inter6;
+  //  [DEFAULT]=nu_inter7=2.5
+  double nu_inter7;
+  //  [DEFAULT]=nu_inter8=2.5
+  double nu_inter8;
+  //  [DEFAULT]=nu_inter9=2.5
+  double nu_inter9;
+  //  [DEFAULT]=nu_inter10=2.5
+  double nu_inter10;
+
+  //  [DEFAULT]=npts_inter=35.0
+  double npts_inter;
+  //  [DEFAULT]=nu_inter0=2.5
+  double nu_inter0;
+  //  [DEFAULT]=rad_core0=50.0 [A]
+  double rad_core0;
+} SphereSLDParameters;
+
+
+static double sphere_sld_kernel(double dp[], double q) {
+  int n = dp[0];
+  int i,j,k;
+
+  double scale = dp[1];
+  double thick_inter_core = dp[2];
+  double sld_core = dp[4];
+  double sld_solv = dp[5];
+  double background = dp[6];
+  double npts = dp[57]; //number of sub_layers in each interface
+  double nsl=npts;//21.0; //nsl = Num_sub_layer:  MUST ODD number in double //no other number works now
+  int n_s;
+
+  double sld_i,sld_f,dz,bes,fun,f,vol,vol_pre,vol_sub,qr,r,contr,f2;
+  double sign,slope=0.0;
+  double pi;
+
+  int* fun_type;
+  double* sld;
+  double* thick_inter;
+  double* thick;
+  double* fun_coef;
+
+  double total_thick=0.0;
+
+  fun_type = (int*)malloc((n+2)*sizeof(int));
+  sld = (double*)malloc((n+2)*sizeof(double));
+  thick_inter = (double*)malloc((n+2)*sizeof(double));
+  thick = (double*)malloc((n+2)*sizeof(double));
+  fun_coef = (double*)malloc((n+2)*sizeof(double));
+
+  fun_type[0] = dp[3];
+  fun_coef[0] = fabs(dp[58]);
+  for (i =1; i<=n; i++){
+    sld[i] = dp[i+6];
+    thick_inter[i]= dp[i+16];
+    thick[i] = dp[i+26];
+    fun_type[i] = dp[i+36];
+    fun_coef[i] = fabs(dp[i+46]);
+    total_thick += thick[i] + thick_inter[i];
+  }
+  sld[0] = sld_core;
+  sld[n+1] = sld_solv;
+  thick[0] = dp[59];
+  thick[n+1] = total_thick/5.0;
+  thick_inter[0] = thick_inter_core;
+  thick_inter[n+1] = 0.0;
+  fun_coef[n+1] = 0.0;
+
+  pi = 4.0*atan(1.0);
+  f = 0.0;
+  r = 0.0;
+  vol = 0.0;
+  vol_pre = 0.0;
+  vol_sub = 0.0;
+  sld_f = sld_core;
+
+  //floor_nsl = floor(nsl/2.0);
+
+  dz = 0.0;
+  // iteration for # of shells + core + solvent
+  for (i=0;i<=n+1; i++){
+    //iteration for N sub-layers
+    //if (fabs(thick[i]) <= 1e-24){
+    //   continue;
+    //}
+    // iteration for flat and interface
+    for (j=0;j<2;j++){
+      // iteration for sub_shells in the interface
+      // starts from #1 sub-layer
+      for (n_s=1;n_s<=nsl; n_s++){
+        // for solvent, it doesn't have an interface
+        if (i==n+1 && j==1)
+          break;
+        // for flat layers
+        if (j==0){
+          dz = thick[i];
+          sld_i = sld[i];
+          slope = 0.0;
+        }
+        // for interfacial sub_shells
+        else{
+          dz = thick_inter[i]/nsl;
+          // find sld_i at the outer boundary of sub-layer #n_s
+          sld_i = intersldfunc(fun_type[i],nsl, n_s, fun_coef[i], sld[i], sld[i+1]);
+          // calculate slope
+          slope= (sld_i -sld_f)/dz;
+        }
+        contr = sld_f-slope*r;
+        // iteration for the left and right boundary of the shells(or sub_shells)
+        for (k=0; k<2; k++){
+          // At r=0, the contribution to I is zero, so skip it.
+          if ( i == 0 && j == 0 && k == 0){
+            continue;
+          }
+          // On the top of slovent there is no interface; skip it.
+          if (i == n+1 && k == 1){
+            continue;
+          }
+          // At the right side (outer) boundary
+          if ( k == 1){
+            sign = 1.0;
+            r += dz;
+          }
+          // At the left side (inner) boundary
+          else{
+            sign = -1.0;
+          }
+          qr = q * r;
+          fun = 0.0;
+          if(qr == 0.0){
+            // sigular point
+            bes = sign * 1.0;
+          }
+          else{
+            // for flat sub-layer
+            bes = sign *  3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
+            // with linear slope
+            if (fabs(slope) > 0.0 ){
+              fun = sign * 3.0 * r * (2.0*qr*sin(qr)-((qr*qr)-2.0)*cos(qr))/(qr * qr * qr * qr);
+            }
+          }
+          // update total volume
+          vol = 4.0 * pi / 3.0 * r * r * r;
+          // we won't do the following volume correction for now.
+          // substrate empty area of volume
+          //if (k == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && fun_type[i]==0){
+          //  vol_sub += (vol_pre - vol);
+          //}
+          f += vol * (bes * contr + fun * slope);
+        }
+        // remember this sld as sld_f
+        sld_f =sld_i;
+        // no sub-layer iteration (n_s loop) for the flat layer
+        if (j==0)
+          break;
+      }
+    }
+  }
+  //vol += vol_sub;
+  f2 = f * f / vol * 1.0e8;
+  f2 *= scale;
+  f2 += background;
+
+  free(fun_type);
+  free(sld);
+  free(thick_inter);
+  free(thick);
+  free(fun_coef);
+
+  return (f2);
 }
 
+
 SphereSLDModel :: SphereSLDModel() {
-	n_shells = Parameter(1.0);
-	scale = Parameter(1.0);
-	thick_inter0 = Parameter(1.0, true);
-	thick_inter0.set_min(0.0);
-	func_inter0 = Parameter(0);
-	sld_core0 = Parameter(2.07e-06);
-	sld_solv = Parameter(1.0e-06);
-    background = Parameter(0.0);
-
-
-    sld_flat1 = Parameter(2.7e-06);
-    sld_flat2 = Parameter(3.5e-06);
-    sld_flat3 = Parameter(4.0e-06);
-    sld_flat4 = Parameter(3.5e-06);
-    sld_flat5 = Parameter(4.0e-06);
-    sld_flat6 = Parameter(3.5e-06);
-    sld_flat7 = Parameter(4.0e-06);
-    sld_flat8 = Parameter(3.5e-06);
-    sld_flat9 = Parameter(4.0e-06);
-    sld_flat10 = Parameter(3.5e-06);
-
-
-    thick_inter1 = Parameter(1.0);
-    thick_inter2 = Parameter(1.0);
-    thick_inter3 = Parameter(1.0);
-    thick_inter4 = Parameter(1.0);
-    thick_inter5 = Parameter(1.0);
-    thick_inter6 = Parameter(1.0);
-    thick_inter7 = Parameter(1.0);
-    thick_inter8 = Parameter(1.0);
-    thick_inter9 = Parameter(1.0);
-    thick_inter10 = Parameter(1.0);
-
-
-    thick_flat1 = Parameter(100.0);
-    thick_flat2 = Parameter(100.0);
-    thick_flat3 = Parameter(100.0);
-    thick_flat4 = Parameter(100.0);
-    thick_flat5 = Parameter(100.0);
-    thick_flat6 = Parameter(100.0);
-    thick_flat7 = Parameter(100.0);
-    thick_flat8 = Parameter(100.0);
-    thick_flat9 = Parameter(100.0);
-    thick_flat10 = Parameter(100.0);
-
-
-    func_inter1 = Parameter(0);
-    func_inter2 = Parameter(0);
-    func_inter3 = Parameter(0);
-    func_inter4 = Parameter(0);
-    func_inter5 = Parameter(0);
-    func_inter6 = Parameter(0);
-    func_inter7 = Parameter(0);
-    func_inter8 = Parameter(0);
-    func_inter9 = Parameter(0);
-    func_inter10 = Parameter(0);
-
-    nu_inter1 = Parameter(2.5);
-    nu_inter2 = Parameter(2.5);
-    nu_inter3 = Parameter(2.5);
-    nu_inter4 = Parameter(2.5);
-    nu_inter5 = Parameter(2.5);
-    nu_inter6 = Parameter(2.5);
-    nu_inter7 = Parameter(2.5);
-    nu_inter8 = Parameter(2.5);
-    nu_inter9 = Parameter(2.5);
-    nu_inter10 = Parameter(2.5);
-
-    npts_inter = Parameter(35.0);
-    nu_inter0 = Parameter(2.5);
-    rad_core0 = Parameter(60.0, true);
-    rad_core0.set_min(0.0);
+  n_shells = Parameter(1.0);
+  scale = Parameter(1.0);
+  thick_inter0 = Parameter(1.0, true);
+  thick_inter0.set_min(0.0);
+  func_inter0 = Parameter(0);
+  sld_core0 = Parameter(2.07e-06);
+  sld_solv = Parameter(1.0e-06);
+  background = Parameter(0.0);
+
+
+  sld_flat1 = Parameter(2.7e-06);
+  sld_flat2 = Parameter(3.5e-06);
+  sld_flat3 = Parameter(4.0e-06);
+  sld_flat4 = Parameter(3.5e-06);
+  sld_flat5 = Parameter(4.0e-06);
+  sld_flat6 = Parameter(3.5e-06);
+  sld_flat7 = Parameter(4.0e-06);
+  sld_flat8 = Parameter(3.5e-06);
+  sld_flat9 = Parameter(4.0e-06);
+  sld_flat10 = Parameter(3.5e-06);
+
+
+  thick_inter1 = Parameter(1.0);
+  thick_inter2 = Parameter(1.0);
+  thick_inter3 = Parameter(1.0);
+  thick_inter4 = Parameter(1.0);
+  thick_inter5 = Parameter(1.0);
+  thick_inter6 = Parameter(1.0);
+  thick_inter7 = Parameter(1.0);
+  thick_inter8 = Parameter(1.0);
+  thick_inter9 = Parameter(1.0);
+  thick_inter10 = Parameter(1.0);
+
+
+  thick_flat1 = Parameter(100.0);
+  thick_flat2 = Parameter(100.0);
+  thick_flat3 = Parameter(100.0);
+  thick_flat4 = Parameter(100.0);
+  thick_flat5 = Parameter(100.0);
+  thick_flat6 = Parameter(100.0);
+  thick_flat7 = Parameter(100.0);
+  thick_flat8 = Parameter(100.0);
+  thick_flat9 = Parameter(100.0);
+  thick_flat10 = Parameter(100.0);
+
+
+  func_inter1 = Parameter(0);
+  func_inter2 = Parameter(0);
+  func_inter3 = Parameter(0);
+  func_inter4 = Parameter(0);
+  func_inter5 = Parameter(0);
+  func_inter6 = Parameter(0);
+  func_inter7 = Parameter(0);
+  func_inter8 = Parameter(0);
+  func_inter9 = Parameter(0);
+  func_inter10 = Parameter(0);
+
+  nu_inter1 = Parameter(2.5);
+  nu_inter2 = Parameter(2.5);
+  nu_inter3 = Parameter(2.5);
+  nu_inter4 = Parameter(2.5);
+  nu_inter5 = Parameter(2.5);
+  nu_inter6 = Parameter(2.5);
+  nu_inter7 = Parameter(2.5);
+  nu_inter8 = Parameter(2.5);
+  nu_inter9 = Parameter(2.5);
+  nu_inter10 = Parameter(2.5);
+
+  npts_inter = Parameter(35.0);
+  nu_inter0 = Parameter(2.5);
+  rad_core0 = Parameter(60.0, true);
+  rad_core0.set_min(0.0);
 }
 
@@ -91,110 +386,110 @@
  */
 double SphereSLDModel :: operator()(double q) {
-	double dp[60];
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = n_shells();
-	dp[1] = scale();
-	dp[2] = thick_inter0();
-	dp[3] = func_inter0();
-	dp[4] = sld_core0();
-	dp[5] = sld_solv();
-	dp[6] = 0.0;
-
-	dp[7] = sld_flat1();
-	dp[8] = sld_flat2();
-	dp[9] = sld_flat3();
-	dp[10] = sld_flat4();
-	dp[11] = sld_flat5();
-	dp[12] = sld_flat6();
-	dp[13] = sld_flat7();
-	dp[14] = sld_flat8();
-	dp[15] = sld_flat9();
-	dp[16] = sld_flat10();
-
-	dp[17] = thick_inter1();
-	dp[18] = thick_inter2();
-	dp[19] = thick_inter3();
-	dp[20] = thick_inter4();
-	dp[21] = thick_inter5();
-	dp[22] = thick_inter6();
-	dp[23] = thick_inter7();
-	dp[24] = thick_inter8();
-	dp[25] = thick_inter9();
-	dp[26] = thick_inter10();
-
-	dp[27] = thick_flat1();
-	dp[28] = thick_flat2();
-	dp[29] = thick_flat3();
-	dp[30] = thick_flat4();
-	dp[31] = thick_flat5();
-	dp[32] = thick_flat6();
-	dp[33] = thick_flat7();
-	dp[34] = thick_flat8();
-	dp[35] = thick_flat9();
-	dp[36] = thick_flat10();
-
-	dp[37] = func_inter1();
-	dp[38] = func_inter2();
-	dp[39] = func_inter3();
-	dp[40] = func_inter4();
-	dp[41] = func_inter5();
-	dp[42] = func_inter6();
-	dp[43] = func_inter7();
-	dp[44] = func_inter8();
-	dp[45] = func_inter9();
-	dp[46] = func_inter10();
-
-	dp[47] = nu_inter1();
-	dp[48] = nu_inter2();
-	dp[49] = nu_inter3();
-	dp[50] = nu_inter4();
-	dp[51] = nu_inter5();
-	dp[52] = nu_inter6();
-	dp[53] = nu_inter7();
-	dp[54] = nu_inter8();
-	dp[55] = nu_inter9();
-	dp[56] = nu_inter10();
-
-
-	dp[57] = npts_inter();
-	dp[58] = nu_inter0();
-	dp[59] = rad_core0();
-
-	// No polydispersion supported in this model.
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_rad_core0;
-	rad_core0.get_weights(weights_rad_core0);
-	vector<WeightPoint> weights_thick_inter0;
-	thick_inter0.get_weights(weights_thick_inter0);
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
-
-	// Loop over core weight points
-	for(size_t i=0; i<weights_rad_core0.size(); i++) {
-		dp[59] = weights_rad_core0[i].value;
-		// Loop over thick_inter0 weight points
-		for(size_t j=0; j<weights_thick_inter0.size(); j++) {
-			dp[2] = weights_thick_inter0[j].value;
-
-			//Un-normalize Sphere by volume
-			sum += weights_rad_core0[i].weight * weights_thick_inter0[j].weight
-				* sphere_sld_kernel(dp,q) * pow((weights_rad_core0[i].value +
-						weights_thick_inter0[j].value),3.0);
-			//Find average volume
-			vol += weights_rad_core0[i].weight * weights_thick_inter0[j].weight
-				* pow((weights_rad_core0[i].value+weights_thick_inter0[j].value),3.0);
-
-			norm += weights_rad_core0[i].weight * weights_thick_inter0[j].weight;
-		}
-	}
-
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
-
-	return sum/norm + background();
+  double dp[60];
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = n_shells();
+  dp[1] = scale();
+  dp[2] = thick_inter0();
+  dp[3] = func_inter0();
+  dp[4] = sld_core0();
+  dp[5] = sld_solv();
+  dp[6] = 0.0;
+
+  dp[7] = sld_flat1();
+  dp[8] = sld_flat2();
+  dp[9] = sld_flat3();
+  dp[10] = sld_flat4();
+  dp[11] = sld_flat5();
+  dp[12] = sld_flat6();
+  dp[13] = sld_flat7();
+  dp[14] = sld_flat8();
+  dp[15] = sld_flat9();
+  dp[16] = sld_flat10();
+
+  dp[17] = thick_inter1();
+  dp[18] = thick_inter2();
+  dp[19] = thick_inter3();
+  dp[20] = thick_inter4();
+  dp[21] = thick_inter5();
+  dp[22] = thick_inter6();
+  dp[23] = thick_inter7();
+  dp[24] = thick_inter8();
+  dp[25] = thick_inter9();
+  dp[26] = thick_inter10();
+
+  dp[27] = thick_flat1();
+  dp[28] = thick_flat2();
+  dp[29] = thick_flat3();
+  dp[30] = thick_flat4();
+  dp[31] = thick_flat5();
+  dp[32] = thick_flat6();
+  dp[33] = thick_flat7();
+  dp[34] = thick_flat8();
+  dp[35] = thick_flat9();
+  dp[36] = thick_flat10();
+
+  dp[37] = func_inter1();
+  dp[38] = func_inter2();
+  dp[39] = func_inter3();
+  dp[40] = func_inter4();
+  dp[41] = func_inter5();
+  dp[42] = func_inter6();
+  dp[43] = func_inter7();
+  dp[44] = func_inter8();
+  dp[45] = func_inter9();
+  dp[46] = func_inter10();
+
+  dp[47] = nu_inter1();
+  dp[48] = nu_inter2();
+  dp[49] = nu_inter3();
+  dp[50] = nu_inter4();
+  dp[51] = nu_inter5();
+  dp[52] = nu_inter6();
+  dp[53] = nu_inter7();
+  dp[54] = nu_inter8();
+  dp[55] = nu_inter9();
+  dp[56] = nu_inter10();
+
+
+  dp[57] = npts_inter();
+  dp[58] = nu_inter0();
+  dp[59] = rad_core0();
+
+  // No polydispersion supported in this model.
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad_core0;
+  rad_core0.get_weights(weights_rad_core0);
+  vector<WeightPoint> weights_thick_inter0;
+  thick_inter0.get_weights(weights_thick_inter0);
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over core weight points
+  for(size_t i=0; i<weights_rad_core0.size(); i++) {
+    dp[59] = weights_rad_core0[i].value;
+    // Loop over thick_inter0 weight points
+    for(size_t j=0; j<weights_thick_inter0.size(); j++) {
+      dp[2] = weights_thick_inter0[j].value;
+
+      //Un-normalize Sphere by volume
+      sum += weights_rad_core0[i].weight * weights_thick_inter0[j].weight
+          * sphere_sld_kernel(dp,q) * pow((weights_rad_core0[i].value +
+              weights_thick_inter0[j].value),3.0);
+      //Find average volume
+      vol += weights_rad_core0[i].weight * weights_thick_inter0[j].weight
+          * pow((weights_rad_core0[i].value+weights_thick_inter0[j].value),3.0);
+
+      norm += weights_rad_core0[i].weight * weights_thick_inter0[j].weight;
+    }
+  }
+
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -206,6 +501,6 @@
  */
 double SphereSLDModel :: 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);
 }
 
@@ -218,5 +513,5 @@
  */
 double SphereSLDModel :: evaluate_rphi(double q, double phi) {
-	return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 
@@ -231,106 +526,106 @@
 // sld of solvent.
 double SphereSLDModel :: calculate_ER() {
-	SphereSLDParameters dp;
-
-	dp.n_shells = n_shells();
-
-	dp.rad_core0 = rad_core0();
-	dp.thick_flat1 = thick_flat1();
-	dp.thick_flat2 = thick_flat2();
-	dp.thick_flat3 = thick_flat3();
-	dp.thick_flat4 = thick_flat4();
-	dp.thick_flat5 = thick_flat5();
-	dp.thick_flat6 = thick_flat6();
-	dp.thick_flat7 = thick_flat7();
-	dp.thick_flat8 = thick_flat8();
-	dp.thick_flat9 = thick_flat9();
-	dp.thick_flat10 = thick_flat10();
-
-	dp.thick_inter0 = thick_inter0();
-	dp.thick_inter1 = thick_inter1();
-	dp.thick_inter2 = thick_inter2();
-	dp.thick_inter3 = thick_inter3();
-	dp.thick_inter4 = thick_inter4();
-	dp.thick_inter5 = thick_inter5();
-	dp.thick_inter6 = thick_inter6();
-	dp.thick_inter7 = thick_inter7();
-	dp.thick_inter8 = thick_inter8();
-	dp.thick_inter9 = thick_inter9();
-	dp.thick_inter10 = thick_inter10();
-
-	double rad_out = 0.0;
-	double out = 0.0;
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_rad_core0;
-	rad_core0.get_weights(weights_rad_core0);
-
-	// Get the dispersion points for the thick 1
-	vector<WeightPoint> weights_thick_inter0;
-	thick_inter0.get_weights(weights_thick_inter0);
-	// Loop over radius weight points
-	for(size_t i=0; i<weights_rad_core0.size(); i++) {
-		dp.rad_core0 = weights_rad_core0[i].value;
-		// Loop over radius weight points
-		for(size_t j=0; j<weights_thick_inter0.size(); j++) {
-			dp.thick_inter0 = weights_thick_inter0[j].value;
-			rad_out = dp.rad_core0 + dp.thick_inter0;
-			if (dp.n_shells > 0)
-				rad_out += dp.thick_flat1 + dp.thick_inter1;
-			if (dp.n_shells > 1)
-				rad_out += dp.thick_flat2 + dp.thick_inter2;
-			if (dp.n_shells > 2)
-				rad_out += dp.thick_flat3 + dp.thick_inter3;
-			if (dp.n_shells > 3)
-				rad_out += dp.thick_flat4 + dp.thick_inter4;
-			if (dp.n_shells > 4)
-				rad_out += dp.thick_flat5 + dp.thick_inter5;
-			if (dp.n_shells > 5)
-				rad_out += dp.thick_flat6 + dp.thick_inter6;
-			if (dp.n_shells > 6)
-				rad_out += dp.thick_flat7 + dp.thick_inter7;
-			if (dp.n_shells > 7)
-				rad_out += dp.thick_flat8 + dp.thick_inter8;
-			if (dp.n_shells > 8)
-				rad_out += dp.thick_flat9 + dp.thick_inter9;
-			if (dp.n_shells > 9)
-				rad_out += dp.thick_flat10 + dp.thick_inter10;
-			sum += weights_rad_core0[i].weight*weights_thick_inter0[j].weight
-				* (rad_out);
-			norm += weights_rad_core0[i].weight*weights_thick_inter0[j].weight;
-		}
-	}
-	if (norm != 0){
-		//return the averaged value
-		out =  sum/norm;}
-	else{
-		//return normal value
-		out = dp.rad_core0 + dp.thick_inter0;
-		if (dp.n_shells > 0)
-			out += dp.thick_flat1 + dp.thick_inter1;
-		if (dp.n_shells > 1)
-			out += dp.thick_flat2 + dp.thick_inter2;
-		if (dp.n_shells > 2)
-			out += dp.thick_flat3 + dp.thick_inter3;
-		if (dp.n_shells > 3)
-			out += dp.thick_flat4 + dp.thick_inter4;
-		if (dp.n_shells > 4)
-			out += dp.thick_flat5 + dp.thick_inter5;
-		if (dp.n_shells > 5)
-			out += dp.thick_flat6 + dp.thick_inter6;
-		if (dp.n_shells > 6)
-			out += dp.thick_flat7 + dp.thick_inter7;
-		if (dp.n_shells > 7)
-			out += dp.thick_flat8 + dp.thick_inter8;
-		if (dp.n_shells > 8)
-			out += dp.thick_flat9 + dp.thick_inter9;
-		if (dp.n_shells > 9)
-			out += dp.thick_flat10 + dp.thick_inter10;
-		}
-
-	return out;
+  SphereSLDParameters dp;
+
+  dp.n_shells = n_shells();
+
+  dp.rad_core0 = rad_core0();
+  dp.thick_flat1 = thick_flat1();
+  dp.thick_flat2 = thick_flat2();
+  dp.thick_flat3 = thick_flat3();
+  dp.thick_flat4 = thick_flat4();
+  dp.thick_flat5 = thick_flat5();
+  dp.thick_flat6 = thick_flat6();
+  dp.thick_flat7 = thick_flat7();
+  dp.thick_flat8 = thick_flat8();
+  dp.thick_flat9 = thick_flat9();
+  dp.thick_flat10 = thick_flat10();
+
+  dp.thick_inter0 = thick_inter0();
+  dp.thick_inter1 = thick_inter1();
+  dp.thick_inter2 = thick_inter2();
+  dp.thick_inter3 = thick_inter3();
+  dp.thick_inter4 = thick_inter4();
+  dp.thick_inter5 = thick_inter5();
+  dp.thick_inter6 = thick_inter6();
+  dp.thick_inter7 = thick_inter7();
+  dp.thick_inter8 = thick_inter8();
+  dp.thick_inter9 = thick_inter9();
+  dp.thick_inter10 = thick_inter10();
+
+  double rad_out = 0.0;
+  double out = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad_core0;
+  rad_core0.get_weights(weights_rad_core0);
+
+  // Get the dispersion points for the thick 1
+  vector<WeightPoint> weights_thick_inter0;
+  thick_inter0.get_weights(weights_thick_inter0);
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad_core0.size(); i++) {
+    dp.rad_core0 = weights_rad_core0[i].value;
+    // Loop over radius weight points
+    for(size_t j=0; j<weights_thick_inter0.size(); j++) {
+      dp.thick_inter0 = weights_thick_inter0[j].value;
+      rad_out = dp.rad_core0 + dp.thick_inter0;
+      if (dp.n_shells > 0)
+        rad_out += dp.thick_flat1 + dp.thick_inter1;
+      if (dp.n_shells > 1)
+        rad_out += dp.thick_flat2 + dp.thick_inter2;
+      if (dp.n_shells > 2)
+        rad_out += dp.thick_flat3 + dp.thick_inter3;
+      if (dp.n_shells > 3)
+        rad_out += dp.thick_flat4 + dp.thick_inter4;
+      if (dp.n_shells > 4)
+        rad_out += dp.thick_flat5 + dp.thick_inter5;
+      if (dp.n_shells > 5)
+        rad_out += dp.thick_flat6 + dp.thick_inter6;
+      if (dp.n_shells > 6)
+        rad_out += dp.thick_flat7 + dp.thick_inter7;
+      if (dp.n_shells > 7)
+        rad_out += dp.thick_flat8 + dp.thick_inter8;
+      if (dp.n_shells > 8)
+        rad_out += dp.thick_flat9 + dp.thick_inter9;
+      if (dp.n_shells > 9)
+        rad_out += dp.thick_flat10 + dp.thick_inter10;
+      sum += weights_rad_core0[i].weight*weights_thick_inter0[j].weight
+          * (rad_out);
+      norm += weights_rad_core0[i].weight*weights_thick_inter0[j].weight;
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    out =  sum/norm;}
+  else{
+    //return normal value
+    out = dp.rad_core0 + dp.thick_inter0;
+    if (dp.n_shells > 0)
+      out += dp.thick_flat1 + dp.thick_inter1;
+    if (dp.n_shells > 1)
+      out += dp.thick_flat2 + dp.thick_inter2;
+    if (dp.n_shells > 2)
+      out += dp.thick_flat3 + dp.thick_inter3;
+    if (dp.n_shells > 3)
+      out += dp.thick_flat4 + dp.thick_inter4;
+    if (dp.n_shells > 4)
+      out += dp.thick_flat5 + dp.thick_inter5;
+    if (dp.n_shells > 5)
+      out += dp.thick_flat6 + dp.thick_inter6;
+    if (dp.n_shells > 6)
+      out += dp.thick_flat7 + dp.thick_inter7;
+    if (dp.n_shells > 7)
+      out += dp.thick_flat8 + dp.thick_inter8;
+    if (dp.n_shells > 8)
+      out += dp.thick_flat9 + dp.thick_inter9;
+    if (dp.n_shells > 9)
+      out += dp.thick_flat10 + dp.thick_inter10;
+  }
+
+  return out;
 
 }
Index: nsmodels/src/python_wrapper/CSphereSLDModel.cpp
===================================================================
--- sansmodels/src/python_wrapper/CSphereSLDModel.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ 	(revision )
@@ -1,844 +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
- */
-
-/** CSphereSLDModel
- *
- * C extension 
- *
- * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
- *          DO NOT MODIFY THIS FILE, MODIFY spheresld.h
- *          AND RE-RUN THE GENERATOR SCRIPT
- *
- */
-#define NO_IMPORT_ARRAY
-#define PY_ARRAY_UNIQUE_SYMBOL PyArray_API_sans
- 
-extern "C" {
-#include <Python.h>
-#include <arrayobject.h>
-#include "structmember.h"
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <time.h>
-#include "spheresld.h"
-}
-
-#include "models.hh"
-#include "dispersion_visitor.hh"
-
-/// Error object for raised exceptions
-static PyObject * CSphereSLDModelError = NULL;
-
-
-// Class definition
-typedef struct {
-    PyObject_HEAD
-    /// Parameters
-    PyObject * params;
-    /// Dispersion parameters
-    PyObject * dispersion;
-    /// Underlying model object
-    SphereSLDModel * model;
-    /// Log for unit testing
-    PyObject * log;
-} CSphereSLDModel;
-
-
-static void
-CSphereSLDModel_dealloc(CSphereSLDModel* self)
-{
-    Py_DECREF(self->params);
-    Py_DECREF(self->dispersion);
-    Py_DECREF(self->log);
-    delete self->model;
-    self->ob_type->tp_free((PyObject*)self);
-    
-
-}
-
-static PyObject *
-CSphereSLDModel_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
-{
-    CSphereSLDModel *self;
-    
-    self = (CSphereSLDModel *)type->tp_alloc(type, 0);
-   
-    return (PyObject *)self;
-}
-
-static int
-CSphereSLDModel_init(CSphereSLDModel *self, PyObject *args, PyObject *kwds)
-{
-    if (self != NULL) {
-    	
-    	// Create parameters
-        self->params = PyDict_New();
-        self->dispersion = PyDict_New();
-        self->model = new SphereSLDModel();
-        
-        // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"thick_flat8",Py_BuildValue("d",100.000000000000));
-        PyDict_SetItemString(self->params,"thick_flat9",Py_BuildValue("d",100.000000000000));
-        PyDict_SetItemString(self->params,"thick_flat2",Py_BuildValue("d",100.000000000000));
-        PyDict_SetItemString(self->params,"thick_flat3",Py_BuildValue("d",100.000000000000));
-        PyDict_SetItemString(self->params,"thick_flat1",Py_BuildValue("d",100.000000000000));
-        PyDict_SetItemString(self->params,"thick_flat6",Py_BuildValue("d",100.000000000000));
-        PyDict_SetItemString(self->params,"thick_flat7",Py_BuildValue("d",100.000000000000));
-        PyDict_SetItemString(self->params,"thick_flat4",Py_BuildValue("d",100.000000000000));
-        PyDict_SetItemString(self->params,"thick_flat5",Py_BuildValue("d",100.000000000000));
-        PyDict_SetItemString(self->params,"thick_inter0",Py_BuildValue("d",50.000000000000));
-        PyDict_SetItemString(self->params,"thick_inter1",Py_BuildValue("d",50.000000000000));
-        PyDict_SetItemString(self->params,"thick_inter2",Py_BuildValue("d",50.000000000000));
-        PyDict_SetItemString(self->params,"thick_inter3",Py_BuildValue("d",50.000000000000));
-        PyDict_SetItemString(self->params,"sld_core0",Py_BuildValue("d",0.000002070000));
-        PyDict_SetItemString(self->params,"thick_inter5",Py_BuildValue("d",50.000000000000));
-        PyDict_SetItemString(self->params,"thick_inter6",Py_BuildValue("d",50.000000000000));
-        PyDict_SetItemString(self->params,"thick_inter7",Py_BuildValue("d",50.000000000000));
-        PyDict_SetItemString(self->params,"thick_inter8",Py_BuildValue("d",50.000000000000));
-        PyDict_SetItemString(self->params,"thick_inter9",Py_BuildValue("d",50.000000000000));
-        PyDict_SetItemString(self->params,"rad_core0",Py_BuildValue("d",50.000000000000));
-        PyDict_SetItemString(self->params,"npts_inter",Py_BuildValue("d",35.000000000000));
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000000000));
-        PyDict_SetItemString(self->params,"thick_inter10",Py_BuildValue("d",50.000000000000));
-        PyDict_SetItemString(self->params,"func_inter8",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"func_inter9",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"sld_solv",Py_BuildValue("d",0.000001000000));
-        PyDict_SetItemString(self->params,"func_inter0",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"func_inter1",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"func_inter2",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"func_inter3",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"thick_flat10",Py_BuildValue("d",100.000000000000));
-        PyDict_SetItemString(self->params,"func_inter5",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"func_inter6",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"func_inter7",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"sld_flat4",Py_BuildValue("d",0.000003500000));
-        PyDict_SetItemString(self->params,"sld_flat5",Py_BuildValue("d",0.000004000000));
-        PyDict_SetItemString(self->params,"sld_flat6",Py_BuildValue("d",0.000003500000));
-        PyDict_SetItemString(self->params,"sld_flat7",Py_BuildValue("d",0.000004000000));
-        PyDict_SetItemString(self->params,"sld_flat1",Py_BuildValue("d",0.000004000000));
-        PyDict_SetItemString(self->params,"sld_flat2",Py_BuildValue("d",0.000003500000));
-        PyDict_SetItemString(self->params,"sld_flat3",Py_BuildValue("d",0.000004000000));
-        PyDict_SetItemString(self->params,"sld_flat8",Py_BuildValue("d",0.000003500000));
-        PyDict_SetItemString(self->params,"sld_flat9",Py_BuildValue("d",0.000004000000));
-        PyDict_SetItemString(self->params,"nu_inter10",Py_BuildValue("d",2.500000000000));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"n_shells",Py_BuildValue("d",1.000000000000));
-        PyDict_SetItemString(self->params,"sld_flat10",Py_BuildValue("d",0.000003500000));
-        PyDict_SetItemString(self->params,"func_inter4",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"nu_inter8",Py_BuildValue("d",2.500000000000));
-        PyDict_SetItemString(self->params,"nu_inter5",Py_BuildValue("d",2.500000000000));
-        PyDict_SetItemString(self->params,"nu_inter4",Py_BuildValue("d",2.500000000000));
-        PyDict_SetItemString(self->params,"nu_inter7",Py_BuildValue("d",2.500000000000));
-        PyDict_SetItemString(self->params,"nu_inter6",Py_BuildValue("d",2.500000000000));
-        PyDict_SetItemString(self->params,"nu_inter1",Py_BuildValue("d",2.500000000000));
-        PyDict_SetItemString(self->params,"nu_inter0",Py_BuildValue("d",2.500000000000));
-        PyDict_SetItemString(self->params,"nu_inter3",Py_BuildValue("d",2.500000000000));
-        PyDict_SetItemString(self->params,"nu_inter2",Py_BuildValue("d",2.500000000000));
-        PyDict_SetItemString(self->params,"thick_inter4",Py_BuildValue("d",50.000000000000));
-        PyDict_SetItemString(self->params,"nu_inter9",Py_BuildValue("d",2.500000000000));
-        PyDict_SetItemString(self->params,"func_inter10",Py_BuildValue("d",0.000000000000));
-        // Initialize dispersion / averaging parameter dict
-        DispersionVisitor* visitor = new DispersionVisitor();
-        PyObject * disp_dict;
-        disp_dict = PyDict_New();
-        self->model->rad_core0.dispersion->accept_as_source(visitor, self->model->rad_core0.dispersion, disp_dict);
-        PyDict_SetItemString(self->dispersion, "rad_core0", disp_dict);
-        disp_dict = PyDict_New();
-        self->model->thick_inter0.dispersion->accept_as_source(visitor, self->model->thick_inter0.dispersion, disp_dict);
-        PyDict_SetItemString(self->dispersion, "thick_inter0", disp_dict);
-
-
-         
-        // Create empty log
-        self->log = PyDict_New();
-        
-        
-
-    }
-    return 0;
-}
-
-static char name_params[] = "params";
-static char def_params[] = "Parameters";
-static char name_dispersion[] = "dispersion";
-static char def_dispersion[] = "Dispersion parameters";
-static char name_log[] = "log";
-static char def_log[] = "Log";
-
-static PyMemberDef CSphereSLDModel_members[] = {
-    {name_params, T_OBJECT, offsetof(CSphereSLDModel, params), 0, def_params},
-	{name_dispersion, T_OBJECT, offsetof(CSphereSLDModel, dispersion), 0, def_dispersion},     
-    {name_log, T_OBJECT, offsetof(CSphereSLDModel, log), 0, def_log},
-    {NULL}  /* Sentinel */
-};
-
-/** Read double from PyObject
-    @param p PyObject
-    @return double
-*/
-double CSphereSLDModel_readDouble(PyObject *p) {
-    if (PyFloat_Check(p)==1) {
-        return (double)(((PyFloatObject *)(p))->ob_fval);
-    } else if (PyInt_Check(p)==1) {
-        return (double)(((PyIntObject *)(p))->ob_ival);
-    } else if (PyLong_Check(p)==1) {
-        return (double)PyLong_AsLong(p);
-    } else {
-        return 0.0;
-    }
-}
-/**
- * Function to call to evaluate model
- * @param args: input numpy array q[] 
- * @return: numpy array object 
- */
- 
-static PyObject *evaluateOneDim(SphereSLDModel* model, PyArrayObject *q){
-    PyArrayObject *result;
-   
-    // Check validity of array q , q must be of dimension 1, an array of double
-    if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE)
-    {
-        //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num;
-        //PyErr_SetString(PyExc_ValueError , message);
-        return NULL;
-    }
-    result = (PyArrayObject *)PyArray_FromDims(q->nd, (int *)(q->dimensions), PyArray_DOUBLE);
-	if (result == NULL) {
-        const char * message= "Could not create result ";
-        PyErr_SetString(PyExc_RuntimeError , message);
-		return NULL;
-	}
-#pragma omp parallel for
-	 for (int i = 0; i < q->dimensions[0]; i++){
-      double q_value  = *(double *)(q->data + i*q->strides[0]);
-      double *result_value = (double *)(result->data + i*result->strides[0]);
-      *result_value =(*model)(q_value);
-	}
-    return PyArray_Return(result); 
- }
-
- /**
- * Function to call to evaluate model
- * @param args: input numpy array  [x[],y[]]
- * @return: numpy array object 
- */
- static PyObject * evaluateTwoDimXY( SphereSLDModel* model, 
-                              PyArrayObject *x, PyArrayObject *y)
- {
-    PyArrayObject *result;
-    int x_len, y_len, dims[1];
-    //check validity of input vectors
-    if (x->nd != 1 || x->descr->type_num != PyArray_DOUBLE
-        || y->nd != 1 || y->descr->type_num != PyArray_DOUBLE
-        || y->dimensions[0] != x->dimensions[0]){
-        const char * message= "evaluateTwoDimXY  expect 2 numpy arrays";
-        PyErr_SetString(PyExc_ValueError , message); 
-        return NULL;
-    }
-   
-	if (PyArray_Check(x) && PyArray_Check(y)) {
-		
-	    x_len = dims[0]= x->dimensions[0];
-        y_len = dims[0]= y->dimensions[0];
-	    
-	    // Make a new double matrix of same dims
-        result=(PyArrayObject *) PyArray_FromDims(1,dims,NPY_DOUBLE);
-        if (result == NULL){
-	    const char * message= "Could not create result ";
-        PyErr_SetString(PyExc_RuntimeError , message);
-	    return NULL;
-	    }
-       
-        /* Do the calculation. */
-#pragma omp parallel for
-        for (int i=0; i< x_len; i++) {
-            double x_value = *(double *)(x->data + i*x->strides[0]);
-  		    double y_value = *(double *)(y->data + i*y->strides[0]);
-  			double *result_value = (double *)(result->data +
-  			      i*result->strides[0]);
-  			*result_value = (*model)(x_value, y_value);
-        }           
-        return PyArray_Return(result); 
-        
-        }else{
-		    PyErr_SetString(CSphereSLDModelError, 
-                   "CSphereSLDModel.evaluateTwoDimXY couldn't run.");
-	        return NULL;
-		}      	
-}
-/**
- *  evalDistribution function evaluate a model function with input vector
- *  @param args: input q as vector or [qx, qy] where qx, qy are vectors
- *
- */ 
-static PyObject * evalDistribution(CSphereSLDModel *self, PyObject *args){
-	PyObject *qx, *qy;
-	PyArrayObject * pars;
-	int npars ,mpars;
-	
-	// Get parameters
-	
-	    // Reader parameter dictionary
-    self->model->thick_flat8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat8") );
-    self->model->thick_flat9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat9") );
-    self->model->thick_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat2") );
-    self->model->thick_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat3") );
-    self->model->thick_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat1") );
-    self->model->thick_flat6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat6") );
-    self->model->thick_flat7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat7") );
-    self->model->thick_flat4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat4") );
-    self->model->thick_flat5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat5") );
-    self->model->thick_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter0") );
-    self->model->thick_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter1") );
-    self->model->thick_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter2") );
-    self->model->thick_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter3") );
-    self->model->sld_core0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_core0") );
-    self->model->thick_inter5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter5") );
-    self->model->thick_inter6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter6") );
-    self->model->thick_inter7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter7") );
-    self->model->thick_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter8") );
-    self->model->thick_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter9") );
-    self->model->rad_core0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rad_core0") );
-    self->model->npts_inter = PyFloat_AsDouble( PyDict_GetItemString(self->params, "npts_inter") );
-    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
-    self->model->thick_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter10") );
-    self->model->func_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter8") );
-    self->model->func_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter9") );
-    self->model->sld_solv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_solv") );
-    self->model->func_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter0") );
-    self->model->func_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter1") );
-    self->model->func_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter2") );
-    self->model->func_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter3") );
-    self->model->thick_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat10") );
-    self->model->func_inter5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter5") );
-    self->model->func_inter6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter6") );
-    self->model->func_inter7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter7") );
-    self->model->sld_flat4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat4") );
-    self->model->sld_flat5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat5") );
-    self->model->sld_flat6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat6") );
-    self->model->sld_flat7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat7") );
-    self->model->sld_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat1") );
-    self->model->sld_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat2") );
-    self->model->sld_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat3") );
-    self->model->sld_flat8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat8") );
-    self->model->sld_flat9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat9") );
-    self->model->nu_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter10") );
-    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
-    self->model->n_shells = PyFloat_AsDouble( PyDict_GetItemString(self->params, "n_shells") );
-    self->model->sld_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat10") );
-    self->model->func_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter4") );
-    self->model->nu_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter8") );
-    self->model->nu_inter5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter5") );
-    self->model->nu_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter4") );
-    self->model->nu_inter7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter7") );
-    self->model->nu_inter6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter6") );
-    self->model->nu_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter1") );
-    self->model->nu_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter0") );
-    self->model->nu_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter3") );
-    self->model->nu_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter2") );
-    self->model->thick_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter4") );
-    self->model->nu_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter9") );
-    self->model->func_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter10") );
-    // Read in dispersion parameters
-    PyObject* disp_dict;
-    DispersionVisitor* visitor = new DispersionVisitor();
-    disp_dict = PyDict_GetItemString(self->dispersion, "rad_core0");
-    self->model->rad_core0.dispersion->accept_as_destination(visitor, self->model->rad_core0.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "thick_inter0");
-    self->model->thick_inter0.dispersion->accept_as_destination(visitor, self->model->thick_inter0.dispersion, disp_dict);
-
-	
-	// Get input and determine whether we have to supply a 1D or 2D return value.
-	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
-	    PyErr_SetString(CSphereSLDModelError, 
-	    	"CSphereSLDModel.evalDistribution expects a q value.");
-		return NULL;
-	}
-    // Check params
-	
-    if(PyArray_Check(pars)==1) {
-		
-	    // Length of list should 1 or 2
-	    npars = pars->nd; 
-	    if(npars==1) {
-	        // input is a numpy array
-	        if (PyArray_Check(pars)) {
-		        return evaluateOneDim(self->model, (PyArrayObject*)pars); 
-		    }
-		}else{
-		    PyErr_SetString(CSphereSLDModelError, 
-                   "CSphereSLDModel.evalDistribution expect numpy array of one dimension.");
-	        return NULL;
-		}
-    }else if( PyList_Check(pars)==1) {
-    	// Length of list should be 2 for I(qx,qy)
-	    mpars = PyList_GET_SIZE(pars); 
-	    if(mpars!=2) {
-	    	PyErr_SetString(CSphereSLDModelError, 
-	    		"CSphereSLDModel.evalDistribution expects a list of dimension 2.");
-	    	return NULL;
-	    }
-	     qx = PyList_GET_ITEM(pars,0);
-	     qy = PyList_GET_ITEM(pars,1);
-	     if (PyArray_Check(qx) && PyArray_Check(qy)) {
-	         return evaluateTwoDimXY(self->model, (PyArrayObject*)qx,
-		           (PyArrayObject*)qy);
-		 }else{
-		    PyErr_SetString(CSphereSLDModelError, 
-                   "CSphereSLDModel.evalDistribution expect 2 numpy arrays in list.");
-	        return NULL;
-	     }
-	}
-	PyErr_SetString(CSphereSLDModelError, 
-                   "CSphereSLDModel.evalDistribution couln't be run.");
-	return NULL;
-	
-}
-
-/**
- * Function to call to evaluate model
- * @param args: input q or [q,phi]
- * @return: function value
- */
-static PyObject * run(CSphereSLDModel *self, PyObject *args) {
-	double q_value, phi_value;
-	PyObject* pars;
-	int npars;
-	
-	// Get parameters
-	
-	    // Reader parameter dictionary
-    self->model->thick_flat8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat8") );
-    self->model->thick_flat9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat9") );
-    self->model->thick_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat2") );
-    self->model->thick_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat3") );
-    self->model->thick_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat1") );
-    self->model->thick_flat6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat6") );
-    self->model->thick_flat7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat7") );
-    self->model->thick_flat4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat4") );
-    self->model->thick_flat5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat5") );
-    self->model->thick_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter0") );
-    self->model->thick_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter1") );
-    self->model->thick_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter2") );
-    self->model->thick_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter3") );
-    self->model->sld_core0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_core0") );
-    self->model->thick_inter5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter5") );
-    self->model->thick_inter6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter6") );
-    self->model->thick_inter7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter7") );
-    self->model->thick_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter8") );
-    self->model->thick_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter9") );
-    self->model->rad_core0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rad_core0") );
-    self->model->npts_inter = PyFloat_AsDouble( PyDict_GetItemString(self->params, "npts_inter") );
-    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
-    self->model->thick_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter10") );
-    self->model->func_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter8") );
-    self->model->func_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter9") );
-    self->model->sld_solv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_solv") );
-    self->model->func_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter0") );
-    self->model->func_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter1") );
-    self->model->func_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter2") );
-    self->model->func_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter3") );
-    self->model->thick_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat10") );
-    self->model->func_inter5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter5") );
-    self->model->func_inter6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter6") );
-    self->model->func_inter7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter7") );
-    self->model->sld_flat4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat4") );
-    self->model->sld_flat5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat5") );
-    self->model->sld_flat6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat6") );
-    self->model->sld_flat7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat7") );
-    self->model->sld_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat1") );
-    self->model->sld_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat2") );
-    self->model->sld_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat3") );
-    self->model->sld_flat8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat8") );
-    self->model->sld_flat9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat9") );
-    self->model->nu_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter10") );
-    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
-    self->model->n_shells = PyFloat_AsDouble( PyDict_GetItemString(self->params, "n_shells") );
-    self->model->sld_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat10") );
-    self->model->func_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter4") );
-    self->model->nu_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter8") );
-    self->model->nu_inter5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter5") );
-    self->model->nu_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter4") );
-    self->model->nu_inter7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter7") );
-    self->model->nu_inter6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter6") );
-    self->model->nu_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter1") );
-    self->model->nu_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter0") );
-    self->model->nu_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter3") );
-    self->model->nu_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter2") );
-    self->model->thick_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter4") );
-    self->model->nu_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter9") );
-    self->model->func_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter10") );
-    // Read in dispersion parameters
-    PyObject* disp_dict;
-    DispersionVisitor* visitor = new DispersionVisitor();
-    disp_dict = PyDict_GetItemString(self->dispersion, "rad_core0");
-    self->model->rad_core0.dispersion->accept_as_destination(visitor, self->model->rad_core0.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "thick_inter0");
-    self->model->thick_inter0.dispersion->accept_as_destination(visitor, self->model->thick_inter0.dispersion, disp_dict);
-
-	
-	// Get input and determine whether we have to supply a 1D or 2D return value.
-	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
-	    PyErr_SetString(CSphereSLDModelError, 
-	    	"CSphereSLDModel.run expects a q value.");
-		return NULL;
-	}
-	  
-	// Check params
-	if( PyList_Check(pars)==1) {
-		
-		// Length of list should be 2 for I(q,phi)
-	    npars = PyList_GET_SIZE(pars); 
-	    if(npars!=2) {
-	    	PyErr_SetString(CSphereSLDModelError, 
-	    		"CSphereSLDModel.run expects a double or a list of dimension 2.");
-	    	return NULL;
-	    }
-	    // We have a vector q, get the q and phi values at which
-	    // to evaluate I(q,phi)
-	    q_value = CSphereSLDModel_readDouble(PyList_GET_ITEM(pars,0));
-	    phi_value = CSphereSLDModel_readDouble(PyList_GET_ITEM(pars,1));
-	    // Skip zero
-	    if (q_value==0) {
-	    	return Py_BuildValue("d",0.0);
-	    }
-		return Py_BuildValue("d",(*(self->model)).evaluate_rphi(q_value,phi_value));
-
-	} else {
-
-		// We have a scalar q, we will evaluate I(q)
-		q_value = CSphereSLDModel_readDouble(pars);		
-		
-		return Py_BuildValue("d",(*(self->model))(q_value));
-	}	
-}
-/**
- * Function to call to calculate_ER
- * @return: effective radius value 
- */
-static PyObject * calculate_ER(CSphereSLDModel *self) {
-
-	// Get parameters
-	
-	    // Reader parameter dictionary
-    self->model->thick_flat8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat8") );
-    self->model->thick_flat9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat9") );
-    self->model->thick_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat2") );
-    self->model->thick_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat3") );
-    self->model->thick_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat1") );
-    self->model->thick_flat6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat6") );
-    self->model->thick_flat7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat7") );
-    self->model->thick_flat4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat4") );
-    self->model->thick_flat5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat5") );
-    self->model->thick_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter0") );
-    self->model->thick_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter1") );
-    self->model->thick_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter2") );
-    self->model->thick_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter3") );
-    self->model->sld_core0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_core0") );
-    self->model->thick_inter5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter5") );
-    self->model->thick_inter6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter6") );
-    self->model->thick_inter7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter7") );
-    self->model->thick_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter8") );
-    self->model->thick_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter9") );
-    self->model->rad_core0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rad_core0") );
-    self->model->npts_inter = PyFloat_AsDouble( PyDict_GetItemString(self->params, "npts_inter") );
-    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
-    self->model->thick_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter10") );
-    self->model->func_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter8") );
-    self->model->func_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter9") );
-    self->model->sld_solv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_solv") );
-    self->model->func_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter0") );
-    self->model->func_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter1") );
-    self->model->func_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter2") );
-    self->model->func_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter3") );
-    self->model->thick_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat10") );
-    self->model->func_inter5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter5") );
-    self->model->func_inter6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter6") );
-    self->model->func_inter7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter7") );
-    self->model->sld_flat4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat4") );
-    self->model->sld_flat5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat5") );
-    self->model->sld_flat6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat6") );
-    self->model->sld_flat7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat7") );
-    self->model->sld_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat1") );
-    self->model->sld_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat2") );
-    self->model->sld_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat3") );
-    self->model->sld_flat8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat8") );
-    self->model->sld_flat9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat9") );
-    self->model->nu_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter10") );
-    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
-    self->model->n_shells = PyFloat_AsDouble( PyDict_GetItemString(self->params, "n_shells") );
-    self->model->sld_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat10") );
-    self->model->func_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter4") );
-    self->model->nu_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter8") );
-    self->model->nu_inter5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter5") );
-    self->model->nu_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter4") );
-    self->model->nu_inter7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter7") );
-    self->model->nu_inter6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter6") );
-    self->model->nu_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter1") );
-    self->model->nu_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter0") );
-    self->model->nu_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter3") );
-    self->model->nu_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter2") );
-    self->model->thick_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter4") );
-    self->model->nu_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter9") );
-    self->model->func_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter10") );
-    // Read in dispersion parameters
-    PyObject* disp_dict;
-    DispersionVisitor* visitor = new DispersionVisitor();
-    disp_dict = PyDict_GetItemString(self->dispersion, "rad_core0");
-    self->model->rad_core0.dispersion->accept_as_destination(visitor, self->model->rad_core0.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "thick_inter0");
-    self->model->thick_inter0.dispersion->accept_as_destination(visitor, self->model->thick_inter0.dispersion, disp_dict);
-
-		
-	return Py_BuildValue("d",(*(self->model)).calculate_ER());
-
-}
-/**
- * Function to call to evaluate model in cartesian coordinates
- * @param args: input q or [qx, qy]]
- * @return: function value
- */
-static PyObject * runXY(CSphereSLDModel *self, PyObject *args) {
-	double qx_value, qy_value;
-	PyObject* pars;
-	int npars;
-	
-	// Get parameters
-	
-	    // Reader parameter dictionary
-    self->model->thick_flat8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat8") );
-    self->model->thick_flat9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat9") );
-    self->model->thick_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat2") );
-    self->model->thick_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat3") );
-    self->model->thick_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat1") );
-    self->model->thick_flat6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat6") );
-    self->model->thick_flat7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat7") );
-    self->model->thick_flat4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat4") );
-    self->model->thick_flat5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat5") );
-    self->model->thick_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter0") );
-    self->model->thick_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter1") );
-    self->model->thick_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter2") );
-    self->model->thick_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter3") );
-    self->model->sld_core0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_core0") );
-    self->model->thick_inter5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter5") );
-    self->model->thick_inter6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter6") );
-    self->model->thick_inter7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter7") );
-    self->model->thick_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter8") );
-    self->model->thick_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter9") );
-    self->model->rad_core0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rad_core0") );
-    self->model->npts_inter = PyFloat_AsDouble( PyDict_GetItemString(self->params, "npts_inter") );
-    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
-    self->model->thick_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter10") );
-    self->model->func_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter8") );
-    self->model->func_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter9") );
-    self->model->sld_solv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_solv") );
-    self->model->func_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter0") );
-    self->model->func_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter1") );
-    self->model->func_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter2") );
-    self->model->func_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter3") );
-    self->model->thick_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat10") );
-    self->model->func_inter5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter5") );
-    self->model->func_inter6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter6") );
-    self->model->func_inter7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter7") );
-    self->model->sld_flat4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat4") );
-    self->model->sld_flat5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat5") );
-    self->model->sld_flat6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat6") );
-    self->model->sld_flat7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat7") );
-    self->model->sld_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat1") );
-    self->model->sld_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat2") );
-    self->model->sld_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat3") );
-    self->model->sld_flat8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat8") );
-    self->model->sld_flat9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat9") );
-    self->model->nu_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter10") );
-    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
-    self->model->n_shells = PyFloat_AsDouble( PyDict_GetItemString(self->params, "n_shells") );
-    self->model->sld_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat10") );
-    self->model->func_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter4") );
-    self->model->nu_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter8") );
-    self->model->nu_inter5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter5") );
-    self->model->nu_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter4") );
-    self->model->nu_inter7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter7") );
-    self->model->nu_inter6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter6") );
-    self->model->nu_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter1") );
-    self->model->nu_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter0") );
-    self->model->nu_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter3") );
-    self->model->nu_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter2") );
-    self->model->thick_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter4") );
-    self->model->nu_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter9") );
-    self->model->func_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter10") );
-    // Read in dispersion parameters
-    PyObject* disp_dict;
-    DispersionVisitor* visitor = new DispersionVisitor();
-    disp_dict = PyDict_GetItemString(self->dispersion, "rad_core0");
-    self->model->rad_core0.dispersion->accept_as_destination(visitor, self->model->rad_core0.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "thick_inter0");
-    self->model->thick_inter0.dispersion->accept_as_destination(visitor, self->model->thick_inter0.dispersion, disp_dict);
-
-	
-	// Get input and determine whether we have to supply a 1D or 2D return value.
-	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
-	    PyErr_SetString(CSphereSLDModelError, 
-	    	"CSphereSLDModel.run expects a q value.");
-		return NULL;
-	}
-	  
-	// Check params
-	if( PyList_Check(pars)==1) {
-		
-		// Length of list should be 2 for I(qx, qy))
-	    npars = PyList_GET_SIZE(pars); 
-	    if(npars!=2) {
-	    	PyErr_SetString(CSphereSLDModelError, 
-	    		"CSphereSLDModel.run expects a double or a list of dimension 2.");
-	    	return NULL;
-	    }
-	    // We have a vector q, get the qx and qy values at which
-	    // to evaluate I(qx,qy)
-	    qx_value = CSphereSLDModel_readDouble(PyList_GET_ITEM(pars,0));
-	    qy_value = CSphereSLDModel_readDouble(PyList_GET_ITEM(pars,1));
-	    return Py_BuildValue("d",(*(self->model))(qx_value,qy_value));
-
-	} else {
-
-		// We have a scalar q, we will evaluate I(q)
-		qx_value = CSphereSLDModel_readDouble(pars);		
-		
-		return Py_BuildValue("d",(*(self->model))(qx_value));
-	}	
-}
-
-static PyObject * reset(CSphereSLDModel *self, PyObject *args) {
-    
-
-    return Py_BuildValue("d",0.0);
-}
-
-static PyObject * set_dispersion(CSphereSLDModel *self, PyObject *args) {
-	PyObject * disp;
-	const char * par_name;
-
-	if ( !PyArg_ParseTuple(args,"sO", &par_name, &disp) ) {
-	    PyErr_SetString(CSphereSLDModelError,
-	    	"CSphereSLDModel.set_dispersion expects a DispersionModel object.");
-		return NULL;
-	}
-	void *temp = PyCObject_AsVoidPtr(disp);
-	DispersionModel * dispersion = static_cast<DispersionModel *>(temp);
-
-
-	// Ugliness necessary to go from python to C
-	    // TODO: refactor this
-    if (!strcmp(par_name, "rad_core0")) {
-        self->model->rad_core0.dispersion = dispersion;
-    } else    if (!strcmp(par_name, "thick_inter0")) {
-        self->model->thick_inter0.dispersion = dispersion;
-    } else {
-	    PyErr_SetString(CSphereSLDModelError,
-	    	"CSphereSLDModel.set_dispersion expects a valid parameter name.");
-		return NULL;
-	}
-
-	DispersionVisitor* visitor = new DispersionVisitor();
-	PyObject * disp_dict = PyDict_New();
-	dispersion->accept_as_source(visitor, dispersion, disp_dict);
-	PyDict_SetItemString(self->dispersion, par_name, disp_dict);
-    return Py_BuildValue("i",1);
-}
-
-
-static PyMethodDef CSphereSLDModel_methods[] = {
-    {"run",      (PyCFunction)run     , METH_VARARGS,
-      "Evaluate the model at a given Q or Q, phi"},
-    {"runXY",      (PyCFunction)runXY     , METH_VARARGS,
-      "Evaluate the model at a given Q or Qx, Qy"},
-    {"calculate_ER",      (PyCFunction)calculate_ER     , METH_VARARGS,
-      "Evaluate the model at a given Q or Q, phi"},
-      
-    {"evalDistribution",  (PyCFunction)evalDistribution , METH_VARARGS,
-      "Evaluate the model at a given Q or Qx, Qy vector "},
-    {"reset",    (PyCFunction)reset   , METH_VARARGS,
-      "Reset pair correlation"},
-    {"set_dispersion",      (PyCFunction)set_dispersion     , METH_VARARGS,
-      "Set the dispersion model for a given parameter"},
-   {NULL}
-};
-
-static PyTypeObject CSphereSLDModelType = {
-    PyObject_HEAD_INIT(NULL)
-    0,                         /*ob_size*/
-    "CSphereSLDModel",             /*tp_name*/
-    sizeof(CSphereSLDModel),             /*tp_basicsize*/
-    0,                         /*tp_itemsize*/
-    (destructor)CSphereSLDModel_dealloc, /*tp_dealloc*/
-    0,                         /*tp_print*/
-    0,                         /*tp_getattr*/
-    0,                         /*tp_setattr*/
-    0,                         /*tp_compare*/
-    0,                         /*tp_repr*/
-    0,                         /*tp_as_number*/
-    0,                         /*tp_as_sequence*/
-    0,                         /*tp_as_mapping*/
-    0,                         /*tp_hash */
-    0,                         /*tp_call*/
-    0,                         /*tp_str*/
-    0,                         /*tp_getattro*/
-    0,                         /*tp_setattro*/
-    0,                         /*tp_as_buffer*/
-    Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
-    "CSphereSLDModel objects",           /* tp_doc */
-    0,		               /* tp_traverse */
-    0,		               /* tp_clear */
-    0,		               /* tp_richcompare */
-    0,		               /* tp_weaklistoffset */
-    0,		               /* tp_iter */
-    0,		               /* tp_iternext */
-    CSphereSLDModel_methods,             /* tp_methods */
-    CSphereSLDModel_members,             /* tp_members */
-    0,                         /* tp_getset */
-    0,                         /* tp_base */
-    0,                         /* tp_dict */
-    0,                         /* tp_descr_get */
-    0,                         /* tp_descr_set */
-    0,                         /* tp_dictoffset */
-    (initproc)CSphereSLDModel_init,      /* tp_init */
-    0,                         /* tp_alloc */
-    CSphereSLDModel_new,                 /* tp_new */
-};
-
-
-//static PyMethodDef module_methods[] = {
-//    {NULL} 
-//};
-
-/**
- * Function used to add the model class to a module
- * @param module: module to add the class to
- */ 
-void addCSphereSLDModel(PyObject *module) {
-	PyObject *d;
-	
-    if (PyType_Ready(&CSphereSLDModelType) < 0)
-        return;
-
-    Py_INCREF(&CSphereSLDModelType);
-    PyModule_AddObject(module, "CSphereSLDModel", (PyObject *)&CSphereSLDModelType);
-    
-    d = PyModule_GetDict(module);
-    static char error_name[] = "CSphereSLDModel.error";
-    CSphereSLDModelError = PyErr_NewException(error_name, NULL, NULL);
-    PyDict_SetItemString(d, "CSphereSLDModelError", CSphereSLDModelError);
-}
-
