[Multi-surface mode]

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;
 	}
 }

src/GenerateB.c

@@ -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,35 @@ 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]));
+					if (cabs(sub.m[sub.N-1]-1)<ROUND_ERR && cabs(ki)<SQRT_RND_ERR) kt=ki;
+					else kt=cSqrtCut(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]));
+						if (cabs(sub.m[sub.N-1]-1)<ROUND_ERR && cabs(ki)<SQRT_RND_ERR) kt=ki;
+						else kt=cSqrtCut(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 of wave transmitted into substrate
 						ktVec[0]=prop_0[0];
 						ktVec[1]=prop_0[1];
@@ -138,15 +138,15 @@ 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);
 				}
 			}
 			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,22 +354,22 @@ 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);
 					}
 					else { // p-polarized
 						crCrossProd(ey,ktVec,eIncTran);
-						tc=FresnelTP(ki,kt,1/msub);
+						tc=FresnelTP(ki,kt,1/sub.m[sub.N-1]);
 					}
 					// 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);
@@ -390,12 +390,12 @@ 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;
+						if (sub.mInf) rcp1=0;
 						else rcp1=FresnelTS(ki,kt);
 					}
 					else { // p-polarized
-						if (msubInf) rcp1=2;
-						else rcp1=msub*FresnelTP(ki,kt,msub);
+						if (sub.mInf) rcp1=2;
+						else rcp1=sub.m[sub.N-1]*FresnelTP(ki,kt,sub.m[sub.N-1]);
 					}
 					// main part
 					for (i=0;i<local_nvoid_Ndip;i++) {
@@ -440,7 +440,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 +457,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 +635,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

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)

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 10       // 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

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,44 @@ 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) {
+		if (sub.mInf) {
 			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;
+		  // since kt is not further needed, we directly calculate cs and cp (equivalent to kt=ki)
+		else if (cabs(sub.m[sub.N-1]-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]));
+			kt=cSqrtCut(sub.m[sub.N-1]*sub.m[sub.N-1] - (nN[0]*nN[0]+nN[1]*nN[1]));
 			cs=FresnelRS(ki,kt);
-			cp=FresnelRP(ki,kt,msub);
+			cp=FresnelRP(ki,kt,sub.m[sub.N-1]);
 		}
-		phSh=imExp(2*WaveNum*hsub*creal(ki)); // assumes real ki
+		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];
+		if (cabs(sub.m[sub.N-1]-1)<ROUND_ERR && cabs(kt)<SQRT_RND_ERR) ki=kt;
+		else ki=cSqrtCut(1 - sub.m[sub.N-1]*sub.m[sub.N-1]*(nF[0]*nF[0]+nF[1]*nF[1]));
 		// 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);
+		cp=FresnelTP(kt,ki,1/sub.m[sub.N-1]);
 		// coefficient comes from  k0->k in definition of F(n) (in denominator)
-		phSh=msub*cexp(I*WaveNum*hsub*(ki-kt));
+		phSh=sub.m[sub.N-1]*cexp(I*WaveNum*sub.hP*(ki-kt));
 	}
 #ifndef SPARSE
 	// prepare values of exponents, along each of the coordinates
@@ -807,7 +807,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 +1009,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;
 	}

src/interaction.c

@@ -937,7 +937,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 +951,7 @@ 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));
+		surfRCn=sub.mInf ? -1 : ((1-sub.m[sub.N - 1] * sub.m[sub.N - 1])/(1+sub.m[sub.N - 1] * sub.m[sub.N - 1]));
 	}
 
 #ifdef USE_SSE3

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