Implementation of point-dipole incident field (issue 149). Should wor…

…k both for free space and in the presence of surface (dipole should be above the surface). Controlled by new beam type: '-beam dipole <x> <y> <z>', dipole is assumed to have unity amplitude with direction, given by '-prop ...' command line option (e_z by default).

To do this the logic in interaction.c/h was significantly modified. Now two types of functions are defined (visible outside) for both direct and reflected interaction, accepting either integers (in units of d, used for filling interaction matrix) or real values (in um, useful for other applications, such as incident dipole field or nearfields). Thus the whole existing machinery of interaction.c/h (all formulations, numerical routines, etc.) can easily be accessed from other parts of the code.
- internally, those functions are usually wrappers which reuse most of the code by calling inline functions.
- some of the interaction formulations (igt_so and so) are incompatible with arbitrary real arguments, since they are tightly linked to existing tables. If used they produce a careful exception (new flag InteractionRealArgs added for that).
- computing reflected Green's tensor (based on Sommerfeld integrals) was split into modules, including the tabulation. Now further optimization efforts can be directed precisely.

Still a lot to do - to test, integrate (add exceptions) to other parts of ADDA, add enhancement rates (radiative and non-radiative) to the output.

Other (related) changes:
- a number of exceptions with respect to combination of different beam types, -prop, and -surf were moved from VariablesInterconnect to InitBeam. Now they are more accurate.
- beam info for Gaussian beam was made more uniform. For centered beams, "Center is in the origin" was changed to "Center position: (0,0,0)".
- new functions cSymMatrVecReal and cReflMatrVecReal in cmplx.h
- version incremented to 1.3b2.
- test suites in tests/2exec/ were updated to include new features.

src/GenerateB.c

@@ -30,6 +30,7 @@
 #include "cmplx.h"
 #include "comm.h"
 #include "io.h"
+#include "interaction.h"
 #include "param.h"
 #include "vars.h"
 // system headers
@@ -89,6 +90,10 @@ void InitBeam(void)
 			if (IFROOT) strcpy(beam_descr,"plane wave");
 			beam_asym=false;
 			if (surface) {
+				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 "
+					"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
@@ -119,17 +124,30 @@ void InitBeam(void)
 				}
 			}
 			return;
+		case B_DIPOLE:
+			vCopy(beam_pars,beam_center_0);
+			if (surface && beam_center_0[2]<=-hsub)
+				PrintErrorHelp("External dipole should be placed strictly above the surface");
+			// in weird scenarios the dipole can be positioned exactly at the origin; reused code from Gaussian beams
+			beam_asym=(beam_center_0[0]!=0 || beam_center_0[1]!=0 || beam_center_0[2]!=0);
+			if (beam_asym) { // 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;
+				if (beam_center_0[2]!=0) symZ=false;
+			}
+			else vInit(beam_center);
+			if (IFROOT) sprintf(beam_descr,"point dipole at "GFORMDEF3V,COMP3V(beam_center_0));
+			return;
 		case B_LMINUS:
 		case B_DAVIS3:
 		case B_BARTON5:
-			if (surface) PrintErrorHelp("Currently, Gaussian incident beam is not supported for '-surf'");
+			if (surface) PrintError("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) {
-				vCopy(beam_pars+1,beam_center_0);
-				// if necessary break the symmetry of the problem
+			vCopy(beam_pars+1,beam_center_0);
+			beam_asym=(beam_Npars==4 && (beam_center_0[0]!=0 || beam_center_0[1]!=0 || beam_center_0[2]!=0));
+			if (beam_asym) { // 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;
 				if (beam_center_0[2]!=0) symZ=false;
@@ -154,10 +172,8 @@ void InitBeam(void)
 						break;
 					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,COMP3V(beam_center_0));
-				else strcat(beam_descr,"\tCenter is in the origin");
+				sprintf(beam_descr+strlen(beam_descr),"\tWidth="GFORMDEF" (confinement factor s="GFORMDEF")\n"
+				                                      "\tCenter position: "GFORMDEF3V,w0,s,COMP3V(beam_center_0));
 			}
 			return;
 		case B_READ:
