diff --git a/src/CalculateE.c b/src/CalculateE.c
index c6c41a87..5a5cc175 100644
--- a/src/CalculateE.c
+++ b/src/CalculateE.c
@@ -107,7 +107,7 @@ static void ComputeMuellerMatrix(double matrix[4][4], const doublecomplex s1,con
 
 	if (surface) {
 		double scale=inc_scale;
-		if (TestBelowDeg(theta)) scale*=creal(msub);
+		if (TestBelowDeg(theta)) scale*=creal(sub.m[sub.N-1]);
 		if (fabs(scale-1)>ROUND_ERR) for (int i=0;i<4;i++) for (int j=0;j<4;j++) matrix[i][j]*=scale;
 	}
 }
diff --git a/src/GenerateB.c b/src/GenerateB.c
index 100bbb03..3b05b5b4 100644
--- a/src/GenerateB.c
+++ b/src/GenerateB.c
@@ -50,8 +50,8 @@ int vorticity;                 // Vorticity of vortex beams (besN for Bessel bea
 // used in param.c
 const char *beam_descr; // string for log file with beam parameters
 /* Propagation (phase) directions of secondary incident beams above (refl) and below (tran) the surface (unit vectors)
- * When msub is complex, one of this doesn't tell the complete story, since the corresponding wave is inhomogeneous,
- * given by the complex wavenumber ktVec
+ * When the reflactive index of the surface is complex, one of this doesn't tell the complete story,
+ * since the corresponding wave is inhomogeneous, given by the complex wavenumber ktVec
  */
 double prIncRefl[3],prIncTran[3];
 
@@ -72,10 +72,10 @@ static doublecomplex besM[4];     // Matrix M defining the generalized Bessel be
  * afterwards. If you need local, intermediate variables, put them into the beginning of the corresponding function.
  * Add descriptive comments, use 'static'.
  */
- 
+
 // EXTERNAL FUNCTIONS
 
-#ifndef NO_FORTRAN 
+#ifndef NO_FORTRAN
 void bjndd_(const int *n,const double *x,double *bj,double *dj,double *fj);
 #endif
 
@@ -89,7 +89,7 @@ void InitBeam(void)
 	/* TO ADD NEW BEAM
 	 * Add here all intermediate variables, which are used only inside this function.
 	 */
-	
+
 	// initialization of global option index for error messages
 	opt=opt_beam;
 	// beam initialization
