Index: nsmodels/src/c_extensions/refl_adv.c
===================================================================
--- sansmodels/src/c_extensions/refl_adv.c	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ 	(revision )
@@ -1,268 +1,0 @@
-// The original code, of which work was not DANSE funded,
-// was provided by J. Cho.
-/**
- * NR model Parratt method
- */
-#include <math.h>
-#include "refl_adv.h"
-#include "libmultifunc/librefl.h"
-#include <stdio.h>
-#include <stdlib.h>
-
-#define lamda 4.62
-
-
-double re_adv_kernel(double dp[], double q) {
-	int n = dp[0];
-	int i,j;
-	double nsl;
-
-	double scale = dp[1];
-	double thick_inter_sub = dp[2];
-	double sld_sub = dp[4];
-	double sld_super = dp[5];
-	double background = dp[6];
-	double npts = dp[69]; //number of sub_layers in each interface
-
-	double total_thick=0.0;
-
-  int n_s;
-  double sld_i,sldim_i,dz,phi,R,ko2;
-  double pi;
-
-  int* fun_type;
-  double* sld;
-  double* sld_im;
-  double* thick_inter;
-  double* thick;
-  double* fun_coef;
-  complex  phi1,alpha,alpha2,kn,fnm,fnp,rn,Xn,nn,nn2,an,nnp1,one,two,n_sub,n_sup,knp1,Xnp1;
-
-  fun_type = (int*)malloc((n+2)*sizeof(int));
-  sld = (double*)malloc((n+2)*sizeof(double));
-  sld_im = (double*)malloc((n+2)*sizeof(double));
-  thick_inter = (double*)malloc((n+2)*sizeof(double));
-  thick = (double*)malloc((n+2)*sizeof(double));
-  fun_coef = (double*)malloc((n+2)*sizeof(double));
-
-  fun_type[0] = dp[3];
-	fun_coef[0] = fabs(dp[70]);
-	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];
-		sld_im[i] = dp[i+46];
-		fun_coef[i] = fabs(dp[i+56]);
-		//printf("type_func2 =%g\n",fun_coef[i]);
-
-		total_thick += thick[i] + thick_inter[i];
-	}
-	sld[0] = sld_sub;
-	sld[n+1] = sld_super;
-	sld_im[0] = fabs(dp[1+66]);
-	sld_im[n+1] = fabs(dp[2+66]);
-	thick[0] = total_thick/5.0;
-	thick[n+1] = total_thick/5.0;
-	thick_inter[0] = thick_inter_sub;
-	thick_inter[n+1] = 0.0;
-	fun_coef[n+1] = 0.0;
-
-	nsl=npts;//21.0; //nsl = Num_sub_layer:  MUST ODD number in double //no other number works now
-
-	pi = 4.0*atan(1.0);
-    one = cassign(1.0,0.0);
-	Xn = cassign(0.0,0.0);
-	two = cassign(0.0,-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]);
-    ko2 = pow(2.0*pi/lamda,2.0);
-
-    phi = asin(lamda*q/(4.0*pi));
-    phi1 = cplx_div(rcmult(phi,one),n_sup);
-    alpha = cplx_mult(n_sup,cplx_cos(phi1));
-	alpha2 = cplx_mult(alpha,alpha);
-
-    nnp1=n_sub;
-    knp1=cplx_sqrt(rcmult(ko2,cplx_sub(cplx_mult(nnp1,nnp1),alpha2)));  //nnp1*ko*sin(phinp1)
-    Xnp1=cassign(0.0,0.0);
-    dz = 0.0;
-    // iteration for # of layers +sub from the top
-    for (i=1;i<=n+1; i++){
-    	//if (fun_coef[i]==0.0)
-    	//	// this condition protects an error in numerical multiplication
-    	//	fun_coef[i] = 1e-14;
-		//iteration for 9 sub-layers
-		for (j=0;j<2;j++){
-			for (n_s=0;n_s<nsl; n_s++){
-				// for flat layer
-				if (j==1){
-					if (i==n+1)
-						break;
-					dz = thick[i];
-					sld_i = sld[i];
-					sldim_i = sld_im[i];
-				}
-				// for interface
-				else{
-					dz = thick_inter[i-1]/nsl;
-					if (sld[i-1] == sld[i]){
-						sld_i = sld[i];
-					}
-					else{
-						sld_i = intersldfunc(fun_type[i-1],nsl, n_s+0.5, fun_coef[i-1], sld[i-1], sld[i]);
-					}
-					if (sld_im[i-1] == sld_im[i]){
-						sldim_i = sld_im[i];
-					}
-					else{
-						sldim_i = intersldfunc(fun_type[i-1],nsl, n_s+0.5, fun_coef[i-1], sld_im[i-1], sld_im[i]);
-					}
-				}
-				nn = cassign(1.0-sld_i*pow(lamda,2.0)/(2.0*pi),pow(lamda,2.0)/(2.0*pi)*sldim_i);
-				nn2=cplx_mult(nn,nn);
-
-				kn=cplx_sqrt(rcmult(ko2,cplx_sub(nn2,alpha2)));        //nn*ko*sin(phin)
-				an=cplx_exp(rcmult(dz,cplx_mult(two,kn)));
-
-				fnm=cplx_sub(kn,knp1);
-				fnp=cplx_add(kn,knp1);
-				rn=cplx_div(fnm,fnp);
-				Xn=cplx_mult(an,cplx_div(cplx_add(rn,Xnp1),cplx_add(one,cplx_mult(rn,Xnp1))));    //Xn=an*((rn+Xnp1*anp1)/(1+rn*Xnp1*anp1))
-
-				Xnp1=Xn;
-				knp1=kn;
-				// no for-loop for flat layer
-				if (j==1)
-					break;
-			}
-		}
-    }
-    R=pow(Xn.re,2.0)+pow(Xn.im,2.0);
-    // This temperarily fixes the total reflection for Rfunction and linear.
-    // ToDo: Show why it happens that Xn.re=0 and Xn.im >1!
-    if (Xn.im == 0.0 || R > 1){
-    	R=1.0;
-    }
-    R *= scale;
-    R += background;
-
-    free(fun_type);
-    free(sld);
-    free(sld_im);
-    free(thick_inter);
-    free(thick);
-    free(fun_coef);
-
-    return R;
-
-}
-
-/**
- * Function to evaluate NR function
- * @param pars: parameters of refl
- * @param q: q-value
- * @return: function value
- */
-
-double refl_adv_analytical_1D(ReflAdvParameters *pars, double q) {
-	double dp[71];
-
-	dp[0] = pars->n_layers;
-	dp[1] = pars->scale;
-	dp[2] = pars->thick_inter0;
-	dp[3] = pars->func_inter0;
-	dp[4] = pars->sld_bottom0;
-	dp[5] = pars->sld_medium;
-	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->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->nu_inter1;
-	dp[58] = pars->nu_inter2;
-	dp[59] = pars->nu_inter3;
-	dp[60] = pars->nu_inter4;
-	dp[61] = pars->nu_inter5;
-	dp[62] = pars->nu_inter6;
-	dp[63] = pars->nu_inter7;
-	dp[64] = pars->nu_inter8;
-	dp[65] = pars->nu_inter9;
-	dp[66] = pars->nu_inter10;
-
-	dp[67] = pars->sldIM_sub0;
-	dp[68] = pars->sldIM_medium;
-	dp[69] = pars->npts_inter;
-	dp[70] = pars->nu_inter0;
-
-	return re_adv_kernel(dp, q);
-}
-
-/**
- * Function to evaluate NR function
- * @param pars: parameters of NR
- * @param q: q-value
- * @return: function value
- */
-double refl_adv_analytical_2D(ReflAdvParameters *pars, double q, double phi) {
-	return refl_adv_analytical_1D(pars,q);
-}
-
-double refl_adv_analytical_2DXY(ReflAdvParameters *pars, double qx, double qy){
-	return refl_adv_analytical_1D(pars,sqrt(qx*qx+qy*qy));
-}
Index: sansmodels/src/c_extensions/refl_adv.h
===================================================================
--- sansmodels/src/c_extensions/refl_adv.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/refl_adv.h	(revision dd60b455f137744c64501646f2107cdb72dab9c1)
@@ -3,201 +3,204 @@
 #if !defined(o_h)
 #define refl_adv_h
