First step towards DDA with particles above surface. Currently, mostl…

…y infrastructure parts have been handled. The only limitation of the scattering problem is that the particle is located fully above the surface and the upper medium is vacuum (the lower has any complex refractive index). Also, -orient, non-plane waves, -init_field wkb, and radiation forces are currently incompatible with surface.

Implemented:
- command line options '-surf ...', '-int_surf ...' and -scat_plane. The latter defines that the scattering should be calculate in the scattering plane, which passes through prop and ez (in laboratory frame). Two difference from the default behavior (-yz): for default prop, it is the xz-plane, and for non-default prop, the angle is counted from ez, not from prop.
- new definitions of scattered angle as -scat_plane for default calculation, and with respect to the laboratory (surface) reference frame for scat_grid and alldir. Such new definitions are now used only for surface, but may become universal ones in the future (issue 170). 
- incident beam due to surface both for propagation from above and from below surface, including evanescent fields and other complex waves.

Still to do:
- Calculation of scattered fields in the presence of surface (simple reflection and transmission).
- Actually plugging in the surface-interaction Green's tensor, both rigorous (Sommerfeld integrals) and approximate (image dipole) ones
- FFT acceleration of surface mode

Also done:
- A few macros were implemented to specify all components of a real or complex vector (or both real and imaginary parts of the complex number) at once in printf-family arguments.
- Special test mode was added to tests/2exec to test new version with trivial surface (of refractive index 1) against the standard implementation. Also added a combination of -orient -prop, and -scat_matr both to all tests. Improved handling of diffs in IncBeam files.
- Several functions have been added to cmplx.h. In particular, s few memcpy calls were replaced by vCopy.
- Handling of scattering directions based on angles theta and phi (alldir and scat_grid) was localized to a separate function SetScatPlane in crosssec.c. Currently, it contains rather complicated logic to handle differences with surface and all special cases.

The following tests were performed:
- In the default mode no changes with the previous version - tested by tests/2exec (including the new surf_standard test mode, where applicable)
- Use of -scat_plane is equivalent to -store_scat_grid with (phi: type=values, values=90)
- The above equivalence also holds for non-trivial -prop with corresponding shift of theta range. For example, -prop 1 0 1 can be compensated by (theta: min=-45, max=315, and double value of N) in scat_params.dat
- When using non-trivial prop with alldir and trivial -surf (like '-surf 1 0 10'), Csca is the same, while g and Csca.g are rotated from beam RF to laboratory RF.
- When using the non-trivial surface the calculation of incident field was tested (both by -store_beam and by comparing calculated Cext, which is supposed to scale accordingly, since other effects of surface are not yet implemented).

src/GenerateB.c

@@ -47,12 +47,22 @@ extern opt_index opt_beam;
 
 // used in crosssec.c
 double beam_center_0[3]; // position of the beam center in laboratory reference frame
+/* complex wave amplitudes of secondary waves (with phase relative to particle center);
+ * they satisfy (e,e)=e_x^2+e_y^2+e_z^2=1, which doesn't necessarily means that (e,e*)=||e||^2=|e_x|^2+|e_y|^2+|e_z|^2=1
+ * the latter condition may be broken for eIncTran of inhomogeneous wave (when msub is complex). In that case
+ * for E=E0*e, ||E||!=|E0| (which is counterintuitive)
+ *
+ * !!! TODO: determine whether they are actually needed in crosssec.c, or make them static here
+ */
+doublecomplex eIncRefl[3],eIncTran[3];
 // used in param.c
 char beam_descr[MAX_MESSAGE2]; // string for log file with beam parameters
 
 // LOCAL VARIABLES
 static double s,s2;            // beam confinement factor and its square
 static double scale_x,scale_z; // multipliers for scaling coordinates
