Index: sansmodels/src/c_extensions/refl.c
===================================================================
--- sansmodels/src/c_extensions/refl.c	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ 	(revision )
@@ -1,215 +1,0 @@
-/**
- *  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>
-
-#define lamda 4.62
-
-
-double re_kernel(double dp[], double q) {
-	int n = dp[0];
-	int i,j;
-
-	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 total_thick=0.0;
-  double nsl=21.0; //nsl = Num_sub_layer:
-  int n_s;
-  double sld_i,dz,phi,R,ko2;
-  double fun;
-  double pi;
-
-  double* sld;
-  double* thick_inter;
-  double* thick;
-  int*fun_type;
-  complex  phi1,alpha,alpha2,kn,fnm,fnp,rn,Xn,nn,nn2,an,nnp1,one,two,n_sub,n_sup,knp1,Xnp1;
-
-  sld = (double*)malloc((n+2)*sizeof(double));
-  thick_inter = (double*)malloc((n+2)*sizeof(double));
-  thick = (double*)malloc((n+2)*sizeof(double));
-  fun_type = (int*)malloc((n+2)*sizeof(int));
-
-  fun_type[0] = dp[3];
-	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];
-		total_thick += thick[i] + thick_inter[i];
-	}
-	sld[0] = sld_sub;
-	sld[n+1] = sld_super;
-
-	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;
-
-	pi = 4.0*atan(1.0);
-  	Xn = cassign(0.0,0.0);
-  	one = cassign(1.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),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);
-
-    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_type[i-1]==1)
-    		fun = 5;
-    	else
-			fun = 0;
-		//iteration for 9 sub-layers
-		for (j=0;j<2;j++){
-			for (n_s=0;n_s<nsl; n_s++){
-				if (j==1){
-					if (i==n+1)
-						break;
-					dz = thick[i];
-					sld_i = sld[i];
-				}
-				else{
-					dz = thick_inter[i-1]/nsl;//nsl;
-					if (sld[i-1] == sld[i]){
-						sld_i = sld[i];
-					}
-					else{
-						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),0.0);
-				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;
-
-				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.0;
-    }
-    R *= scale;
-    R += background;
-
-    free(sld);
-    free(thick_inter);
-    free(thick);
-    free(fun_type);
-
-    return R;
-
-}
-/**
- * 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[47];
-
-	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;
-
-	return re_kernel(dp, q);
-}
-
-/**
- * Function to evaluate NR function
- * @param pars: parameters of NR
- * @param q: q-value
- * @return: function value
- */
-double refl_analytical_2D(ReflParameters *pars, double q, double phi) {
-	return refl_analytical_1D(pars,q);
-}
-
-double refl_analytical_2DXY(ReflParameters *pars, double qx, double qy){
-	return refl_analytical_1D(pars,sqrt(qx*qx+qy*qy));
-}
Index: sansmodels/src/c_extensions/refl.h
===================================================================
--- sansmodels/src/c_extensions/refl.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/refl.h	(revision 2d1b7008c4fcc2791722e1581f47a3559d64f820)
@@ -3,4 +3,5 @@
 #if !defined(o_h)
 #define refl_h