+#include "parameters.hh"
 
 /**
  * Structure definition for sphere parameters
  */
- //[PYTHONCLASS] = ReflAdvModel
- //[DISP_PARAMS] = thick_inter0
- //[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
- //				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_layers;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 layers
-	//  [DEFAULT]=n_layers=1
-	int n_layers;
-    /// 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_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)]
-	double sld_medium;
-	/// 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 [A]
-    double thick_inter1;
-    //  [DEFAULT]=thick_inter2=50 [A]
-    double thick_inter2;
-    //  [DEFAULT]=thick_inter3=50 [A]
-    double thick_inter3;
-    //  [DEFAULT]=thick_inter4=50 [A]
-    double thick_inter4;
-    //  [DEFAULT]=thick_inter5=50 [A]
-    double thick_inter5;
-    //  [DEFAULT]=thick_inter6=50 [A]
-    double thick_inter6;
-    //  [DEFAULT]=thick_inter7=50 [A]
-    double thick_inter7;
-    //  [DEFAULT]=thick_inter8=50 [A]
-    double thick_inter8;
-    //  [DEFAULT]=thick_inter9=50 [A]
-    double thick_inter9;
-    //  [DEFAULT]=thick_inter10=50 [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]=sldIM_flat1=0 [1/A^(2)]
-    double sldIM_flat1;
-    //  [DEFAULT]=sldIM_flat2=0 [1/A^(2)]
-    double sldIM_flat2;
-    //  [DEFAULT]=sldIM_flat3=0 [1/A^(2)]
-    double sldIM_flat3;
-    //  [DEFAULT]=sldIM_flat4=0 [1/A^(2)]
-    double sldIM_flat4;
-    //  [DEFAULT]=sldIM_flat5=0 [1/A^(2)]
-    double sldIM_flat5;
-    //  [DEFAULT]=sldIM_flat6=0 [1/A^(2)]
-    double sldIM_flat6;
-    //  [DEFAULT]=sldIM_flat7=0 [1/A^(2)]
-    double sldIM_flat7;
-    //  [DEFAULT]=sldIM_flat8=0 [1/A^(2)]
-    double sldIM_flat8;
-    //  [DEFAULT]=sldIM_flat9=0 [1/A^(2)]
-    double sldIM_flat9;
-    //  [DEFAULT]=sldIM_flat10=0 [1/A^(2)]
-    double sldIM_flat10;
-
-    //  [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]=sldIM_sub0=0
-    double sldIM_sub0;
-    //  [DEFAULT]=sldIM_medium=0
-    double sldIM_medium;
-    //  [DEFAULT]=npts_inter=21.0
-    double npts_inter;
-    //  [DEFAULT]=nu_inter0=2.5
-    double nu_inter0;
-} ReflAdvParameters;
-
-double re_adv_kernel(double dq[], double q);
-
-/// 1D scattering function
-double refl_adv_analytical_1D(ReflAdvParameters *pars, double q);
-
-/// 2D scattering function
-double refl_adv_analytical_2D(ReflAdvParameters *pars, double q, double phi);
-double refl_adv_analytical_2DXY(ReflAdvParameters *pars, double qx, double qy);
+//[PYTHONCLASS] = ReflAdvModel
+//[DISP_PARAMS] = thick_inter0
+//[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
+//				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_layers;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>
+
+class ReflAdvModel{
+public:
+  // Model parameters
+  /// number of layers
+  //  [DEFAULT]=n_layers=1
+  Parameter n_layers;
+  /// Scale factor
+  //  [DEFAULT]=scale= 1.0
+  Parameter scale;
+  /// thick_inter0 [A]
+  //  [DEFAULT]=thick_inter0=50.0 [A]
+  Parameter thick_inter0;
+  /// func_inter0
+  //  [DEFAULT]=func_inter0= 0
+  Parameter func_inter0;
+  /// sld_bottom0 [1/A^(2)]
+  //  [DEFAULT]=sld_bottom0= 2.07e-6 [1/A^(2)]
+  Parameter sld_bottom0;
+  /// sld_medium [1/A^(2)]
+  //  [DEFAULT]=sld_medium= 1.0e-6 [1/A^(2)]
+  Parameter sld_medium;
+  /// Background
+  //  [DEFAULT]=background=0
+  Parameter background;
+
+  //  [DEFAULT]=sld_flat1=4.0e-06 [1/A^(2)]
+  Parameter sld_flat1;
+  //  [DEFAULT]=sld_flat2=3.5e-06 [1/A^(2)]
+  Parameter sld_flat2;
+  //  [DEFAULT]=sld_flat3=4.0e-06 [1/A^(2)]
+  Parameter sld_flat3;
+  //  [DEFAULT]=sld_flat4=3.5e-06 [1/A^(2)]
+  Parameter sld_flat4;
+  //  [DEFAULT]=sld_flat5=4.0e-06 [1/A^(2)]
+  Parameter sld_flat5;
+  //  [DEFAULT]=sld_flat6=3.5e-06 [1/A^(2)]
+  Parameter sld_flat6;
+  //  [DEFAULT]=sld_flat7=4.0e-06 [1/A^(2)]
+  Parameter sld_flat7;
+  //  [DEFAULT]=sld_flat8=3.5e-06 [1/A^(2)]
+  Parameter sld_flat8;
+  //  [DEFAULT]=sld_flat9=4.0e-06 [1/A^(2)]
+  Parameter sld_flat9;
+  //  [DEFAULT]=sld_flat10=3.5e-06 [1/A^(2)]
+  Parameter sld_flat10;
+
+  //  [DEFAULT]=thick_inter1=50 [A]
+  Parameter thick_inter1;
+  //  [DEFAULT]=thick_inter2=50 [A]
+  Parameter thick_inter2;
+  //  [DEFAULT]=thick_inter3=50 [A]
+  Parameter thick_inter3;
+  //  [DEFAULT]=thick_inter4=50 [A]
+  Parameter thick_inter4;
+  //  [DEFAULT]=thick_inter5=50 [A]
+  Parameter thick_inter5;
+  //  [DEFAULT]=thick_inter6=50 [A]
+  Parameter thick_inter6;
+  //  [DEFAULT]=thick_inter7=50 [A]
+  Parameter thick_inter7;
+  //  [DEFAULT]=thick_inter8=50 [A]
+  Parameter thick_inter8;
+  //  [DEFAULT]=thick_inter9=50 [A]
+  Parameter thick_inter9;
+  //  [DEFAULT]=thick_inter10=50 [A]
+  Parameter thick_inter10;
+
+  //  [DEFAULT]=thick_flat1=100 [A]
+  Parameter thick_flat1;
+  //  [DEFAULT]=thick_flat2=100 [A]
+  Parameter thick_flat2;
+  //  [DEFAULT]=thick_flat3=100 [A]
+  Parameter thick_flat3;
+  //  [DEFAULT]=thick_flat4=100 [A]
+  Parameter thick_flat4;
+  //  [DEFAULT]=thick_flat5=100 [A]
+  Parameter thick_flat5;
+  //  [DEFAULT]=thick_flat6=100 [A]
+  Parameter thick_flat6;
+  //  [DEFAULT]=thick_flat7=100 [A]
+  Parameter thick_flat7;
+  //  [DEFAULT]=thick_flat8=100 [A]
+  Parameter thick_flat8;
+  //  [DEFAULT]=thick_flat9=100 [A]
+  Parameter thick_flat9;
+  //  [DEFAULT]=thick_flat10=100 [A]
+  Parameter thick_flat10;
+
+  //  [DEFAULT]=func_inter1=0
+  Parameter func_inter1;
+  //  [DEFAULT]=func_inter2=0
+  Parameter func_inter2;
+  //  [DEFAULT]=func_inter3=0
+  Parameter func_inter3;
+  //  [DEFAULT]=func_inter4=0
+  Parameter func_inter4;
+  //  [DEFAULT]=func_inter5=0
+  Parameter func_inter5;
+  //  [DEFAULT]=func_inter6=0
+  Parameter func_inter6;
+  //  [DEFAULT]=func_inter7=0
+  Parameter func_inter7;
+  //  [DEFAULT]=func_inter8=0
+  Parameter func_inter8;
+  //  [DEFAULT]=func_inter9=0
+  Parameter func_inter9;
+  //  [DEFAULT]=func_inter10=0
+  Parameter func_inter10;
+
+  //  [DEFAULT]=sldIM_flat1=0 [1/A^(2)]
+  Parameter sldIM_flat1;
+  //  [DEFAULT]=sldIM_flat2=0 [1/A^(2)]
+  Parameter sldIM_flat2;
+  //  [DEFAULT]=sldIM_flat3=0 [1/A^(2)]
+  Parameter sldIM_flat3;
+  //  [DEFAULT]=sldIM_flat4=0 [1/A^(2)]
+  Parameter sldIM_flat4;
+  //  [DEFAULT]=sldIM_flat5=0 [1/A^(2)]
+  Parameter sldIM_flat5;
+  //  [DEFAULT]=sldIM_flat6=0 [1/A^(2)]
+  Parameter sldIM_flat6;
+  //  [DEFAULT]=sldIM_flat7=0 [1/A^(2)]
+  Parameter sldIM_flat7;
+  //  [DEFAULT]=sldIM_flat8=0 [1/A^(2)]
+  Parameter sldIM_flat8;
+  //  [DEFAULT]=sldIM_flat9=0 [1/A^(2)]
+  Parameter sldIM_flat9;
+  //  [DEFAULT]=sldIM_flat10=0 [1/A^(2)]
+  Parameter sldIM_flat10;
+
+  //  [DEFAULT]=nu_inter1=2.5
+  Parameter nu_inter1;
+  //  [DEFAULT]=nu_inter2=2.5
+  Parameter nu_inter2;
+  //  [DEFAULT]=nu_inter3=2.5
+  Parameter nu_inter3;
+  //  [DEFAULT]=nu_inter4=2.5
+  Parameter nu_inter4;
+  //  [DEFAULT]=nu_inter5=2.5
+  Parameter nu_inter5;
+  //  [DEFAULT]=nu_inter6=2.5
+  Parameter nu_inter6;
+  //  [DEFAULT]=nu_inter7=2.5
+  Parameter nu_inter7;
+  //  [DEFAULT]=nu_inter8=2.5
+  Parameter nu_inter8;
+  //  [DEFAULT]=nu_inter9=2.5
+  Parameter nu_inter9;
+  //  [DEFAULT]=nu_inter10=2.5
+  Parameter nu_inter10;
+
+
+  //  [DEFAULT]=sldIM_sub0=0
+  Parameter sldIM_sub0;
+  //  [DEFAULT]=sldIM_medium=0
+  Parameter sldIM_medium;
+  //  [DEFAULT]=npts_inter=21.0
+  Parameter npts_inter;
+  //  [DEFAULT]=nu_inter0=2.5
+  Parameter nu_inter0;
+
+  // Constructor
+  ReflAdvModel();
+
+  // Operators to get I(Q)
+  double operator()(double q);
+  double operator()(double qx, double qy);
+  double calculate_ER();
+  double evaluate_rphi(double q, double phi);
+};
 
 #endif