@@ -100,35 +100,38 @@ void InitBeam(void)
 		case B_PLANE:
 			if (IFROOT) beam_descr="plane wave";
 			if (surface) {
-				/* here we assume that prop_0 will not change further (e.g., by rotation of particle), 
+				/* here we assume that prop_0 will not change further (e.g., by rotation of particle),
 				 * i.e. prop=prop_0 in GenerateBeam() below
-				 */				
+				 */
 				if (prop_0[2]==0) PrintError("Ambiguous setting of beam propagating along the surface. Please specify "
 					"the incident direction to have (arbitrary) small positive or negative z-component");
-				if (msubInf && prop_0[2]>0) PrintError("Perfectly reflecting surface ('-surf ... inf') is incompatible "
+				if (sub.mInf && prop_0[2]>0) PrintError("Perfectly reflecting surface ('-surf ... inf') is incompatible "
 					"with incident direction from below (including the default one)");
 				// Here we set ki,kt,ktVec and propagation directions prIncRefl,prIncTran
 				if (prop_0[2]>0) { // beam comes from the substrate (below)
-					// here msub should always be defined
-					inc_scale=1/creal(msub);
-					ki=msub*prop_0[2];
-					/* Special case for msub near 1 to remove discontinuities for near-grazing incidence. The details
+					// here sub.m[sub.N-1] should always be defined
+					inc_scale=1/creal(sub.m[sub.N-1]);
+					ki=sub.m[sub.N-1]*prop_0[2];
+					/* Special case for sub.m[sub.N-1] near 1 to remove discontinuities for near-grazing incidence. The details
 					 * are discussed in CalcFieldSurf() in crosssec.c.
 					 */
-					if (cabs(msub-1)<ROUND_ERR && cabs(ki)<SQRT_RND_ERR) kt=ki;
-					else kt=cSqrtCut(1 - msub*msub*(prop_0[0]*prop_0[0]+prop_0[1]*prop_0[1]));
+                    kt = CalculateKt(
+                            ki,
+                            sub.m[sub.N-1],
+                            1,
+                            sub.m[sub.N-1]*sub.m[sub.N-1]*(prop_0[0]*prop_0[0]+prop_0[1]*prop_0[1])
+                            );
 					// determine propagation direction and full wavevector of wave transmitted into substrate
-					ktVec[0]=msub*prop_0[0];
-					ktVec[1]=msub*prop_0[1];
+					ktVec[0]=sub.m[sub.N-1]*prop_0[0];
+					ktVec[1]=sub.m[sub.N-1]*prop_0[1];
 					ktVec[2]=kt;
 				}
 				else if (prop_0[2]<0) { // beam comes from above the substrate
 					inc_scale=1;
 					ki=-prop_0[2]; // always real
-					if (!msubInf) {
+					if (!sub.mInf) {
 						// same special case as above
-						if (cabs(msub-1)<ROUND_ERR && cabs(ki)<SQRT_RND_ERR) kt=ki;
-						else kt=cSqrtCut(msub*msub - (prop_0[0]*prop_0[0]+prop_0[1]*prop_0[1]));
+                        kt = CalculateKt(ki, 1, sub.m[sub.N-1], prop_0[0]*prop_0[0]+prop_0[1]*prop_0[1]);
 						// determine propagation direction of wave transmitted into substrate
 						ktVec[0]=prop_0[0];
 						ktVec[1]=prop_0[1];
@@ -138,7 +141,7 @@ void InitBeam(void)
 				else LogError(ONE_POS,"Ambiguous setting of beam propagating along the surface. Please specify the"
 					"incident direction to have (arbitrary) small positive or negative z-component");
 				vRefl(prop_0,prIncRefl);
-				if (!msubInf) {
+				if (!sub.mInf) {
 					vReal(ktVec,prIncTran);
 					vNormalize(prIncTran);
 				}
@@ -146,7 +149,7 @@ void InitBeam(void)
 			return;
 		case B_DIPOLE:
 			if (surface) {
-				if (beam_center_0[2]<=-hsub)
+				if (beam_center_0[2]<=-sub.hP)
 					PrintErrorHelp("External dipole should be placed strictly above the surface");
 				inc_scale=1; // but scaling of Mueller matrix is weird anyway
 			}
@@ -354,25 +357,32 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 			if (surface) {
 				/* With respect to normalization we use here the same assumption as in the free space - the origin is in
 				 * the particle center, and amplitude of incoming plane wave is equal to 1. Then irradiance of the beam
-				 * coming from below is c*Re(msub)/(8pi), different from that coming from above.
+				 * coming from below is c*Re(sub.m[sub.N-1])/(8pi), different from that coming from above.
 				 * Original incident (incoming) beam propagating from the vacuum (above) is Exp(i*k*r.a), while - from
-				 * the substrate (below) is Exp(i*k*msub*r.a). We assume that the incoming beam is homogeneous in its
+				 * the substrate (below) is Exp(i*k*sub.m[sub.N-1]*r.a). We assume that the incoming beam is homogeneous in its
 				 * original medium.
 				 */
-				double hbeam=hsub+beam_center[2]; // height of beam center above the surface 
+				double hbeam=sub.hP+beam_center[2]; // height of beam center above the surface
 				if (prop[2]>0) { // beam comes from the substrate (below)
 					doublecomplex tc; // transmission coefficients
-					//  determine amplitude of the transmitted wave; here msub is always defined
+                    //  determine amplitude of the reflected and transmitted waves; here sub.m[sub.N-1] is always defined
 					if (which==INCPOL_Y) { // s-polarized
 						cvBuildRe(ex,eIncTran);
-						tc=FresnelTS(ki,kt);
+						SubstrateFresnel(sub,WaveNum,true,
+										 sub.m[sub.N-1]*sub.m[sub.N-1]*(prop_0[0]*prop_0[0]+prop_0[1]*prop_0[1]),ki,&tc,
+										 NULL,NULL,NULL, NULL);
 					}
 					else { // p-polarized
 						crCrossProd(ey,ktVec,eIncTran);
-						tc=FresnelTP(ki,kt,1/msub);
+						SubstrateFresnel(sub,WaveNum,true,
+										 sub.m[sub.N-1]*sub.m[sub.N-1]*(prop_0[0]*prop_0[0]+prop_0[1]*prop_0[1]),ki,
+										 NULL,NULL,&tc,NULL,NULL);
 					}
+					double hbeam_eff = hbeam;
+					for (int k = 0; k < sub.N - 1; ++k)
+						hbeam_eff += sub.h[k];
 					// phase shift due to the beam center relative to the origin and surface
-					cvMultScal_cmplx(tc*cexp(I*WaveNum*((kt-ki)*hbeam-crDotProd(ktVec,beam_center))),eIncTran,eIncTran);
+					cvMultScal_cmplx(tc*cexp(I*WaveNum*((kt*hbeam -ki*hbeam_eff)-crDotProd(ktVec,beam_center))),eIncTran,eIncTran);
 					// main part
 					for (i=0;i<local_nvoid_Ndip;i++) {
 						j=3*i;
@@ -390,13 +400,14 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 					// determine reflection coefficient + 1
 					doublecomplex rcp1;
 					if (which==INCPOL_Y) { // s-polarized
-						if (msubInf) rcp1=0;
-						else rcp1=FresnelTS(ki,kt);
+						SubstrateFresnel(sub,WaveNum,false,
+										 prop_0[0]*prop_0[0]+prop_0[1]*prop_0[1],ki,NULL,&rcp1,NULL,NULL, NULL);
 					}
 					else { // p-polarized
-						if (msubInf) rcp1=2;
-						else rcp1=msub*FresnelTP(ki,kt,msub);
+						SubstrateFresnel(sub,WaveNum,false,
+								prop_0[0]*prop_0[0]+prop_0[1]*prop_0[1],ki,NULL,NULL,NULL,&rcp1, NULL);
 					}
+					rcp1 = rcp1 + 1;
 					// main part
 					for (i=0;i<local_nvoid_Ndip;i++) {
 						j=3*i;
@@ -440,7 +451,7 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 				(*InterTerm_real)(r1,gt);
 				cSymMatrVecReal(gt,dip_p,b+j);
 				if (surface) { // add reflected field
-					r1[2]=DipoleCoord[j+2]+beam_center[2]+2*hsub;
+					r1[2]=DipoleCoord[j+2]+beam_center[2]+2*sub.hP;
 					(*ReflTerm_real)(r1,gt);
 					cReflMatrVecReal(gt,dip_p,v1);
 					cvAdd(v1,b+j,b+j);
@@ -457,7 +468,7 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 			C0dipole=2*FOUR_PI_OVER_THREE*temp*temp;
 			if (surface) {
 				r1[0]=r1[1]=0;
-				r1[2]=2*(beam_center[2]+hsub);
+				r1[2]=2*(beam_center[2]+sub.hP);
 				(*ReflTerm_real)(r1,gt);
 				double tmp;
 				/* the following expression uses that dip_p is real and a specific (anti-)symmetry of the gt
@@ -635,7 +646,7 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 				cvMultScal_cmplx(ctemp,v1,b+j);
 			}
 			return;
-#endif // !NO_FORTRAN		
+#endif // !NO_FORTRAN
 		case B_READ:
 			if (which==INCPOL_Y) fname=beam_fnameY;
 			else fname=beam_fnameX; // which==INCPOL_X
diff --git a/src/calculator.c b/src/calculator.c
index 6a2daf3a..cfeb9f8e 100644
--- a/src/calculator.c
+++ b/src/calculator.c
@@ -60,7 +60,8 @@ double * restrict muel_alpha; // mueller matrix for different values of alpha
 
 // used in crosssec.c
 doublecomplex * restrict E_ad; // complex field E, calculated for alldir
-double * restrict E2_alldir; // square of E (scaled with msub, so ~ Poynting vector or dC/dOmega), calculated for alldir
+double * restrict E2_alldir; /* square of E (scaled with the refractive index of the last layer, so ~ Poynting vector or
+                                dC/dOmega), calculated for alldir */
 doublecomplex cc[MAX_NMAT][3]; // couple constants
 #ifndef SPARSE
 doublecomplex * restrict expsX,* restrict expsY,* restrict expsZ; // arrays of exponents along 3 axes (for calc_field)
diff --git a/src/cmplx.c b/src/cmplx.c
index 774425f7..e856b258 100644
--- a/src/cmplx.c
+++ b/src/cmplx.c
@@ -239,3 +239,234 @@ void imExp_arr(const doublecomplex arg,const int size,doublecomplex *c)
 		}
 	}
 }
+
+static inline void set_mutual_S_coef(
+		const doublecomplex ki,
+		const doublecomplex km,
+		doublecomplex* t_im,
+		doublecomplex* r_mi
+		)
+/* A helper function for a following SubstrateFresnel function.
+ * Calculates s-polarization fresnel coefficients required for all possible combinations of output coefficients in a
+ * two-layer substrate case
+ */
+{
+	*t_im = FresnelTS(ki, km);
+	*r_mi = FresnelRS(km, ki);
+}
+
+static inline void set_mutual_P_coef(
+		const doublecomplex ki,
+		const doublecomplex km,
+		const doublecomplex m_i,
+		const doublecomplex m_m,
+		doublecomplex* t_im,
+		doublecomplex* r_mi
+		)
+/* A helper function for the following SubstrateFresnel function.
+ * Calculates p-polarization fresnel coefficients required for all possible combinations of output coefficients in a
+ * two-layer substrate case
+ */
+{
+	*t_im = FresnelTP(ki, km, m_m/m_i);
+	*r_mi = FresnelRP(km, ki, m_i/m_m);
+}
+
+static inline void set_reflection_S_coef(
+		const doublecomplex ki,
+		const doublecomplex km,
+		doublecomplex* t_mi,
+		doublecomplex* r_im
+		)
+/* A helper function for the following SubstrateFresnel function.
+ * Calculates fresnel coefficients required for the output s-polarization reflection coefficient in a two-layer
+ * substrate case
+ */
+{
+	*t_mi = FresnelTS(km, ki);
+	*r_im = FresnelRS(ki, km);
+}
+
+static inline void set_reflection_P_coef(
+		const doublecomplex ki,
+		const doublecomplex km,
+		const doublecomplex m_i,
+		const doublecomplex m_m,
+		doublecomplex* t_mi,
+		doublecomplex* r_im
+		)
+/* A helper function for the following SubstrateFresnel function.
+ * Calculates fresnel coefficients required for the output p-polarization reflection coefficient in a two-layer
+ * substrate case
+ */
+{
+	*t_mi = FresnelTP(km, ki, m_i / m_m);
+	*r_im = FresnelRP(ki, km, m_m / m_i);
+}
+
+static inline doublecomplex get_two_layer_t(
+		const doublecomplex t_im,
+		const doublecomplex t_mt,
+		const doublecomplex r_mi,
+		const doublecomplex r_mt,
+		const doublecomplex eL
+		)
+/* A helper function for the following SubstrateFresnel function.
+ * Calculates transmission coefficient in a two-layer substrate case
+ */
+{
+	return t_mt * t_im / (1 - r_mi * r_mt * eL * eL) * eL;
+}
+
+static inline doublecomplex get_two_layer_r(
+		const doublecomplex t_im,
+		const doublecomplex t_mi,
+		const doublecomplex r_im,
+		const doublecomplex r_mi,
+		const doublecomplex r_mt,
+		const doublecomplex eL
+		)
+/* A helper function for the following SubstrateFresnel function.
+ * Calculates reflection coefficient in a two-layer substrate case
+ */
+{
+	return r_im + t_mi * r_mt * t_im * eL * eL / (1 - r_mi * r_mt * eL * eL);
+}
+void SubstrateFresnel(
+		const struct Substrate sub,
+		const double wave_num,
+		const bool is_positive_z_direction,
+		const doublecomplex sqr_long_k,
+		const doublecomplex ki,
+		doublecomplex * const ts_out,
+		doublecomplex * const rs_out,
+		doublecomplex * const tp_out,
+		doublecomplex * const rp_out,
+		doublecomplex * const kt_out
+		)
+{
+	if (sub.N > 2)
+		PrintError("substrates with more then 2 layers are not supported yet");
+
+	if (sub.N == 1) {
+		if (sub.mInf) {
+			if (ts_out!=NULL) *ts_out = 0;
+			if (tp_out!=NULL) *tp_out = 0;
+			if (rs_out!=NULL) *rs_out = -1;
+			if (rp_out!=NULL) *rp_out = 1;
+			return;
+		}
+		doublecomplex mi = is_positive_z_direction ? sub.m[0] : 1;
+		doublecomplex mt = is_positive_z_direction ? 1 : sub.m[0];
+		doublecomplex kt = CalculateKt(ki, mi, mt, sqr_long_k);
+		if (kt_out!=NULL) *kt_out = kt;
+		if (ts_out!=NULL) *ts_out = FresnelTS(ki, kt);
+		if (tp_out!=NULL) *tp_out = FresnelTP(ki, kt, mt/mi);
+		if (rs_out!=NULL) *rs_out = FresnelRS(ki, kt);
+		if (rp_out!=NULL) *rp_out = FresnelRP(ki, kt, mt/mi);
+		return;
+	}
+
+	// Two layer case:
+
+	doublecomplex m_i, m_m, m_t;
+	doublecomplex km, kt;
+	doublecomplex eL;
+	doublecomplex t_im, t_mi, t_mt, r_mt, r_mi, r_im;
+
+	m_m = sub.m[0];
+	if (is_positive_z_direction) {
+		m_i = sub.m[1];
+		m_t = 1;
+	}
+	else {
+		m_i = 1;
+		if (!sub.mInf) m_t = sub.m[1];
+	}
+	km = CalculateKt(ki, m_i, m_m, sqr_long_k);
+	if (!sub.mInf) {
+		kt = CalculateKt(km, m_m, m_t, sqr_long_k);
+		if (kt_out!=NULL) *kt_out = kt;
+	}
+	eL = cexp(I * wave_num * sub.h[0] * km);
+
+	switch (((ts_out!=NULL)<<2)|((rs_out!=NULL)<<1)|sub.mInf) { // encoding configuration for output coefficients
+		case 7: // equals 0b111: set ts and rs (sub.mInf is true)
+			set_mutual_S_coef(ki, km, &t_im, &r_mi);
+			r_mt = -1;
+			set_reflection_S_coef(ki, km, &t_mi, &r_im);
+			*ts_out = 0;
+			*rs_out = get_two_layer_r(t_im, t_mi, r_im, r_mi, r_mt, eL);
+			break;
+		case 6: // equals 0b110: set ts and rs (sub.mInf is false)
+			set_mutual_S_coef(ki, km, &t_im, &r_mi);
+			r_mt = FresnelRS(km, kt);
+			set_reflection_S_coef(ki, km, &t_mi, &r_im);
+			t_mt = FresnelTS(km, kt);
+			*ts_out = get_two_layer_t(t_im, t_mt, r_mi, r_mt, eL);
+			*rs_out = get_two_layer_r(t_im, t_mi, r_im, r_mi, r_mt, eL);
+			break;
+		case 5: // equals 0b101: set ts (sub.mInf is true)
+			*ts_out = 0;
+		case 4: // equals 0b100: set ts (sub.mInf is false)
+			set_mutual_S_coef(ki, km, &t_im, &r_mi);
+			r_mt = FresnelRS(km, kt);
+			t_mt = FresnelTS(km, kt);
+			*ts_out = get_two_layer_t(t_im, t_mt, r_mi, r_mt, eL);
+			break;
+		case 3: // equals 0b011 set rs (sub.mInf is true)
+			set_mutual_S_coef(ki, km, &t_im, &r_mi);
+			r_mt = -1;
+			set_reflection_S_coef(ki, km, &t_mi, &r_im);
+			*rs_out = get_two_layer_r(t_im, t_mi, r_im, r_mi, r_mt, eL);
+			break;
+		case 2: // equals 0b010: set rs (sub.mInf is false)
+			set_mutual_S_coef(ki, km, &t_im, &r_mi);
+			r_mt = FresnelRS(km, kt);
+			set_reflection_S_coef(ki, km, &t_mi, &r_im);
+			*rs_out = get_two_layer_r(t_im, t_mi, r_im, r_mi, r_mt, eL);
+			break;
+		default:
+			break;
+	}
+
+	switch (((tp_out!=NULL)<<2)|((rp_out!=NULL)<<1)|sub.mInf) {
+		case 7:
+			set_mutual_P_coef(ki, km, m_i, m_m, &t_im, &r_mi);
+			r_mt = 1;
+			set_reflection_P_coef(ki, km, m_i, m_m, &t_mi, &r_im);
+			*tp_out = 0;
+			*rp_out = get_two_layer_r(t_im, t_mi, r_im, r_mi, r_mt, eL);
+			break;
+		case 6:
+			set_mutual_P_coef(ki, km, m_i, m_m, &t_im, &r_mi);
+			r_mt = FresnelRP(km, kt, m_t/m_m);
+			set_reflection_P_coef(ki, km, m_i, m_m, &t_mi, &r_im);
+			t_mt = FresnelTP(km, kt, m_t/m_m);
+			*tp_out = get_two_layer_t(t_im, t_mt, r_mi, r_mt, eL);
+			*rp_out = get_two_layer_r(t_im, t_mi, r_im, r_mi, r_mt, eL);
+			break;
+		case 5:
+			*tp_out = 0;
+		case 4:
+			set_mutual_P_coef(ki, km, m_i, m_m, &t_im, &r_mi);
+			r_mt = FresnelRP(km, kt, m_t/m_m);
+			t_mt = FresnelTP(km, kt, m_t/m_m);
+			*tp_out = get_two_layer_t(t_im, t_mt, r_mi, r_mt, eL);
+			break;
+		case 3:
+			set_mutual_P_coef(ki, km, m_i, m_m, &t_im, &r_mi);
+			r_mt = 1;
+			set_reflection_P_coef(ki, km, m_i, m_m, &t_mi, &r_im);
+			*rp_out = get_two_layer_r(t_im, t_mi, r_im, r_mi, r_mt, eL);
+			break;
+		case 2:
+			set_mutual_P_coef(ki, km, m_i, m_m, &t_im, &r_mi);
+			r_mt = FresnelRP(km, kt, m_t/m_m);
+			set_reflection_P_coef(ki, km, m_i, m_m, &t_mi, &r_im);
+			*rp_out = get_two_layer_r(t_im, t_mi, r_im, r_mi, r_mt, eL);
+			break;
+		default:
+			break;
+	}
+}
diff --git a/src/cmplx.h b/src/cmplx.h
index e6a988aa..a8bab70c 100644
--- a/src/cmplx.h
+++ b/src/cmplx.h
@@ -235,6 +235,14 @@ static inline double cvNorm2(const doublecomplex a[static 3])
 	return cAbs2(a[0]) + cAbs2(a[1]) + cAbs2(a[2]);
 }
 
+//======================================================================================================================
+
+static inline double cdNorm2(const doublecomplex a[static 6])
+// square of the norm of a complex dyad[6]
+{
+	return cAbs2(a[0]) + cAbs2(a[1]) + cAbs2(a[2]) + cAbs2(a[3]) + cAbs2(a[4]) + cAbs2(a[5]);
+}
+
 
 //======================================================================================================================
 
@@ -252,8 +260,8 @@ static inline double cDotProd_Im(const doublecomplex a[static 3],const doublecom
  */
 {
 	return ( cimag(a[0])*creal(b[0]) - creal(a[0])*cimag(b[0])
-	       + cimag(a[1])*creal(b[1]) - creal(a[1])*cimag(b[1])
-	       + cimag(a[2])*creal(b[2]) - creal(a[2])*cimag(b[2]) );
+		   + cimag(a[1])*creal(b[1]) - creal(a[1])*cimag(b[1])
+		   + cimag(a[2])*creal(b[2]) - creal(a[2])*cimag(b[2]) );
 }
 
 //======================================================================================================================
@@ -276,6 +284,19 @@ static inline void cvAdd(const doublecomplex a[static 3],const doublecomplex b[s
 
 //======================================================================================================================
 
+static inline void cdAdd(const doublecomplex a[static 6],const doublecomplex b[static 6],doublecomplex c[static 6])
+// add two complex dyad[6]; c=a+b;
+{
+    c[0] = a[0] + b[0];
+    c[1] = a[1] + b[1];
+    c[2] = a[2] + b[2];
+    c[3] = a[3] + b[3];
+    c[4] = a[4] + b[4];
+    c[5] = a[5] + b[5];
+}
+
+//======================================================================================================================
+
 static inline void cvSubtr(const doublecomplex a[static 3],const doublecomplex b[static 3],doublecomplex c[static 3])
 // add two complex vector[3]; c=a-b;
 {
@@ -556,7 +577,7 @@ static inline double QuadForm(const double matr[static 6],const double vec[stati
 // value of a quadratic form matr (symmetric matrix stored as a vector of size 6) over a vector vec;
 {
 	return ( vec[0]*vec[0]*matr[0] + vec[1]*vec[1]*matr[2] + vec[2]*vec[2]*matr[5]
-	       + 2*(vec[0]*vec[1]*matr[1] + vec[0]*vec[2]*matr[3] + vec[1]*vec[2]*matr[4]) );
+		   + 2*(vec[0]*vec[1]*matr[1] + vec[0]*vec[2]*matr[3] + vec[1]*vec[2]*matr[4]) );
 }
 
 //======================================================================================================================
@@ -677,6 +698,39 @@ static inline doublecomplex FresnelTP(const doublecomplex ki,const doublecomplex
 	return 2*mr*ki/(mr*mr*ki+kt);
 }
 
+static inline doublecomplex CalculateKt(
+		const doublecomplex ki,
+		const doublecomplex mi,
+		const doublecomplex mt,
+		const doublecomplex sqr_long_k /* square of the longitudinal component of kVec. The value is the same for all
+										* transitions due to the law of refraction:
+										* mi*sin(ai)==mt*sin(at)==sqrt(sqr_long_k)
+										*/
+		)
+{
+	if (cabs(mt-mi)<ROUND_ERR && cabs(ki)<SQRT_RND_ERR)
+		return ki;
+	else
+		return cSqrtCut(mt*mt - sqr_long_k);
+}
+
+void SubstrateFresnel(
+		struct Substrate sub,
+		double wave_num,
+		bool is_positive_z_direction, // indicates the direction of incidence field
+		doublecomplex sqr_long_k, // a square of longitudinal component of the k vector
+		doublecomplex ki,
+		// pointers to output values. All unnecessary pointers should be passed as NULL
+		doublecomplex * ts_out,
+		doublecomplex * rs_out,
+		doublecomplex * tp_out,
+		doublecomplex * rp_out,
+		doublecomplex * kt_out
+);
+/* A function to calculate the substrate transmission and/or reflection coefficients. kt may be stored in kt_out for
+ * optimization purposes
+ */
+
 #ifdef USE_SSE3
 
 //======================================================================================================================
diff --git a/src/const.h b/src/const.h
index f13bad94..bd48a896 100644
--- a/src/const.h
+++ b/src/const.h
@@ -108,6 +108,7 @@ the compilation may fail or produce wrong results. If you still want to try, ena
 #define MAX_NMAT         15   // maximum number of different refractive indices (<256)
 #define MAX_N_SH_PARMS   25   // maximum number of shape parameters
 #define MAX_N_BEAM_PARMS 10   // maximum number of beam parameters
+#define MAX_N_LAYERS     2    // maximum number of layers in multi-surface mode
 
 // sizes of filenames and other strings
 /* There is MAX_PATH constant that equals 260 on Windows. However, even this OS allows ways to override this limit. On
diff --git a/src/crosssec.c b/src/crosssec.c
index b496c226..440f3e73 100644
--- a/src/crosssec.c
+++ b/src/crosssec.c
@@ -628,12 +628,12 @@ static void CalcFieldSurf(doublecomplex ebuff[static restrict 3], // where to wr
 	cvInit(sumN);
 	if (above) cvInit(sumF); //additional storage for directly propagated scattering
 
-	/* There is an inherent discontinuity for msub approaching 1 and scattering angle 90 degrees (nF[2]=0). The problem
+	/* There is an inherent discontinuity for sub.m[sub.N-1] approaching 1 and scattering angle 90 degrees (nF[2]=0). The problem
 	 * is that for m=1+-0, but |m-1|>>(nF[2])^2, ki<<kt<<1 => rs=rp=-1
 	 * while for m=1 (exactly) the limit of nF[2]->0 results in kt=ki => rs=rp=0
 	 * Therefore, below is a certain logic, which behaves in an intuitively expected way, for common special cases.
-	 * However, it is still not expected to be continuous for fine-changing parameters (like msub approaching 1).
-	 * In particular, the jump occurs when msub crosses 1+-ROUND_ERR boundary.
+	 * However, it is still not expected to be continuous for fine-changing parameters (like sub.m[sub.N-1] approaching 1).
+	 * In particular, the jump occurs when sub.m[sub.N-1] crosses 1+-ROUND_ERR boundary.
 	 * Still, the discontinuity should apply only to scattering at exactly 90 degrees, but not to, e.g., integral
 	 * quantities, like Csca (if sufficient large number of integration points is chosen).
 	 */
@@ -644,44 +644,34 @@ static void CalcFieldSurf(doublecomplex ebuff[static restrict 3], // where to wr
 		 * a particle at a planar interface," J. Opt. Soc. Am. A 30, 2519-2525 (2013) for theoretical discussion of
 		 * this fact.
 		 */
-		if (fabs(nF[2])<ROUND_ERR && cabs(msub-1)>ROUND_ERR) {
+		if (fabs(nF[2])<ROUND_ERR && cabs(sub.m[sub.N-1]-1)>ROUND_ERR) {
 			cvInit(ebuff);
 			return;
 		}
 		cvBuildRe(nF,nN);
 		nN[2]*=-1;
 		ki=nF[2]; // real
-		if (msubInf) {
-			cs=-1;
-			cp=1;
-		}
-		// since kt is not further needed, we directly calculate cs and cp (equivalent to kt=ki)
-		else if (cabs(msub-1)<ROUND_ERR && cabs(ki)<SQRT_RND_ERR) cs=cp=0;
-		else { // no special treatment here, since other cases, including 90deg-scattering, are taken care above.
-			kt=cSqrtCut(msub*msub - (nN[0]*nN[0]+nN[1]*nN[1]));
-			cs=FresnelRS(ki,kt);
-			cp=FresnelRP(ki,kt,msub);
-		}
-		phSh=imExp(2*WaveNum*hsub*creal(ki)); // assumes real ki
+		SubstrateFresnel(sub,WaveNum,false,nN[0]*nN[0]+nN[1]*nN[1],ki,NULL,&cs,NULL,&cp,NULL);
+		phSh=imExp(2*WaveNum*sub.hP*creal(ki)); // assumes real ki
 	}
 	else { // transmission; here nF[2] is negative
 		// formulae correspond to plane wave incoming from below, but with change ki<->kt
-		if (msubInf) { // no transmission for perfectly reflecting substrate => zero result
+		if (sub.mInf) { // no transmission for perfectly reflecting substrate => zero result
 			cvInit(ebuff);
 			return;
 		}
-		kt=-msub*nF[2];
-		if (cabs(msub-1)<ROUND_ERR && cabs(kt)<SQRT_RND_ERR) ki=kt;
-		else ki=cSqrtCut(1 - msub*msub*(nF[0]*nF[0]+nF[1]*nF[1]));
+		kt=-sub.m[sub.N-1]*nF[2];
+		SubstrateFresnel(sub,WaveNum,true,sub.m[sub.N-1]*sub.m[sub.N-1]*(nF[0]*nF[0]+nF[1]*nF[1]),kt,&cs,NULL,&cp,NULL,
+						 &ki);
 		// here nN may be complex, but normalized to n.n=1
-		nN[0]=msub*nF[0];
-		nN[1]=msub*nF[1];
+		nN[0]=sub.m[sub.N-1]*nF[0];
+		nN[1]=sub.m[sub.N-1]*nF[1];
 		nN[2]=-ki;
-		// these formulae works fine for ki=kt (even very small), and ki=kt=0 is impossible here
-		cs=FresnelTS(kt,ki);
-		cp=FresnelTP(kt,ki,1/msub);
 		// coefficient comes from  k0->k in definition of F(n) (in denominator)
-		phSh=msub*cexp(I*WaveNum*hsub*(ki-kt));
+		double hP_eff = sub.hP;
+		for (int k = 0; k < sub.N - 1; ++k)
+			hP_eff += sub.h[k];
+		phSh=sub.m[sub.N-1]*cexp(I*WaveNum * (sub.hP*ki - hP_eff*kt));
 	}
 #ifndef SPARSE
 	// prepare values of exponents, along each of the coordinates
@@ -807,7 +797,7 @@ double ExtCross(const double * restrict incPol)
 	double sum;
 	size_t i;
 
-	// this can be considered a legacy case, which works only for the simplest plane way centered at the particle 
+	// this can be considered a legacy case, which works only for the simplest plane way centered at the particle
 	if (beamtype==B_PLANE && !surface && !beam_asym) {
 		CalcField (ebuff,prop);
 		sum=crDotProd_Re(ebuff,incPol); // incPol is real, so no conjugate is needed
@@ -1009,13 +999,13 @@ void CalcAlldir(void)
 	Accumulate(E_ad,cmplx_type,2*npoints,&Timing_EFieldADComm);
 	// calculate square of the field
 	for (point=0;point<npoints;point++) E2_alldir[point] = cAbs2(E_ad[2*point]) + cAbs2(E_ad[2*point+1]);
-	/* when below surface we scale E2 by Re(1/msub) in accordance with formula for the Poynting vector (and factor of
+	/* when below surface we scale E2 by Re(1/sub.m[sub.N-1]) in accordance with formula for the Poynting vector (and factor of
 	 * k_sca^2). After that Csca (and g) computed using the standard formula should correctly describe the energy and
-	 * momentum balance for any (even complex) msub (since the scattered wave is homogeneous at far-field), but doesn't
+	 * momentum balance for any (even complex) sub.m[sub.N-1] (since the scattered wave is homogeneous at far-field), but doesn't
 	 * include energy or momentum absorbed (obtained) by the medium.
 	 */
-	if (surface && !msubInf) {
-		double scale=creal(1/msub); // == Re(msub)/|msub|^2
+	if (surface && !sub.mInf) {
+		double scale=creal(1/sub.m[sub.N-1]); // == Re(sub.m[sub.N-1])/|sub.m[sub.N-1]|^2
 		for (i=0;i<theta_int.N;++i) if (TestBelowDeg(theta_int.val[i]))
 			for (j=0,point=phi_int.N*i;j<phi_int.N;j++,point++) E2_alldir[point]*=scale;
 	}
diff --git a/src/interaction.c b/src/interaction.c
index 9b822158..f91d2205 100644
--- a/src/interaction.c
+++ b/src/interaction.c
@@ -73,7 +73,7 @@ int local_Nz_Rm; // number of local layers in Rmatrix, not greater than 2*boxZ-1
 
 // KroneckerDelta[mu,nu] - can serve both as multiplier, and as bool
 static const double dmunu[6] = {1.0, 0.0, 0.0, 1.0, 0.0, 1.0};
-static doublecomplex surfRCn; // reflection coefficient for normal incidence
+static doublecomplex surfRCn, surfBcn, surfAcn; // reflection coefficient for normal incidence
 static bool XlessY; // whether boxX is not larger than boxY (used for SomTable)
 static doublecomplex * restrict somTable; // table of Sommerfeld integrals
 static size_t * restrict somIndex; // array for indexing somTable (in the xy-plane)
@@ -650,13 +650,13 @@ static inline void UnitsGridToCoordShift(const int i,const int j,const int k,dou
 //=====================================================================================================================
 
 static inline void ReflTerm_core(const double kr,const double kr2,const double invr3,const double qmunu[static 6],
-	doublecomplex *expval,doublecomplex result[static 6])
+	doublecomplex *expval,doublecomplex scale_coef,doublecomplex result[static 6])
 // Core routine that calculates the reflection interaction between probe dipole and image of source dipole
 {
 	// this is a modification of InterTerm_core, multiplying by refl. coef. and inverting z-components
 	const double t1=(3-kr2), t2=-3*kr, t3=(kr2-1);
 	*expval=invr3*imExp(kr);
-	doublecomplex scale=surfRCn*(*expval);
+	doublecomplex scale=scale_coef*(*expval);
 
 #define INTERACT_DIAG(ind) { result[ind] = ((t1*qmunu[ind]+t3) + I*(kr+t2*qmunu[ind]))*scale; }
 #define INTERACT_NONDIAG(ind) { result[ind] = (t1+I*t2)*qmunu[ind]*scale; }
@@ -676,7 +676,7 @@ static inline void ReflTerm_core(const double kr,const double kr2,const double i
 }
 //=====================================================================================================================
 
-static inline void ReflTerm_img(double qvec[static 3],doublecomplex result[static 6])
+static inline void ReflTerm_img(double qvec[static 3],const doublecomplex scale_coef,doublecomplex result[static 6])
 /* Reflection term using the image-dipole approximation;
  * qvec is the real distance between probe and image of source, so its z-component is the sum of heights of source and
  * probe points above the surface. result is for produced output
@@ -688,10 +688,54 @@ static inline void ReflTerm_img(double qvec[static 3],doublecomplex result[stati
 	doublecomplex expval; // exp(ikR)/|R|^3
 
 	InterParams(qvec,qmunu,&rr,&invr3,&kr,&kr2);
-	ReflTerm_core(kr,kr2,invr3,qmunu,&expval,result);
+	ReflTerm_core(kr,kr2,invr3,qmunu,&expval,scale_coef,result);
+}
+
+void ReflTerm_img_add_series(const double qvec_in[static 3], doublecomplex result[static restrict 6]) {
+	double qvec[3];
+	doublecomplex term_buffer[6];
+	doublecomplex term_scale_coef = surfBcn;
+	doublecomplex series_buffer[6] = {0,0,0,0,0,0};
+	double z_layer_shift;
+	vCopy(qvec_in, qvec);
+
+	z_layer_shift = sub.h[0] * 2;
+	qvec[2] += z_layer_shift;
+	ReflTerm_img(qvec, term_scale_coef, term_buffer);
+	cdAdd(term_buffer, series_buffer, series_buffer);
+	double error = cdNorm2(series_buffer) + cdNorm2(result);
+	const double error_stop = error * ROUND_ERR * ROUND_ERR;
+	term_scale_coef *= surfAcn;
+
+	while (error > error_stop) {
+		vCopy(qvec_in, qvec);
+		z_layer_shift += sub.h[0] * 2;
+		qvec[2] += z_layer_shift;
+		ReflTerm_img(qvec, term_scale_coef, term_buffer);
+		cdAdd(term_buffer, series_buffer, series_buffer);
+		error = cdNorm2(term_buffer);
+		term_scale_coef *= surfAcn;
+	}
+	cdAdd(series_buffer,result,result);
+}
+
+void ReflTerm_img_int(const int i,const int j,const int k,doublecomplex result[static restrict 6]) {
+	double qvec[3];
+	UnitsGridToCoordShift(i,j,k,qvec);
+	ReflTerm_img(qvec, surfRCn, result);
+	if (sub.N == 2){
+		UnitsGridToCoordShift(i,j,k,qvec);
+		ReflTerm_img_add_series(qvec, result);
+	}
 }
 
-WRAPPERS_REFL(ReflTerm_img)
+void ReflTerm_img_real(const double qvec_in[static restrict 3],doublecomplex result[static restrict 6]) {
+	double qvec[3];
+	vCopy(qvec_in, qvec);
+	ReflTerm_img(qvec, surfRCn, result);
+	if (sub.N == 2)
+		ReflTerm_img_add_series(qvec_in, result);
+}
 
 //=====================================================================================================================
 
@@ -822,7 +866,7 @@ void ReflTerm_som_int(const int i,const int j,const int k,doublecomplex result[s
 	doublecomplex expval; // exp(ikR)/|R|^3
 	UnitsGridToCoordShift(i,j,k,qvec);
 	InterParams(qvec,qmunu,&rr,&invr3,&kr,&kr2);
-	ReflTerm_core(kr,kr2,invr3,qmunu,&expval,result);
+	ReflTerm_core(kr,kr2,invr3,qmunu,&expval,surfRCn,result);
 
 	// second, Sommerfeld integral part
 	// compute table index
@@ -859,7 +903,7 @@ void ReflTerm_som_real(const double qvec[static restrict 3],doublecomplex result
 	double qv[3]; // copy of qvec
 	vCopy(qvec,qv);
 	InterParams(qv,qmunu,&rr,&invr3,&kr,&kr2);
-	ReflTerm_core(kr,kr2,invr3,qmunu,&expval,result);
+	ReflTerm_core(kr,kr2,invr3,qmunu,&expval,surfRCn,result);
 
 	// second, Sommerfeld integral part
 	double x=qvec[0];
@@ -937,7 +981,7 @@ void InitInteraction(void)
 			case GR_IMG:  SET_FUNC_POINTERS(ReflTerm,img); break;
 			case GR_SOM:
 				SET_FUNC_POINTERS(ReflTerm,som);
-				if (!prognosis) som_init(msub*msub);
+				if (!prognosis) som_init(sub.m[sub.N - 1] * sub.m[sub.N - 1]);
 				CalcSomTable();
 				break;
 			/* TO ADD NEW REFLECTION FORMULATION
@@ -951,7 +995,15 @@ void InitInteraction(void)
 			default: LogError(ONE_POS, "Invalid reflection term calculation method: %d",(int)ReflRelation);
 				// no break
 		}
-		surfRCn=msubInf ? -1 : ((1-msub*msub)/(1+msub*msub));
+		if (sub.N == 1)
+			surfRCn=sub.mInf ? -1 : ((1-sub.m[0] * sub.m[0])/(1+sub.m[0] * sub.m[0]));
+		else if (sub.N == 2) {
+			doublecomplex e2 = sub.m[0]*sub.m[0];
+			surfRCn=(1-e2)/(1+e2);
+			doublecomplex surfRCn_mt = sub.mInf ? -1 : (e2 - sub.m[1]*sub.m[1])/(e2 + sub.m[1]*sub.m[1]);
+			surfBcn=4*e2/(1+e2)/(1+e2)*surfRCn_mt;
+			surfAcn=-surfRCn*surfRCn_mt;
+		}
 	}
 
 #ifdef USE_SSE3
diff --git a/src/make_particle.c b/src/make_particle.c
index 2c152d18..84696344 100644
--- a/src/make_particle.c
+++ b/src/make_particle.c
@@ -2427,14 +2427,14 @@ void MakeParticle(void)
 		if (minZco>DipoleCoord[i3+2]) minZco=DipoleCoord[i3+2]; // crude way to find the minimum on the way
 	}
 	/* test that particle is wholly above the substrate; strictly speaking, we test dipole centers to be above the
-	 * substrate - hsub+minZco>0, while the geometric boundary of the particle may still intersect with the substrate.
+	 * substrate - sub.hP+minZco>0, while the geometric boundary of the particle may still intersect with the substrate.
 	 * However, the current test is sufficient to ensure that corresponding routines to calculate reflected Green's
 	 * tensor do not fail. And accuracy of the DDA itself is anyway questionable when some of the dipoles are very close
 	 * to the substrate (whether they cross it or not).
 	 */
-	if (surface && hsub<=-minZco) LogError(ALL_POS,"The particle must be entirely above the substrate. There exist a "
+	if (surface && sub.hP<=-minZco) LogError(ALL_POS,"The particle must be entirely above the substrate. There exist a "
 		"dipole with z="GFORMDEF" (relative to the center), making specified height of the center ("GFORMDEF") too "
-		"small",minZco,hsub);
+		"small",minZco,sub.hP);
 	// save geometry
 	if (save_geom)
 #ifndef SPARSE
@@ -2459,10 +2459,10 @@ void MakeParticle(void)
 	box_origin_unif[1]=-dsY*cY;
 #ifndef SPARSE
 	box_origin_unif[2]=dsZ*(local_z0_unif-cZ);
-	if (surface) ZsumShift=2*(hsub+(local_z0-cZ)*dsZ);
+	if (surface) ZsumShift=2*(sub.hP+(local_z0-cZ)*dsZ);
 #else
 	box_origin_unif[2]=-dsZ*cZ;
-	if (surface) ZsumShift=2*(hsub-cZ*dsZ);
+	if (surface) ZsumShift=2*(sub.hP-cZ*dsZ);
 #	ifdef PARALLEL
 	AllGather(NULL,position_full,int3_type,NULL);
 #	endif
diff --git a/src/param.c b/src/param.c
index f286cc1f..1a0df5e9 100644
--- a/src/param.c
+++ b/src/param.c
@@ -686,12 +686,15 @@ static struct opt_struct options[]={
 #endif
 	{PAR(store_int_field),"","Save internal fields to a file",0,NULL},
 	{PAR(store_scat_grid),"","Calculate Mueller matrix for a grid of scattering angles and save it to a file.",0,NULL},
-	{PAR(surf),"<h> {<mre> <mim>|inf}","Specifies that scatterer is located above the plane surface, parallel to the "
-		"xy-plane. <h> specifies the height of particle center above the surface (along the z-axis, in um). Particle "
-		"must be entirely above the substrate. Following argument(s) specify the refractive index of the substrate "
-		"(below the surface), assuming that the vacuum is above the surface. It is done either by two values (real and "
-		"imaginary parts of the complex value) or as effectively infinite 'inf' which corresponds to perfectly"
-		"reflective surface. The latter implies certain simplifications during calculations.",UNDEF,NULL},
+	{PAR(surf),"<h> {<m1Re> <m1Im> [<h1> <m2Re> <m2Im> ... <h_n-1> {<mnRe> <mnIm>|inf}]|inf}",
+		"Specifies that scatterer is located (entirely) above the n-layered plane substrate (n>=1), parallel to the "
+		"xy-plane. <h> is the height of particle center above the substrate (its top interface) along the z-axis "
+		"(in um), must be larger than corresponding distance from the center to the lowermost dipole. Additional "
+		"arguments <h1> ... <h_n-1> specify each layer thickness (from top to bottom, also in um) - the last (n-th) "
+		"layer is always semi-infinite (but can be a vacuum one). Refractive indices (real and imaginary parts: "
+		"<m1Re>, <m1Im>, ... <mnRe>, <mnIm>) are specified for each layer, assuming the vacuum is above the substrate. "
+		"The last refractive index can also be set with single 'inf' value, corresponding to a perfect electric "
+		"conductor as a bottom layer. The latter implies certain simplifications during calculations.",UNDEF,NULL},
 	{PAR(sym),"{auto|no|enf}","Automatically determine particle symmetries ('auto'), do not take them into account "
 		"('no'), or enforce them ('enf').\n"
 		"Default: auto",1,NULL},
@@ -1613,18 +1616,42 @@ PARSE_FUNC(store_scat_grid)
 PARSE_FUNC(surf)
 {
 	double mre,mim;
-	if (Narg!=2 && Narg!=3) NargError(Narg,"2 or 3");
-	ScanDoubleError(argv[1],&hsub);
-	TestPositive(hsub,"height above surface");
-	if (strcmp(argv[2],"inf")==0) {
-		if (Narg>2) PrintErrorHelp("Additional arguments are not allowed for '-surf ... inf'");
-		msubInf=true;
+	int i,arg_pos;
+
+	if (Narg%3 == 1) NargError(Narg,"3n or 3n-1, if the last argument is 'inf',");
+	sub.N=(Narg+1)/3;
+	if (sub.N>MAX_N_LAYERS) PrintErrorHelp("Too many layers in the substrate (%d), maximum %d are supported. You may "
+		"increase parameter MAX_N_LAYERS in const.h and recompile.",sub.N,MAX_N_LAYERS);
+	/* In principle, everywhere in the command line ADDA accepts inf instead of float arguments (that's how sscanf
+	 * works). ADDA may even work fine in some cases, but that is definitely not supported and is left for user's
+	 * responsibility. Still, we do not explicitly test for or forbid it.
+	 * By contrast, here we test all arguments against inf, since it has legitimate usage here, but accidentally
+	 * removing or adding one variable before it may lead to nonsense simulation (and wasted computational resources).
+	 * So we do our best to produce error message as early as possible. Still, something like 'inf1' will still cause
+	 * strange bevavior.
+	 */
+	sub.mInf=false;
+	for (i=1;i<=Narg;i++) if (strcmp(argv[i],"inf")==0) {
+		if (i<Narg) PrintErrorHelp("Specification of 'inf' as the layer refractive index must be the last argument to "
+			"'-surf ...'");
+		if (Narg%3 != 2) NargError(Narg,"if the last argument is 'inf', 3n-1");
+		sub.mInf=true;
 	}
-	else {
-		if (Narg!=3) NargError(Narg,"total 3 or second argument 'inf'");
-		ScanDoubleError(argv[2],&mre);
-		ScanDoubleError(argv[3],&mim);
-		msub = mre + I*mim;
+	for (i=0;i<sub.N;i++) {
+		arg_pos=3*i+1;
+		if (i==0) { // first h is named differently
+			ScanDoubleError(argv[arg_pos],&sub.hP);
+			TestPositive(sub.hP,"height above substrate");
+		}
+		else {
+			ScanDoubleError(argv[arg_pos],sub.h+i-1);
+			TestPositive(sub.h[i-1],"layer thickness");
+		}
+		if (!(i==sub.N-1 && sub.mInf)) { // exception for infinite last refractive index (processed above)
+			ScanDoubleError(argv[arg_pos+1],&mre);
+			ScanDoubleError(argv[arg_pos+2],&mim);
+			sub.m[i] = mre + I*mim;
+		}
 	}
 	surface = true;
 	symZ=false;
@@ -2029,7 +2056,7 @@ void InitVariables(void)
 	InitField=IF_AUTO;
 	recalc_resid=false;
 	surface=false;
-	msubInf=false;
+	sub.mInf=false;
 	ReflRelation=(enum refl)UNDEF;
 	// sometimes the following two are left uninitialized
 	beam_fnameX=NULL;
@@ -2229,8 +2256,8 @@ void VariablesInterconnect(void)
 		if (orient_used) PrintError("Currently '-orient' and '-surf' can not be used together");
 		if (calc_mat_force) PrintError("Currently calculation of radiation forces is incompatible with '-surf'");
 		if (InitField==IF_WKB) PrintError("'-init_field wkb' and '-surf' can not be used together");
-		if (ReflRelation==(enum refl)UNDEF) ReflRelation = msubInf ? GR_IMG : GR_SOM;
-		else if (msubInf && ReflRelation!=GR_IMG) PrintError("For perfectly reflecting surface interaction is always "
+		if (ReflRelation==(enum refl)UNDEF) ReflRelation = sub.mInf ? GR_IMG : GR_SOM;
+		else if (sub.mInf && ReflRelation!=GR_IMG) PrintError("For perfectly reflecting surface interaction is always "
 			"computed through an image dipole. So this case is incompatible with other options to '-int_surf ...'");
 		/* TO ADD NEW REFLECTION FORMULATION
 		 * Take a look at the above logic, and revise if the new formulation is not fully consistent with it
@@ -2477,9 +2504,9 @@ void PrintInfo(void)
 			}
 		}
 		if (surface) {
-			if (msubInf) fprintf(logfile,"Particle is placed near the perfectly reflecting substrate\n");
-			else fprintf(logfile,"Particle is placed near the substrate with refractive index "CFORM"\n",REIM(msub));
-			fprintf(logfile,"  height of the particle center: "GFORMDEF"\n",hsub);
+			if (sub.mInf) fprintf(logfile,"Particle is placed near the perfectly reflecting substrate\n");
+			else fprintf(logfile,"Particle is placed near the substrate with refractive index "CFORM"\n",REIM(sub.m[sub.N-1]));
+			fprintf(logfile,"  height of the particle center: "GFORMDEF"\n",sub.hP);
 		}
 		if (rectDip) fprintf(logfile,"Dipole size: "GFORMDEF"x"GFORMDEF"x"GFORMDEF" (relative ratios "GFORMDEF":"
 			GFORMDEF":"GFORMDEF")\n",dsX,dsY,dsZ,rectScaleX,rectScaleY,rectScaleZ);
@@ -2505,7 +2532,7 @@ void PrintInfo(void)
 		if (surface && beamtype!=B_DIPOLE) { // include surface-specific vectors
 			fprintf(logfile,"Reflected propagation vector: "GFORMDEF3V"\n",COMP3V(prIncRefl));
 			fprintf(logfile,"Transmitted propagation vector: ");
-			if (msubInf) fprintf (logfile,"none\n");
+			if (sub.mInf) fprintf (logfile,"none\n");
 			else {
 				fprintf(logfile,GFORMDEF3V,COMP3V(prIncTran));
 				if (prIncTran[2]==0) fprintf(logfile," (evanescent)\n");
diff --git a/src/types.h b/src/types.h
index 43125e14..452336ee 100644
--- a/src/types.h
+++ b/src/types.h
@@ -55,4 +55,13 @@ typedef struct	        // integration parameters
 	angle_set phi;      // values of phi
 } scat_grid_angles;
 
+struct Substrate
+{
+	bool mInf;                     // whether the reflactive index of the last layer is inf
+	int N;                         // total number of layers
+	double hP;                     // distance between the substrate (top intefrace) and the particle center
+	doublecomplex m[MAX_N_LAYERS]; // refractive indices of the layers
+	double h[MAX_N_LAYERS-1];      // thicknesses of the layers (the last one is always semi-infinite)
+};
+
 #endif // __types_h
diff --git a/src/vars.c b/src/vars.c
index 3ab83bf2..d1349758 100644
--- a/src/vars.c
+++ b/src/vars.c
@@ -70,7 +70,7 @@ double box_origin_unif[3];    /* coordinates of the center of the first dipole i
                                  uniform distribution of non-void dipoles among all processors) */
 
 // 3D vectors (in laboratory reference frame)
-double prop_0[3];             // incident direction 
+double prop_0[3];             // incident direction
 double beam_center_0[3];      // coordinates of the beam center
 
 // file info
@@ -129,10 +129,8 @@ TIME_TYPE Timing_EField,      // time for calculating scattered fields
 // related to a nearby surface
 bool surface;           // whether nearby surface is present
 enum refl ReflRelation; // method to calculate reflected Green's tensor
-doublecomplex msub;     // complex refractive index of the substrate
-double inc_scale;       // scale to account for irradiance of the incident beam - 1/Re(msub)
-bool msubInf;           // whether msub is infinite (perfectly reflecting surface)
-double hsub;            // height of particle center above surface
+double inc_scale;       // scale to account for irradiance of the incident beam - 1/Re(sub.m[sub.N-1])
+struct Substrate sub;   //structure with surface related parameters
 
 #ifndef SPARSE // These variables are exclusive to the FFT mode
 
diff --git a/src/vars.h b/src/vars.h
index a7ebd51e..b350394e 100644
--- a/src/vars.h
+++ b/src/vars.h
@@ -81,10 +81,10 @@ extern size_t local_Ndip,local_nvoid_Ndip,local_nRows,local_nvoid_d0,local_nvoid
 extern TIME_TYPE Timing_EField,Timing_FileIO,Timing_Integration,tstart_main;
 
 // related to a nearby surface
-extern bool surface,msubInf;
+extern struct Substrate sub;
+extern bool surface;
 extern enum refl ReflRelation;
-extern doublecomplex msub;
-extern double inc_scale,hsub;
+extern double inc_scale;
 
 #ifndef SPARSE // These variables are exclusive to the FFT mode
 
diff --git a/tests/equiv/surface.sh b/tests/equiv/surface.sh
new file mode 100644
index 00000000..f9648880
--- /dev/null
+++ b/tests/equiv/surface.sh
@@ -0,0 +1,47 @@
+#!/bin/bash
+# Copyright (C) ADDA contributors
+# GNU General Public License version 3
+
+# Variables
+ADDA='../../src/seq/adda'
+ADDAREF='../../win64/adda'
+
+atol=14
+rtol=8
+
+function numdiff {
+  diff $1 $2 | awk -f ../2exec/diff_numeric.awk -v abs_tol=1e-$atol -v rel_tol=1e-$rtol -  >&2
+}
+
+function mydiff {
+  tmp=""
+  if ! numdiff $1 $2; then
+    echo -e "DIFF above is between files '$1' and '$2'" >&2
+    echo "TEST $3 FAILED."
+  fi
+
+}
+
+$ADDAREF -surf 5 1.3 0.1 -store_beam -int_surf img -dir dirref >outref
+$ADDA -surf 5 1.3 0.1 -store_beam -int_surf img -dir dir1 >out1
+$ADDA -surf 4 1 0 1 1.3 0.1 -store_beam -int_surf img -dir dir2 >out2
+$ADDA -surf 5 1.3 0.1 1 1.3 0.1 -store_beam -int_surf img -dir dir3 >out3
+mydiff dirref/IncBeam-Y dir1/IncBeam-Y 1
+mydiff dir1/IncBeam-Y dir2/IncBeam-Y 2
+mydiff dir1/IncBeam-Y dir3/IncBeam-Y 3
+mydiff dirref/CrossSec-Y dir1/CrossSec-Y 4
+mydiff dir1/CrossSec-Y dir2/CrossSec-Y 5
+mydiff dir1/CrossSec-Y dir3/CrossSec-Y 6
+rm -rf dirref outref dir1 dir2 dir3 out1 out2 out3
+
+$ADDAREF -surf 5 1.3 0.1 -store_beam -prop 0 0 -1 -int_surf img -dir dirref >outref
+$ADDA -surf 5 1.3 0.1 -store_beam -prop 0 0 -1 -int_surf img -dir dir1 >out1
+$ADDA -surf 4 1 0 1 1.3 0.1 -store_beam -prop 0 0 -1 -int_surf img -dir dir2 >out2
+$ADDA -surf 5 1.3 0.1 1 1.3 0.1 -store_beam -prop 0 0 -1 -int_surf img -dir dir3 >out3
+mydiff dirref/IncBeam-Y dir1/IncBeam-Y 11
+mydiff dir1/IncBeam-Y dir2/IncBeam-Y 12
+mydiff dir1/IncBeam-Y dir3/IncBeam-Y 13
+mydiff dirref/CrossSec-Y dir1/CrossSec-Y 14
+mydiff dir1/CrossSec-Y dir2/CrossSec-Y 15
+mydiff dir1/CrossSec-Y dir3/CrossSec-Y 16
+rm -rf dirref outref dir1 dir2 dir3 out1 out2 out3
\ No newline at end of file
diff --git a/tests/pytests/__init__.py b/tests/pytests/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/pytests/wrappers/__init__.py b/tests/pytests/wrappers/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/pytests/wrappers/setup.py b/tests/pytests/wrappers/setup.py
new file mode 100644
index 00000000..a25402aa
--- /dev/null
+++ b/tests/pytests/wrappers/setup.py
@@ -0,0 +1,15 @@
+from distutils.core import setup, Extension
+import os
+
+
+source_dir = os.path.abspath('./../../../src')
+
+if __name__ == "__main__":
+    setup(
+        name="ADDAPyWrappers",
+        version="1.0.0",
+        description="Python interface for some ADDA functions",
+        ext_modules=[Extension(
+            "ADDAWrappers", ["wrappers.c", os.path.join(source_dir, "cmplx.c")],
+        )]
+    )
diff --git a/tests/pytests/wrappers/wrappers.c b/tests/pytests/wrappers/wrappers.c
new file mode 100644
index 00000000..0e92a582
--- /dev/null
+++ b/tests/pytests/wrappers/wrappers.c
@@ -0,0 +1,113 @@
+#define PY_SSIZE_T_CLEAN
+#include <Python.h>
+#ifndef OVERRIDE_STDC_TEST
+#define OVERRIDE_STDC_TEST
+#endif
+#include "../../../src/cmplx.h"
+
+#define FRESNEL_DECL(POLAR_TYPE) {            \
+	"Fresnel"#POLAR_TYPE,                     \
+	Py_Fresnel##POLAR_TYPE,                   \
+	METH_VARARGS,                             \
+	"Python interface for Fresnel"#POLAR_TYPE \
+	}
+
+#define FRESNEL_ARGNUM_S 2
+#define FRESNEL_ARGNUM_P 3
+#define FRESNEL_ARRAY_SEQ_S(a) a[0], a[1]
+#define FRESNEL_ARRAY_SEQ_P(a) a[0], a[1], a[2]
+#define FRESNEL_FORMAT_S "DD"
+#define FRESNEL_FORMAT_P "DDD"
+
+#define MAKE_FRESNEL(T_OR_R, S_OR_P)                                                     \
+    static PyObject * Py_Fresnel##T_OR_R##S_OR_P(PyObject * self, PyObject * args) {     \
+		doublecomplex c_arg[FRESNEL_ARGNUM_##S_OR_P];                                    \
+		Py_complex py_arg[FRESNEL_ARGNUM_##S_OR_P];                                      \
+		if (!PyArg_ParseTuple(args, FRESNEL_FORMAT_##S_OR_P,                             \
+			                  FRESNEL_ARRAY_SEQ_##S_OR_P(&py_arg))                       \
+		)                                                                                \
+			return NULL;                                                                 \
+		for (int i = 0; i < FRESNEL_ARGNUM_##S_OR_P; ++i)                                \
+			c_arg[i] = py_arg[i].real + I * py_arg[i].imag;                              \
+        doublecomplex res = Fresnel##T_OR_R##S_OR_P(FRESNEL_ARRAY_SEQ_##S_OR_P(c_arg));  \
+		Py_complex py_res = {creal(res), cimag(res)};                                    \
+		return PyComplex_FromCComplex(py_res);                                           \
+	}
+
+MAKE_FRESNEL(T, S);
+MAKE_FRESNEL(R, S);
+MAKE_FRESNEL(T, P);
+MAKE_FRESNEL(R, P);
+
+void PrintError(const char * restrict fmt, ... ) {
+	printf(fmt);
+};
+
+#TODO THIS WRAPPER IS INCOMPLETE
+static PyObject * Py_SubstrateFresnel(PyObject * self, PyObject * args) {
+	PyObject  py_sub;
+	struct Substrate c_sub;
+	Py_complex py_wave_num;
+	doublecomplex c_wave_num;
+	bool is_z_positive;
+	Py_complex py_sqr_long_k;
+	doublecomplex c_sqr_long_k;
+	Py_complex py_ki;
+	doublecomplex c_ki;
+	PyObject py_coef_set;
+	doublecomplex c_coef[4];
+	doublecomplex * c_coef_ptr[4];
+	for (int i = 0; i < 4; ++i)
+		c_coef_ptr[i] = NULL;
+
+	if (!PyArg_ParseTuple(args, "ODpDDO", &py_sub, &py_wave_num, &is_z_positive, &py_sqr_long_k, &py_ki, &py_coef_set))
+		return NULL;
+	c_wave_num = py_wave_num.real + I * py_wave_num.imag;
+	c_sqr_long_k = py_sqr_long_k.real + I * py_sqr_long_k.imag;
+	c_ki = py_ki.real + I * py_ki.imag;
+	PyObject * tmp_obj;
+	tmp_obj = PyDict_GetItemString(&py_sub, "mInf");
+	if (!PyArg_ParseTuple(tmp_obj, "p", &c_sub.mInf))
+		return NULL;
+	tmp_obj = PyDict_GetItemString(&py_sub, "N");
+	if (!PyArg_ParseTuple(tmp_obj, "i", &c_sub.N))
+		return NULL;
+	tmp_obj = PyDict_GetItemString(&py_sub, "hP");
+	if (!PyArg_ParseTuple(tmp_obj, "d", &c_sub.hP))
+		return NULL;
+	Py_complex m[2];
+	tmp_obj = PyDict_GetItemString(&py_sub, "m");
+	if (!PyArg_ParseTuple(tmp_obj, "(DD)", m))
+		return NULL;
+	c_sub.m[0] = m[0].real + I * m[0].imag;
+	c_sub.m[1] = m[0].real + I * m[1].imag;
+	tmp_obj = PyDict_GetItemString(&py_sub, "h");
+	if (!PyArg_ParseTuple(tmp_obj, "(dd)", c_sub.h))
+		return NULL;
+	SubstrateFresnel(c_sub, c_wave_num, is_z_positive, c_sqr_long_k, c_ki,
+					 &c_coef[0], &c_coef[1], &c_coef[2], &c_coef[3], NULL);
+	Py_complex py_res = {creal(c_coef[0]), cimag(c_coef[0])};
+		return PyComplex_FromCComplex(py_res);
+}
+
+static PyMethodDef ADDAWrappersMethods[] = {
+		FRESNEL_DECL(TS),
+		FRESNEL_DECL(RS),
+		FRESNEL_DECL(TP),
+		FRESNEL_DECL(RP),
+		{"SubstrateFresnel", Py_SubstrateFresnel, METH_VARARGS, "Python interface for SubstrateFresnel"},
+		{NULL, NULL, 0, NULL}
+};
+
+static struct PyModuleDef ADDAWrappersModule = {
+		PyModuleDef_HEAD_INIT,
+		"Fresnel",
+		"Python interface for fresnel coefficients",
+		-1,
+		ADDAWrappersMethods
+};
+
+PyMODINIT_FUNC PyInit_ADDAWrappers(void) {
+	return PyModule_Create(&ADDAWrappersModule);
+}
+
diff --git a/tests/unittests/fresnel/CMakeLists.txt b/tests/unittests/fresnel/CMakeLists.txt
new file mode 100644
index 00000000..2e068008
--- /dev/null
+++ b/tests/unittests/fresnel/CMakeLists.txt
@@ -0,0 +1,13 @@
+project(TestFresnel)
+cmake_minimum_required(VERSION 3.0)
+set (CMAKE_C_STANDARD 99)
+set (ADDA_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/../../../src/)
+set (SOURCE_FILES TestFresnel.c ${ADDA_SOURCE_DIR}/cmplx.c)
+find_library(CHECK_LIBRARY REQUIRED
+        NAMES Check check)
+set(CHECK_LIBRARIES "${CHECK_LIBRARY}")
+include_directories(${CHECK_INCLUDE_DIRS})
+MARK_AS_ADVANCED( CHECK_INCLUDE_DIR CHECK_LIBRARIES )
+
+add_executable(test_fresnel ${SOURCE_FILES})
+target_link_libraries(test_fresnel ${CHECK_LIBRARIES})
\ No newline at end of file
diff --git a/tests/unittests/fresnel/TestFresnel.c b/tests/unittests/fresnel/TestFresnel.c
new file mode 100644
index 00000000..2b3ce1e9
--- /dev/null
+++ b/tests/unittests/fresnel/TestFresnel.c
@@ -0,0 +1,183 @@
+#include <check.h>
+#include <stdio.h>
+#include "../../../src/cmplx.h"
+
+#define ck_assert_doublecomplex_eq_tol(a, b, tol)                                         \
+	do {                                                                                  \
+		ck_assert_double_eq_tol(creal((doublecomplex) a), creal((doublecomplex) b), tol); \
+		ck_assert_double_eq_tol(cimag((doublecomplex) a), cimag((doublecomplex) b), tol); \
+	} while (false)
+
+void PrintError(const char * restrict fmt, ... ) {
+	printf(fmt);
+};
+
+START_TEST(test_FresnelS) {
+	double tol = 1e-16;
+	doublecomplex t = FresnelTS(1.3 + I * 0.1, 1.3 + I * 0.1);
+	ck_assert_doublecomplex_eq_tol(t, 1, tol);
+	doublecomplex r = FresnelRS(1.3 + I * 0.1, 1.3 + I * 0.1);
+	ck_assert_doublecomplex_eq_tol(r, 0, tol);
+
+}
+END_TEST
+
+START_TEST(test_FresnelP) {
+	double tol = 1e-16;
+	doublecomplex t = FresnelTP(1.3 + I * 0.1, 1.3 + I * 0.1, 1);
+	ck_assert_doublecomplex_eq_tol(t, 1, tol);
+	doublecomplex r = FresnelRP(1.3 + I * 0.1, 1.3 + I * 0.1, 1);
+	ck_assert_doublecomplex_eq_tol(t, 1, tol);
+
+}
+END_TEST
+
+void check_all_cases_SubstrateFresnel (const struct Substrate sub, const double wave_num, const bool z_positive,
+		                               const doublecomplex sqr_long_k, const doublecomplex ki, 
+									   const doublecomplex ts_ref, const doublecomplex rs_ref, 
+		                               const doublecomplex tp_ref, const doublecomplex rp_ref, doublecomplex kt_ref,
+		                               const double tol) {
+#define ARG_SEQ_FRESNEL sub, wave_num, z_positive, sqr_long_k, ki
+#define CHECK_FRESNEL_HELPER(c) ck_assert_doublecomplex_eq_tol(c, c##_ref, tol)
+#define REINIT_COEF ts = 2, rs = 2, tp = 2, rp = 2, kt = 2
+	
+	doublecomplex ts, rs, tp, rp, kt;
+	REINIT_COEF;
+	
+	SubstrateFresnel(ARG_SEQ_FRESNEL, &ts, NULL, NULL, NULL, &kt);
+	CHECK_FRESNEL_HELPER(ts);
+	CHECK_FRESNEL_HELPER(kt);
+	REINIT_COEF;
+	
+	SubstrateFresnel(ARG_SEQ_FRESNEL, NULL, &rs, NULL, NULL, NULL);
+	CHECK_FRESNEL_HELPER(rs);
+	REINIT_COEF;
+	SubstrateFresnel(ARG_SEQ_FRESNEL, NULL, NULL, &tp, NULL, NULL);
+	CHECK_FRESNEL_HELPER(tp);
+	REINIT_COEF;
+
+	SubstrateFresnel(ARG_SEQ_FRESNEL, NULL, NULL, NULL, &rp, NULL);
+	CHECK_FRESNEL_HELPER(rp);
+	REINIT_COEF;
+
+	SubstrateFresnel(ARG_SEQ_FRESNEL, &ts, &rs, NULL, NULL, NULL);
+	CHECK_FRESNEL_HELPER(ts);
+	CHECK_FRESNEL_HELPER(rs);
+	REINIT_COEF;
+
+	SubstrateFresnel(ARG_SEQ_FRESNEL, NULL, NULL, &tp, &rp, NULL);
+	CHECK_FRESNEL_HELPER(tp);
+	CHECK_FRESNEL_HELPER(rp);
+	REINIT_COEF;
+
+	SubstrateFresnel(ARG_SEQ_FRESNEL, &ts, &rs, &tp, &rp, &kt);
+	CHECK_FRESNEL_HELPER(ts);
+	CHECK_FRESNEL_HELPER(rs);
+	CHECK_FRESNEL_HELPER(tp);
+	CHECK_FRESNEL_HELPER(rp);
+	CHECK_FRESNEL_HELPER(kt);
+	REINIT_COEF;
+	
+	
+#undef ARG_SEQ_FRESNEL
+#undef CHECK_FRESNEL_HELPER
+#undef REINIT_COEF
+}
+
+START_TEST(test_SubstrateFresnel_1) {
+	struct Substrate sub;
+	sub.mInf = false;
+	sub.m[0] = 1.3;
+	sub.m[1] = 1.3;
+	sub.h[0] = 2;
+	sub.N = 2;
+	double wave_num = 1;
+	doublecomplex ki = 1;
+	doublecomplex kt_ref = CalculateKt(ki, 1, sub.m[1], 0);
+	doublecomplex eL = cexp(I * wave_num * sub.h[0] * kt_ref);
+	doublecomplex ts_ref = FresnelTS(ki, kt_ref) * eL;
+	doublecomplex tp_ref = FresnelTP(ki, kt_ref, 1.3) * eL;
+	doublecomplex rs_ref = FresnelRS(ki, kt_ref);
+	doublecomplex rp_ref = FresnelRP(ki, kt_ref, 1.3);
+	check_all_cases_SubstrateFresnel(sub, wave_num, false, 0, ki,
+                                     ts_ref, rs_ref, tp_ref, rp_ref, kt_ref, 1e-14);
+
+}
+END_TEST
+
+START_TEST(test_SubstrateFresnel_2) {
+    struct Substrate sub;
+    sub.mInf = false;
+    sub.m[0] = 1;
+    sub.m[1] = 1.3;
+    sub.h[0] = 2;
+    sub.N = 2;
+    double wave_num = 1;
+    doublecomplex ki = 1;
+    doublecomplex eL = cexp(I * wave_num * sub.h[0] * ki);
+    doublecomplex kt_ref = CalculateKt(ki, 1, sub.m[1], 0);
+    doublecomplex ts_ref = FresnelTS(ki, kt_ref) * eL;
+    doublecomplex tp_ref = FresnelTP(ki, kt_ref, 1.3) * eL;
+    doublecomplex rs_ref = FresnelRS(ki, kt_ref) * eL * eL;
+    doublecomplex rp_ref = FresnelRP(ki, kt_ref, 1.3) * eL * eL;
+    check_all_cases_SubstrateFresnel(sub, wave_num, false, 0, ki,
+                                     ts_ref, rs_ref, tp_ref, rp_ref, kt_ref, 1e-14);
+
+}
+END_TEST
+
+START_TEST(test_SubstrateFresnel_3) {
+    struct Substrate sub;
+    sub.mInf = false;
+    sub.m[0] = 1.3;
+    sub.m[1] = 1.3;
+    sub.h[0] = 2;
+    sub.N = 2;
+    double wave_num = 1;
+    doublecomplex ki = 1.3;
+    doublecomplex eL = cexp(I * wave_num * sub.h[0] * ki);
+    doublecomplex kt_ref = CalculateKt(ki, 1.3, 1, 0);
+    doublecomplex ts_ref = FresnelTS(ki, kt_ref) * eL;
+    doublecomplex tp_ref = FresnelTP(ki, kt_ref, 1/1.3) * eL;
+    doublecomplex rs_ref = FresnelRS(ki, kt_ref) * eL * eL;
+    doublecomplex rp_ref = FresnelRP(ki, kt_ref, 1/1.3) * eL * eL;
+    check_all_cases_SubstrateFresnel(sub, 1, true, 0, ki,
+                                     ts_ref, rs_ref, tp_ref, rp_ref, kt_ref, 1e-14);
+
+}
+END_TEST
+
+Suite * fresnel_suite(void)
+{
+	Suite *s;
+	TCase *tc_Fresnel;
+
+	s = suite_create("Fresnel");
+
+	tc_Fresnel = tcase_create("Fresnel");
+
+
+	tcase_add_test(tc_Fresnel, test_FresnelS);
+	tcase_add_test(tc_Fresnel, test_FresnelP);
+	tcase_add_test(tc_Fresnel, test_SubstrateFresnel_1);
+    tcase_add_test(tc_Fresnel, test_SubstrateFresnel_2);
+    tcase_add_test(tc_Fresnel, test_SubstrateFresnel_3);
+	suite_add_tcase(s, tc_Fresnel);
+
+	return s;
+}
+
+int main(void)
+{
+	int number_failed;
+	Suite *s;
+	SRunner *sr;
+
+	s = fresnel_suite();
+	sr = srunner_create(s);
+
+	srunner_run_all(sr, CK_NORMAL);
+	number_failed = srunner_ntests_failed(sr);
+	srunner_free(sr);
+	return (number_failed == 0) ? 0 : 1;
+}
\ No newline at end of file