+static doublecomplex ki,kt;   // abs of normal components of k_inc/k0, and ktran/k0
+static doublecomplex ktVec[3]; // k_tran/k0, abs of its normal component
 /* TO ADD NEW BEAM
  * Add here all internal variables (beam parameters), which you initialize in InitBeam() and use in GenerateB()
  * afterwards. If you need local, intermediate variables, put them into the beginning of the corresponding function.
@@ -76,16 +86,45 @@ void InitBeam(void)
 		case B_PLANE:
 			if (IFROOT) strcpy(beam_descr,"plane wave");
 			beam_asym=false;
+			if (surface) {
+				// Here we set ki,kt,ktVec and propagation directions prIncRefl,prIncTran
+				doublecomplex q2; // square of relative component of k_inc/k0 along the surface
+				if (prop_0[2]>0) { // beam comes from the substrate (below)
+					ki=msub*prop_0[2];
+					q2=msub*msub-ki*ki;
+					kt=cSqrtCut(1-q2);
+					// determine propagation direction and full wavevector of wave transmitted into substrate
+					ktVec[0]=msub*prop_0[0];
+					ktVec[1]=msub*prop_0[1];
+					ktVec[2]=kt;
+				}
+				else if (prop_0[2]<0) { // beam comes from above the substrate
+					vRefl(prop_0,prIncRefl);
+					ki=-prop_0[2];
+					q2=1-ki*ki;
+					kt=cSqrtCut(msub*msub-q2);
+					// determine propagation direction of wave transmitted into substrate
+					ktVec[0]=prop_0[0];
+					ktVec[1]=prop_0[1];
+					ktVec[2]=-kt;
+				}
+				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);
+				vReal(ktVec,prIncTran);
+				vNormalize(prIncTran);
+			}
 			return;
 		case B_LMINUS:
 		case B_DAVIS3:
 		case B_BARTON5:
+			if (surface) PrintErrorHelp("Currently, Gaussian incident beam is not supported for '-surf'");
 			// initialize parameters
 			w0=beam_pars[0];
 			TestPositive(w0,"beam width");
 			beam_asym=(beam_Npars==4 && (beam_pars[1]!=0 || beam_pars[2]!=0 || beam_pars[3]!=0));
 			if (beam_asym) {
-				memcpy(beam_center_0,beam_pars+1,3*sizeof(double));
+				vCopy(beam_pars+1,beam_center_0);
 				// if necessary break the symmetry of the problem
 				if (beam_center_0[0]!=0) symX=symR=false;
 				if (beam_center_0[1]!=0) symY=symR=false;
@@ -112,8 +151,8 @@ void InitBeam(void)
 					default: break;
 				}
 				sprintf(beam_descr+strlen(beam_descr),"\tWidth="GFORMDEF" (confinement factor s="GFORMDEF")\n",w0,s);
-				if (beam_asym) sprintf(beam_descr+strlen(beam_descr),"\tCenter position: "GFORMDEF3V,beam_center_0[0],
-					beam_center_0[1],beam_center_0[2]);
+				if (beam_asym)
+					sprintf(beam_descr+strlen(beam_descr),"\tCenter position: "GFORMDEF3V,COMP3V(beam_center_0));
 				else strcat(beam_descr,"\tCenter is in the origin");
 			}
 			return;
@@ -177,12 +216,71 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 	}
 	else { // which==INCPOL_X
 		ex=incPolX;
-		memcpy(ey,incPolY,3*sizeof(double));
+		vCopy(incPolY,ey);
 	}
 
 	switch (beamtype) {
-		case B_PLANE: // plane is separate to be fast
-			for (i=0;i<local_nvoid_Ndip;i++) {
+		case B_PLANE: // plane is separate to be fast (for non-surface)
+			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 beam irradiance is equal to that of a unity-amplitude field in the vacuum
+				 * (i.e. 1/8pi in CGS). Thus, original incident beam propagating from the vacuum (above) is
+				 * Exp(i*k*r.a), while - from the substrate (below) is Exp(i*k*msub*r.a)/sqrt(Re(msub)). We assume that
+				 * the incident beam is homogeneous in its original medium.
+				 */
+				doublecomplex rc,tc; // reflection and transmission coefficients
+				if (prop[2]>0) { // beam comes from the substrate (below)
+					//  determine amplitude of the reflected and transmitted waves
+					if (which==INCPOL_Y) { // s-polarized
+						cvBuildRe(ex,eIncRefl);
+						cvBuildRe(ex,eIncTran);
+						rc=FresnelRS(ki,kt);
+						tc=FresnelTS(ki,kt);
+					}
+					else { // p-polarized
+						vInvRefl_cr(ex,eIncRefl);
+						crCrossProd(ey,ktVec,eIncTran);
+						rc=FresnelRP(ki,kt,1/msub);
+						tc=FresnelTP(ki,kt,1/msub);
+					}
+					// phase shift due to the origin at height hsub
+					cvMultScal_cmplx(rc*cexp(-2*I*WaveNum*ki*hsub)/sqrt(creal(msub)),eIncRefl,eIncRefl);
+					cvMultScal_cmplx(tc*cexp(I*WaveNum*(kt-ki)*hsub)/sqrt(creal(msub)),eIncTran,eIncTran);
+					// main part
+					for (i=0;i<local_nvoid_Ndip;i++) {
+						j=3*i;
+						// b[i] = eIncTran*exp(ik*kt.r)
+						cvMultScal_cmplx(cexp(I*WaveNum*crDotProd(ktVec,DipoleCoord+j)),eIncTran,b+j);
+					}
+				}
+				else if (prop[2]<0) { // beam comes from above the substrate
+					// determine amplitude of the reflected and transmitted waves
+					if (which==INCPOL_Y) { // s-polarized
+						cvBuildRe(ex,eIncRefl);
+						cvBuildRe(ex,eIncTran);
+						rc=FresnelRS(ki,kt);
+						tc=FresnelTS(ki,kt);
+					}
+					else { // p-polarized
+						vInvRefl_cr(ex,eIncRefl);
+						crCrossProd(ey,ktVec,eIncTran);
+						cvMultScal_cmplx(1/msub,eIncTran,eIncTran); // normalize eIncTran by ||ktVec||=msub
+						rc=FresnelRP(ki,kt,msub);
+						tc=FresnelTP(ki,kt,msub);
+					}
+					// phase shift due to the origin at height hsub
+					cvMultScal_cmplx(rc*imExp(2*WaveNum*ki*hsub),eIncRefl,eIncRefl);
+					cvMultScal_cmplx(tc*cexp(I*WaveNum*(ki-kt)*hsub),eIncTran,eIncTran);
+					// main part
+					for (i=0;i<local_nvoid_Ndip;i++) {
+						j=3*i;
+						// b[i] = ex*exp(ik*r.a) + eIncRefl*exp(ik*prIncRefl.r)
+						cvMultScal_RVec(imExp(WaveNum*DotProd(DipoleCoord+j,prop)),ex,b+j);
+						cvLinComb1_cmplx(eIncRefl,b+j,imExp(WaveNum*DotProd(DipoleCoord+j,prIncRefl)),b+j);
+					}
+				}
+			}
+			else for (i=0;i<local_nvoid_Ndip;i++) { // standard (non-surface) plane wave
 				j=3*i;
 				ctemp=imExp(WaveNum*DotProd(DipoleCoord+j,prop)); // ctemp=exp(ik*r.a)
 				cvMultScal_RVec(ctemp,ex,b+j); // b[i]=ctemp*ex

src/calculator.c

@@ -52,8 +52,9 @@ extern size_t TotalEval;
 // used in CalculateE.c
 double * restrict muel_phi; // used to store values of Mueller matrix for different phi (to integrate)
 double * restrict muel_phi_buf; // additional for integrating with different multipliers
-	// scattered E (for scattering in one plane) for two incident polarizations
+	// scattered E (for scattering in the default scattering plane) for two incident polarizations
 doublecomplex * restrict EplaneX, * restrict EplaneY;
+doublecomplex * restrict EyzplX, * restrict EyzplY; // same for scattering in yz-plane
 double dtheta_deg,dtheta_rad; // delta theta in degrees and radians
 doublecomplex * restrict ampl_alphaX,* restrict ampl_alphaY; // amplitude matrix for different values of alpha
 double * restrict muel_alpha; // mueller matrix for different values of alpha
@@ -184,12 +185,11 @@ static void CoupleConstant(doublecomplex *mrel,const enum incpol which,doublecom
 		}
 	}
 	if (asym || anisotropy) {
-		if (!orient_avg && IFROOT) PrintBoth(logfile, "CoupleConstant:"CFORM3V"\n",creal(coup_con[0]),
-			cimag(coup_con[0]),creal(coup_con[1]),cimag(coup_con[1]),creal(coup_con[2]),cimag(coup_con[2]));
+		if (!orient_avg && IFROOT) PrintBoth(logfile, "CoupleConstant:"CFORM3V"\n",REIM3V(coup_con));
 	}
 	else {
 		coup_con[2]=coup_con[1]=coup_con[0];
-		if (!orient_avg && IFROOT) PrintBoth(logfile,"CoupleConstant:"CFORM"\n",creal(coup_con[0]),cimag(coup_con[0]));
+		if (!orient_avg && IFROOT) PrintBoth(logfile,"CoupleConstant:"CFORM"\n",REIM(coup_con[0]));
 	}
 	memcpy(res,coup_con,3*sizeof(doublecomplex));
 }