Index: sansmodels/src/c_models/models.hh
===================================================================
--- sansmodels/src/c_models/models.hh	(revision 2d1b7008c4fcc2791722e1581f47a3559d64f820)
+++ sansmodels/src/c_models/models.hh	(revision dd60b455f137744c64501646f2107cdb72dab9c1)
@@ -27,95 +27,5 @@
 using namespace std;
 
-class ReflAdvModel{
-public:
-	// Model parameters
-	Parameter n_layers;
-	Parameter scale;
-	Parameter thick_inter0;
-	Parameter func_inter0;
-	Parameter sld_bottom0;
-	Parameter sld_medium;
-	Parameter background;
-
-	Parameter sld_flat1;
-	Parameter sld_flat2;
-	Parameter sld_flat3;
-	Parameter sld_flat4;
-	Parameter sld_flat5;
-	Parameter sld_flat6;
-	Parameter sld_flat7;
-	Parameter sld_flat8;
-	Parameter sld_flat9;
-	Parameter sld_flat10;
-
-	Parameter thick_inter1;
-	Parameter thick_inter2;
-	Parameter thick_inter3;
-	Parameter thick_inter4;
-	Parameter thick_inter5;
-	Parameter thick_inter6;
-	Parameter thick_inter7;
-	Parameter thick_inter8;
-	Parameter thick_inter9;
-	Parameter thick_inter10;
-
-	Parameter thick_flat1;
-	Parameter thick_flat2;
-	Parameter thick_flat3;
-	Parameter thick_flat4;
-	Parameter thick_flat5;
-	Parameter thick_flat6;
-	Parameter thick_flat7;
-	Parameter thick_flat8;
-	Parameter thick_flat9;
-	Parameter thick_flat10;
-
-	Parameter func_inter1;
-	Parameter func_inter2;
-	Parameter func_inter3;
-	Parameter func_inter4;
-	Parameter func_inter5;
-	Parameter func_inter6;
-	Parameter func_inter7;
-	Parameter func_inter8;
-	Parameter func_inter9;
-	Parameter func_inter10;
-
-	Parameter sldIM_flat1;
-	Parameter sldIM_flat2;
-	Parameter sldIM_flat3;
-	Parameter sldIM_flat4;
-	Parameter sldIM_flat5;
-	Parameter sldIM_flat6;
-	Parameter sldIM_flat7;
-	Parameter sldIM_flat8;
-	Parameter sldIM_flat9;
-	Parameter sldIM_flat10;
-
-	Parameter nu_inter1;
-	Parameter nu_inter2;
-	Parameter nu_inter3;
-	Parameter nu_inter4;
-	Parameter nu_inter5;
-	Parameter nu_inter6;
-	Parameter nu_inter7;
-	Parameter nu_inter8;
-	Parameter nu_inter9;
-	Parameter nu_inter10;
-
-	Parameter sldIM_sub0;
-	Parameter sldIM_medium;
-	Parameter npts_inter;
-	Parameter nu_inter0;
-
-	// Constructor
-	ReflAdvModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
+
 
 
Index: sansmodels/src/c_models/refl_adv.cpp
===================================================================
--- sansmodels/src/c_models/refl_adv.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/refl_adv.cpp	(revision dd60b455f137744c64501646f2107cdb72dab9c1)
@@ -1,96 +1,249 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
+#include <stdlib.h>
+#include "refl_adv.h"
+
+extern "C" {
+#include "libmultifunc/librefl.h"
+}
+
 using namespace std;
 
-extern "C" {
-	#include "refl_adv.h"
-}
-
+#define lamda 4.62
+
+
+double re_adv_kernel(double dp[], double q) {
+  int n = dp[0];
+  int i,j;
+  double nsl;
+
+  double scale = dp[1];
+  double thick_inter_sub = dp[2];
+  double sld_sub = dp[4];
+  double sld_super = dp[5];
+  double background = dp[6];
+  double npts = dp[69]; //number of sub_layers in each interface
+
+  double total_thick=0.0;
+
+  int n_s;
+  double sld_i,sldim_i,dz,phi,R,ko2;
+  double pi;
+
+  int* fun_type;
+  double* sld;
+  double* sld_im;
+  double* thick_inter;
+  double* thick;
+  double* fun_coef;
+  complex  phi1,alpha,alpha2,kn,fnm,fnp,rn,Xn,nn,nn2,an,nnp1,one,two,n_sub,n_sup,knp1,Xnp1;
+
+  fun_type = (int*)malloc((n+2)*sizeof(int));
+  sld = (double*)malloc((n+2)*sizeof(double));
+  sld_im = (double*)malloc((n+2)*sizeof(double));
+  thick_inter = (double*)malloc((n+2)*sizeof(double));
+  thick = (double*)malloc((n+2)*sizeof(double));
+  fun_coef = (double*)malloc((n+2)*sizeof(double));
+
+  fun_type[0] = dp[3];
+  fun_coef[0] = fabs(dp[70]);
+  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];
+    sld_im[i] = dp[i+46];
+    fun_coef[i] = fabs(dp[i+56]);
+    //printf("type_func2 =%g\n",fun_coef[i]);
+
+    total_thick += thick[i] + thick_inter[i];
+  }
+  sld[0] = sld_sub;
+  sld[n+1] = sld_super;
+  sld_im[0] = fabs(dp[1+66]);
+  sld_im[n+1] = fabs(dp[2+66]);
+  thick[0] = total_thick/5.0;
+  thick[n+1] = total_thick/5.0;
+  thick_inter[0] = thick_inter_sub;
+  thick_inter[n+1] = 0.0;
+  fun_coef[n+1] = 0.0;
+
+  nsl=npts;//21.0; //nsl = Num_sub_layer:  MUST ODD number in double //no other number works now
+
+  pi = 4.0*atan(1.0);
+  one = cassign(1.0,0.0);
+  Xn = cassign(0.0,0.0);
+  two = cassign(0.0,-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]);
+  ko2 = pow(2.0*pi/lamda,2.0);
+
+  phi = asin(lamda*q/(4.0*pi));
+  phi1 = cplx_div(rcmult(phi,one),n_sup);
+  alpha = cplx_mult(n_sup,cplx_cos(phi1));
+  alpha2 = cplx_mult(alpha,alpha);
+
+  nnp1=n_sub;
+  knp1=cplx_sqrt(rcmult(ko2,cplx_sub(cplx_mult(nnp1,nnp1),alpha2)));  //nnp1*ko*sin(phinp1)
+  Xnp1=cassign(0.0,0.0);
+  dz = 0.0;
+  // iteration for # of layers +sub from the top
+  for (i=1;i<=n+1; i++){
+    //if (fun_coef[i]==0.0)
+    //  // this condition protects an error in numerical multiplication
+    //  fun_coef[i] = 1e-14;
+    //iteration for 9 sub-layers
+    for (j=0;j<2;j++){
+      for (n_s=0;n_s<nsl; n_s++){
+        // for flat layer
+        if (j==1){
+          if (i==n+1)
+            break;
+          dz = thick[i];
+          sld_i = sld[i];
+          sldim_i = sld_im[i];
+        }
+        // for interface
+        else{
+          dz = thick_inter[i-1]/nsl;
+          if (sld[i-1] == sld[i]){
+            sld_i = sld[i];
+          }
+          else{
+            sld_i = intersldfunc(fun_type[i-1],nsl, n_s+0.5, fun_coef[i-1], sld[i-1], sld[i]);
+          }
+          if (sld_im[i-1] == sld_im[i]){
+            sldim_i = sld_im[i];
+          }
+          else{
+            sldim_i = intersldfunc(fun_type[i-1],nsl, n_s+0.5, fun_coef[i-1], sld_im[i-1], sld_im[i]);
+          }
+        }
+        nn = cassign(1.0-sld_i*pow(lamda,2.0)/(2.0*pi),pow(lamda,2.0)/(2.0*pi)*sldim_i);
+        nn2=cplx_mult(nn,nn);
+
+        kn=cplx_sqrt(rcmult(ko2,cplx_sub(nn2,alpha2)));        //nn*ko*sin(phin)
+        an=cplx_exp(rcmult(dz,cplx_mult(two,kn)));
+
+        fnm=cplx_sub(kn,knp1);
+        fnp=cplx_add(kn,knp1);
+        rn=cplx_div(fnm,fnp);
+        Xn=cplx_mult(an,cplx_div(cplx_add(rn,Xnp1),cplx_add(one,cplx_mult(rn,Xnp1))));    //Xn=an*((rn+Xnp1*anp1)/(1+rn*Xnp1*anp1))
+
+        Xnp1=Xn;
+        knp1=kn;
+        // no for-loop for flat layer
+        if (j==1)
+          break;
+      }
+    }
+  }
+  R=pow(Xn.re,2.0)+pow(Xn.im,2.0);
+  // This temperarily fixes the total reflection for Rfunction and linear.
+  // ToDo: Show why it happens that Xn.re=0 and Xn.im >1!
+  if (Xn.im == 0.0 || R > 1){
+    R=1.0;
+  }
+  R *= scale;
+  R += background;
+
+  free(fun_type);
+  free(sld);
+  free(sld_im);
+  free(thick_inter);
+  free(thick);
+  free(fun_coef);
+
+  return R;
+
+}
 ReflAdvModel :: ReflAdvModel() {
-	n_layers = Parameter(1.0);
-	scale = Parameter(1.0);
-	thick_inter0 = Parameter(1.0);
-	func_inter0 = Parameter(0);
-	sld_bottom0 = Parameter(2.07e-06);
-	sld_medium = 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(15.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);
-
-    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);
-
-    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);
-
-    sldIM_sub0 = Parameter(0.0);
-    sldIM_medium = Parameter(0.0);
-    npts_inter = Parameter(21.0);
-    nu_inter0 = Parameter(2.5);
+  n_layers = Parameter(1.0);
+  scale = Parameter(1.0);
+  thick_inter0 = Parameter(1.0);
+  func_inter0 = Parameter(0);
+  sld_bottom0 = Parameter(2.07e-06);
+  sld_medium = 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(15.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);
+
+  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);
+
+  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);
+
+  sldIM_sub0 = Parameter(0.0);
+  sldIM_medium = Parameter(0.0);
+  npts_inter = Parameter(21.0);
+  nu_inter0 = Parameter(2.5);
 }
 