@@ -186,6 +202,8 @@ void InitBeam(void)
 	 *    beam_asym - whether beam center does not coincide with the reference frame origin. If it is set to true, then
 	 *                set also beam_center_0 - 3D radius-vector of beam center in the laboratory reference frame (it
 	 *                will be then automatically transformed to particle reference frame, if required).
+	 * 5) Consider the case of surface (substrate near the particle). If the new beam type is incompatible with it, add
+	 *    an explicit exception, like "if (surface) PrintErrorHelp(...);".
 	 * All other auxiliary variables, which are used in beam generation (GenerateB(), see below), should be defined in
 	 * the beginning of this file. If you need temporary local variables (which are used only in this part of the code),
 	 * define them in the beginning of this function.
@@ -200,7 +218,7 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 {
 	size_t i,j;
 	doublecomplex psi0,Q,Q2;
-	doublecomplex v1[3],v2[3],v3[3];
+	doublecomplex v1[3],v2[3],v3[3],gt[6];
 	double ro2,ro4;
 	double x,y,z,x2_s,xy_s;
 	doublecomplex t1,t2,t3,t4,t5,t6,t7,t8,ctemp;
@@ -296,6 +314,20 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 				cvMultScal_RVec(ctemp,ex,b+j); // b[i]=ctemp*ex
 			}
 			return;
+		case B_DIPOLE:
+			for (i=0;i<local_nvoid_Ndip;i++) { // here we explicitly use that dipole moment (prop) is real
+				j=3*i;
+				LinComb(DipoleCoord+j,beam_center,1,-1,r1);
+				(*InterTerm_real)(r1,gt);
+				cSymMatrVecReal(gt,prop,b+j);
+				if (surface) { // add reflected field
+					r1[2]=DipoleCoord[j+2]+beam_center[2]+2*hsub;
+					(*ReflTerm_real)(r1,gt);
+					cReflMatrVecReal(gt,prop,v1);
+					cvAdd(v1,b+j,b+j);
+				}
+			}
+			return;
 		case B_LMINUS:
 		case B_DAVIS3:
 		case B_BARTON5:
@@ -372,6 +404,8 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 	 * 4) 'ey' – complementary unity vector of polarization (orthogonal to both 'prop' and 'ex');
 	 * 5) 'beam_center' – beam center in the particle reference frame (automatically calculated from 'beam_center_0'
 	 *                    defined in InitBeam).
+	 * If the new beam type is compatible with '-surf', include here the corresponding code. For that you will need
+	 * the variables, related to surface - see vars.c after "// related to a nearby surface".
 	 * If you need temporary local variables (which are used only in this part of the code), either use 't1'-'t8' or
 	 * define your own (with more informative names) in the beginning of this function.
 	 */

src/cmplx.h