@@ -263,10 +263,7 @@ static void calculate_one_orientation(double * restrict res)
 	D("MuellerMatrix finished");
 	if (IFROOT && orient_avg) {
 		tstart=GET_TIME();
-		/* it is more logical to use store_mueller in the following test, but for orientation averaging these two flags
-		 * are identical
-		 */
-		if (yzplane) printf("\nError of alpha integration (Mueller) is "GFORMDEF"\n",
+		if (store_mueller) printf("\nError of alpha integration (Mueller) is "GFORMDEF"\n",
 			Romberg1D(parms_alpha,block_theta,muel_alpha,res+2));
 		memcpy(res,muel_alpha-2,2*sizeof(double));
 		D("Integration over alpha completed on root");
@@ -369,6 +366,16 @@ static void AllocateEverything(void)
 	MALLOC_VECTOR(expsZ,complex,local_Nz_unif,ALL);
 #endif // !SPARSE
 	if (yzplane) {
+		tmp=2*(double)nTheta;
+		if (!prognosis) {
+			CheckOverflow(2*tmp,ONE_POS_FUNC);
+			temp_int=tmp;
+			MALLOC_VECTOR(EyzplX,complex,temp_int,ALL);
+			MALLOC_VECTOR(EyzplY,complex,temp_int,ALL);
+		}
+		memory+=2*tmp*sizeof(doublecomplex);
+	}
+	if (scat_plane) {
 		tmp=2*(double)nTheta;
 		if (!prognosis) {
 			CheckOverflow(2*tmp,ONE_POS_FUNC);
@@ -419,7 +426,7 @@ static void AllocateEverything(void)
 		if (!prognosis) {
 			// this covers these 2 and next 2 malloc calls
 			CheckOverflow(8*tmp+2,ONE_POS_FUNC);
-			if (yzplane) {
+			if (store_mueller) {
 				temp_int=tmp;
 				MALLOC_VECTOR(ampl_alphaX,complex,temp_int,ONE);
 				MALLOC_VECTOR(ampl_alphaY,complex,temp_int,ONE);
@@ -517,6 +524,10 @@ static void FreeEverything(void)
 	 * in AllocateEverything() above.
 	 */
 	if (yzplane) {
+		Free_cVector(EyzplX);
+		Free_cVector(EyzplY);
+	}
+	if (scat_plane) {
 		Free_cVector(EplaneX);
 		Free_cVector(EplaneY);
 	}
@@ -542,7 +553,7 @@ static void FreeEverything(void)
 
 	if (orient_avg) {
 		if (IFROOT) {
-			if (yzplane) {
+			if (store_mueller) {
 				Free_cVector(ampl_alphaX);
 				Free_cVector(ampl_alphaY);
 			}