@@ -101,91 +254,91 @@
  */
 double ReflAdvModel :: operator()(double q) {
-	double dp[71];
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = n_layers();
-	dp[1] = scale();
-	dp[2] = thick_inter0();
-	dp[3] = func_inter0();
-	dp[4] = sld_bottom0();
-	dp[5] = sld_medium();
-	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] = 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] = nu_inter1();
-	dp[58] = nu_inter2();
-	dp[59] = nu_inter3();
-	dp[60] = nu_inter4();
-	dp[61] = nu_inter5();
-	dp[62] = nu_inter6();
-	dp[63] = nu_inter7();
-	dp[64] = nu_inter8();
-	dp[65] = nu_inter9();
-	dp[66] = nu_inter10();
-
-	dp[67] = sldIM_sub0();
-	dp[68] = sldIM_medium();
-	dp[69] = npts_inter();
-	dp[70] = nu_inter0();
-	// Get the dispersion points for the radius
-	//vector<WeightPoint> weights_thick;
-	//thick_inter0.get_weights(weights_thick_inter0);
-
-
-	return re_adv_kernel(dp,q);
+  double dp[71];
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = n_layers();
+  dp[1] = scale();
+  dp[2] = thick_inter0();
+  dp[3] = func_inter0();
+  dp[4] = sld_bottom0();
+  dp[5] = sld_medium();
+  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] = 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] = nu_inter1();
+  dp[58] = nu_inter2();
+  dp[59] = nu_inter3();
+  dp[60] = nu_inter4();
+  dp[61] = nu_inter5();
+  dp[62] = nu_inter6();
+  dp[63] = nu_inter7();
+  dp[64] = nu_inter8();
+  dp[65] = nu_inter9();
+  dp[66] = nu_inter10();
+
+  dp[67] = sldIM_sub0();
+  dp[68] = sldIM_medium();
+  dp[69] = npts_inter();
+  dp[70] = nu_inter0();
+  // Get the dispersion points for the radius
+  //vector<WeightPoint> weights_thick;
+  //thick_inter0.get_weights(weights_thick_inter0);
+
+
+  return re_adv_kernel(dp,q);
 }
 