@@ -308,6 +308,17 @@ static inline void cSymMatrVec(const doublecomplex matr[static 6],const doubleco
 
 //======================================================================================================================
 
+static inline void cSymMatrVecReal(const doublecomplex matr[static 6],const double vec[static 3],
+	doublecomplex res[static 3])
+// multiplication of complex symmetric matrix[6] by  vec[3]; res=matr.vec; !!! vec and res must not alias
+{
+	res[0] = matr[0]*vec[0] + matr[1]*vec[1] + matr[2]*vec[2];
+	res[1] = matr[1]*vec[0] + matr[3]*vec[1] + matr[4]*vec[2];
+	res[2] = matr[2]*vec[0] + matr[4]*vec[1] + matr[5]*vec[2];
+}
+
+//======================================================================================================================
+
 static inline void cReflMatrVec(const doublecomplex matr[static 6],const doublecomplex vec[static 3],
 	doublecomplex res[static 3])
 /* multiplication of matrix[6] by complex vec[3]; res=matr.vec; passed components are the same as for symmetric matrix:
@@ -320,6 +331,20 @@ static inline void cReflMatrVec(const doublecomplex matr[static 6],const doublec
 	res[2] = - matr[2]*vec[0] - matr[4]*vec[1] + matr[5]*vec[2];
 }
 
+//======================================================================================================================
+
+static inline void cReflMatrVecReal(const doublecomplex matr[static 6],const double vec[static 3],
+	doublecomplex res[static 3])
+/* multiplication of matrix[6] by complex vec[3]; res=matr.vec; passed components are the same as for symmetric matrix:
+ * 11,12,13,22,23,33, but the matrix has the following symmetry - M21=M12, M31=-M13, M32=-M23
+ * !!! vec and res must not alias
+ */
+{
+	res[0] = matr[0]*vec[0] + matr[1]*vec[1] + matr[2]*vec[2];
+	res[1] = matr[1]*vec[0] + matr[3]*vec[1] + matr[4]*vec[2];
+	res[2] = - matr[2]*vec[0] - matr[4]*vec[1] + matr[5]*vec[2];
+}
+
 //======================================================================================================================
 // operations on real vectors
 
@@ -420,6 +445,14 @@ static inline double DotProd(const double a[static 3],const double b[static 3])
 
 //======================================================================================================================
 
+static inline double vNorm(const double a[static 3])
+// norm of a real vector[3]
+{
+	return sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);
+}
+
+//======================================================================================================================
+
 static inline void CrossProd(const double a[static 3],const double b[static 3],double c[static 3])
 // cross product of two real vectors; c = a x b; !!! vectors must not alias
 {

src/const.h

@@ -18,7 +18,7 @@
 #define __const_h
 
 // version number (string)
-#define ADDA_VERSION "1.3b1"
+#define ADDA_VERSION "1.3b2"
 
 /* ADDA uses certain C99 extensions, which are widely supported by GNU and Intel compilers. However, they may be not
  * completely supported by e.g. Microsoft Visual Studio compiler. Therefore, we check the version of the standard here
@@ -256,6 +256,7 @@ enum incpol {
 enum beam { // beam types
 	B_BARTON5, // 5th order description of the Gaussian beam
 	B_DAVIS3,  // 3rd order description of the Gaussian beam
+	B_DIPOLE,  // field of a point dipole
 	B_LMINUS,  // 1st order description of the Gaussian beam
 	B_PLANE,   // infinite plane wave
 	B_READ     // read from file

src/fft.c

@@ -803,7 +803,7 @@ static void InitRmatrix(const double invNgrid)
 	// fill Rmatrix with values of reflected Green's tensor
 	for(k=0;k<local_Nz_Rm;k++) for (j=jstartR;j<boxY;j++) for (i=1-boxX;i<boxX;i++) {
 			index=NDCOMP*Index2matrix(i,j,k,R2sizeY);
-			(*CalcReflTerm)(i,j,k,Rmatrix+index);
+			(*ReflTerm_int)(i,j,k,Rmatrix+index);
 	} // end of i,j,k loop
 	if (IFROOT) printf("Fourier transform of Rmatrix");
 	for(Rcomp=0;Rcomp<NDCOMP;Rcomp++) { // main cycle over components of Rmatrix
@@ -1105,7 +1105,7 @@ void InitDmatrix(void)
 			 * 1) complicate the loops to remove the zero element in the beginning (move tests to the upper level)
 			 * 2) call the function with zero - it will produce NaN. Then set this element to zero after the loop.
 			 */
-			if (i!=0 || j!=0 || kcor!=0) (*CalcInterTerm)(i,j,kcor,Dmatrix+index);
+			if (i!=0 || j!=0 || kcor!=0) (*InterTerm_int)(i,j,kcor,Dmatrix+index);
 		}
 	} // end of i,j,k loop
 	if (IFROOT) printf("Fourier transform of Dmatrix");