Index: sansmodels/src/sans/models/c_extensions/refl_adv.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/refl_adv.c	(revision 96656e3ad6cdb9f8256e2118c60bac29635136af)
+++ sansmodels/src/sans/models/c_extensions/refl_adv.c	(revision 96656e3ad6cdb9f8256e2118c60bac29635136af)
@@ -0,0 +1,245 @@
+// 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,fun_type[n+2];
+
+	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 sld[n+2],sld_im[n+2],thick_inter[n+2],thick[n+2],fun_coef[n+2],total_thick;
+	fun_type[0] = dp[3];
+	fun_coef[0] = fabs(dp[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;
+
+	double nsl=npts;//21.0; //nsl = Num_sub_layer:  MUST ODD number in double //no other number works now
+	int n_s;
+
+    double sld_i,sldim_i,dz,phi,R,ko2;
+    double sign,erfunc;
+    double pi;
+	complex  inv_n,phi1,alpha,alpha2,kn,fnm,fnp,rn,Xn,nn,nn2,an,nnp1,one,zero,two,n_sub,n_sup,knp1,Xnp1;
+	pi = 4.0*atan(1.0);
+    one = cassign(1.0,0.0);
+	//zero = 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 = cdiv(rcmult(phi,one),n_sup);
+    alpha = cmult(n_sup,ccos(phi1));
+	alpha2 = cmult(alpha,alpha);
+
+    nnp1=n_sub;
+    knp1=csqrt(rcmult(ko2,csub(cmult(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=cmult(nn,nn);
+
+				kn=csqrt(rcmult(ko2,csub(nn2,alpha2)));        //nn*ko*sin(phin)
+				an=cexp(rcmult(dz,cmult(two,kn)));
+
+				fnm=csub(kn,knp1);
+				fnp=cadd(kn,knp1);
+				rn=cdiv(fnm,fnp);
+				Xn=cmult(an,cdiv(cadd(rn,Xnp1),cadd(one,cmult(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;
+
+    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/sans/models/c_extensions/refl_adv.h
===================================================================
--- sansmodels/src/sans/models/c_extensions/refl_adv.h	(revision 96656e3ad6cdb9f8256e2118c60bac29635136af)
+++ sansmodels/src/sans/models/c_extensions/refl_adv.h	(revision 96656e3ad6cdb9f8256e2118c60bac29635136af)
@@ -0,0 +1,203 @@
+// The original code, of which work was not DANSE funded,
+// was provided by J. Cho.
+#if !defined(o_h)
+#define refl_adv_h
+
+/**
+ * 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);
+
+#endif