@@ -197,10 +350,10 @@
  */
 double ReflAdvModel :: operator()(double qx, double qy) {
-	// For 2D set qy as q, ignoring qx.
-	double q = qy;//sqrt(qx*qx + qy*qy);
-	if (q < 0.0){
-		return 0.0;
-	}
-	return (*this).operator()(q);
+  // For 2D set qy as q, ignoring qx.
+  double q = qy;//sqrt(qx*qx + qy*qy);
+  if (q < 0.0){
+    return 0.0;
+  }
+  return (*this).operator()(q);
 }
 
@@ -213,5 +366,5 @@
  */
 double ReflAdvModel :: evaluate_rphi(double q, double phi) {
-	return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 
@@ -221,5 +374,5 @@
  */
 double ReflAdvModel :: calculate_ER() {
-	//NOT implemented yet!!!
-	return 0.0;
-}
+  //NOT implemented yet!!!
+  return 0.0;
+}
Index: nsmodels/src/python_wrapper/CReflAdvModel.cpp
===================================================================
--- sansmodels/src/python_wrapper/CReflAdvModel.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ 	(revision )
@@ -1,886 +1,0 @@
-/**
-	This software was developed by the University of Tennessee as part of the
-	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
-	project funded by the US National Science Foundation.
-
-	If you use DANSE applications to do scientific research that leads to
-	publication, we ask that you acknowledge the use of the software with the
-	following sentence:
-
-	"This work benefited from DANSE software developed under NSF award DMR-0520547."
-
-	copyright 2008, University of Tennessee
- */
-
-/** CReflAdvModel
- *
- * C extension 
- *
- * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
- *          DO NOT MODIFY THIS FILE, MODIFY refl_adv.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 "refl_adv.h"
-}
-
-#include "models.hh"
-#include "dispersion_visitor.hh"
-
-/// Error object for raised exceptions
-static PyObject * CReflAdvModelError = NULL;
-
-
-// Class definition
-typedef struct {
-    PyObject_HEAD
-    /// Parameters
-    PyObject * params;
-    /// Dispersion parameters
-    PyObject * dispersion;
-    /// Underlying model object
-    ReflAdvModel * model;
-    /// Log for unit testing
-    PyObject * log;
-} CReflAdvModel;
-
-
-static void
-CReflAdvModel_dealloc(CReflAdvModel* 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 *
-CReflAdvModel_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
-{
-    CReflAdvModel *self;
-    
-    self = (CReflAdvModel *)type->tp_alloc(type, 0);
-   
-    return (PyObject *)self;
-}
-
-static int
-CReflAdvModel_init(CReflAdvModel *self, PyObject *args, PyObject *kwds)
-{
-    if (self != NULL) {
-    	
-    	// Create parameters
-        self->params = PyDict_New();
-        self->dispersion = PyDict_New();
-        self->model = new ReflAdvModel();
-        
-        // 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",100.000000000000));
-        PyDict_SetItemString(self->params,"thick_flat3",Py_BuildValue("d",100.000000000000));
-        PyDict_SetItemString(self->params,"n_layers",Py_BuildValue("d",1.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,"sld_bottom0",Py_BuildValue("d",0.000002070000));
-        PyDict_SetItemString(self->params,"thick_inter3",Py_BuildValue("d",50.000000000000));
-        PyDict_SetItemString(self->params,"thick_inter4",Py_BuildValue("d",50.000000000000));
-        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,"npts_inter",Py_BuildValue("d",21.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",50.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));
-        PyDict_SetItemString(self->params,"func_inter0",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"sld_medium",Py_BuildValue("d",0.000001000000));
-        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,"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,"nu_inter10",Py_BuildValue("d",2.500000000000));
-        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,"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",50.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,"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 char name_params[] = "params";
-static char def_params[] = "Parameters";
-static char name_dispersion[] = "dispersion";
-static char def_dispersion[] = "Dispersion parameters";
-static char name_log[] = "log";
-static char def_log[] = "Log";
-
-static PyMemberDef CReflAdvModel_members[] = {
-    {name_params, T_OBJECT, offsetof(CReflAdvModel, params), 0, def_params},
-	{name_dispersion, T_OBJECT, offsetof(CReflAdvModel, dispersion), 0, def_dispersion},     
-    {name_log, T_OBJECT, offsetof(CReflAdvModel, log), 0, def_log},
-    {NULL}  /* Sentinel */
-};
-
-/** Read double from PyObject
-    @param p PyObject
-    @return double
-*/
-double CReflAdvModel_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(ReflAdvModel* model, PyArrayObject *q){
-    PyArrayObject *result;
-   
-    // Check validity of array q , q must be of dimension 1, an array of double
-    if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE)
-    {
-        //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num;
-        //PyErr_SetString(PyExc_ValueError , message);
-        return NULL;
-    }
-    result = (PyArrayObject *)PyArray_FromDims(q->nd, (int *)(q->dimensions), PyArray_DOUBLE);
-	if (result == NULL) {
-        const char * message= "Could not create result ";
-        PyErr_SetString(PyExc_RuntimeError , message);
-		return NULL;
-	}
-#pragma omp parallel for
-	 for (int i = 0; i < q->dimensions[0]; i++){
-      double q_value  = *(double *)(q->data + i*q->strides[0]);
-      double *result_value = (double *)(result->data + i*result->strides[0]);
-      *result_value =(*model)(q_value);
-	}
-    return PyArray_Return(result); 
- }
-
- /**
- * Function to call to evaluate model
- * @param args: input numpy array  [x[],y[]]
- * @return: numpy array object 
- */
- static PyObject * evaluateTwoDimXY( ReflAdvModel* model, 
-                              PyArrayObject *x, PyArrayObject *y)
- {
-    PyArrayObject *result;
-    int x_len, y_len, dims[1];
-    //check validity of input vectors
-    if (x->nd != 1 || x->descr->type_num != PyArray_DOUBLE
-        || y->nd != 1 || y->descr->type_num != PyArray_DOUBLE
-        || y->dimensions[0] != x->dimensions[0]){
-        const char * message= "evaluateTwoDimXY  expect 2 numpy arrays";
-        PyErr_SetString(PyExc_ValueError , message); 
-        return NULL;
-    }
-   
-	if (PyArray_Check(x) && PyArray_Check(y)) {
-		
-	    x_len = dims[0]= x->dimensions[0];
-        y_len = dims[0]= y->dimensions[0];
-	    
-	    // Make a new double matrix of same dims
-        result=(PyArrayObject *) PyArray_FromDims(1,dims,NPY_DOUBLE);
-        if (result == NULL){
-	    const char * message= "Could not create result ";
-        PyErr_SetString(PyExc_RuntimeError , message);
-	    return NULL;
-	    }
-       
-        /* Do the calculation. */
-#pragma omp parallel for
-        for (int i=0; i< x_len; i++) {
-            double x_value = *(double *)(x->data + i*x->strides[0]);
-  		    double y_value = *(double *)(y->data + i*y->strides[0]);
-  			double *result_value = (double *)(result->data +
-  			      i*result->strides[0]);
-  			*result_value = (*model)(x_value, y_value);
-        }           
-        return PyArray_Return(result); 
-        
-        }else{
-		    PyErr_SetString(CReflAdvModelError, 
-                   "CReflAdvModel.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(CReflAdvModel *self, PyObject *args){
-	PyObject *qx, *qy;
-	PyArrayObject * pars;
-	int npars ,mpars;
-	
-	// Get parameters
-	
-	    // 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") );
-    self->model->n_layers = PyFloat_AsDouble( PyDict_GetItemString(self->params, "n_layers") );
-    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->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") );
-    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->npts_inter = PyFloat_AsDouble( PyDict_GetItemString(self->params, "npts_inter") );
-    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") );
-    self->model->func_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter0") );
-    self->model->sld_medium = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_medium") );
-    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->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->nu_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter10") );
-    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->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->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->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(CReflAdvModelError, 
-	    	"CReflAdvModel.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(CReflAdvModelError, 
-                   "CReflAdvModel.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(CReflAdvModelError, 
-	    		"CReflAdvModel.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(CReflAdvModelError, 
-                   "CReflAdvModel.evalDistribution expect 2 numpy arrays in list.");
-	        return NULL;
-	     }
-	}
-	PyErr_SetString(CReflAdvModelError, 
-                   "CReflAdvModel.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(CReflAdvModel *self, PyObject *args) {
-	double q_value, phi_value;
-	PyObject* pars;
-	int npars;
-	
-	// Get parameters
-	
-	    // 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") );
-    self->model->n_layers = PyFloat_AsDouble( PyDict_GetItemString(self->params, "n_layers") );
-    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->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") );
-    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->npts_inter = PyFloat_AsDouble( PyDict_GetItemString(self->params, "npts_inter") );
-    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") );
-    self->model->func_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter0") );
-    self->model->sld_medium = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_medium") );
-    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->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->nu_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter10") );
-    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->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->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->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(CReflAdvModelError, 
-	    	"CReflAdvModel.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(CReflAdvModelError, 
-	    		"CReflAdvModel.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 = CReflAdvModel_readDouble(PyList_GET_ITEM(pars,0));
-	    phi_value = CReflAdvModel_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 = CReflAdvModel_readDouble(pars);		
-		
-		return Py_BuildValue("d",(*(self->model))(q_value));
-	}	
-}
-/**
- * Function to call to calculate_ER
- * @return: effective radius value 
- */
-static PyObject * calculate_ER(CReflAdvModel *self) {
-
-	// Get parameters
-	
-	    // 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") );
-    self->model->n_layers = PyFloat_AsDouble( PyDict_GetItemString(self->params, "n_layers") );
-    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->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") );
-    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->npts_inter = PyFloat_AsDouble( PyDict_GetItemString(self->params, "npts_inter") );
-    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") );
-    self->model->func_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter0") );
-    self->model->sld_medium = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_medium") );
-    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->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->nu_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter10") );
-    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->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->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->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(CReflAdvModel *self, PyObject *args) {
-	double qx_value, qy_value;
-	PyObject* pars;
-	int npars;
-	
-	// Get parameters
-	
-	    // 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") );
-    self->model->n_layers = PyFloat_AsDouble( PyDict_GetItemString(self->params, "n_layers") );
-    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->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") );
-    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->npts_inter = PyFloat_AsDouble( PyDict_GetItemString(self->params, "npts_inter") );
-    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") );
-    self->model->func_inter0 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "func_inter0") );
-    self->model->sld_medium = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_medium") );
-    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->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->nu_inter10 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nu_inter10") );
-    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->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->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->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(CReflAdvModelError, 
-	    	"CReflAdvModel.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(CReflAdvModelError, 
-	    		"CReflAdvModel.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 = CReflAdvModel_readDouble(PyList_GET_ITEM(pars,0));
-	    qy_value = CReflAdvModel_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 = CReflAdvModel_readDouble(pars);		
-		
-		return Py_BuildValue("d",(*(self->model))(qx_value));
-	}	
-}
-
-static PyObject * reset(CReflAdvModel *self, PyObject *args) {
-    
-
-    return Py_BuildValue("d",0.0);
-}
-
-static PyObject * set_dispersion(CReflAdvModel *self, PyObject *args) {
-	PyObject * disp;
-	const char * par_name;
-
-	if ( !PyArg_ParseTuple(args,"sO", &par_name, &disp) ) {
-	    PyErr_SetString(CReflAdvModelError,
-	    	"CReflAdvModel.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(CReflAdvModelError,
-	    	"CReflAdvModel.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 CReflAdvModel_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 CReflAdvModelType = {
-    PyObject_HEAD_INIT(NULL)
-    0,                         /*ob_size*/
-    "CReflAdvModel",             /*tp_name*/
-    sizeof(CReflAdvModel),             /*tp_basicsize*/
-    0,                         /*tp_itemsize*/
-    (destructor)CReflAdvModel_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*/
-    "CReflAdvModel objects",           /* tp_doc */
-    0,		               /* tp_traverse */
-    0,		               /* tp_clear */
-    0,		               /* tp_richcompare */
-    0,		               /* tp_weaklistoffset */
-    0,		               /* tp_iter */
-    0,		               /* tp_iternext */
-    CReflAdvModel_methods,             /* tp_methods */
-    CReflAdvModel_members,             /* tp_members */
-    0,                         /* tp_getset */
-    0,                         /* tp_base */
-    0,                         /* tp_dict */
-    0,                         /* tp_descr_get */
-    0,                         /* tp_descr_set */
-    0,                         /* tp_dictoffset */
-    (initproc)CReflAdvModel_init,      /* tp_init */
-    0,                         /* tp_alloc */
-    CReflAdvModel_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 addCReflAdvModel(PyObject *module) {
-	PyObject *d;
-	
-    if (PyType_Ready(&CReflAdvModelType) < 0)
-        return;
-
-    Py_INCREF(&CReflAdvModelType);
-    PyModule_AddObject(module, "CReflAdvModel", (PyObject *)&CReflAdvModelType);
-    
-    d = PyModule_GetDict(module);
-    static char error_name[] = "CReflAdvModel.error";
-    CReflAdvModelError = PyErr_NewException(error_name, NULL, NULL);
-    PyDict_SetItemString(d, "CReflAdvModelError", CReflAdvModelError);
-}
-
