Index: sansmodels/src/sans/models/c_extensions/libmultifunc/librefl.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/libmultifunc/librefl.c	(revision 0164899a491d20467a714225f7ea955c852265f3)
+++ sansmodels/src/sans/models/c_extensions/libmultifunc/librefl.c	(revision 0164899a491d20467a714225f7ea955c852265f3)
@@ -0,0 +1,338 @@
+// The original code, of which work was not DANSE funded,
+// was provided by J. Cho.
+// And modified to fit sansmodels/sansview: JC
+
+#include <math.h>
+#include "libmultifunc/librefl.h"
+#include <stdio.h>
+#include <stdlib.h>
+
+
+complex cassign(real, imag)
+	double real, imag;
+{
+	complex x;
+	x.re = real;
+	x.im = imag;
+	return x;
+}
+
+
+complex cadd(x,y)
+	complex x,y;
+{
+	complex z;
+	z.re = x.re + y.re;
+	z.im = x.im + y.im;
+	return z;
+}
+
+complex rcmult(x,y)
+	double x;
+    complex y;
+{
+	complex z;
+	z.re = x*y.re;
+	z.im = x*y.im;
+	return z;
+}
+
+complex csub(x,y)
+	complex x,y;
+{
+	complex z;
+	z.re = x.re - y.re;
+	z.im = x.im - y.im;
+	return z;
+}
+
+
+complex cmult(x,y)
+	complex x,y;
+{
+	complex z;
+	z.re = x.re*y.re - x.im*y.im;
+	z.im = x.re*y.im + x.im*y.re;
+	return z;
+}
+
+complex cdiv(x,y)
+	complex x,y;
+{
+	complex z;
+	z.re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im);
+	z.im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im);
+	return z;
+}
+
+complex cexp(b)
+	complex b;
+{
+	complex z;
+	double br,bi;
+	br=b.re;
+	bi=b.im;
+	z.re = exp(br)*cos(bi);
+	z.im = exp(br)*sin(bi);
+	return z;
+}
+
+
+complex csqrt(z)    //see Schaum`s Math Handbook p. 22, 6.6 and 6.10
+	complex z;
+{
+	complex c;
+	double zr,zi,x,y,r,w;
+
+	zr=z.re;
+	zi=z.im;
+
+	if (zr==0.0 && zi==0.0)
+	{
+    c.re=0.0;
+	c.im=0.0;
+    return c;
+	}
+	else
+	{
+		x=fabs(zr);
+		y=fabs(zi);
+		if (x>y)
+		{
+			r=y/x;
+			w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
+		}
+		else
+		{
+			r=x/y;
+			w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
+		}
+		if (zr >=0.0)
+		{
+			c.re=w;
+			c.im=zi/(2.0*w);
+		}
+		else
+		{
+			c.im=(zi >= 0) ? w : -w;
+			c.re=zi/(2.0*c.im);
+		}
+		return c;
+	}
+}
+
+complex ccos(b)
+	complex b;
+{
+	complex neg,negb,zero,two,z,i,bi,negbi;
+	zero = cassign(0.0,0.0);
+	two = cassign(2.0,0.0);
+	i = cassign(0.0,1.0);
+	bi = cmult(b,i);
+	negbi = csub(zero,bi);
+	z = cdiv(cadd(cexp(bi),cexp(negbi)),two);
+	return z;
+}
+
+// normalized and modified erf
+//   |
+// 1 +                __  - - - -
+//   |             _
+//	 |            _
+//   |        __
+// 0 + - - -
+//   |-------------+------------+--
+//   0           center       n_sub    --->
+//                                     ind
+//
+// n_sub = total no. of bins(or sublayers)
+// ind = x position: 0 to max
+// nu = max x to integration
+double err_mod_func(n_sub, ind, nu)
+	double n_sub,nu, ind;
+{
+	if (nu == 0.0)
+		nu = 1e-14;
+	if (n_sub == 0.0)
+		n_sub = 1.0;
+
+	double center, func;
+	//ind = (n_sub-1.0)/2.0-1.0 +ind;
+	center = n_sub/2.0;
+	// transform it so that min(ind) = 0
+	ind -= center;
+	// normalize by max limit
+	ind /= center;
+	// divide by sqrt(2) to get Gaussian func
+	nu /= sqrt(2.0);
+	ind *= nu;
+	// re-scale and normalize it so that max(erf)=1, min(erf)=0
+	func = erf(ind)/erf(nu)/2.0;
+	// shift it by +0.5 in y-direction so that min(erf) = 0
+	func += 0.5;
+
+	return func;
+}
+double linearfunc(n_sub, ind, nu)
+	double n_sub,nu, ind;
+{
+	if (n_sub == 0.0)
+		n_sub = 1.0;
+
+	double bin_size, func;
+	bin_size = 1.0/n_sub;  //size of each sub-layer
+	// rescale
+	ind *= bin_size;
+	func = ind;
+
+	return func;
+}
+// use the right hand side from the center of power func
+double power_r(n_sub, ind, nu)
+	double n_sub, nu, ind;
+{
+	if (nu == 0.0)
+		nu = 1e-14;
+	if (n_sub == 0.0)
+		n_sub = 1.0;
+
+	double bin_size,func;
+	bin_size = 1.0/n_sub;  //size of each sub-layer
+	// rescale
+	ind *= bin_size;
+	func = pow(ind, nu);
+
+	return func;
+}
+// use the left hand side from the center of power func
+double power_l(n_sub, ind, nu)
+	double n_sub, nu, ind;
+{
+	if (nu == 0.0)
+		nu = 1e-14;
+	if (n_sub == 0.0)
+		n_sub = 1.0;
+
+	double bin_size, func;
+	bin_size = 1.0/n_sub;  //size of each sub-layer
+	// rescale
+	ind *= bin_size;
+	func = 1.0-pow((1.0-ind),nu);
+
+	return func;
+}
+// use 1-exp func from x=0 to x=1
+double exp_r(n_sub, ind, nu)
+	double n_sub, nu, ind;
+{
+	if (nu == 0.0)
+		nu = 1e-14;
+	if (n_sub == 0.0)
+		n_sub = 1.0;
+
+	double bin_size, func;
+	bin_size = 1.0/n_sub;  //size of each sub-layer
+	// rescale
+	ind *= bin_size;
+	// modify func so that func(0) =0 and func(max)=1
+	func = 1.0-exp(-nu*ind);
+	// normalize by its max
+	func /= (1.0-exp(-nu));
+
+	return func;
+}
+
+// use the left hand side mirror image of exp func
+double exp_l(n_sub, ind, nu)
+	double n_sub, nu, ind;
+{
+	if (nu == 0.0)
+		nu = 1e-14;
+	if (n_sub == 0.0)
+		n_sub = 1.0;
+
+	double bin_size, func;
+	bin_size = 1.0/n_sub;  //size of each sub-layer
+	// rescale
+	ind *= bin_size;
+	// modify func
+	func = exp(-nu*(1.0-ind))-exp(-nu);
+	// normalize by its max
+	func /= (1.0-exp(-nu));
+
+	return func;
+}
+
+// To select function called
+// At nu = 0 (singular point), call line function
+double intersldfunc(fun_type, n_sub, i, nu, sld_l, sld_r)
+	double n_sub, nu, sld_l, sld_r,i;
+	int fun_type;
+{
+	double sld_i, func;
+	// this condition protects an error from the singular point
+	if (nu == 0.0){
+		nu = 1e-13;
+	}
+	// select func
+	switch(fun_type){
+		case 1 :
+			func = power_r(n_sub, i, nu);
+			break;
+		case 2 :
+			func = power_l(n_sub, i, nu);
+			break;
+		case 3 :
+			func = exp_r(n_sub, i, nu);
+			break;
+		case 4 :
+			func = exp_l(n_sub, i, nu);
+			break;
+		case 5 :
+			func = linearfunc(n_sub, i, nu);
+			break;
+		default:
+			func = err_mod_func(n_sub, i, nu);
+			break;
+	}
+	// compute sld
+	if (sld_r>sld_l){
+		sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
+	}
+	else if (sld_r<sld_l){
+		func = 1.0-func;
+		sld_i = (sld_l-sld_r)*func+sld_r; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
+	}
+	else{
+		sld_i = sld_r;
+	}
+	return sld_i;
+}
+
+
+// used by refl.c
+double interfunc(fun_type, n_sub, i, sld_l, sld_r)
+	double n_sub, sld_l, sld_r, i;
+	int fun_type;
+{
+	double sld_i, func;
+	switch(fun_type){
+		case 0 :
+			func = err_mod_func(n_sub, i, 2.5);
+			break;
+		default:
+			func = linearfunc(n_sub, i, 1.0);
+			break;
+	}
+	if (sld_r>sld_l){
+		sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
+	}
+	else if (sld_r<sld_l){
+		func = 1.0-func;
+		sld_i = (sld_l-sld_r)*func+sld_r; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
+	}
+	else{
+		sld_i = sld_r;
+	}
+	return sld_i;
+}
Index: sansmodels/src/sans/models/c_extensions/libmultifunc/librefl.h
===================================================================
--- sansmodels/src/sans/models/c_extensions/libmultifunc/librefl.h	(revision 0164899a491d20467a714225f7ea955c852265f3)
+++ sansmodels/src/sans/models/c_extensions/libmultifunc/librefl.h	(revision 0164899a491d20467a714225f7ea955c852265f3)
@@ -0,0 +1,37 @@
+#if !defined(o_h)
+#define librefl_h
+
+typedef struct {
+	double re;
+	double im;
+} complex;
+
+typedef struct {
+	complex a;
+	complex b;
+	complex c;
+	complex d;
+} matrix;
+
+complex cassign(double real, double imag);
+
+complex cadd(complex x,complex y);
+
+complex rcmult(double x,complex y);
+
+complex csub(complex x,complex y);
+
+complex cmult(complex x,complex y);
+
+complex cdiv(complex x,complex y);
+
+complex cexp(complex b);
+
+complex csqrt(complex z);
+
+complex ccos(complex b);
+
+double intersldfunc(int fun_type, double n_sub, double i, double nu, double sld_l, double sld_r);
+double interfunc(int fun_type, double n_sub, double i, double sld_l, double sld_r);
+double linePq(int fun_type, double thick, double sld_in, double sld_out,double r, double q);
+#endif
Index: sansmodels/src/sans/models/c_extensions/onion.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/onion.c	(revision 339ce675656e3ab89eff7b84d8292a2e0ffe6109)
+++ sansmodels/src/sans/models/c_extensions/onion.c	(revision 0164899a491d20467a714225f7ea955c852265f3)
@@ -1,4 +1,4 @@
 /**
- * Scattering model for a sphere
+ * Scattering model for a onion
  */
 
