Index: src/sas/models/c_extension/c_models/refl.cpp
===================================================================
--- src/sas/models/c_extension/c_models/refl.cpp	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
+++ src/sas/models/c_extension/c_models/refl.cpp	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
@@ -0,0 +1,303 @@
+
+#include <math.h>
+#include "parameters.hh"
+#include <stdio.h>
+#include <stdlib.h>
+#include "refl.h"
+using namespace std;
+
+extern "C" {
+#include "libmultifunc/librefl.h"
+}
+
+#define lamda 4.62
+
+static 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;
+
+}
+ReflModel :: ReflModel() {
+  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(3.0e-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);
+  thick_inter2 = Parameter(1);
+  thick_inter3 = Parameter(1);
+  thick_inter4 = Parameter(1);
+  thick_inter5 = Parameter(1);
+  thick_inter6 = Parameter(1);
+  thick_inter7 = Parameter(1);
+  thick_inter8 = Parameter(1);
+  thick_inter9 = Parameter(1);
+  thick_inter10 = Parameter(1);
+
+
+  thick_flat1 = Parameter(15);
+  thick_flat2 = Parameter(100);
+  thick_flat3 = Parameter(100);
+  thick_flat4 = Parameter(100);
+  thick_flat5 = Parameter(100);
+  thick_flat6 = Parameter(100);
+  thick_flat7 = Parameter(100);
+  thick_flat8 = Parameter(100);
+  thick_flat9 = Parameter(100);
+  thick_flat10 = Parameter(100);
+
+
+  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);
+
+}
+
+/**
+ * Function to evaluate 1D NR function
+ * @param q: q-value
+ * @return: function value
+ */
+double ReflModel :: operator()(double q) {
+  double dp[47];
+  // 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();
+
+  // Get the dispersion points for the radius
+  //vector<WeightPoint> weights_thick;
+  //thick_inter0.get_weights(weights_thick_inter0);
+
+
+  return re_kernel(dp,q);
+}
+
+/**
+ * Function to evaluate 2D NR function
+ * @param q_x: value of Q along x
+ * @param q_y: value of Q along y
+ * @return: function value
+ */
+double ReflModel :: 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);
+}
+
+/**
+ * Function to evaluate 2D NR function
+ * @param pars: parameters of the sphere
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double ReflModel :: evaluate_rphi(double q, double phi) {
+  return (*this).operator()(q);
+}
+
+/**
+ * Function to calculate effective radius
+ * @return: effective radius value
+ */
+double ReflModel :: calculate_ER() {
+  //NOT implemented yet!!!
+  return 0.0;
+}
+double ReflModel :: calculate_VR() {
+  return 1.0;
+}
