Changes in tests

src/GenerateB.c

@@ -49,6 +49,11 @@ double C0dipole,C0dipole_refl; // inherent cross sections of exciting dipole (in
 int vorticity;                 // Vorticity of vortex beams (besN for Bessel beams)
 // 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
+ */
+double prIncRefl[3],prIncTran[3];
 
 // LOCAL VARIABLES
 static double s,s2;            // beam confinement factor and its square
@@ -119,7 +124,6 @@ void InitBeam(void)
 				}
 				else if (prop_0[2]<0) { // beam comes from above the substrate
 					inc_scale=1;
-					vRefl(prop_0,prIncRefl);
 					ki=-prop_0[2]; // always real
 					if (!msubInf) {
 						// same special case as above
@@ -179,7 +183,7 @@ void InitBeam(void)
 					default: LogError(ONE_POS,"Incompatibility error in GenerateB");
 				}
 				beam_descr=dyn_sprintf("Gaussian beam (%s)\n"
-				                       "\tWidth="GFORMDEF" (confinement factor s="GFORMDEF")\n",tmp_str,w0,s);
+				                       "\tWidth="GFORMDEF" (confinement factor s="GFORMDEF")",tmp_str,w0,s);
 			}
 			return;
 #ifndef NO_FORTRAN
@@ -260,7 +264,7 @@ void InitBeam(void)
 			symR=symX=symY=symZ=false;
 			// beam info
 			if (IFROOT) beam_descr=dyn_sprintf("Bessel beam (%s)\n"
-				                               "\tOrder: %d, half-cone angle: "GFORMDEF" deg\n",
+				                               "\tOrder: %d, half-cone angle: "GFORMDEF" deg",
 				                               tmp_str,besN,beam_pars[1]);
 			return;
 #endif // !NO_FORTRAN
@@ -287,7 +291,10 @@ void InitBeam(void)
 	 *    false here. Do not set any of them to true, as they can be set to false by other factors.
 	 *    symX, symY, symZ - symmetries of reflection over planes YZ, XZ, XY respectively.
 	 *    symR - symmetry of rotation for 90 degrees over the Z axis
-	 * 4) initialize beam_descr - descriptive string, which will appear in log file.
+	 * 4) initialize the following:
+	 *    beam_descr - descriptive string, which will appear in log file  (should NOT end with \n).
+	 *    vorticity - (only for vortex beam) integer value, how many turns the phase experience, when one makes a full
+	 *                turn around the beam axis.
 	 * 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(...);". Otherwise, you also need to define inc_scale.
 	 * All other auxiliary variables, which are used in beam generation (GenerateB(), see below), should be defined in
@@ -311,14 +318,14 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 	const double *ex; // coordinate axis of the beam reference frame
 	double ey[3];
 	double r1[3];
-	/* complex wave amplitudes of secondary waves (with phase relative to particle center);
+	/* complex wave amplitudes of transmitted wave (with phase relative to beam center);
 	 * The transmitted wave can be inhomogeneous wave (when msub is complex), then eIncTran (e) is normalized
 	 * counter-intuitively. Before multiplying by tc, it satisfies (e,e)=1!=||e||^2. This normalization is consistent
 	 * with used formulae for transmission coefficients. So this transmission coefficient is not (generally) equal to
 	 * the ratio of amplitudes of the electric fields. In particular, when E=E0*e, ||E||!=|E0|*||e||, where
 	 * ||e||^2=(e,e*)=|e_x|^2+|e_y|^2+|e_z|^2=1
 	 */
-	doublecomplex eIncRefl[3],eIncTran[3];
+	doublecomplex eIncTran[3];
 #ifndef NO_FORTRAN
 	// for Bessel beams
 	int n1,q;
@@ -352,10 +359,10 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 				 * the substrate (below) is Exp(i*k*msub*r.a). We assume that the incoming beam is homogeneous in its
 				 * original medium.
 				 */