@@ -7,5 +7,5 @@
 #include <stdio.h>
 #include <stdlib.h>
-
+// some details can be found in sld_cal.c
 double so_kernel(double dp[], double q) {
 	int n = dp[0];
@@ -107,7 +107,7 @@
 
                         vol = 4.0 * pi / 3.0 * r * r * r;
-                        if (j == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && A[i]==0.0){
-                        	vol_sub += (vol_pre - vol);
-                        }
+                        //if (j == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && A[i]==0.0){
+                        //	vol_sub += (vol_pre - vol);
+                        //}
                         f += vol * (contr * (fun) + (sld_in[i]-slope[i]) * bes);
                         }
@@ -148,7 +148,7 @@
                                 }
                             vol = 4.0 * pi / 3.0 * r * r * r;
-                            if (j == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && fun_type[i]==0){
-                            	vol_sub += (vol_pre - vol);
-                            }
+                            //if (j == 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[i]);
                             }
@@ -157,5 +157,5 @@
 
         }
-    vol += vol_sub;
+    //vol += vol_sub;
     f2 = f * f / vol * 1.0e8;
 	f2 *= scale;
Index: sansmodels/src/sans/models/c_extensions/refl.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/refl.c	(revision 339ce675656e3ab89eff7b84d8292a2e0ffe6109)
+++ sansmodels/src/sans/models/c_extensions/refl.c	(revision 0164899a491d20467a714225f7ea955c852265f3)
@@ -1,8 +1,10 @@
 /**
- * Scattering model for a sphere
+ *  model for NR
  */
-
+// The original code, of which work was not DANSE funded,
+// was provided by J. Cho.
 #include <math.h>
 #include "refl.h"
+#include "libmultifunc/librefl.h"
 #include <stdio.h>
 #include <stdlib.h>
@@ -10,251 +12,4 @@
 #define lamda 4.62
 
