From 3372bc25d03711aa2e61540fa483c4ca5d886639 Mon Sep 17 00:00:00 2001
From: stefanyagl <42984899+stefanyagl@users.noreply.github.com>
Date: Sun, 9 Jun 2019 00:08:30 +0700
Subject: [PATCH 01/46] Add files via upload
Test version of Bessel beams implementation in ADDA
---
GenerateB.c | 529 +++++++++++
bess.c | 352 +++++++
bess.h | 27 +
const.h | 427 +++++++++
param.c | 2568 +++++++++++++++++++++++++++++++++++++++++++++++++++
5 files changed, 3903 insertions(+)
create mode 100644 GenerateB.c
create mode 100644 bess.c
create mode 100644 bess.h
create mode 100644 const.h
create mode 100644 param.c
diff --git a/GenerateB.c b/GenerateB.c
new file mode 100644
index 00000000..ebd6f92c
--- /dev/null
+++ b/GenerateB.c
@@ -0,0 +1,529 @@
+/* File: GenerateB.c
+ * $Date:: $
+ * Descr: generate a incident beam
+ *
+ * Lminus beam is based on: G. Gouesbet, B. Maheu, G. Grehan, "Light scattering from a sphere arbitrary located
+ * in a Gaussian beam, using a Bromwhich formulation", J.Opt.Soc.Am.A 5,1427-1443 (1988).
+ * Eq.(22) - complex conjugate.
+ *
+ * Davis beam is based on: L. W. Davis, "Theory of electromagnetic beams," Phys.Rev.A 19, 1177-1179 (1979).
+ * Eqs.(15a),(15b) - complex conjugate; in (15a) "Q" changed to "Q^2" (typo).
+ *
+ * Barton beam is based on: J. P. Barton and D. R. Alexander, "Fifth-order corrected electromagnetic-field
+ * components for a fundamental Gaussian-beam," J.Appl.Phys. 66,2800-2802 (1989).
+ * Eqs.(25)-(28) - complex conjugate.
+ *
+ * Copyright (C) 2006-2014 ADDA contributors
+ * This file is part of ADDA.
+ *
+ * ADDA 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.
+ *
+ * ADDA is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License along with ADDA. If not, see
+ * .
+ */
+#include "const.h" // keep this first
+// project headers
+#include "cmplx.h"
+#include "comm.h"
+#include "io.h"
+#include "interaction.h"
+#include "param.h"
+#include "vars.h"
+#include
+#include
+#include "bess.h"
+
+// SEMI-GLOBAL VARIABLES
+
+// defined and initialized in param.c
+extern const int beam_Npars;
+extern const double beam_pars[];
+extern const char *beam_fnameY;
+extern const char *beam_fnameX;
+extern const opt_index opt_beam;
+
+// used in CalculateE.c
+double C0dipole,C0dipole_refl; // inherent cross sections of exciting dipole (in free space and addition due to surface)
+
+// 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);
+ * 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
+ *
+ * !!! 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
+static double p0; // amplitude of the incident dipole moment
+static int n0; // Bessel beam order
+static double alpha0; // half-cone angle
+/* 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.
+ * Add descriptive comments, use 'static'.
+ */
+
+//======================================================================================================================
+
+void InitBeam(void)
+// initialize beam; produce description string
+{
+ double w0; // beam width
+ /* 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
+ switch (beamtype) {
+ case B_PLANE:
+ 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
+ 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
+ * are discussed in CalcFieldSurf() in crosssec.c.
+ */
+ if (cabs(msub-1)0) { // beam comes from the substrate (below)
+ // determine amplitude of the reflected and transmitted waves; here msub is always defined
+ 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),eIncRefl,eIncRefl);
+ cvMultScal_cmplx(tc*cexp(I*WaveNum*(kt-ki)*hsub),eIncTran,eIncTran);
+ // main part
+ for (i=0;i.
+ */
+#include
+#include
+#include
+
+double BESSJ0 (double X) {
+ const double
+ P1=1.0, P2=-0.1098628627E-2, P3=0.2734510407E-4,
+ P4=-0.2073370639E-5, P5= 0.2093887211E-6,
+ Q1=-0.1562499995E-1, Q2= 0.1430488765E-3, Q3=-0.6911147651E-5,
+ Q4= 0.7621095161E-6, Q5=-0.9349451520E-7,
+ R1= 57568490574.0, R2=-13362590354.0, R3=651619640.7,
+ R4=-11214424.18, R5= 77392.33017, R6=-184.9052456,
+ S1= 57568490411.0, S2=1029532985.0, S3=9494680.718,
+ S4= 59272.64853, S5=267.8532712, S6=1.0;
+ double
+ AX,FR,FS,Z,FP,FQ,XX,Y, TMP;
+ if (X==0.0) return 1.0;
+ AX = fabs(X);
+ if (AX < 8.0) {
+ Y = X*X;
+ FR = R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))));
+ FS = S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))));
+ TMP = FR/FS;
+ }
+ else {
+ Z = 8./AX;
+ Y = Z*Z;
+ XX = AX-0.785398164;
+ FP = P1+Y*(P2+Y*(P3+Y*(P4+Y*P5)));
+ FQ = Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)));
+ TMP = sqrt(0.636619772/AX)*(FP*cos(XX)-Z*FQ*sin(XX));
+ }
+return TMP;
+}
+
+double Sign(double X, double Y) {
+ if (Y<0.0) return (-fabs(X));
+ else return (fabs(X));
+}
+
+double BESSJ1 (double X) {
+ const double
+ P1=1.0, P2=0.183105E-2, P3=-0.3516396496E-4, P4=0.2457520174E-5,
+ P5=-0.240337019E-6, P6=0.636619772,
+ Q1= 0.04687499995, Q2=-0.2002690873E-3, Q3=0.8449199096E-5,
+ Q4=-0.88228987E-6, Q5= 0.105787412E-6,
+ R1= 72362614232.0, R2=-7895059235.0, R3=242396853.1,
+ R4=-2972611.439, R5=15704.48260, R6=-30.16036606,
+ S1=144725228442.0, S2=2300535178.0, S3=18583304.74,
+ S4=99447.43394, S5=376.9991397, S6=1.0;
+
+ double AX,FR,FS,Y,Z,FP,FQ,XX, TMP;
+
+ AX = fabs(X);
+ if (AX < 8.0) {
+ Y = X*X;
+ FR = R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))));
+ FS = S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))));
+ TMP = X*(FR/FS);
+ }
+ else {
+ Z = 8.0/AX;
+ Y = Z*Z;
+ XX = AX-2.35619491;
+ FP = P1+Y*(P2+Y*(P3+Y*(P4+Y*P5)));
+ FQ = Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)));
+ TMP = sqrt(P6/AX)*(cos(XX)*FP-Z*sin(XX)*FQ)*Sign(S6,X);
+ }
+ return TMP;
+}
+
+void BESSJCS (int N2, double X, double ARRJ[5]) {
+ const IACC = 40;
+ const double BIGNO = 1e10, BIGNI = 1e-10;
+
+ double TOX,BJM,BJ,BJP,SUM,TMP,BJM1,BJ1,BJP1,SUM1,TMP1;
+ int J, JSUM, JSUM1, M, K, N, sgn, sgn1;
+
+ N = (N2 - 2);
+ if (X == 0.0) {
+ switch(N2) {
+ case -2:
+ ARRJ[0] = 0.0;
+ ARRJ[1] = 0.0;
+ ARRJ[2] = 0.0;
+ ARRJ[3] = 0.0;
+ ARRJ[4] = 1.0;
+ break;
+ case -1:
+ ARRJ[0] = 0.0;
+ ARRJ[1] = 0.0;
+ ARRJ[2] = 0.0;
+ ARRJ[3] = 1.0;
+ ARRJ[4] = 0.0;
+ break;
+ case 0:
+ ARRJ[0] = 0.0;
+ ARRJ[1] = 0.0;
+ ARRJ[2] = 1.0;
+ ARRJ[3] = 0.0;
+ ARRJ[4] = 0.0;
+ break;
+ case 1:
+ ARRJ[0] = 0.0;
+ ARRJ[1] = 1.0;
+ ARRJ[2] = 0.0;
+ ARRJ[3] = 0.0;
+ ARRJ[4] = 0.0;
+ break;
+ case 2:
+ ARRJ[0] = 1.0;
+ ARRJ[1] = 0.0;
+ ARRJ[2] = 0.0;
+ ARRJ[3] = 0.0;
+ ARRJ[4] = 0.0;
+ break;
+ default:
+ ARRJ[0] = 0.0;
+ ARRJ[1] = 0.0;
+ ARRJ[2] = 0.0;
+ ARRJ[3] = 0.0;
+ ARRJ[4] = 0.0;
+ break;
+ }
+ }
+ else {
+ if ((N == -5)||(N == -4)||(N == -3)||(N == -2)||(N == -1)||(N == 0)||(N == 1)) {
+ switch(N) {
+ case -5:
+ ARRJ[4] = -BESSJ1(X);
+ ARRJ[3] = -(2.0/X)*ARRJ[4]-BESSJ0(X);
+ ARRJ[2] = -(4.0/X)*ARRJ[3]-ARRJ[4];
+ ARRJ[1] = -(6.0/X)*ARRJ[2]-ARRJ[3];
+ ARRJ[0] = -(8.0/X)*ARRJ[1]-ARRJ[2];
+ break;
+ case -4:
+ ARRJ[4] = BESSJ0(X);
+ ARRJ[3] = -BESSJ1(X);
+ ARRJ[2] = -(2.0/X)*ARRJ[3]-ARRJ[4];
+ ARRJ[1] = -(4.0/X)*ARRJ[2]-ARRJ[3];
+ ARRJ[0] = -(6.0/X)*ARRJ[1]-ARRJ[2];
+ break;
+ case -3:
+ ARRJ[4] = BESSJ1(X);
+ ARRJ[3] = BESSJ0(X);
+ ARRJ[2] = -ARRJ[4];;
+ ARRJ[1] = -(2.0/X)*ARRJ[2]-ARRJ[3];
+ ARRJ[0] = -(4.0/X)*ARRJ[1]-ARRJ[2];
+ break;
+ case -2:
+ ARRJ[3] = BESSJ1(X);
+ ARRJ[2] = BESSJ0(X);
+ ARRJ[1] = -ARRJ[4];;
+ ARRJ[1] = -(2.0/X)*ARRJ[2]-ARRJ[3];
+ ARRJ[0] = -(4.0/X)*ARRJ[1]-ARRJ[2];
+ break;
+ case -1:
+ ARRJ[0] = -BESSJ1(X);
+ ARRJ[1] = BESSJ0(X);
+ ARRJ[2] = -ARRJ[0];
+ ARRJ[3] = (2.0/X)*ARRJ[2]-ARRJ[1];
+ ARRJ[4] = (4.0/X)*ARRJ[3]-ARRJ[2];
+ break;
+ case 0:
+ ARRJ[0] = BESSJ0(X);
+ ARRJ[1] = BESSJ1(X);
+ ARRJ[2] = (2.0/X)*ARRJ[1]-ARRJ[0];
+ ARRJ[3] = (4.0/X)*ARRJ[2]-ARRJ[1];
+ ARRJ[4] = (6.0/X)*ARRJ[3]-ARRJ[2];
+ break;
+ case 1:
+ ARRJ[0] = BESSJ1(X);
+ ARRJ[1] = (2.0/X)*ARRJ[0]-BESSJ0(X);
+ ARRJ[2] = (4.0/X)*ARRJ[1]-ARRJ[0];
+ ARRJ[3] = (6.0/X)*ARRJ[2]-ARRJ[1];
+ ARRJ[4] = (8.0/X)*ARRJ[3]-ARRJ[2];
+ break;
+ }
+ }
+ else {
+ TOX = 2.0/X;
+ sgn = 1;
+ sgn1 = 1;
+ N = abs(N2-2);
+ if ((N2-2 <0)&((N2-2)%2 != 0)) sgn = -1;
+ if ((N2-1 <0)&((N2-1)%2 != 0)) sgn1 = -1;
+ if (X > 1.0*N) {
+ BJM = BESSJ0(X);
+ BJ = BESSJ1(X);
+
+ for (J=1; J0; J--) {
+ BJM1 = J*TOX*BJ1-BJP1;
+ BJP1 = BJ1;
+ BJ1 = BJM1;
+ if (fabs(BJ1) > BIGNO) {
+ BJ1 = BJ1*BIGNI;
+ BJP1 = BJP1*BIGNI;
+ TMP1 = TMP1*BIGNI;
+ SUM1 = SUM1*BIGNI;
+ }
+ if (JSUM1 != 0) SUM1 += BJ1;
+ JSUM1 = 1-JSUM1;
+ if (J == N+1) TMP1 = BJP1;
+ if (J <= M) {
+ BJM = J*TOX*BJ-BJP;
+ BJP = BJ;
+ BJ = BJM;
+ if (fabs(BJ) > BIGNO) {
+ BJ = BJ*BIGNI;
+ BJP = BJP*BIGNI;
+ TMP = TMP*BIGNI;
+ SUM = SUM*BIGNI;
+ }
+ if (JSUM != 0) SUM += BJ;
+ JSUM = 1-JSUM;
+ if (J == N) TMP = BJP;
+ }
+ }
+
+ ARRJ[0] = sgn*TMP/(2.0*SUM-BJ);
+ ARRJ[1] = sgn1*TMP1/(2.0*SUM1-BJ1);
+ ARRJ[2] = (N+2)*(2.0/X)*ARRJ[1]-ARRJ[0];
+ ARRJ[3] = (N+3)*(2.0/X)*ARRJ[2]-ARRJ[1];
+ ARRJ[4] = (N+4)*(2.0/X)*ARRJ[3]-ARRJ[2];
+
+ }
+ }
+ }
+}
+
+void BESSJLP (int N, double X, double ARRJ[2]) {
+ const IACC = 40;
+ const double BIGNO = 1e10, BIGNI = 1e-10;
+
+ double TOX,BJM,BJ,BJP,SUM,TMP,BJM1,BJ1,BJP1,SUM1,TMP1;
+ int J, JSUM, JSUM1, M, K;
+ if (X == 0.0) {
+ if (N == 0) { ARRJ[0] = 1.0; ARRJ[1] = 0.0;}
+ else { ARRJ[0] = 0.0; ARRJ[1] = 0.0;}
+ }
+ else {
+ if ((N == 0)||(N == 1)) {
+ if (N == 0) { ARRJ[0] = BESSJ0(X); ARRJ[1] = BESSJ1(X);}
+ if (N == 1) { ARRJ[0] = BESSJ1(X); ARRJ[1] = (2.0/X)*ARRJ[0]-BESSJ0(X);}
+ }
+ else {
+ TOX = 2.0/X;
+ if (X > 1.0*N) {
+ BJM = BESSJ0(X);
+ BJ = BESSJ1(X);
+
+ for (J=1; J0; J--) {
+ BJM1 = J*TOX*BJ1-BJP1;
+ BJP1 = BJ1;
+ BJ1 = BJM1;
+ if (fabs(BJ1) > BIGNO) {
+ BJ1 = BJ1*BIGNI;
+ BJP1 = BJP1*BIGNI;
+ TMP1 = TMP1*BIGNI;
+ SUM1 = SUM1*BIGNI;
+ }
+ if (JSUM1 != 0) SUM1 += BJ1;
+ JSUM1 = 1-JSUM1;
+ if (J == N+1) TMP1 = BJP1;
+ if (J <= M) {
+ BJM = J*TOX*BJ-BJP;
+ BJP = BJ;
+ BJ = BJM;
+ if (fabs(BJ) > BIGNO) {
+ BJ = BJ*BIGNI;
+ BJP = BJP*BIGNI;
+ TMP = TMP*BIGNI;
+ SUM = SUM*BIGNI;
+ }
+ if (JSUM != 0) SUM += BJ;
+ JSUM = 1-JSUM;
+ if (J == N) TMP = BJP;
+ }
+ }
+
+ ARRJ[0] = TMP/(2.0*SUM-BJ);
+ ARRJ[1] = TMP1/(2.0*SUM1-BJ1);
+ }
+ }
+ }
+}
diff --git a/bess.h b/bess.h
new file mode 100644
index 00000000..71e6048d
--- /dev/null
+++ b/bess.h
@@ -0,0 +1,27 @@
+/* File: bess.h
+ * $Date:: $
+ * Descr:
+ *
+ * Copyright (C) 2006,2008-2013
+ * This file is part of ADDA.
+ *
+ * ADDA 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.
+ *
+ * ADDA is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License along with ADDA. If not, see
+ * .
+ */
+
+#include
+#include
+#include
+
+double BESSJ0 (double X);
+double Sign(double X, double Y);
+double BESSJ1 (double X);
+
+void BESSJCS (int N2, double X, double ARRJ[5]);
+void BESSJLP (int N, double X, double ARRJ[2]);
diff --git a/const.h b/const.h
new file mode 100644
index 00000000..c1fcc441
--- /dev/null
+++ b/const.h
@@ -0,0 +1,427 @@
+/* File: const.h
+ * $Date:: $
+ * Descr: all the constants used by ADDA code, including enum constants, also defines some useful macros
+ *
+ * Copyright (C) 2006-2014 ADDA contributors
+ * This file is part of ADDA.
+ *
+ * ADDA 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.
+ *
+ * ADDA is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License along with ADDA. If not, see
+ * .
+ */
+#ifndef __const_h
+#define __const_h
+
+// version number (string)
+#define ADDA_VERSION "1.4.0-alpha"
+
+/* 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
+ * and produce a strong warning, if it is not satisfied. The list of C99 features, used by ADDA, include (but may be not
+ * limited to): stdbool.h, snprintf, %z argument in printf, '//' comments, restricted pointers, variadic macros
+*/
+# if !defined(OVERRIDE_STDC_TEST) && (!defined(__STDC_VERSION__) || (__STDC_VERSION__ < 199901L))
+# error "Support for C99 standard (at least many of its parts) is strongly recommended for compilation. Otherwise \
+ the compilation will may fail or produce wrong results. If you still want to try, you may enable an override in \
+ the Makefile."
+#endif
+
+/* The following is to ensure that mingw64 with "-std=c99" will use c99-compliant printf-family functions. For some
+ * (philosophical) reasons mingw64 developers have not implemented this behavior as the default one. So we need to set
+ * it manually. This macro should be defined before any system includes, hence inclusion of "const.h" should be the
+ * first one in all sources. This is also convenient for testing c99 standard above. However, there is no simple way
+ * then to test for MinGW64 at this point, since, e.g., __MINGW64_VERSION_STR is defined by the system header (not by
+ * the compiler itself). So the code executes always. Not the most reliable way, but seems the only way to keep the
+ * following definition in one place.
+ */
+#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L)
+# define __USE_MINGW_ANSI_STDIO 1
+#endif
+
+// basic constants
+#define UNDEF -1 // should be used only for variables, which are naturally non-negative
+ // denotes that shape accepts filename arguments; used in definitions of options
+#define FNAME_ARG -2 // single filename
+#define FNAME_ARG_2 -3 // two filenames
+#define FNAME_ARG_1_2 -4 // 1 or 2 filenames
+ // macro to test for occurrence of one of FNAME_ARG
+#define IS_FNAME_ARG(A) (((A)==FNAME_ARG) || ((A)==FNAME_ARG_2) || ((A)==FNAME_ARG_1_2))
+
+// simple functions
+#define MIN(A,B) (((A) > (B)) ? (B) : (A))
+#define MAX(A,B) (((A) < (B)) ? (B) : (A))
+#define MAXIMIZE(A,B) {if ((A)<(B)) (A)=(B);}
+#define IS_EVEN(A) (((A)%2) == 0)
+#define DIV_CEILING(A,B) (((A)%(B)==0) ? (A)/(B) : ((A)/(B))+1 ) // valid only for nonnegative A and B
+#define LENGTH(A) ((int)(sizeof(A)/sizeof(A[0]))) // length of any array (converted to int)
+#define STRINGIFY(A) #A
+#define GREATER_EQ2(a1,a2,b1,b2) ( (a1)>(b1) || ( (a1)==(b1) && (a2)>=(b2) )) // a1.a2>=b1.b2
+
+// parallel definitions
+#ifdef ADDA_MPI
+# define PARALLEL
+#endif
+
+/* ringid of root processor. Using ADDA_ROOT!=0 should work, however it was not thoroughly tested. Hence do not change
+ * without necessity.
+ */
+#define ADDA_ROOT 0
+
+// math constants rounded for 32 decimals; C99 standard specifies that they are encoded as double
+#define PI 3.1415926535897932384626433832795
+#define TWO_PI 6.283185307179586476925286766559
+#define FOUR_PI 12.566370614359172953850573533118
+#define EIGHT_PI 25.132741228718345907701147066236
+#define FOUR_PI_OVER_THREE 4.1887902047863909846168578443727
+#define PI_OVER_TWO 1.5707963267948966192313216916398
+#define PI_OVER_FOUR 0.78539816339744830961566084581988
+#define PI_OVER_SIX 0.52359877559829887307710723054658
+#define INV_PI 0.31830988618379067153776752674503
+#define TWO_OVER_PI 0.63661977236758134307553505349006
+#define THREE_OVER_FOUR_PI 0.23873241463784300365332564505877
+#define SIX_OVER_PI 1.9098593171027440292266051604702
+#define ONE_THIRD 0.33333333333333333333333333333333
+#define PI_OVER_180 0.017453292519943295769236907684886
+#define INV_PI_180 57.295779513082320876798154814105
+#define SQRT_PI 1.7724538509055160272981674833411
+#define TWO_OVER_SQRT_PI 1.1283791670955125738961589031215
+#define SQRT2 1.4142135623730950488016887242097
+#define SQRT3 1.7320508075688772935274463415059
+#define SQRT1_2 0.70710678118654752440084436210485
+#define SQRT1_2PI 0.39894228040143267793994605993438
+#define SQRT2_9PI 0.26596152026762178529329737328959
+#define EULER 0.57721566490153286060651209008241
+#define FULL_ANGLE 360.0
+#define MICRO 1E-6
+
+// sets the maximum box size; otherwise 'position' should be changed
+#define BOX_MAX USHRT_MAX
+
+// sizes of some arrays
+#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
+
+// 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
+ * POSIX this constant may be much larger and have even less reliability. So the values below are conservatively high,
+ * but are not guaranteed to suffice. However, the functions in the code are designed to survive if this buffer won't
+ * suffice.
+ */
+#define MAX_DIRNAME 300 // maximum length of dirname; increase THIS if any errors appear
+#define MAX_FNAME_SH 100 // maximum length of filename (used for known names)
+#define MAX_TMP_FNAME_SH 15 // maximum length of names of temporary files (short)
+#define MAX_SYSTEM_CALL 10 // maximum string length of system call (itself)
+#define MAX_WORD 10 // maximum length of a short word
+#define MAX_LINE 100 // maximum length of a line
+#define BUF_LINE 300 // size of buffer for reading lines (longer lines are handled robustly)
+#define MAX_PARAGRAPH 600 // maximum length of a paragraph (few lines)
+
+// derived sizes
+ // maximum string to create directory
+#define MAX_DIRSYS (MAX_DIRNAME + MAX_SYSTEM_CALL)
+ // maximum length of filename (including directory name)
+#define MAX_FNAME (MAX_DIRNAME + MAX_FNAME_SH)
+ // maximum length of temporary filename (including directory name)
+#define MAX_TMP_FNAME (MAX_DIRNAME + MAX_TMP_FNAME_SH)
+ // maximum message that may include a filename (for PrintError)
+#define MAX_MESSAGE (MAX_FNAME + MAX_PARAGRAPH)
+ // maximum message that may include 2 filenames (for LogError)
+#define MAX_MESSAGE2 (2*MAX_FNAME + MAX_PARAGRAPH)
+
+// widths of terminal used for output
+#define DEF_TERM_WIDTH 80 // default
+#define MIN_TERM_WIDTH 20 // ADDA never takes value less than that from environmental variables
+
+// formats for outputs of float values
+#define EFORM "%.10E" // fixed width
+#define GFORM "%.10g" // variable width (showing significant digits)
+#define GFORMDEF "%g" // default output for non-precise values
+#define GFORM_DEBUG "%.2g" // for debug and error output
+#define CFORM "%.10g%+.10gi" // for complex numbers; may be defined in terms of GFORM
+ // derived formats; starting "" is to avoid redundant syntax errors in Eclipse
+#define GFORM3V "("GFORM","GFORM","GFORM")"
+#define GFORM3L ""GFORM" "GFORM" "GFORM
+#define GFORM6L ""GFORM" "GFORM" "GFORM" "GFORM" "GFORM" "GFORM
+#define GFORM7L ""GFORM" "GFORM" "GFORM" "GFORM" "GFORM" "GFORM" "GFORM
+#define GFORM10L ""GFORM" "GFORM" "GFORM" "GFORM" "GFORM" "GFORM" "GFORM" "GFORM" "GFORM" "GFORM
+#define GFORMDEF3V "("GFORMDEF","GFORMDEF","GFORMDEF")"
+#define CFORM3V "("CFORM","CFORM","CFORM")"
+ // macros to shorten printing of all vector components
+#define COMP3V(a) (a)[0],(a)[1],(a)[2]
+#define COMP16V(a) (a)[0],(a)[1],(a)[2],(a)[3],(a)[4],(a)[5],(a)[6],(a)[7],(a)[8],(a)[9],(a)[10],(a)[11],(a)[12],\
+ (a)[13],(a)[14],(a)[15]
+
+enum sh { // shape types
+ SH_AXISYMMETRIC, // axisymmetric
+ SH_BICOATED, // two coated spheres
+ SH_BIELLIPSOID, // two general ellipsoids
+ SH_BISPHERE, // two spheres
+ SH_BOX, // box (may be rectangular)
+ SH_CAPSULE, // capsule
+ SH_CHEBYSHEV, // Chebyshev particle (axisymmetric)
+ SH_COATED, // coated sphere
+ SH_CYLINDER, // cylinder
+ SH_EGG, // egg
+ SH_ELLIPSOID, // general ellipsoid
+ SH_LINE, // line with width of one dipole
+ SH_PLATE, // plate
+ SH_PRISM, // right rectangular prism
+ SH_RBC, // Red Blood Cell
+ SH_READ, // read from file
+ //SH_SDISK_ROT, // disc cut of a sphere -- not operational
+ SH_SPHERE, // sphere
+ SH_SPHEREBOX // sphere in a box
+ /* TO ADD NEW SHAPE
+ * add an identifier starting with 'SH_' and a descriptive comment to this list in alphabetical order.
+ */
+};
+
+enum pol { // which way to calculate coupleconstant
+ POL_CLDR, // Corrected Lattice Dispersion Relation
+ POL_CM, // Clausius-Mossotti
+ POL_DGF, // Digitized Green's Function (second order approximation of LAK)
+ POL_FCD, // Filtered Coupled Dipoles
+ POL_IGT_SO, // Second order approximation to Green's tensor integrated over a cube
+ POL_LAK, // Exact result of IGT for sphere
+ POL_LDR, // Lattice Dispersion Relation
+ POL_NLOC, // non-local extension (Gaussian dipole-density, formula based on lattice sums)
+ POL_NLOC_AV, // same as NLOC, but based on averaging of Gaussian over the dipole volume
+ POL_RRC, // Radiative Reaction correction
+ POL_SO // Second Order formulation
+ /* TO ADD NEW POLARIZABILITY FORMULATION
+ * add an identifier starting with 'POL_' and a descriptive comment to this list in the alphabetical order.
+ */
+};
+// in alphabetical order
+
+enum scat { // how to calculate scattering quantities
+ SQ_DRAINE, // classical, as Draine
+ SQ_FINDIP, /* Same as Draine, but with correction of radiation energy of a _finite_ dipole when calculating
+ absorption cross section */
+ SQ_IGT_SO, // Integration of Green's tensor (second order in kd approximation)
+ SQ_SO // Second Order formulation
+};
+// in alphabetical order
+
+enum inter { // how to calculate interaction term
+ G_FCD, // Filtered Green's tensor (Filtered Coupled Dipoles)
+ G_FCD_ST, // quasi-static version of FCD
+ G_IGT, // (direct) integration of Green's tensor
+ G_IGT_SO, // approximate integration of Green's tensor (based on ideas of SO)
+ G_NLOC, // non-local extension (interaction of Gaussian dipole-densities)
+ G_NLOC_AV, // same as NLOC, but based on averaging of Gaussian over the dipole volume
+ G_POINT_DIP, // as point dipoles
+ G_SO // Second Order formulation
+ /* TO ADD NEW INTERACTION FORMULATION
+ * add an identifier starting with 'G_' and a descriptive comment to this list in the alphabetical order.
+ */
+};
+enum refl { // how to calculate interaction of dipoles through the nearby surface (reflected G)
+ GR_IMG, // approximate expression based on a single image dipole
+ GR_SOM // direct evaluation of Sommerfeld integrals
+ /* TO ADD NEW REFLECTION FORMULATION
+ * add an identifier starting with 'GR_' and a descriptive comment to this list in the alphabetical order.
+ */
+};// in alphabetical order
+
+// ldr constants
+#define LDR_B1 1.8915316
+#define LDR_B2 -0.1648469
+#define LDR_B3 1.7700004
+
+// 2nd_order constants; derived from c1=ln(5+3^(3/2))-ln(2)/2-pi/4 and c2=pi/6
+#define SO_B1 1.5867182426530356710958782335228 // 4c1/3
+#define SO_B2 0.13488017286410948123541594310740 // c1/3 - c2/2
+#define SO_B3 0.11895825700597042937085940122438 // (5/2)c1 - c2
+
+// other constants for polarizability
+#define DGF_B1 1.6119919540164696407169668466385 // (4pi/3)^(1/3)
+#define LAK_C 0.62035049089940001666800681204778 // (4pi/3)^(-1/3)
+
+// two boundaries for separation between G_SO 'close', 'median', and 'far'
+#define G_BOUND_CLOSE 1 // k*R^2/d < GB_CLOSE => 'close'
+#define G_BOUND_MEDIAN 1 // k*R < GB_MEDIAN => 'median'
+
+enum iter { // iterative methods
+ IT_BCGS2, // Enhanced Bi-Conjugate Gradient Stabilized (2)
+ IT_BICG_CS, // Bi-Conjugate Gradient for Complex-Symmetric matrices
+ IT_BICGSTAB, // Bi-Conjugate Gradient Stabilized
+ IT_CGNR, // Conjugate Gradient for Normalized equations minimizing Residual norm
+ IT_CSYM, // Algorithm CSYM
+ IT_QMR_CS, // Quasi-minimal residual for Complex-Symmetric matrices
+ IT_QMR_CS_2 // 2-term QMR (better roundoff properties)
+ /* TO ADD NEW ITERATIVE SOLVER
+ * add an identifier starting with 'IT_' and a descriptive comment to this list in the alphabetical order.
+ */
+};
+
+enum Eftype { // type of E field calculation
+ CE_NORMAL, // normal
+ CE_PARPER // use symmetry to calculate both incident polarizations from one calculation of internal fields
+};
+
+enum incpol {
+ INCPOL_Y, // y-polarization, assumed to be used first
+ INCPOL_X // x-polarization
+};
+// path and size of tables
+#define TAB_PATH "tables/"
+#define TAB_FNAME(a) "t"#a"f.dat" // a is a number, e.g. TAB_FNAME(2) -> "t2f.dat"
+#define TAB_SIZE 142
+#define TAB_RMAX 10
+
+enum beam { // beam types
+ B_BARTON5, // 5th order description of the Gaussian beam
+ B_BESSELCS, // circularly symmetric Bessel beam derived using ASR
+ B_BESSELLP, // linearly polarized Bessel beam derived using Hertz vector potentials
+ 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
+ /* TO ADD NEW BEAM
+ * add an identifier starting with 'B_' and a descriptive comment to this list in alphabetical order.
+ */
+};
+
+enum scatgrid { // types of scattering grid
+ SG_GRID, // grid of angles
+ SG_PAIRS // set of independent pairs
+};
+enum angleset {// types of angles set
+ AS_RANGE, // range with uniformly spaced points
+ AS_VALUES // any set of values
+};
+
+// types of phi_integr (should be different one-bit numbers)
+#define PHI_UNITY 1 // just integrate
+#define PHI_COS2 2 // integrate with cos(2*phi)
+#define PHI_SIN2 4 // integrate with sin(2*phi)
+#define PHI_COS4 8 // integrate with cos(4*phi)
+#define PHI_SIN4 16 // integrate with sin(4*phi)
+
+enum sym { // ways to treat particle symmetries
+ SYM_AUTO, // automatic
+ SYM_NO, // do not take into account
+ SYM_ENF // enforce
+};
+
+enum chpoint { // types of checkpoint (to save)
+ CHP_NONE, // do not save checkpoint
+ CHP_NORMAL, // save checkpoint if not finished in time and exit
+ CHP_REGULAR, // save checkpoints in regular time intervals (until finished or halted)
+ CHP_ALWAYS /* save checkpoint if either simulation is finished or time elapsed and calculate all scattering
+ quantities */
+};
+
+enum init_field { // how to calculate initial field to be used in the iterative solver
+ IF_AUTO, // automatically choose from ZERO or INC (based on lower residual value)
+ IF_ZERO, // zero
+ IF_INC, // equal to incident field
+ IF_READ, // read from file
+ IF_WKB // from WKB approximation (incident field corrected for phase shift in the particle)
+};
+
+// return values for functions
+#define CHP_EXIT -2 // exit after saving checkpoint
+
+// default values; other are specified in InitVariables (param.c)
+#define DEF_GRID (16*jagged)
+#define MIN_AUTO_GRID 16 // minimum grid, when set from default dpl
+
+// numbers less than this value (compared to unity) are considered to be zero (approximately 10*DBL_EPSILON)
+#define ROUND_ERR 1E-15
+#define SQRT_RND_ERR 3E-8 // sqrt(ROUND_ERR)
+
+// output and input file and directory names (can only be changed at compile time)
+#define F_EXPCOUNT "ExpCount"
+#define F_EXPCOUNT_LCK F_EXPCOUNT ".lck"
+#define F_CS "CrossSec"
+#define F_FRP "RadForce"
+#define F_INTFLD "IntField"
+#define F_DIPPOL "DipPol"
+#define F_BEAM "IncBeam"
+#define F_GRANS "granules"
+ // suffixes
+#define F_XSUF "-X"
+#define F_YSUF "-Y"
+ // logs
+#define F_LOG "log"
+#define F_LOG_ERR "logerr.%d" // ringid as argument
+#define F_LOG_ORAVG "log_orient_avg"
+#define F_LOG_INT_CSCA "log_int_Csca"
+#define F_LOG_INT_ASYM "log_int_asym"
+ // log suffixes
+#define F_LOG_X "_x"
+#define F_LOG_Y "_y"
+#define F_LOG_Z "_z"
+ // Mueller files
+#define F_MUEL "mueller"
+#define F_MUEL_SG "mueller_scatgrid"
+#define F_MUEL_INT "mueller_integr"
+#define F_MUEL_C2 "mueller_integr_c2"
+#define F_MUEL_S2 "mueller_integr_s2"
+#define F_MUEL_C4 "mueller_integr_c4"
+#define F_MUEL_S4 "mueller_integr_s4"
+ // files for amplitude matrix
+#define F_AMPL "ampl"
+#define F_AMPL_SG "ampl_scatgrid"
+ // temporary files; used in printf with ringid as argument
+#define F_FRP_TMP "rf%d.tmp"
+#define F_BEAM_TMP "b%d.tmp"
+#define F_INTFLD_TMP "f%d.tmp"
+#define F_DIPPOL_TMP "p%d.tmp"
+#define F_GEOM_TMP "g%d.tmp"
+ // checkpoint files
+#define F_CHP_LOG "chp.log"
+#define F_CHP "chp.%d" // ringid as argument
+
+// default file and directory names; can be changed by command line options
+#define FD_ALLDIR_PARMS "alldir_params.dat"
+#define FD_AVG_PARMS "avg_params.dat"
+#define FD_SCAT_PARMS "scat_params.dat"
+#define FD_CHP_DIR "chpoint"
+
+/* number of components of D. Really, it can't be easily changed, but using constant instead of 6 adds more
+ * understanding for reader
+ */
+#define NDCOMP 6
+
+// shape formats; numbers should be nonnegative
+enum shform {
+ SF_DDSCAT6, // DDSCAT 6 format (FRMFIL), produced by calltarget
+ SF_DDSCAT7, // DDSCAT 7 format (FRMFIL), produced by calltarget
+ SF_TEXT, // ADDA text format for one-domain particles
+ SF_TEXT_EXT // ADDA text format for multi-domain particles
+ /* TO ADD NEW FORMAT OF SHAPE FILE
+ * add an identifier starting with 'SF_' and a descriptive comment to this list in alphabetical order.
+ */
+};
+
+#define POSIT __FILE__,__LINE__ // position of the error in source code
+
+enum enwho { // who is calling
+ ALL, // each processor may report this error
+ ONE // only root processor reports an error
+};
+
+// derived; for simplicity
+#define ALL_POS ALL,POSIT
+#define ONE_POS ONE,POSIT
+#define ALL_POS_FUNC ALL_POS,__func__
+#define ONE_POS_FUNC ONE_POS,__func__
+
+enum ec { // error codes
+ EC_OK, // no error, corresponds to zero exit code
+ EC_ERROR, // error
+ EC_WARN, // warning
+ EC_INFO // slight warning, that does not interfere at all with normal execution
+};
+
+#endif // __const_h
diff --git a/param.c b/param.c
new file mode 100644
index 00000000..68903f0d
--- /dev/null
+++ b/param.c
@@ -0,0 +1,2568 @@
+/* File: param.c
+ * $Date:: $
+ * Descr: initialization, parsing and handling of input parameters; also printout general information; contains file
+ * locking routines
+ *
+ * Copyright (C) 2006-2014 ADDA contributors
+ * This file is part of ADDA.
+ *
+ * ADDA 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.
+ *
+ * ADDA is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License along with ADDA. If not, see
+ * .
+ */
+#include "const.h" // keep this first
+#include "param.h" // corresponding header
+// project headers
+#include "cmplx.h"
+#include "comm.h"
+#include "crosssec.h"
+#include "debug.h"
+#include "fft.h"
+#include "function.h"
+#include "io.h"
+#include "oclcore.h"
+#include "os.h"
+#include "parbas.h"
+#include "vars.h"
+// system headers
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#ifdef CLFFT_AMD
+/* One can also include clFFT.h (the only recommended public header), which can be redundant, but more portable.
+ * However, version macros are not documented anyway (in the manual), so there seem to be no perfectly portable way to
+ * obtain them.
+ */
+# include
+#endif
+
+#ifdef OCL_BLAS
+/* In contrast to clFFT, here including main header (clBLAS.h) is not an option, since it doesn't include the
+ * following header.
+ */
+# include
+#endif
+
+#ifndef NO_GITHASH
+# include "githash.h" // for GITHASH, this file is automatically created during compilation
+#endif
+
+// definitions for file locking
+#ifdef USE_LOCK
+# ifdef WINDOWS
+# define FILEHANDLE HANDLE
+# elif defined(POSIX)
+# include
+# include
+# ifdef LOCK_FOR_NFS
+# include // for error handling of fcntl call
+# endif
+# define FILEHANDLE int
+# else
+# error "Unknown operation system. Creation of lock files is not supported."
+# endif
+# define LOCK_WAIT 1 // in seconds
+# define MAX_LOCK_WAIT_CYCLES 60
+# define ONLY_FOR_LOCK
+#else
+# define FILEHANDLE int
+# define ONLY_FOR_LOCK ATT_UNUSED
+#endif
+
+// GLOBAL VARIABLES
+
+opt_index opt; // main option index; it is also defined as extern in param.h
+
+// SEMI-GLOBAL VARIABLES
+
+// defined and initialized in crosssec.c
+extern const char avg_string[];
+// defined and initialized in GenerateB.c
+extern const char beam_descr[];
+// defined and initialized in make_particle.c
+extern const bool volcor_used;
+extern const char *sh_form_str1,*sh_form_str2;
+#ifndef SPARSE
+extern const int gr_N;
+extern const double gr_vf_real;
+#endif //SPARSE
+extern const size_t mat_count[];
+
+// used in CalculateE.c
+bool store_int_field; // save full internal fields to text file
+bool store_dip_pol; // save dipole polarizations to text file
+bool store_beam; // save incident beam to file
+bool store_scat_grid; // Store the scattered field for grid of angles
+bool calc_Cext; // Calculate the extinction cross-section - always do
+bool calc_Cabs; // Calculate the absorption cross-section - always do
+bool calc_Csca; // Calculate the scattering cross-section by integration
+bool calc_vec; // Calculate the unnormalized asymmetry-parameter
+bool calc_asym; // Calculate the asymmetry-parameter
+bool calc_mat_force; // Calculate the scattering force by matrix-evaluation
+bool store_force; // Write radiation pressure per dipole to file
+bool store_ampl; // Write amplitude matrix to file
+int phi_int_type; // type of phi integration (each bit determines whether to calculate with different multipliers)
+// used in calculator.c
+bool avg_inc_pol; // whether to average CC over incident polarization
+double polNlocRp; // Gaussian width for non-local polarizability
+const char *alldir_parms; // name of file with alldir parameters
+const char *scat_grid_parms; // name of file with parameters of scattering grid
+// used in crosssec.c
+double incPolX_0[3],incPolY_0[3]; // initial incident polarizations (in lab RF)
+enum scat ScatRelation; // type of formulae for scattering quantities
+// used in GenerateB.c
+int beam_Npars;
+double beam_pars[MAX_N_BEAM_PARMS]; // beam parameters
+opt_index opt_beam; // option index of beam option used
+const char *beam_fnameY; // names of files, defining the beam (for two polarizations)
+const char *beam_fnameX;
+// used in interaction.c
+double igt_lim; // limit (threshold) for integration in IGT
+double igt_eps; // relative error of integration in IGT
+double nloc_Rp; // Gaussian width for non-local interaction
+bool InteractionRealArgs; // whether interaction (or reflection) routines can be called with real arguments
+// used in io.c
+char logfname[MAX_FNAME]=""; // name of logfile
+// used in iterative.c
+double iter_eps; // relative error to reach
+enum init_field InitField; // how to calculate initial field for the iterative solver
+const char *infi_fnameY; // names of files, defining the initial field (for two polarizations)
+const char *infi_fnameX;
+bool recalc_resid; // whether to recalculate residual at the end of iterative solver
+enum chpoint chp_type; // type of checkpoint (to save)
+time_t chp_time; // time of checkpoint (in sec)
+char const *chp_dir; // directory name to save/load checkpoint
+// used in make_particle.c
+enum sh shape; // particle shape definition
+int sh_Npars; // number of shape parameters
+double sh_pars[MAX_N_SH_PARMS]; // storage for shape parameters
+double sizeX; // size of particle along x-axis
+double dpl; // number of dipoles per lambda (wavelength)
+double lambda; // incident wavelength (in vacuum)
+int jagged; // size of big dipoles, used to construct a particle
+const char *shape_fname; // name of file, defining the shape
+const char *save_geom_fname; // geometry file name to save dipole configuration
+const char *shapename; // name of the used shape
+bool volcor; // whether to use volume correction
+bool save_geom; // whether to save dipole configuration in .geom file
+opt_index opt_sh; // option index of shape option used
+double gr_vf; // granules volume fraction
+double gr_d; // granules diameter
+int gr_mat; // domain number to granulate
+double a_eq; // volume-equivalent radius of the particle
+enum shform sg_format; // format for saving geometry files
+bool store_grans; // whether to save granule positions to file
+
+// LOCAL VARIABLES
+
+#define GFORM_RI_DIRNAME "%.4g" // format for refractive index in directory name
+
+static const char *run_name; // first part of the dir name ('run' or 'test')
+static const char *avg_parms; // name of file with orientation averaging parameters
+static const char *exename; // name of executable (adda, adda.exe, adda_mpi,...)
+static int Nmat_given; // number of refractive indices given in the command line
+static enum sym sym_type; // how to treat particle symmetries
+/* The following '..._used' flags are, in principle, redundant, since the structure 'options' contains the same flags.
+ * However, the latter can't be easily addressed by the option name (a search over the whole options is required).
+ * When thinking about adding a new one, first consider using UNDEF machinery instead
+ */
+static bool prop_used; // whether '-prop ...' was used in the command line
+static bool orient_used; // whether '-orient ...' was used in the command line
+static bool yz_used; // whether '-yz ...' was used in the command line
+static bool scat_plane_used; // whether '-scat_plane ...' was used in the command line
+
+/* TO ADD NEW COMMAND LINE OPTION
+ * If you need new variables or flags to implement effect of the new command line option, define them here. If a
+ * variable is used only in this source file, define it as local (static). If it is used in one another file, define
+ * them as semi-global, i.e. define them here under corresponding 'used in ...' line and put 'extern' declaration in
+ * corresponding source file under 'defined and initialized in param.c' line. If it is used in this and two or more
+ * files, define it as global, by putting definition in vars.c and 'extern' declaration in vars.h.
+ */
+
+// structure definitions
+struct subopt_struct {
+ const char *name; // name of option
+ const char *usage; // how to use (argument list)
+ const char *help; // help string
+ const int narg; /* possible number of arguments; UNDEF -> should not be checked; may contain also some special
+ negative codes, like FNAME_ARG. Currently it is assumed that all arguments are float (if no
+ special argument is specified), but it can be changed if will become a limitation. */
+ const int type; // type of suboption
+};
+struct opt_struct {
+ const char *name; // name of option
+ void (*func)(int Narg,char **argv); // pointer to a function, parsing this parameter
+ bool used; // flag to indicate, if the option was already used
+ const char *usage; // how to use (argument list)
+ const char *help; // help string
+ const int narg; // possible number of arguments; UNDEF -> should not be checked
+ const struct subopt_struct *sub; // suboptions
+};
+// const string for usage of ADDA
+static const char exeusage[]="[- [] [-]...]]";
+
+/* initializations of suboptions (comments to elements of subopt_struct are above); Contrary to 'options', suboptions
+ * are defined as null-terminated array, because they may be referenced not directly by their names but also as
+ * options[i].sub. In the latter case macro LENGTH can't be used to estimate the length of the array. So getting this
+ * length is at least nontrivial (e.g. can be done with some intermediate variables). So using NULL-termination seems
+ * to be the easiest.
+ */
+static const struct subopt_struct beam_opt[]={
+ {"barton5"," []","5th order approximation of the Gaussian beam (by Barton). The beam width is "
+ "obligatory and x, y, z coordinates of the center of the beam (in laboratory reference frame) are optional "
+ "(zero, by default). All arguments are in um. This is recommended option for simulation of the Gaussian beam.",
+ UNDEF,B_BARTON5},
+ {"besselCS","[]","Circularly symmetric Bessel beam of integer order derived "
+ "using ASR method. The tilt angle is measured from the z axis. Arguments x, y, z are coordinates of the "
+ "center of the beam (in laboratory reference frame). Coordinate arguments are in um. All arguments are "
+ "optional (zero, by default).",UNDEF,B_BESSELCS},
+ {"besselLP","[]","Linearly polarized Bessel beam of integer order derived using "
+ "Hertz vector potentials (m- type). The tilt angle is measured from the z axis. Arguments x, y, z are coordinates of "
+ "the center of the beam (in laboratory reference frame). Coordinate arguments are in um. All arguments are "
+ "optional (zero, by default).",UNDEF,B_BESSELLP},
+ {"davis3"," []","3rd order approximation of the Gaussian beam (by Davis). The beam width is "
+ "obligatory and x, y, z coordinates of the center of the beam (in laboratory reference frame) are optional "
+ "(zero, by default). All arguments are in um.",UNDEF,B_DAVIS3},
+ {"dipole","","Field of a unit point dipole placed at x, y, z coordinates (in laboratory reference "
+ "frame). All arguments are in um. Orientation of the dipole is determined by -prop command line option."
+ "Implies '-scat_matr none'. If '-surf' is used, dipole position should be above the surface.",3,B_DIPOLE},
+ {"lminus"," []","Simplest approximation of the Gaussian beam. The beam width is obligatory and "
+ "x, y, z coordinates of the center of the beam (in laboratory reference frame) are optional (zero, by"
+ " default). All arguments are in um.",UNDEF,B_LMINUS},
+ {"plane","","Infinite plane wave",0,B_PLANE},
+ {"read"," []","Defined by separate files, which names are given as arguments. Normally two "
+ "files are required for Y- and X-polarizations respectively, but a single filename is sufficient if only "
+ "Y-polarization is used (e.g. due to symmetry). Incident field should be specified in a particle reference "
+ "frame in the same format as used by '-store_beam'.",FNAME_ARG_1_2,B_READ},
+ /* TO ADD NEW BEAM
+ * add a row to this list in alphabetical order. It contains: beam name (used in command line), usage string, help
+ * string, possible number of float parameters, beam identifier (defined inside 'enum beam' in const.h). Usage and
+ * help string are shown either when -h option is invoked or error is found in input command line. Usage string is
+ * one line giving a list of possible arguments or argument combinations. Do not include beam name in it, use
+ * to denote argument, [...] for optional arguments, and {...|...|...} for multiple options of an
+ * argument. Help string should contain general description of the beam type and its arguments. Instead of number
+ * of parameters UNDEF can be used (if beam can accept variable number of parameters, then check it explicitly in
+ * function PARSE_FUNC(beam) below) or one of FNAME_ARG family (see explanation in const.h) if beam accepts a
+ * filenames as arguments. Number of parameters should not be greater than MAX_N_BEAM_PARMS (defined in const.h).
+ */
+ {NULL,NULL,NULL,0,0}
+};
+static const struct subopt_struct shape_opt[]={
+#ifndef SPARSE
+ {"axisymmetric","","Axisymmetric homogeneous shape, defined by its contour in ro-z plane of the "
+ "cylindrical coordinate system. Its symmetry axis coincides with the z-axis, and the contour is read from "
+ "file.",FNAME_ARG,SH_AXISYMMETRIC},
+ {"bicoated","","Two identical concentric coated spheres with outer diameter d (first domain), "
+ "inner diameter d_in, and center-to-center distance R_cc (along the z-axis). It describes both separate and "
+ "sintered coated spheres. In the latter case sintering is considered symmetrically for cores and shells.",
+ 2,SH_BICOATED},
+ {"biellipsoid","","Two general ellipsoids in default orientations with "
+ "centers on the z-axis, touching each other. Their semi-axes are x1,y1,z1 (lower one, first domain) and "
+ "x2,y2,z2 (upper one, second domain).",5,SH_BIELLIPSOID},
+ {"bisphere"," ","Two identical spheres with diameter d and center-to-center distance R_cc (along the "
+ "z-axis). It describe both separate and sintered spheres.",1,SH_BISPHERE},
+ {"box","[]","Homogeneous cube (if no arguments are given) or a rectangular parallelepiped with edges "
+ "x,y,z.",UNDEF,SH_BOX},
+ {"capsule","","Homogeneous capsule (cylinder with half-spherical end caps) with cylinder height h and "
+ "diameter d (its axis of symmetry coincides with the z-axis).",1,SH_CAPSULE},
+ {"chebyshev","","Axisymmetric Chebyshev particle of amplitude eps and order n, r=r_0[1+eps*cos(n*theta)]. "
+ "eps is a real number, such that |eps|<=1, while n is a natural number",2,SH_CHEBYSHEV},
+ {"coated"," []","Sphere with a spherical inclusion; outer sphere has a diameter d (first "
+ "domain). The included sphere has a diameter d_in (optional position of the center: x,y,z).",UNDEF,SH_COATED},
+ {"cylinder","","Homogeneous cylinder with height (length) h and diameter d (its axis of symmetry coincides "
+ "with the z-axis).",1,SH_CYLINDER},
+ {"egg","","Axisymmetric egg shape given by a^2=r^2+nu*r*z-(1-eps)z^2, where 'a' is scaling factor. "
+ "Parameters must satisfy 0","Homogeneous general ellipsoid with semi-axes x,y,z",2,SH_ELLIPSOID},
+ {"line","","Line along the x-axis with the width of one dipole",0,SH_LINE},
+ {"plate", "","Homogeneous plate (cylinder with rounded side) with cylinder height h and full diameter d (i.e. "
+ "diameter of the constituent cylinder is d-h). Its axis of symmetry coincides with the z-axis.",1,SH_PLATE},
+ {"prism","","Homogeneous right prism with height (length along the z-axis) h based on a regular polygon "
+ "with n sides of size 'a'. The polygon is oriented so that positive x-axis is a middle perpendicular for one "
+ "of its sides. Dx is size of the polygon along the x-axis, equal to 2Ri and Rc+Ri for even and odd n "
+ "respectively. Rc=a/[2sin(pi/n)] and Ri=Rc*cos(pi/n) are radii of circumscribed and inscribed circles "
+ "respectively.",2,SH_PRISM},
+ {"rbc","","Red Blood Cell, an axisymmetric (over z-axis) biconcave homogeneous particle, which is "
+ "characterized by diameter d, maximum and minimum width h, b, and diameter at the position of the maximum "
+ "width c. The surface is described by ro^4+2S*ro^2*z^2+z^4+P*ro^2+Q*z^2+R=0, ro^2=x^2+y^2, P,Q,R,S are "
+ "determined by the described parameters.",3,SH_RBC},
+#endif // !SPARSE
+ {"read","","Read a particle geometry from file ",FNAME_ARG,SH_READ},
+#ifndef SPARSE
+ {"sphere","","Homogeneous sphere",0,SH_SPHERE},
+ {"spherebox","","Sphere (diameter d_sph) in a cube (size Dx, first domain)",1,SH_SPHEREBOX},
+#endif // !SPARSE
+ /* TO ADD NEW SHAPE
+ * add a row to this list in alphabetical order. It contains: shape name (used in command line), usage string, help
+ * string, possible number of float parameters, shape identifier (defined inside 'enum sh' in const.h). Usage and
+ * help string are shown either when -h option is invoked or error is found in input command line. Usage string is
+ * one line giving a list of possible arguments or argument combinations. Do not include shape name in it, use
+ * to denote argument, [...] for optional arguments, and {...|...|...} for multiple options of an
+ * argument. Help string should contain general description of the shape and its arguments. Instead of number of
+ * parameters UNDEF can be used (if shape can accept variable number of parameters, then check it explicitly in
+ * function PARSE_FUNC(shape) below) or FNAME_ARG (if the shape accepts a single string argument with file name).
+ * Number of parameters should not be greater than MAX_N_SH_PARMS (defined in const.h). It is recommended to use
+ * dimensionless shape parameters, e.g. aspect ratios.
+ */
+ {NULL,NULL,NULL,0,0}
+};
+/* TO ADD NEW COMMAND LINE OPTION
+ * If a new option requires separate description of suboptions, add a static structure here (similar to already existing
+ * ones). It should contain a number of rows corresponding to different suboptions (see definition of subopt_struct
+ * above for details) and terminate with a NULL row.
+ */
+
+// EXTERNAL FUNCTIONS
+
+// GenerateB.c
+void InitBeam(void);
+
+//======================================================================================================================
+/* Prototypes of parsing functions; definitions are given below. defines are for conciseness. Since we use one common
+ * prototype style for all parsing functions and many of the latter do not actually use the passed parameters, these
+ * parameters have 'unused' attribute to eliminate spurious warnings.
+ */
+#define PARSE_NAME(a) parse_##a
+#define PARSE_FUNC(a) static void PARSE_NAME(a)(int Narg ATT_UNUSED,char **argv ATT_UNUSED)
+#define PAR(a) #a,PARSE_NAME(a),false
+PARSE_FUNC(alldir_inp);
+PARSE_FUNC(anisotr);
+PARSE_FUNC(asym);
+PARSE_FUNC(beam);
+PARSE_FUNC(chp_dir);
+PARSE_FUNC(chp_load);
+PARSE_FUNC(chp_type);
+PARSE_FUNC(chpoint);
+PARSE_FUNC(Cpr);
+PARSE_FUNC(Csca);
+PARSE_FUNC(dir);
+PARSE_FUNC(dpl);
+PARSE_FUNC(eps);
+PARSE_FUNC(eq_rad);
+#ifdef OPENCL
+PARSE_FUNC(gpu);
+#endif
+#ifndef SPARSE
+PARSE_FUNC(granul);
+#endif
+PARSE_FUNC(grid);
+PARSE_FUNC(h) ATT_NORETURN;
+PARSE_FUNC(init_field);
+PARSE_FUNC(int);
+PARSE_FUNC(int_surf);
+PARSE_FUNC(iter);
+PARSE_FUNC(jagged);
+PARSE_FUNC(lambda);
+PARSE_FUNC(m);
+PARSE_FUNC(maxiter);
+PARSE_FUNC(no_reduced_fft);
+PARSE_FUNC(no_vol_cor);
+PARSE_FUNC(ntheta);
+PARSE_FUNC(opt);
+PARSE_FUNC(orient);
+PARSE_FUNC(phi_integr);
+PARSE_FUNC(pol);
+PARSE_FUNC(prognosis);
+PARSE_FUNC(prop);
+PARSE_FUNC(recalc_resid);
+PARSE_FUNC(rect_dip);
+#ifndef SPARSE
+PARSE_FUNC(save_geom);
+#endif
+PARSE_FUNC(scat);
+PARSE_FUNC(scat_grid_inp);
+PARSE_FUNC(scat_matr);
+PARSE_FUNC(scat_plane);
+#ifndef SPARSE
+PARSE_FUNC(sg_format);
+#endif
+PARSE_FUNC(shape);
+PARSE_FUNC(size);
+PARSE_FUNC(store_beam);
+PARSE_FUNC(store_dip_pol);
+PARSE_FUNC(store_force);
+#ifndef SPARSE
+PARSE_FUNC(store_grans);
+#endif
+PARSE_FUNC(store_int_field);
+PARSE_FUNC(store_scat_grid);
+PARSE_FUNC(surf);
+PARSE_FUNC(sym);
+PARSE_FUNC(test);
+PARSE_FUNC(V) ATT_NORETURN;
+PARSE_FUNC(vec);
+PARSE_FUNC(yz);
+/* TO ADD NEW COMMAND LINE OPTION
+ * add a function prototype to this list. Add a line 'PARSE_FUNC(option_name);' in alphabetical order. It will be
+ * expanded automatically using defines specified above.
+ */
+
+static struct opt_struct options[]={
+ {PAR(alldir_inp),"","Specifies a file with parameters of the grid of scattering angles for calculating "
+ "integral scattering quantities.\n"
+ "Default: "FD_ALLDIR_PARMS,1,NULL},
+ {PAR(anisotr),"","Specifies that refractive index is anisotropic (its tensor is limited to be diagonal in particle "
+ "reference frame). '-m' then accepts 6 arguments per each domain. Can not be used with CLDR polarizability and "
+ "all SO formulations.",0,NULL},
+ {PAR(asym),"","Calculate the asymmetry vector. Implies '-Csca' and '-vec'",0,NULL},
+ {PAR(beam)," []","Sets the incident beam, either predefined or 'read' from file. All parameters of "
+ "predefined beam types (if present) are floats.\n"
+ "Default: plane",UNDEF,beam_opt},
+ {PAR(chp_dir),"","Sets directory for the checkpoint (both for saving and loading).\n"
+ "Default: "FD_CHP_DIR,1,NULL},
+ {PAR(chp_load),"","Restart a simulation from a checkpoint",0,NULL},
+ {PAR(chp_type),"{normal|regular|always}",
+ "Sets type of the checkpoint. All types, except 'always', require '-chpoint'.\n"
+ "Default: normal",1,NULL},
+ {PAR(chpoint),"