-				doublecomplex rc,tc; // reflection and transmission coefficients
 				double hbeam=hsub+beam_center[2]; // height of beam center above the surface 
 				if (prop[2]>0) { // beam comes from the substrate (below)
-					//  determine amplitude of the reflected and transmitted waves; here msub is always defined
+					doublecomplex tc; // transmission coefficients
+					//  determine amplitude of the transmitted wave; here msub is always defined
 					if (which==INCPOL_Y) { // s-polarized
 						cvBuildRe(ex,eIncTran);
 						tc=FresnelTS(ki,kt);
@@ -374,40 +381,52 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 					}
 				}
 				else if (prop[2]<0) { // beam comes from above the substrate
-					// determine amplitude of the reflected and transmitted waves
+					/* The following code takes extra care to be stable (not to lose precision) for grazing incidence.
+					 * While it seems more complicated for general incidence, it requires only a few more complex
+					 * arithmetic operations that are negligible compared to two complex exponents. For grazing
+					 * incidence, it is not only stable but may also be faster than the straightforward computation,
+					 * since one of the complex exponents is computed from a very small argument
+					 */
+					// determine reflection coefficient + 1
+					doublecomplex rcp1;
 					if (which==INCPOL_Y) { // s-polarized
-						cvBuildRe(ex,eIncRefl);
-						if (msubInf) rc=-1;
-						else {
-							cvBuildRe(ex,eIncTran);
-							rc=FresnelRS(ki,kt);
-						}
+						if (msubInf) rcp1=0;
+						else rcp1=FresnelTS(ki,kt);
 					}
 					else { // p-polarized
-						vInvRefl_cr(ex,eIncRefl);
-						if (msubInf) rc=1;
-						else {
-							crCrossProd(ey,ktVec,eIncTran);
-							cvMultScal_cmplx(1/msub,eIncTran,eIncTran); // normalize eIncTran by ||ktVec||=msub
-							rc=FresnelRP(ki,kt,msub);
-						}
+						if (msubInf) rcp1=2;
+						else rcp1=msub*FresnelTP(ki,kt,msub);
 					}
-					// phase shift due to the beam center relative to the surface
-					cvMultScal_cmplx(rc*imExp(2*WaveNum*creal(ki)*hbeam),eIncRefl,eIncRefl); // assumes real ki
 					// main part
 					for (i=0;i<local_nvoid_Ndip;i++) {
 						j=3*i;
-						// together with above phase shift, beam_center affects only the common phase factor
-						LinComb(DipoleCoord+j,beam_center,1,-1,r1); 
-						// b[i] = ex*exp(ik*r.a) + eIncRefl*exp(ik*prIncRefl.r)
-						cvMultScal_RVec(imExp(WaveNum*DotProd(r1,prop)),ex,b+j);
-						cvLinComb1_cmplx(eIncRefl,b+j,imExp(WaveNum*DotProd(r1,prIncRefl)),b+j);
+						vSubtr(DipoleCoord+j,beam_center,r1); // beam_center affects only the common phase factor
+						/* b[i] = ex*exp(ik*r.a) + eIncRefl*exp[ik(prIncRefl.r-2az*hbeam)],
+						 * where prIncRefl is prop with reflected z-component, while eIncRefl=rc*ex for s-polarization
+						 * and (additionally) with reflected transverse (x,y) components for p-polarization
+						 */
+						ctemp=imExp(WaveNum*DotProd(r1,prop)); // exp(ik*r.a)
+						t1=imExpM1(-2*WaveNum*prop[2]*(r1[2]+hbeam));
+						/* t2=exp(ik*r.a){1+rc*exp[-2i*kz(z+hbeam)]}, but the terms in parentheses (rc and exp(...)) are
+						 * replaced by their differences with -1 and 1, respectively (eliminating 1 inside)
+						 */
+						t2=ctemp*(rcp1*(t1+1)-t1);
+						cvMultScal_RVec(t2,ex,b+j);
+						if (which==INCPOL_X) {
+							/* here we add (eIncRefl-ex*rc)*exp[ik(prIncRefl.r-2az*hbeam)] (non-zero only for
+							 * p-polarization), expressing rc*exp(...) as t2-ctemp
+							 */
+							t3=2*(t2-ctemp);
+							b[j]-=t3*ex[0];
+							b[j+1]-=t3*ex[1];
+						}
 					}
 				}
 			}
 			else for (i=0;i<local_nvoid_Ndip;i++) { // standard (non-surface) plane wave
 				j=3*i;
-				LinComb(DipoleCoord+j,beam_center,1,-1,r1); // this affects only the common phase factor
+				// can be replaced by complex multiplication (by precomputed phase factor), does not seem beneficial
+				vSubtr(DipoleCoord+j,beam_center,r1);
 				ctemp=imExp(WaveNum*DotProd(r1,prop)); // ctemp=exp(ik*r.a)
 				cvMultScal_RVec(ctemp,ex,b+j); // b[i]=ctemp*ex
 			}