-typedef struct {
-	double re;
-	double im;
-} complex;
-
-typedef struct {
-	complex a;
-	complex b;
-	complex c;
-	complex d;
-} matrix;
-
-complex cassign(real, imag)
-	double real, imag;
-{
-	complex x;
-	x.re = real;
-	x.im = imag;
-	return x;
-}
-
-
-complex cadd(x,y)
-	complex x,y;
-{
-	complex z;
-	z.re = x.re + y.re;
-	z.im = x.im + y.im;
-	return z;
-}
-
-complex rcmult(x,y)
-	double x;
-    complex y;
-{
-	complex z;
-	z.re = x*y.re;
-	z.im = x*y.im;
-	return z;
-}
-
-complex csub(x,y)
-	complex x,y;
-{
-	complex z;
-	z.re = x.re - y.re;
-	z.im = x.im - y.im;
-	return z;
-}
-
-
-complex cmult(x,y)
-	complex x,y;
-{
-	complex z;
-	z.re = x.re*y.re - x.im*y.im;
-	z.im = x.re*y.im + x.im*y.re;
-	return z;
-}
-
-complex cdiv(x,y)
-	complex x,y;
-{
-	complex z;
-	z.re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im);
-	z.im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im);
-	return z;
-}
-
-complex cexp(b)
-	complex b;
-{
-	complex z;
-	double br,bi;
-	br=b.re;
-	bi=b.im;
-	z.re = exp(br)*cos(bi);
-	z.im = exp(br)*sin(bi);
-	return z;
-}
-
-
-complex csqrt(z)   /* see Schaum`s Math Handbook p. 22, 6.6 and 6.10 */
-	complex z;
-{
-	complex c;
-	double zr,zi,x,y,r,w;
-
-	zr=z.re;
-	zi=z.im;
-
-	if (zr==0.0 && zi==0.0)
-	{
-    c.re=0.0;
-	c.im=0.0;
-    return c;
-	}
-	else
-	{
-		x=fabs(zr);
-		y=fabs(zi);
-		if (x>y)
-		{
-			r=y/x;
-			w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
-		}
-		else
-		{
-			r=x/y;
-			w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
-		}
-		if (zr >=0.0)
-		{
-			c.re=w;
-			c.im=zi/(2.0*w);
-		}
-		else
-		{
-			c.im=(zi >= 0) ? w : -w;
-			c.re=zi/(2.0*c.im);
-		}
-		return c;
-	}
-}
-
-complex ccos(b)
-	complex b;
-{
-	complex neg,negb,zero,two,z,i,bi,negbi;
-	zero = cassign(0.0,0.0);
-	two = cassign(2.0,0.0);
-	i = cassign(0.0,1.0);
-	bi = cmult(b,i);
-	negbi = csub(zero,bi);
-	z = cdiv(cadd(cexp(bi),cexp(negbi)),two);
-	return z;
-}
-
-double errfunc(n_sub, i)
-	double n_sub;
-	int i;
-{
-	double bin_size, ind, func;
-	ind = i;
-	// i range = [ -4..4], x range = [ -2.5..2.5]
-	bin_size = n_sub/2.0/2.5;  //size of each sub-layer
-	// rescale erf so that 0 < erf < 1 in -2.5 <= x <= 2.5
-	func = (erf(ind/bin_size)/2.0+0.5);
-	return func;
-}
-
-double linefunc(n_sub, i)
-	double n_sub;
-	int i;
-{
-	double bin_size, ind, func;
-	ind = i + 0.5;
-	// i range = [ -4..4], x range = [ -2.5..2.5]
-	bin_size = 1.0/n_sub;  //size of each sub-layer
-	// rescale erf so that 0 < erf < 1 in -2.5 <= x <= 2.5
-	func = ((ind + floor(n_sub/2.0))*bin_size);
-	return func;
-}
-
-double parabolic_r(n_sub, i)
-	double n_sub;
-	int i;
-{
-	double bin_size, ind, func;
-	ind = i + 0.5;
-	// i range = [ -4..4], x range = [ 0..1]
-	bin_size = 1.0/n_sub;  //size of each sub-layer; n_sub = 0 is a singular point (error)
-	func = ((ind + floor(n_sub/2.0))*bin_size)*((ind + floor(n_sub/2.0))*bin_size);
-	return func;
-}
-
-double parabolic_l(n_sub, i)
-	double n_sub;
-	int i;
-{
-	double bin_size,ind, func;
-	ind = i + 0.5;
-	bin_size = 1.0/n_sub;  //size of each sub-layer; n_sub = 0 is a singular point (error)
-	func =1.0-(((ind + floor(n_sub/2.0))*bin_size) - 1.0) *(((ind + floor(n_sub/2.0))*bin_size) - 1.0);
-	return func;
-}
-
-double cubic_r(n_sub, i)
-	double n_sub;
-	int i;
-{
-	double bin_size,ind, func;
-	ind = i + 0.5;
-	// i range = [ -4..4], x range = [ 0..1]
-	bin_size = 1.0/n_sub;  //size of each sub-layer; n_sub = 0 is a singular point (error)
-	func = ((ind+ floor(n_sub/2.0))*bin_size)*((ind + floor(n_sub/2.0))*bin_size)*((ind + floor(n_sub/2.0))*bin_size);
-	return func;
-}
-
-double cubic_l(n_sub, i)
-	double n_sub;
-	int i;
-{
-	double bin_size,ind, func;
-	ind = i + 0.5;
-	bin_size = 1.0/n_sub;  //size of each sub-layer; n_sub = 0 is a singular point (error)
-	func = 1.0+(((ind + floor(n_sub/2.0)))*bin_size - 1.0)*(((ind + floor(n_sub/2.0)))*bin_size - 1.0)*(((ind + floor(n_sub/2.0)))*bin_size - 1.0);
-	return func;
-}
-
-double interfunc(fun_type, n_sub, i, sld_l, sld_r)
-	double n_sub, sld_l, sld_r;
-	int fun_type, i;
-{
-	double sld_i, func;
-	switch(fun_type){
-		case 1 :
-			func = linefunc(n_sub, i);
-			break;
-		case 2 :
-			func = parabolic_r(n_sub, i);
-			break;
-		case 3 :
-			func = parabolic_l(n_sub, i);
-			break;
-		case 4 :
-			func = cubic_r(n_sub, i);
-			break;
-		case 5 :
-			func = cubic_l(n_sub, i);
-			break;
-		default:
-			func = errfunc(n_sub, i);
-			break;
-	}
-	if (sld_r>sld_l){
-		sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
-	}
-	else if (sld_r<sld_l){
-		func = 1.0-func;
-		sld_i = (sld_l-sld_r)*func+sld_r; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
-	}
-	else{
-		sld_i = sld_r;
-	}
-	return sld_i;
-}
 
 double re_kernel(double dp[], double q) {
@@ -268,5 +23,5 @@
 	double background = dp[6];
 
-	double sld[n+2],sld_im[n+2],thick_inter[n+2],thick[n+2],total_thick;
+	double sld[n+2],thick_inter[n+2],thick[n+2],total_thick;
 	fun_type[0] = dp[3];
 	for (i =1; i<=n; i++){
@@ -275,12 +30,9 @@
 		thick[i] = dp[i+26];
 		fun_type[i] = dp[i+36];
-		sld_im[i] = dp[i+46];
-
 		total_thick += thick[i] + thick_inter[i];
 	}
 	sld[0] = sld_sub;
 	sld[n+1] = sld_super;
-	sld_im[0] = fabs(dp[0+56]);
-	sld_im[n+1] = fabs(dp[1+56]);
+
 	thick[0] = total_thick/5.0;
 	thick[n+1] = total_thick/5.0;
@@ -288,8 +40,8 @@
 	thick_inter[n+1] = 0.0;
 
-	double nsl=21.0; //nsl = Num_sub_layer:  MUST ODD number in double //no other number works now
-	int n_s, floor_nsl;
+	double nsl=21.0; //nsl = Num_sub_layer:
+	int n_s;
     double sld_i,sldim_i,dz,phi,R,ko2;
-    double sign,erfunc;
+    double sign,erfunc, fun;
     double pi;
 	complex  inv_n,phi1,alpha,alpha2,kn,fnm,fnp,rn,Xn,nn,nn2,an,nnp1,one,zero,two,n_sub,n_sup,knp1,Xnp1;
@@ -299,10 +51,8 @@
 	two= cassign(0.0,-2.0);
 
-	floor_nsl = floor(nsl/2.0);
-
 	//Checking if floor is available.
 	//no imaginary sld inputs in this function yet
-    n_sub=cassign(1.0-sld_sub*pow(lamda,2.0)/(2.0*pi),pow(lamda,2.0)/(2.0*pi)*sld_im[0]);
-    n_sup=cassign(1.0-sld_super*pow(lamda,2.0)/(2.0*pi),pow(lamda,2.0)/(2.0*pi)*sld_im[n+1]);
+    n_sub=cassign(1.0-sld_sub*pow(lamda,2.0)/(2.0*pi),0.0);
+    n_sup=cassign(1.0-sld_super*pow(lamda,2.0)/(2.0*pi),0.0);
     ko2 = pow(2.0*pi/lamda,2.0);
 
@@ -318,7 +68,11 @@
     // iteration for # of layers +sub from the top
     for (i=1;i<=n+1; i++){
+    	if (fun_type[i-1]==1)
+    		fun = 5;
+    	else
+			fun = 0;
 		//iteration for 9 sub-layers
 		for (j=0;j<2;j++){
-			for (n_s=-floor_nsl;n_s<=floor_nsl; n_s++){
+			for (n_s=0;n_s<nsl; n_s++){
 				if (j==1){
 					if (i==n+1)
@@ -326,22 +80,15 @@
 					dz = thick[i];
 					sld_i = sld[i];
-					sldim_i = sld_im[i];
 				}
 				else{
-					dz = thick_inter[i-1]/nsl;
+					dz = thick_inter[i-1]/nsl;//nsl;
 					if (sld[i-1] == sld[i]){
 						sld_i = sld[i];
 					}
 					else{
-						sld_i = interfunc(fun_type[i-1],nsl, n_s, sld[i-1], sld[i]);
-					}
-					if (sld_im[i-1] == sld_im[i]){
-						sldim_i = sld_im[i];
-					}
-					else{
-						sldim_i = interfunc(fun_type[i-1],nsl, n_s, sld_im[i-1], sld_im[i]);
+						sld_i = intersldfunc(fun,nsl, n_s+0.5, 2.5, sld[i-1], sld[i]);
 					}
 				}
-				nn = cassign(1.0-sld_i*pow(lamda,2.0)/(2.0*pi),pow(lamda,2.0)/(2.0*pi)*sldim_i);
+				nn = cassign(1.0-sld_i*pow(lamda,2.0)/(2.0*pi),0.0);
 				nn2=cmult(nn,nn);
 
@@ -375,11 +122,11 @@
 }
 /**
- * Function to evaluate 1D scattering function
- * @param pars: parameters of the sphere
+ * Function to evaluate NR function
+ * @param pars: parameters
  * @param q: q-value
  * @return: function value
  */
 double refl_analytical_1D(ReflParameters *pars, double q) {
-	double dp[59];
+	double dp[47];
 
 	dp[0] = pars->n_layers;
@@ -387,5 +134,5 @@
 	dp[2] = pars->thick_inter0;
 	dp[3] = pars->func_inter0;
-	dp[4] = pars->sld_sub0;
+	dp[4] = pars->sld_bottom0;
 	dp[5] = pars->sld_medium;
 	dp[6] = pars->background;
@@ -435,24 +182,10 @@
 	dp[46] = pars->func_inter10;
 
-	dp[47] = pars->sldIM_flat1;
-	dp[48] = pars->sldIM_flat2;
-	dp[49] = pars->sldIM_flat3;
-	dp[50] = pars->sldIM_flat4;
-	dp[51] = pars->sldIM_flat5;
-	dp[52] = pars->sldIM_flat6;
-	dp[53] = pars->sldIM_flat7;
-	dp[54] = pars->sldIM_flat8;
-	dp[55] = pars->sldIM_flat9;
-	dp[56] = pars->sldIM_flat10;
-
-	dp[57] = pars->sldIM_sub0;
-	dp[58] = pars->sldIM_medium;
-
 	return re_kernel(dp, q);
 }
 
 /**
- * Function to evaluate 2D scattering function
- * @param pars: parameters of the sphere
+ * Function to evaluate NR function
+ * @param pars: parameters of NR
  * @param q: q-value
  * @return: function value
Index: sansmodels/src/sans/models/c_extensions/refl.h
===================================================================
--- sansmodels/src/sans/models/c_extensions/refl.h	(revision 339ce675656e3ab89eff7b84d8292a2e0ffe6109)
+++ sansmodels/src/sans/models/c_extensions/refl.h	(revision 0164899a491d20467a714225f7ea955c852265f3)
@@ -1,2 +1,4 @@
+// The original code, of which work was not DANSE funded,
+// was provided by J. Cho.
 #if !defined(o_h)
 #define refl_h
@@ -7,27 +9,21 @@
  //[PYTHONCLASS] = ReflModel
  //[DISP_PARAMS] = thick_inter0
- //[DESCRIPTION] =<text>Form factor of mutishells normalized by the volume. Here each shell is described
- //						 by an exponential function;
- //							I)
- //							For A_shell != 0,
- //							f(r) = B*exp(A_shell*(r-r_in)/thick_shell)+C
- //                         where
- //								B=(sld_out-sld_in)/(exp(A_shell)-1)
- //								C=sld_in-B.
- //							Note that in the above case,
- //								the function becomes a linear function
- //								as A_shell --> 0+ or 0-.
- //							II)
- //							For the exact point of A_shell == 0,
- //							f(r) = sld_in ,i.e., it crosses over flat function
- //							Note that the 'sld_out' becaomes NULL in this case.
- //
- //				background:background,
- //				rad_core: radius of sphere(core)
- //				thick_shell#:the thickness of the shell#
- //				sld_core: the SLD of the sphere
- //				sld_solv: the SLD of the solvent
- //				sld_shell: the SLD of the shell#
- //				A_shell#: the coefficient in the exponential function
+ //[DESCRIPTION] =<text>Calculate neutron reflectivity using the Parratt iterative formula
+ //				Parameters:
+ //				background:background
+ //				scale: scale factor
+ //				sld_bottom0: the SLD of the substrate
+ //				sld_medium: the SLD of the incident medium
+ //					or superstrate
+ //				sld_flatN: the SLD of the flat region of
+ //					the N'th layer
+ //				thick_flatN: the thickness of the flat
+ //					region of the N'th layer
+ //				func_interN: the function used to describe
+ //					the interface of the N'th layer
+ //				thick_interN: the thickness of the interface
+ //					of the N'th layer
+ //				Note: the layer number starts to increase
+ //					from the bottom (substrate) to the top.
  //		</text>
  //[FIXED]=  <text></text>
@@ -48,7 +44,7 @@
 	//  [DEFAULT]=func_inter0= 0
 	double func_inter0;
-	///	sld_sub0 [1/A^(2)]
-	//  [DEFAULT]=sld_sub0= 2.07e-6 [1/A^(2)]
-	double sld_sub0;
+	///	sld_bottom0 [1/A^(2)]
+	//  [DEFAULT]=sld_bottom0= 2.07e-6 [1/A^(2)]
+	double sld_bottom0;
 	///	sld_medium [1/A^(2)]
 	//  [DEFAULT]=sld_medium= 1.0e-6 [1/A^(2)]
@@ -142,28 +138,4 @@
     double func_inter10;
 
-    //  [DEFAULT]=sldIM_flat1=0
-    double sldIM_flat1;
-    //  [DEFAULT]=sldIM_flat2=0
-    double sldIM_flat2;
-    //  [DEFAULT]=sldIM_flat3=0
-    double sldIM_flat3;
-    //  [DEFAULT]=sldIM_flat4=0
-    double sldIM_flat4;
-    //  [DEFAULT]=sldIM_flat5=0
-    double sldIM_flat5;
-    //  [DEFAULT]=sldIM_flat6=0
-    double sldIM_flat6;
-    //  [DEFAULT]=sldIM_flat7=0
-    double sldIM_flat7;
-    //  [DEFAULT]=sldIM_flat8=0
-    double sldIM_flat8;
-    //  [DEFAULT]=sldIM_flat9=0
-    double sldIM_flat9;
-    //  [DEFAULT]=sldIM_flat10=0
-    double sldIM_flat10;
-    //  [DEFAULT]=sldIM_sub0=0
-    double sldIM_sub0;
-    //  [DEFAULT]=sldIM_medium=0
-    double sldIM_medium;
 
 } ReflParameters;
Index: sansmodels/src/sans/models/c_extensions/sld_cal.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/sld_cal.c	(revision 0164899a491d20467a714225f7ea955c852265f3)
+++ sansmodels/src/sans/models/c_extensions/sld_cal.c	(revision 0164899a491d20467a714225f7ea955c852265f3)
@@ -0,0 +1,46 @@
+/**
+ * SLD Calculator
+ * @author: Jae Hie Cho / UTK
+ */
+
+#include "sld_cal.h"
+#include "libmultifunc/librefl.h"
+#include <math.h>
+
+
+/**
+ * Function to calculate sld
+ * @param pars: parameters
+ * @param x: independent param-value
+ * @return: sld value
+ */
+double sld_cal_analytical_1D(SLDCalParameters *pars, double x) {
+	double fun, nsl, n_s, fun_coef, sld_l, sld_r, sld_out;
+
+	fun = pars->fun_type;
+	nsl = pars->npts_inter;
+	n_s = pars->shell_num;
+	fun_coef = pars->nu_inter;
+	sld_l = pars-> sld_left;
+	sld_r = pars-> sld_right;
+
+	int fun_type = floor(fun);
+
+	sld_out = intersldfunc(fun_type, nsl, n_s, fun_coef, sld_l, sld_r);
+
+	return sld_out;
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters
+ * @param q: q-value
+ * @return: function value
+ */
+double sld_cal_analytical_2D(SLDCalParameters *pars, double x, double phi) {
+	return sld_cal_analytical_1D(pars,x);
+}
+
+double sld_cal_analytical_2DXY(SLDCalParameters *pars, double xx, double xy){
+	return sld_cal_analytical_1D(pars,sqrt(xx*xx+xy*xy));
+}
Index: sansmodels/src/sans/models/c_extensions/sld_cal.h
===================================================================
--- sansmodels/src/sans/models/c_extensions/sld_cal.h	(revision 0164899a491d20467a714225f7ea955c852265f3)
+++ sansmodels/src/sans/models/c_extensions/sld_cal.h	(revision 0164899a491d20467a714225f7ea955c852265f3)
@@ -0,0 +1,50 @@
+#if !defined(sld_cal_h)
+#define sld_cal_h
+
+/**
+ * To calculate the sld value
+* [PYTHONCLASS] = SLDCalFunc
+* [DISP_PARAMS] = npts_inter
+  [DESCRIPTION] =<text>To calculate sld values
+ 				</text>
+	[FIXED]= <text>
+			</text>
+
+ **/
+typedef struct {
+    ///	fun_type
+    //  [DEFAULT]=fun_type=0
+    double fun_type;
+
+    /// npts_inter
+    //  [DEFAULT]=npts_inter= 21
+    double npts_inter;
+
+    /// shell_num
+    //  [DEFAULT]=shell_num= 0
+    double shell_num;
+
+    /// nu_inter
+    //  [DEFAULT]=nu_inter= 2.5
+    double nu_inter;
+
+    /// sld_left [1/A^(2)]
+    //  [DEFAULT]=sld_left= 0 [1/A^(2)]
+    double sld_left;
+
+    /// sld_right [1/A^(2)]
+    //  [DEFAULT]=sld_right= 0 [1/A^(2)]
+    double sld_right;
+
+} SLDCalParameters;
+
+
+
+/// 1D function
+double sld_cal_analytical_1D(SLDCalParameters *pars, double q);
+
+/// 2D function
+double sld_cal_analytical_2D(SLDCalParameters *pars, double q, double phi);
+double sld_cal_analytical_2DXY(SLDCalParameters *pars, double qx, double qy);
+
+#endif
Index: sansmodels/src/sans/models/c_extensions/spheresld.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/spheresld.c	(revision 0164899a491d20467a714225f7ea955c852265f3)
+++ sansmodels/src/sans/models/c_extensions/spheresld.c	(revision 0164899a491d20467a714225f7ea955c852265f3)
@@ -0,0 +1,241 @@
+/**
+ * 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,fun_type[n+2];
+
+	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 sld[n+2],thick_inter[n+2],thick[n+2],fun_coef[n+2],total_thick;
+	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;
+
+	double nsl=npts;//21.0; //nsl = Num_sub_layer:  MUST ODD number in double //no other number works now
+	int n_s;
+	int floor_nsl;
+
+    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;
+
+    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;
+    double r0 = 0.0, thick_inter_f;
+
+	//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;
+    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/sans/models/c_extensions/spheresld.h
===================================================================
--- sansmodels/src/sans/models/c_extensions/spheresld.h	(revision 0164899a491d20467a714225f7ea955c852265f3)
+++ sansmodels/src/sans/models/c_extensions/spheresld.h	(revision 0164899a491d20467a714225f7ea955c852265f3)
@@ -0,0 +1,177 @@
+#if !defined(o_h)
+#define sphere_sld_h
+
+/**
+ * Structure definition for sphere parameters
+ */
+ //[PYTHONCLASS] = SphereSLDModel
+ //[DISP_PARAMS] = thick_inter0
+ //[DESCRIPTION] =<text>Calculate neutron reflectivity using the Parratt iterative formula
+ //				Parameters:
+ //				background:background
+ //				scale: scale factor
+ //				sld_core0: the SLD of the substrate
+ //				sld_solv: the SLD of the incident medium
+ //					or superstrate
+ //				sld_flatN: the SLD of the flat region of
+ //					the N'th layer
+ //				thick_flatN: the thickness of the flat
+ //					region of the N'th layer
+ //				func_interN: the function used to describe
+ //					the interface of the N'th layer
+ //				nu_interN: the coefficient for the func_interN
+ //				thick_interN: the thickness of the interface
+ //					of the N'th layer
+ //				Note: the layer number starts to increase
+ //					from the bottom (substrate) to the top.
+ //		</text>
+ //[FIXED]=  <text></text>
+ //[NON_FITTABLE_PARAMS]= <text>n_shells;func_inter0;func_inter1;func_inter2;func_inter3;func_inter4;func_inter5;func_inter5;func_inter7;func_inter8;func_inter9;func_inter10 </text>
+ //[ORIENTATION_PARAMS]= <text> </text>
+
+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;
+
+double sphere_sld_kernel(double dq[], double q);
+
+/// 1D scattering function
+double sphere_sld_analytical_1D(SphereSLDParameters *pars, double q);
+
+/// 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/sans/models/c_models/CReflModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CReflModel.cpp	(revision 339ce675656e3ab89eff7b84d8292a2e0ffe6109)
+++ sansmodels/src/sans/models/c_models/CReflModel.cpp	(revision 0164899a491d20467a714225f7ea955c852265f3)
@@ -90,10 +90,6 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"sldIM_flat4",Py_BuildValue("d",0.000000000000));
         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,"sldIM_flat7",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"sldIM_flat10",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"sldIM_flat5",Py_BuildValue("d",0.000000000000));
         PyDict_SetItemString(self->params,"thick_flat1",Py_BuildValue("d",10.000000000000));
         PyDict_SetItemString(self->params,"thick_flat3",Py_BuildValue("d",100.000000000000));
@@ -103,7 +99,7 @@
         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,"sld_sub0",Py_BuildValue("d",0.000002070000));
+        PyDict_SetItemString(self->params,"thick_inter0",Py_BuildValue("d",1.000000000000));
         PyDict_SetItemString(self->params,"thick_inter1",Py_BuildValue("d",1.000000000000));
-        PyDict_SetItemString(self->params,"thick_inter2",Py_BuildValue("d",1.000000000000));
+        PyDict_SetItemString(self->params,"sld_bottom0",Py_BuildValue("d",0.000002070000));
         PyDict_SetItemString(self->params,"thick_inter3",Py_BuildValue("d",1.000000000000));
         PyDict_SetItemString(self->params,"thick_inter4",Py_BuildValue("d",1.000000000000));
@@ -113,10 +109,7 @@
         PyDict_SetItemString(self->params,"thick_inter8",Py_BuildValue("d",1.000000000000));
         PyDict_SetItemString(self->params,"thick_inter9",Py_BuildValue("d",1.000000000000));
-        PyDict_SetItemString(self->params,"sldIM_flat1",Py_BuildValue("d",0.000000000000));
         PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000000000));