+#include "parameters.hh"
 
 /**
@@ -31,122 +32,123 @@
  //[ORIENTATION_PARAMS]= <text> </text>
 
-typedef struct {
-	/// number of layers
-	//  [DEFAULT]=n_layers=1
-	int n_layers;
+class ReflModel{
+public:
+  // Model parameters
+  /// number of layers
+  //  [DEFAULT]=n_layers=1
+  Parameter n_layers;
     /// Scale factor
     //  [DEFAULT]=scale= 1.0
-	double scale;
-    ///	thick_inter0 [A]
+  Parameter scale;
+    /// thick_inter0 [A]
     //  [DEFAULT]=thick_inter0=1.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;
+  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)]
-    double sld_flat1;
+    Parameter sld_flat1;
     //  [DEFAULT]=sld_flat2=3.5e-06 [1/A^(2)]
-    double sld_flat2;
+    Parameter sld_flat2;
     //  [DEFAULT]=sld_flat3=4.0e-06 [1/A^(2)]
-    double sld_flat3;
+    Parameter sld_flat3;
     //  [DEFAULT]=sld_flat4=3.5e-06 [1/A^(2)]
-    double sld_flat4;
+    Parameter sld_flat4;
     //  [DEFAULT]=sld_flat5=4.0e-06 [1/A^(2)]
-    double sld_flat5;
+    Parameter sld_flat5;
     //  [DEFAULT]=sld_flat6=3.5e-06 [1/A^(2)]
-    double sld_flat6;
+    Parameter sld_flat6;
     //  [DEFAULT]=sld_flat7=4.0e-06 [1/A^(2)]
-    double sld_flat7;
+    Parameter sld_flat7;
     //  [DEFAULT]=sld_flat8=3.5e-06 [1/A^(2)]
-    double sld_flat8;
+    Parameter sld_flat8;
     //  [DEFAULT]=sld_flat9=4.0e-06 [1/A^(2)]
-    double sld_flat9;
+    Parameter sld_flat9;
     //  [DEFAULT]=sld_flat10=3.5e-06 [1/A^(2)]
-    double sld_flat10;
+    Parameter sld_flat10;
 
     //  [DEFAULT]=thick_inter1=1 [A]
-    double thick_inter1;
+    Parameter thick_inter1;
     //  [DEFAULT]=thick_inter2=1 [A]
-    double thick_inter2;
+    Parameter thick_inter2;
     //  [DEFAULT]=thick_inter3=1 [A]
-    double thick_inter3;
+    Parameter thick_inter3;
     //  [DEFAULT]=thick_inter4=1 [A]
-    double thick_inter4;
+    Parameter thick_inter4;
     //  [DEFAULT]=thick_inter5=1 [A]
-    double thick_inter5;
+    Parameter thick_inter5;
     //  [DEFAULT]=thick_inter6=1 [A]
-    double thick_inter6;
+    Parameter thick_inter6;
     //  [DEFAULT]=thick_inter7=1 [A]
-    double thick_inter7;
+    Parameter thick_inter7;
     //  [DEFAULT]=thick_inter8=1 [A]
-    double thick_inter8;
+    Parameter thick_inter8;
     //  [DEFAULT]=thick_inter9=1 [A]
-    double thick_inter9;
+    Parameter thick_inter9;
     //  [DEFAULT]=thick_inter10=1 [A]
-    double thick_inter10;
+    Parameter thick_inter10;
 
     //  [DEFAULT]=thick_flat1=10 [A]
-    double thick_flat1;
+    Parameter thick_flat1;
     //  [DEFAULT]=thick_flat2=100 [A]
-    double thick_flat2;
+    Parameter thick_flat2;
     //  [DEFAULT]=thick_flat3=100 [A]
-    double thick_flat3;
+    Parameter thick_flat3;
     //  [DEFAULT]=thick_flat4=100 [A]
-    double thick_flat4;
+    Parameter thick_flat4;
     //  [DEFAULT]=thick_flat5=100 [A]
-    double thick_flat5;
+    Parameter thick_flat5;
     //  [DEFAULT]=thick_flat6=100 [A]
-    double thick_flat6;
+    Parameter thick_flat6;
     //  [DEFAULT]=thick_flat7=100 [A]
-    double thick_flat7;
+    Parameter thick_flat7;
     //  [DEFAULT]=thick_flat8=100 [A]
-    double thick_flat8;
+    Parameter thick_flat8;
     //  [DEFAULT]=thick_flat9=100 [A]
-    double thick_flat9;
+    Parameter thick_flat9;
     //  [DEFAULT]=thick_flat10=100 [A]
-    double thick_flat10;
+    Parameter thick_flat10;
 
     //  [DEFAULT]=func_inter1=0
-    double func_inter1;
+    Parameter func_inter1;
     //  [DEFAULT]=func_inter2=0
-    double func_inter2;
+    Parameter func_inter2;
     //  [DEFAULT]=func_inter3=0
-    double func_inter3;
+    Parameter func_inter3;
     //  [DEFAULT]=func_inter4=0
-    double func_inter4;
+    Parameter func_inter4;
     //  [DEFAULT]=func_inter5=0
-    double func_inter5;
+    Parameter func_inter5;
     //  [DEFAULT]=func_inter6=0
-    double func_inter6;
+    Parameter func_inter6;
     //  [DEFAULT]=func_inter7=0
-    double func_inter7;
+    Parameter func_inter7;
     //  [DEFAULT]=func_inter8=0
-    double func_inter8;
+    Parameter func_inter8;
     //  [DEFAULT]=func_inter9=0
-    double func_inter9;
+    Parameter func_inter9;
     //  [DEFAULT]=func_inter10=0
-    double func_inter10;
+    Parameter func_inter10;
 
 
-} ReflParameters;
+  // Constructor
+  ReflModel();
 
-double re_kernel(double dq[], double q);
-
-/// 1D scattering function
-double refl_analytical_1D(ReflParameters *pars, double q);
-
-/// 2D scattering function
-double refl_analytical_2D(ReflParameters *pars, double q, double phi);
-double refl_analytical_2DXY(ReflParameters *pars, double qx, double qy);
+  // 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