@@ -417,7 +436,7 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 			vMultScal(p0,prop,dip_p);
 			for (i=0;i<local_nvoid_Ndip;i++) { // here we explicitly use that dip_p is real
 				j=3*i;
-				LinComb(DipoleCoord+j,beam_center,1,-1,r1);
+				vSubtr(DipoleCoord+j,beam_center,r1);
 				(*InterTerm_real)(r1,gt);
 				cSymMatrVecReal(gt,dip_p,b+j);
 				if (surface) { // add reflected field
@@ -459,7 +478,7 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 			for (i=0;i<local_nvoid_Ndip;i++) {
 				j=3*i;
 				// set relative coordinates (in beam's coordinate system)
-				LinComb(DipoleCoord+j,beam_center,1,-1,r1);
+				vSubtr(DipoleCoord+j,beam_center,r1);
 				x=DotProd(r1,ex)*scale_x;
 				y=DotProd(r1,ey)*scale_x;
 				z=DotProd(r1,prop)*scale_z;
@@ -532,7 +551,7 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 			vort=(which==INCPOL_Y) ? cpow(I,besN) : 1;
 			for (i=0;i<local_nvoid_Ndip;i++) {
 				j=3*i;
-				LinComb(DipoleCoord+j,beam_center,1,-1,r1);
+				vSubtr(DipoleCoord+j,beam_center,r1);
 				x=DotProd(r1,ex);
 				y=DotProd(r1,ey);
 				z=DotProd(r1,prop);

src/cmplx.h

@@ -101,6 +101,39 @@ static inline doublecomplex imExp(const double arg)
 #endif
 }
 
+//======================================================================================================================
+
+static inline doublecomplex imExpM1(const double arg)
+/* exp(i*arg) - 1 (should be used for small argument to avoid precision loss
+ * We employ special code only for small arguments, ignoring the case when arg is close to 2piN. The latter can,
+ * in principle, be handled by preliminary range reduction as in imExpTable. We do not implement it here, because such
+ * case is a "coincidence" - may happen for a single dipole (or a plane of dipoles), while the case of small arg may
+ * happen for all dipoles. In the latter case loss of precision affects all computed quantities.
+ * The used expression through tan(arg/2) follows from general expression for cexpm1 below.
+ */
+{
+	if (fabs(arg)<1) {
+		double t=tan(0.5*arg);
+		return -2*t/(I+t); // Alternatively, (I-t)*2t/(1+t^2), but we leave the optimization to compiler
+	}
+	else return imExp(arg)-1;
+}
+
+//======================================================================================================================
+
+static inline doublecomplex cExpM1(const doublecomplex a)
+/* Complex analogue of expm1 function ( exp(a) - 1 ), should be used for small arguments to avoid precision loss
+ * The algorithm is a simplified version of the one published in Section 17.7 of Beebe N.H.F., The Mathematical-Function
+ * Computation Handbook: Programming Using the MathCW Portable Software Library. Springer; 2017.
+ *  */
+{
+	if (fabs(creal(a))+fabs(cimag(a))<1) { // uses faster L1-norm instead of L2-norm
+		doublecomplex t=ctanh(0.5*a);
+		return 2*t/(1-t);
+	}
+	else return cexp(a)-1;
+}
+
 //======================================================================================================================
 // operations on complex vectors
 

src/param.c

@@ -95,6 +95,7 @@ bool emulLinebuf;  // whether to emulate line buffering of stdout, defined as ex
 extern const char *avg_string;
 // defined and initialized in GenerateB.c
 extern const char *beam_descr;
+extern const double prIncRefl[3],prIncTran[3];
 // defined and initialized in make_particle.c
 extern const bool volcor_used;
 extern const char *sh_form_str1,*sh_form_str2;
@@ -1641,7 +1642,7 @@ PARSE_FUNC(test)
 }
 PARSE_FUNC(V)
 {
-	char copyright[]="\n\nCopyright (C) 2006-2021 ADDA contributors\n"
+	char copyright[]="\n\nCopyright (C) 2006-2022 ADDA contributors\n"
 		"This program is free software; you can redistribute it and/or modify it under the terms of the GNU General "
 		"Public License as published by the Free Software Foundation; either version 3 of the License, or (at your "
 		"option) any later version.\n\n"

src/vars.c

@@ -133,11 +133,6 @@ 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
-/* Propagation (phase) directions of secondary incident beams above (A) and below (B) 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
- */
-double prIncRefl[3],prIncTran[3];
 
 #ifndef SPARSE // These variables are exclusive to the FFT mode
 

src/vars.h

@@ -84,7 +84,7 @@ extern TIME_TYPE Timing_EField,Timing_FileIO,Timing_Integration,tstart_main;
 extern bool surface,msubInf;
 extern enum refl ReflRelation;
 extern doublecomplex msub;
-extern double inc_scale,hsub,prIncRefl[3],prIncTran[3];
+extern double inc_scale,hsub;
 
 #ifndef SPARSE // These variables are exclusive to the FFT mode
 

tests/2exec/comp2exec

@@ -97,8 +97,10 @@ ADDASPAOCL="./../../src/ocl/adda_spa_ocl" # this combination is not yet supporte
 # Currently formulated to compare 1.5* with 1.4.0 (more specific exceptions are added in parsing of command lines)
 # These variables must be updated whenever new features are added to ADDA
 
-CMDIGNORE="coated2|superellipsoid|-h anisotr|-h int|-h pol|-h rect_dip|-h scat|-int igt_so"
+CMDIGNORE="coated2|superellipsoid|-h anisotr|-h int|-h pol|-h rect_dip|-h scat|-int igt_so|-h beam|-beam bessel|1e-16"
 OLDIGNORE="^ADDA v\.|^Copyright \(C\) 2006-20|^  coated2|^  superellipsoid|^  -int|^  -pol|^  -scat"
+OLDIGNORE="$OLDIGNORE|^Incident beam center position:|^Incident beam: point dipole|^  -beam_center|Center position:"
+OLDIGNORE="$OLDIGNORE|^parameters of predefined shapes"
 
 # Path to reference binaries, examples: "." or "./../../win64". 
 # This may work out of the box on Windows, but requires changes for other systems

tests/2exec/suite

@@ -75,14 +75,14 @@ all -h beam
 all -h beam plane
 all -beam plane ;mgn;
 all -h beam lminus
-all -beam lminus 2 1 2 3 ;mgn;
+all -beam lminus 2 -beam_center 1 2 3 ;mgn;
 all -h beam davis3
-all -beam davis3 2 1 2 3 ;mgn;
+all -beam davis3 2 -beam_center 1 2 3 ;mgn;
 all -beam davis3 2 ;mg9n;
 all -h beam barton5
-all -beam barton5 2 1 2 3 ;mgn;
+all -beam barton5 2 -beam_center 1 2 3 ;mgn;
 all -beam barton5 2 ;mg9n;
-all -beam dipole 3 2 1 ;p; ;mgn;
+all -beam dipole -beam_center 3 2 1 ;p; ;mgn;
 all -h beam read
 all -beam read IncBeam-Y IncBeam-X ;se; ;mgn;
 
@@ -376,8 +376,8 @@ all -surf 4 2 0 ;mgn;
 all -surf 4 3 4 ;p; ;mgn;
 all -surf 4 3 4 -prop 1 2 -3 ;se; ;mgn;
 all -surf 4 inf -prop 1 2 -3 ;se; ;mgn;
-all -surf 4 2 1 -beam dipole 3 2 1 ;p; ;mgn;
-all -surf 4 inf -beam dipole 3 2 1 ;p; ;mgn;
+all -surf 4 2 1 -beam dipole -beam_center 3 2 1 ;p; ;mgn;
+all -surf 4 inf -beam dipole -beam_center 3 2 1 ;p; ;mgn;
 all -surf 4 2 0 -no_reduced_fft ;mgn;
 !mpi,!mpi_seq -surf 4 2 0 -iter cgnr -no_reduced_fft ;mgn;
 

tests/2exec/suite_rd

@@ -78,15 +78,15 @@ all -h beam
 all -h beam plane
 all -beam plane ;mgn;
 all -h beam lminus
-all -beam lminus 2 1 2 3 ;mgn;
+all -beam lminus 2 -beam_center 1 2 3 ;mgn;
 all -h beam davis3
-all -beam davis3 2 1 2 3 ;mgn;
+all -beam davis3 2 -beam_center 1 2 3 ;mgn;
 all -beam davis3 2 ;mg9n;
 all -h beam barton5
-all -beam barton5 2 1 2 3 ;mgn;
+all -beam barton5 2 -beam_center 1 2 3 ;mgn;
 all -beam barton5 2 ;mg9n;
-!RD_TRICKY -beam dipole 3 2 1 ;p; ;mgn;
-&RD_TRICKY -beam dipole 3 2 1 ;smg12n;
+!RD_TRICKY -beam dipole -beam_center 3 2 1 ;p; ;mgn;
+&RD_TRICKY -beam dipole -beam_center 3 2 1 ;smg12n;
 all -h beam read
 #all -beam read IncBeam-Y IncBeam-X ;se; ;mgn;
 
@@ -340,10 +340,10 @@ all -surf 4 2 0 ;mgn;
 !RD_TRICKY -surf 4 3 4 -prop 1 2 -3 ;se; ;mgn;
 !RD_TRICKY -surf 4 inf -prop 1 2 -3 ;se; ;mgn;
 &RD_TRICKY -surf 4 inf -prop 1 2 -3 ;se; ;m; ;n; -size 5 -grid 12
-!RD_TRICKY -surf 4 2 1 -beam dipole 3 2 1 ;p; ;mgn;
-&RD_TRICKY -surf 4 2 1 -beam dipole 3 2 1 ;p; ;smg12n;
-!RD_TRICKY -surf 4 inf -beam dipole 3 2 1 ;p; ;mgn;
-&RD_TRICKY -surf 4 2 1 -beam dipole 3 2 1 ;p; ;smg12n;
+!RD_TRICKY -surf 4 2 1 -beam dipole -beam_center 3 2 1 ;p; ;mgn;
+&RD_TRICKY -surf 4 2 1 -beam dipole -beam_center 3 2 1 ;p; ;smg12n;
+!RD_TRICKY -surf 4 inf -beam dipole -beam_center 3 2 1 ;p; ;mgn;
+&RD_TRICKY -surf 4 2 1 -beam dipole -beam_center 3 2 1 ;p; ;smg12n;
 all -surf 4 2 0 -no_reduced_fft ;mgn;
 !mpi,!mpi_seq -surf 4 2 0 -iter cgnr -no_reduced_fft ;mgn; 
 

tests/2exec/suite_sparse

@@ -56,14 +56,14 @@ all -h beam
 all -h beam plane
 all -beam plane ;mgn;
 all -h beam lminus
-all -beam lminus 2 1 2 3 ;mgn;
+all -beam lminus 2 -beam_center 1 2 3 ;mgn;
 all -h beam davis3
-all -beam davis3 2 1 2 3 ;mgn;
+all -beam davis3 2 -beam_center 1 2 3 ;mgn;
 all -beam davis3 2 ;mg9n;
 all -h beam barton5
-all -beam barton5 2 1 2 3 ;mgn;
+all -beam barton5 2 -beam_center 1 2 3 ;mgn;
 all -beam barton5 2 ;mg9n;
-all -beam dipole 3 2 1 ;p; ;mgn;
+all -beam dipole -beam_center 3 2 1 ;p; ;mgn;
 all -h beam read
 all -beam read IncBeam-Y IncBeam-X ;se; ;mn;
 
@@ -356,8 +356,8 @@ all -surf 4 2 0 ;mgn;
 all -surf 4 3 4 ;p; ;ss; ;mn;
 all -surf 4 3 4 -prop 1 2 -3 ;se; ;mn;
 all -surf 4 inf -prop 1 2 -3 ;se; ;mn;
-all -surf 4 2 1 -beam dipole 3 2 1 ;p; ;ss; ;mn;
-all -surf 4 inf -beam dipole 3 2 1 ;p; ;ss; ;mn;
+all -surf 4 2 1 -beam dipole -beam_center 3 2 1 ;p; ;ss; ;mn;
+all -surf 4 inf -beam dipole -beam_center 3 2 1 ;p; ;ss; ;mn;
 all -surf 4 2 0 -no_reduced_fft ;mgn;
 !mpi,!mpi_seq -surf 4 2 0 -iter cgnr -no_reduced_fft ;mgn;
 

tests/2exec/suite_surf

@@ -78,14 +78,14 @@ all -h beam
 all -h beam plane
 all -beam plane ;mgn;
 all -h beam lminus
-#all -beam lminus 2 1 2 3 ;mgn;
+#all -beam lminus 2 -beam_center 1 2 3 ;mgn;
 all -h beam davis3
-#all -beam davis3 2 1 2 3 ;mgn;
+#all -beam davis3 2 -beam_center 1 2 3 ;mgn;
 #all -beam davis3 2 ;mg9n;
 all -h beam barton5
-#all -beam barton5 2 1 2 3 ;mgn;
+#all -beam barton5 2 -beam_center 1 2 3 ;mgn;
 #all -beam barton5 2 ;mg9n;
-all -beam dipole 3 2 1 ;p; ;mgn;
+all -beam dipole -beam_center 3 2 1 ;p; ;mgn;
 all -h beam read
 all -beam read IncBeam-Y IncBeam-X ;se; ;mgn;