-        PyDict_SetItemString(self->params,"sldIM_flat2",Py_BuildValue("d",0.000000000000));
         PyDict_SetItemString(self->params,"func_inter9",Py_BuildValue("d",0.000000000000));
         PyDict_SetItemString(self->params,"thick_inter10",Py_BuildValue("d",1.000000000000));
-        PyDict_SetItemString(self->params,"sldIM_flat3",Py_BuildValue("d",0.000000000000));
         PyDict_SetItemString(self->params,"func_inter8",Py_BuildValue("d",0.000000000000));
         PyDict_SetItemString(self->params,"thick_flat2",Py_BuildValue("d",100.000000000000));
@@ -136,16 +129,11 @@
         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,"sldIM_flat8",Py_BuildValue("d",0.000000000000));
         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,"sldIM_flat9",Py_BuildValue("d",0.000000000000));
         PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000000000));
         PyDict_SetItemString(self->params,"func_inter1",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"sldIM_sub0",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"sldIM_medium",Py_BuildValue("d",0.000000000000));
         PyDict_SetItemString(self->params,"sld_flat10",Py_BuildValue("d",0.000003500000));
-        PyDict_SetItemString(self->params,"thick_inter0",Py_BuildValue("d",1.000000000000));
-        PyDict_SetItemString(self->params,"sldIM_flat6",Py_BuildValue("d",0.000000000000));
         PyDict_SetItemString(self->params,"func_inter4",Py_BuildValue("d",0.000000000000));
+        PyDict_SetItemString(self->params,"thick_inter2",Py_BuildValue("d",1.000000000000));
         PyDict_SetItemString(self->params,"func_inter10",Py_BuildValue("d",0.000000000000));
         // Initialize dispersion / averaging parameter dict
@@ -284,10 +272,6 @@
 	
 	    // Reader parameter dictionary
-    self->model->sldIM_flat4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat4") );
     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->sldIM_flat7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat7") );
-    self->model->sldIM_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat10") );
-    self->model->sldIM_flat5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat5") );
     self->model->thick_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat1") );
     self->model->thick_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat3") );
@@ -297,7 +281,7 @@
     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->sld_sub0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_sub0") );
+    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->sld_bottom0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_bottom0") );
     self->model->thick_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter3") );
     self->model->thick_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter4") );
@@ -307,10 +291,7 @@
     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->sldIM_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat1") );
     self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
-    self->model->sldIM_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat2") );
     self->model->func_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter9") );
     self->model->thick_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter10") );
-    self->model->sldIM_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat3") );
     self->model->func_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter8") );
     self->model->thick_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat2") );
@@ -330,16 +311,11 @@
     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->sldIM_flat8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat8") );
     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->sldIM_flat9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat9") );
     self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
     self->model->func_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter1") );
-    self->model->sldIM_sub0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_sub0") );
-    self->model->sldIM_medium = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_medium") );
     self->model->sld_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat10") );
-    self->model->thick_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter0") );
-    self->model->sldIM_flat6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat6") );
     self->model->func_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter4") );
+    self->model->thick_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter2") );
     self->model->func_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter10") );
     // Read in dispersion parameters
@@ -410,10 +386,6 @@
 	
 	    // Reader parameter dictionary
-    self->model->sldIM_flat4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat4") );
     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->sldIM_flat7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat7") );
-    self->model->sldIM_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat10") );
-    self->model->sldIM_flat5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat5") );
     self->model->thick_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat1") );
     self->model->thick_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat3") );
@@ -423,7 +395,7 @@
     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->sld_sub0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_sub0") );
+    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->sld_bottom0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_bottom0") );
     self->model->thick_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter3") );
     self->model->thick_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter4") );
@@ -433,10 +405,7 @@
     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->sldIM_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat1") );
     self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
-    self->model->sldIM_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat2") );
     self->model->func_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter9") );
     self->model->thick_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter10") );
-    self->model->sldIM_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat3") );
     self->model->func_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter8") );
     self->model->thick_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat2") );
@@ -456,16 +425,11 @@
     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->sldIM_flat8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat8") );
     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->sldIM_flat9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat9") );
     self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
     self->model->func_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter1") );
-    self->model->sldIM_sub0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_sub0") );
-    self->model->sldIM_medium = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_medium") );
     self->model->sld_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat10") );
-    self->model->thick_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter0") );
-    self->model->sldIM_flat6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat6") );
     self->model->func_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter4") );
+    self->model->thick_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter2") );
     self->model->func_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter10") );
     // Read in dispersion parameters
@@ -523,10 +487,6 @@
 	
 	    // Reader parameter dictionary
-    self->model->sldIM_flat4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat4") );
     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->sldIM_flat7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat7") );
-    self->model->sldIM_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat10") );
-    self->model->sldIM_flat5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat5") );
     self->model->thick_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat1") );
     self->model->thick_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat3") );
@@ -536,7 +496,7 @@
     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->sld_sub0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_sub0") );
+    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->sld_bottom0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_bottom0") );
     self->model->thick_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter3") );
     self->model->thick_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter4") );
@@ -546,10 +506,7 @@
     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->sldIM_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat1") );
     self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
-    self->model->sldIM_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat2") );
     self->model->func_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter9") );
     self->model->thick_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter10") );
-    self->model->sldIM_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat3") );
     self->model->func_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter8") );
     self->model->thick_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat2") );
@@ -569,16 +526,11 @@
     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->sldIM_flat8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat8") );
     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->sldIM_flat9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat9") );
     self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
     self->model->func_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter1") );
-    self->model->sldIM_sub0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_sub0") );
-    self->model->sldIM_medium = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_medium") );
     self->model->sld_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat10") );
-    self->model->thick_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter0") );
-    self->model->sldIM_flat6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat6") );
     self->model->func_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter4") );
+    self->model->thick_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter2") );
     self->model->func_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter10") );
     // Read in dispersion parameters
@@ -605,10 +557,6 @@
 	
 	    // Reader parameter dictionary
-    self->model->sldIM_flat4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat4") );
     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->sldIM_flat7 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat7") );
-    self->model->sldIM_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat10") );
-    self->model->sldIM_flat5 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat5") );
     self->model->thick_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat1") );
     self->model->thick_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat3") );
@@ -618,7 +566,7 @@
     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->sld_sub0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_sub0") );
+    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->sld_bottom0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_bottom0") );
     self->model->thick_inter3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter3") );
     self->model->thick_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter4") );
@@ -628,10 +576,7 @@
     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->sldIM_flat1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat1") );
     self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
-    self->model->sldIM_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat2") );
     self->model->func_inter9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter9") );
     self->model->thick_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter10") );
-    self->model->sldIM_flat3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat3") );
     self->model->func_inter8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter8") );
     self->model->thick_flat2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_flat2") );
@@ -651,16 +596,11 @@
     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->sldIM_flat8 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat8") );
     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->sldIM_flat9 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat9") );
     self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
     self->model->func_inter1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter1") );
-    self->model->sldIM_sub0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_sub0") );
-    self->model->sldIM_medium = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_medium") );
     self->model->sld_flat10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_flat10") );
-    self->model->thick_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter0") );
-    self->model->sldIM_flat6 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldIM_flat6") );
     self->model->func_inter4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter4") );
+    self->model->thick_inter2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_inter2") );
     self->model->func_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter10") );
     // Read in dispersion parameters
Index: sansmodels/src/sans/models/c_models/CSphereSLDModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CSphereSLDModel.cpp	(revision 0164899a491d20467a714225f7ea955c852265f3)
+++ sansmodels/src/sans/models/c_models/CSphereSLDModel.cpp	(revision 0164899a491d20467a714225f7ea955c852265f3)
@@ -0,0 +1,828 @@
+/**
+	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->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 PyMemberDef CSphereSLDModel_members[] = {
+    {"params", T_OBJECT, offsetof(CSphereSLDModel, params), 0,
+     "Parameters"},
+	{"dispersion", T_OBJECT, offsetof(CSphereSLDModel, dispersion), 0,
+	  "Dispersion parameters"},     
+    {"log", T_OBJECT, offsetof(CSphereSLDModel, log), 0,
+     "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;
+	}
+	 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 i,j, 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. */
+        for ( 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, "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, "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) {
+
+	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, "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, "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, "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);
+    CSphereSLDModelError = PyErr_NewException("CSphereSLDModel.error", NULL, NULL);
+    PyDict_SetItemString(d, "CSphereSLDModelError", CSphereSLDModelError);
+}
+
Index: sansmodels/src/sans/models/c_models/refl.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/refl.cpp	(revision 339ce675656e3ab89eff7b84d8292a2e0ffe6109)
+++ sansmodels/src/sans/models/c_models/refl.cpp	(revision 0164899a491d20467a714225f7ea955c852265f3)
@@ -15,5 +15,5 @@
 	thick_inter0 = Parameter(1.0);
 	func_inter0 = Parameter(0);
-	sld_sub0 = Parameter(2.07e-06);
+	sld_bottom0 = Parameter(2.07e-06);
 	sld_medium = Parameter(1.0e-06);
     background = Parameter(0.0);
@@ -67,17 +67,4 @@
     func_inter10 = Parameter(0);
 
-    sldIM_flat1 = Parameter(0);
-    sldIM_flat2 = Parameter(0);
-    sldIM_flat3 = Parameter(0);
-    sldIM_flat4 = Parameter(0);
-    sldIM_flat5 = Parameter(0);
-    sldIM_flat6 = Parameter(0);
-    sldIM_flat7 = Parameter(0);
-    sldIM_flat8 = Parameter(0);
-    sldIM_flat9 = Parameter(0);
-    sldIM_flat10 = Parameter(0);
-
-    sldIM_sub0 = Parameter(0);
-    sldIM_medium = Parameter(0);
 }
 
@@ -88,5 +75,5 @@
  */
 double ReflModel :: operator()(double q) {
-	double dp[59];
+	double dp[47];
 	// Fill parameter array for IGOR library
 	// Add the background after averaging
@@ -95,5 +82,5 @@
 	dp[2] = thick_inter0();
 	dp[3] = func_inter0();
-	dp[4] = sld_sub0();
+	dp[4] = sld_bottom0();
 	dp[5] = sld_medium();
 	dp[6] = background();
@@ -143,18 +130,4 @@
 	dp[46] = func_inter10();
 
-	dp[47] = sldIM_flat1();
-	dp[48] = sldIM_flat2();
-	dp[49] = sldIM_flat3();
-	dp[50] = sldIM_flat4();
-	dp[51] = sldIM_flat5();
-	dp[52] = sldIM_flat6();
-	dp[53] = sldIM_flat7();
-	dp[54] = sldIM_flat8();
-	dp[55] = sldIM_flat9();
-	dp[56] = sldIM_flat10();
-
-	dp[57] = sldIM_sub0();
-	dp[58] = sldIM_medium();
-
 	// Get the dispersion points for the radius
 	//vector<WeightPoint> weights_thick;
Index: sansmodels/src/sans/models/c_models/sld_cal.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/sld_cal.cpp	(revision 0164899a491d20467a714225f7ea955c852265f3)
+++ sansmodels/src/sans/models/c_models/sld_cal.cpp	(revision 0164899a491d20467a714225f7ea955c852265f3)
@@ -0,0 +1,71 @@
+/**
+ * Scattering model classes
+ * The classes use the IGOR library found in
+ *   sansmodels/c_extensions/libmultifunc/librefl.h
+ */
+
+#include <math.h>
+#include "models.hh"
+#include "parameters.hh"
+#include <stdio.h>
+using namespace std;
+
+extern "C" {
+	#include "sld_cal.h"
+}
+
+SLDCalFunc :: SLDCalFunc() {
+	fun_type     = Parameter(0);
+	npts_inter   = Parameter(21);
+	shell_num	 = Parameter(1);
+	nu_inter	 = Parameter(2.5);
+	sld_left 	 = Parameter(1.0e-06);
+	sld_right	 = Parameter(2.0e-06);
+}
+
+/**
+ * Function to evaluate 1D scattering function
+ * @param q: q-value
+ * @return: function value
+ */
+double SLDCalFunc :: operator()(double q) {
+	SLDCalParameters dp;
+
+	dp.fun_type = fun_type();
+	dp.npts_inter = npts_inter();
+	dp.shell_num = shell_num();
+	dp.nu_inter = nu_inter();
+	dp.sld_left = sld_left();
+	dp.sld_right = sld_right();
+
+	return sld_cal_analytical_1D(&dp, q);
+}
+
+/**
+ * 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 SLDCalFunc :: operator()(double qx, double qy) {
+	double q = sqrt(qx*qx + qy*qy);
+		return (*this).operator()(q);
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double SLDCalFunc :: evaluate_rphi(double q, double phi) {
+	return (*this).operator()(q);
+}
+/**
+ * Function to calculate effective radius
+ * @return: effective radius value
+ */
+double SLDCalFunc :: calculate_ER() {
+//NOT implemented yet!!!
+}
Index: sansmodels/src/sans/models/c_models/spheresld.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/spheresld.cpp	(revision 0164899a491d20467a714225f7ea955c852265f3)
+++ sansmodels/src/sans/models/c_models/spheresld.cpp	(revision 0164899a491d20467a714225f7ea955c852265f3)
@@ -0,0 +1,251 @@
+
+#include <math.h>
+#include "models.hh"
+#include "parameters.hh"
+#include <stdio.h>
+using namespace std;
+
+extern "C" {
+	#include "spheresld.h"
+}
+
+SphereSLDModel :: SphereSLDModel() {
+	n_shells = Parameter(1.0);
+	scale = Parameter(1.0);
+	thick_inter0 = Parameter(1.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);
+}
+
+/**
+ * Function to evaluate 1D SphereSLD function
+ * @param q: q-value
+ * @return: function value
+ */
+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] = background();
+
+	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.
+
+	return sphere_sld_kernel(dp,q);
+}
+
+/**
+ * Function to evaluate 2D SphereSLD function
+ * @param q_x: value of Q along x
+ * @param q_y: value of Q along y
+ * @return: function value
+ */
+double SphereSLDModel :: operator()(double qx, double qy) {
+	double q = sqrt(qx*qx + qy*qy);
+	return (*this).operator()(q);
+}
+
+/**
+ * Function to evaluate SphereSLD function
+ * @param pars: parameters of the SphereSLD
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double SphereSLDModel :: evaluate_rphi(double q, double phi) {
+	return (*this).operator()(q);
+}
+
+/**
+ * Function to calculate TOTAL radius
+ * ToDo: Find What is the effective radius for this model.
+ * @return: effective radius value
+ */
+// No polydispersion supported in this model.
+// Calculate max radius assumming max_radius = effective radius
+// Note that this max radius is not affected by sld of layer, sld of interface, or
+// 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 = 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;
+
+	return rad_out;
+
+}
