From a334ba60ee23526a7eb570d9eb405db7a3f7d9d7 Mon Sep 17 00:00:00 2001 From: yurkin Date: Tue, 20 Aug 2013 12:05:59 +0000 Subject: [PATCH] First step towards DDA with particles above surface. Currently, mostly 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/CalculateE.c | 218 ++++++++++++++++++++------- src/GenerateB.c | 110 +++++++++++++- src/calculator.c | 31 ++-- src/cmplx.h | 157 ++++++++++++++++++++ src/const.h | 9 +- src/crosssec.c | 115 +++++++++++---- src/crosssec.h | 2 + src/debug.c | 5 +- src/debug.h | 2 +- src/interaction.c | 7 +- src/make_particle.c | 12 +- src/param.c | 158 +++++++++++++++----- src/prec_time.c | 2 +- src/timing.c | 2 +- src/vars.c | 47 ++++-- src/vars.h | 15 +- tests/2exec/comp2exec | 16 +- tests/2exec/suite | 4 +- tests/2exec/suite_sparse | 1 + tests/2exec/suite_surf | 311 +++++++++++++++++++++++++++++++++++++++ 20 files changed, 1064 insertions(+), 160 deletions(-) create mode 100644 tests/2exec/suite_surf diff --git a/src/CalculateE.c b/src/CalculateE.c index 64d213f3..5ea1b94b 100644 --- a/src/CalculateE.c +++ b/src/CalculateE.c @@ -39,15 +39,16 @@ // defined and initialized in calculator.c extern double * restrict muel_phi,* restrict muel_phi_buf; -extern doublecomplex * restrict EplaneX, * restrict EplaneY; +extern doublecomplex * restrict EplaneX, * restrict EplaneY, * restrict EyzplX, * restrict EyzplY; extern const double dtheta_deg,dtheta_rad; extern doublecomplex * restrict ampl_alphaX,* restrict ampl_alphaY; extern double * restrict muel_alpha; // defined and initialized in crosssec.c extern const Parms_1D phi_sg; +extern const double ezLab[3],exSP[3]; // defined and initialized in param.c extern const bool store_int_field,store_dip_pol,store_beam,store_scat_grid,calc_Cext,calc_Cabs, -calc_Csca,calc_vec,calc_asym,calc_mat_force,store_force,store_mueller,store_ampl; + calc_Csca,calc_vec,calc_asym,calc_mat_force,store_force,store_ampl; extern const int phi_int_type; // defined and initialized in timing.c extern TIME_TYPE Timing_EPlane,Timing_EPlaneComm,Timing_IntField,Timing_IntFieldOne,Timing_ScatQuan,Timing_IncBeam; @@ -65,6 +66,8 @@ extern size_t TotalEFieldPlane; #define AMPL_FORMAT EFORM" "EFORM" "EFORM" "EFORM" "EFORM" "EFORM" "EFORM" "EFORM #define ANGLE_FORMAT "%.2f" #define RMSE_FORMAT "%.3E" +#define COMP44M(a) (a)[0][0],(a)[0][1],(a)[0][2],(a)[0][3],(a)[1][0],(a)[1][1],(a)[1][2],(a)[1][3],(a)[2][0],\ + (a)[2][1],(a)[2][2],(a)[2][3],(a)[3][0],(a)[3][1],(a)[3][2],(a)[3][3] // EXTERNAL FUNCTIONS @@ -138,9 +141,7 @@ static inline void PrintToIntegrFile(const int type,FILE * restrict file,double err=Romberg1D(phi_sg,16,muel_buf,matrix[0]); } if (err>*maxerr) *maxerr=err; - fprintf(file,ANGLE_FORMAT" "MUEL_FORMAT" "RMSE_FORMAT"\n",theta,matrix[0][0],matrix[0][1],matrix[0][2], - matrix[0][3],matrix[1][0],matrix[1][1],matrix[1][2],matrix[1][3],matrix[2][0],matrix[2][1],matrix[2][2], - matrix[2][3],matrix[3][0],matrix[3][1],matrix[3][2],matrix[3][3],err); + fprintf(file,ANGLE_FORMAT" "MUEL_FORMAT" "RMSE_FORMAT"\n",theta,COMP44M(matrix),err); } } @@ -191,11 +192,18 @@ void MuellerMatrix(void) si=sin(alph); for (i=0;i X - memcpy(incPolpar,incPolY,3*sizeof(double)); // par <=> Y + else { // special case of the above for alpha=0 + vCopy(incPolX,incPolper); // per <=> X + vCopy(incPolY,incPolpar); // par <=> Y } for(orient=0;orient -IncPolX; incPolX -> IncPolY + * plane. Then IncPolY (par) -> -IncPolX; incPolX (per) -> IncPolY */ if (which==INCPOL_Y) choice=INCPOL_X; else choice=INCPOL_Y; // which==INCPOL_X - incPol=incPolper; - incPolper=incPolpar; - incPolpar=incPol; - vInvSign(incPolpar); + vCopy(incPolper,incPol); + vCopy(incPolpar,incPolper); + vMultScal(-1,incPol,incPolpar); } // initialize Eplane if (orient_avg) { @@ -411,8 +439,8 @@ static void CalcEplane(const enum incpol which,const enum Eftype type) else Eplane=ampl_alphaX + 2*nTheta*k_or; // choice==INCPOL_X } else { - if (choice==INCPOL_Y) Eplane=EplaneY; - else Eplane=EplaneX; // choice==INCPOL_X + if (choice==INCPOL_Y) Eplane=EyzplY; + else Eplane=EyzplX; // choice==INCPOL_X } for (i=0;i -Y + vCopy(incPolX,incPolpar); // par <=> X + } + + for(orient=0;orient0) { // 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;i0) { // 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 #endif +// useful macro for printing complex numbers and matrices +#define REIM(a) creal(a),cimag(a) +#define REIM3V(a) REIM((a)[0]),REIM((a)[1]),REIM((a)[2]) + //====================================================================================================================== // operations on complex numbers @@ -41,6 +45,18 @@ static inline double cAbs2(const doublecomplex a) //====================================================================================================================== +static inline doublecomplex cSqrtCut(const doublecomplex a) +// square root of complex number, with explicit handling of branch cut (not to depend on sign of zero of imaginary part) +{ + if (cimag(a)==0) { + if (creal(a)>=0) return sqrt(a); + else return I*sqrt(-a); + } + else return csqrt(a); +} + +//====================================================================================================================== + static inline doublecomplex imExp(const double arg) // exponent of imaginary argument Exp(i*arg); optimization is performed by compiler // this may be faster than using generic cexp, since imaginary type is not supported by all compilers @@ -94,6 +110,48 @@ static inline void vConj(const doublecomplex a[static 3],doublecomplex b[static //====================================================================================================================== +static inline void vReal(const doublecomplex a[static 3],double b[static 3]) +// takes real part of the complex vector; b=Re(a) +{ + b[0]=creal(a[0]); + b[1]=creal(a[1]); + b[2]=creal(a[2]); +} + +//====================================================================================================================== + +static inline void cvBuildRe(const double a[static 3],doublecomplex b[static 3]) +// builds complex vector from real part; b=a + i*0 +{ + b[0]=a[0]; + b[1]=a[1]; + b[2]=a[2]; +} + +//====================================================================================================================== + +static inline void vInvRefl_cr(const double a[static 3],doublecomplex b[static 3]) +/* reflects real vector with respect to the xy-plane and then inverts it, equivalent to reflection about the z-axis + * result is stored into the complex vector + */ +{ + b[0]=-a[0]; + b[1]=-a[1]; + b[2]=a[2]; +} + +//====================================================================================================================== + +static inline void crCrossProd(const double a[static 3],const doublecomplex b[static 3],doublecomplex c[static 3]) +// cross product of real and complex vector; c = a x b +{ + c[0] = a[1]*b[2] - a[2]*b[1]; + c[1] = a[2]*b[0] - a[0]*b[2]; + c[2] = a[0]*b[1] - a[1]*b[0]; +} + +//====================================================================================================================== + static inline void cvMultScal(const double a,const doublecomplex b[static 3],doublecomplex c[static 3]) // multiplication of vector[3] by real scalar; c=ab; { @@ -241,6 +299,17 @@ static inline void cSymMatrVec(const doublecomplex matr[static 6],const doubleco //====================================================================================================================== // operations on real vectors +static inline void vCopy(const double a[static 3],double b[static 3]) +// copies one vector into another; b=a +{ + // can be rewritten through memcpy, but compiler should be able to do it itself if needed + b[0]=a[0]; + b[1]=a[1]; + b[2]=a[2]; +} + +//====================================================================================================================== + static inline void vAdd(const double a[static 3],const double b[static 3],double c[static 3]) // adds two real vectors; c=a+b { @@ -270,6 +339,16 @@ static inline void vInvSign(double a[static 3]) //====================================================================================================================== +static inline void vRefl(const double inc[static 3],double ref[static 3]) +// reflects the incident vector 'inc' with respect to the xy-plane (inverts z-component) +{ + ref[0]=inc[0]; + ref[1]=inc[1]; + ref[2]=-inc[2]; +} + +//====================================================================================================================== + static inline void vMultScal(const double a,const double b[static 3],double c[static 3]) // multiplication of real vector by scalar; c=a*b; { @@ -298,6 +377,28 @@ static inline double DotProd(const double a[static 3],const double b[static 3]) //====================================================================================================================== +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 +{ + c[0] = a[1]*b[2] - a[2]*b[1]; + c[1] = a[2]*b[0] - a[0]*b[2]; + c[2] = a[0]*b[1] - a[1]*b[0]; +} + +//====================================================================================================================== + +static inline void vNormalize(double a[static 3]) +// normalize real vector to have unit norm +{ + double c; + c=1/sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]); + a[0]*=c; + a[1]*=c; + a[2]*=c; +} + +//====================================================================================================================== + static inline void LinComb(const double a[static 3],const double b[static 3],const double c1,const double c2, double c[static 3]) // linear combination of real vectors[3]; c=c1*a+c2*b @@ -338,6 +439,7 @@ static inline double QuadForm(const double matr[static 6],const double vec[stati } //====================================================================================================================== + static inline void MatrVec(double matr[static 3][3],const double vec[static 3],double res[static 3]) // multiplication of matrix[3][3] by vec[3] (all real); res=matr.vec; { @@ -348,6 +450,16 @@ static inline void MatrVec(double matr[static 3][3],const double vec[static 3],d //====================================================================================================================== +static inline void MatrColumn(double matr[static 3][3],const int ind,double vec[static 3]) +// get ind's column of matrix[3][3] and store it into vec[3] (all real, ind starts from zero); vec=matr[.][ind]; +{ + vec[0]=matr[0][ind]; + vec[1]=matr[1][ind]; + vec[2]=matr[2][ind]; +} + +//====================================================================================================================== + static inline void Permutate(double vec[static 3],const int ord[static 3]) // permutate double vector vec using permutation ord { @@ -389,6 +501,51 @@ static inline double Rad2Deg(const double rad) return (INV_PI_180*rad); } +//====================================================================================================================== +// functions used for substrate + +//====================================================================================================================== + +static inline doublecomplex FresnelRS(const doublecomplex ki,const doublecomplex kt) +/* reflection coefficient for s-polarized wave (E perpendicular to the main plane), + * ki,kt are normal (positive) components of incident and transmitted wavevector (with arbitrary mutual scaling) + */ +{ + return (ki-kt)/(ki+kt); +} + +//====================================================================================================================== + +static inline doublecomplex FresnelTS(const doublecomplex ki,const doublecomplex kt) +/* transmission coefficient for s-polarized wave (E perpendicular to the main plane), + * ki,kt are normal (positive) components of incident and transmitted wavevector (with arbitrary mutual scaling) + */ +{ + return 2*ki/(ki+kt); +} + +//====================================================================================================================== + +static inline doublecomplex FresnelRP(const doublecomplex ki,const doublecomplex kt,const doublecomplex mr) +/* reflection coefficient for p-polarized wave (E parallel to the main plane), + * ki,kt are normal (positive) components of incident and transmitted wavevector (with arbitrary mutual scaling) + * mr is the ratio of refractive indices (mt/mi) + */ +{ + return (mr*mr*ki-kt)/(mr*mr*ki+kt); +} + +//====================================================================================================================== + +static inline doublecomplex FresnelTP(const doublecomplex ki,const doublecomplex kt,const doublecomplex mr) +/* transmission coefficient for p-polarized wave (E parallel to the main plane), + * ki,kt are normal (positive) components of incident and transmitted wavevector (with arbitrary mutual scaling) + * mr is the ratio of refractive indices (mt/mi) + */ +{ + return 2*mr*ki/(mr*mr*ki+kt); +} + #ifdef USE_SSE3 //====================================================================================================================== diff --git a/src/const.h b/src/const.h index 7b3f5097..953f586c 100644 --- a/src/const.h +++ b/src/const.h @@ -142,6 +142,10 @@ #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 @@ -198,7 +202,10 @@ enum inter { // how to calculate interaction term G_POINT_DIP, // as point dipoles G_SO // Second Order formulation }; -// in 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 +};// in alphabetical order // ldr constants #define LDR_B1 1.8915316 diff --git a/src/crosssec.c b/src/crosssec.c index 9fd6ccf5..e6faa48b 100644 --- a/src/crosssec.c +++ b/src/crosssec.c @@ -44,8 +44,9 @@ extern doublecomplex * restrict expsX,* restrict expsY,* restrict expsZ; #endif // defined and initialized in GenerateB.c extern const double beam_center_0[3]; +//extern doublecomplex eIncRefl[3],eIncTran[3]; // defined and initialized in param.c -extern const double prop_0[3],incPolX_0[3],incPolY_0[3]; +extern const double incPolX_0[3],incPolY_0[3]; extern const enum scat ScatRelation; // defined and initialized in timing.c extern TIME_TYPE Timing_EFieldAD,Timing_EFieldADComm,Timing_EFieldSG,Timing_EFieldSGComm, @@ -53,6 +54,8 @@ Timing_ScatQuanComm; // used in CalculateE.c Parms_1D phi_sg; +double ezLab[3]; // basis vector ez of laboratory RF transformed into the RF of particle +double exSP[3]; // second vector, which determines the scattering plane passing through prop and ez (in particle RF) // used in calculator.c Parms_1D parms_alpha; // parameters of integration over alpha Parms_1D parms[2]; // parameters for integration over theta,phi or beta,gamma @@ -64,8 +67,9 @@ bool full_al_range; // whether full range of alpha angle is used // LOCAL VARIABLES -double dCabs; // difference between Cabs calculated by 'dr' and 'fin' formulations -bool dCabs_ready=false; // whether dCabs is already calculated +static double dCabs; // difference between Cabs calculated by 'dr' and 'fin' formulations +static bool dCabs_ready=false; // whether dCabs is already calculated +static double exLab[3],eyLab[3]; // basis vectors of laboratory RF transformed into the RF of particle //====================================================================================================================== @@ -86,6 +90,7 @@ void InitRotation (void) double ca,sa,cb,sb,cg,sg; double beta_matr[3][3]; double alph,bet,gam; // in radians + double ex[3]; // second vector of scattering plane in the laboratory RF // initialization of angle values in radians alph=Deg2Rad(alph_deg); @@ -108,10 +113,28 @@ void InitRotation (void) beta_matr[2][0]=ca*sb; beta_matr[2][1]=sa*sb; beta_matr[2][2]=cb; - // rotation of incident field + // rotation of incident field and basis vectors of laboratory RF MatrVec(beta_matr,prop_0,prop); MatrVec(beta_matr,incPolY_0,incPolY); MatrVec(beta_matr,incPolX_0,incPolX); + MatrColumn(beta_matr,0,exLab); + MatrColumn(beta_matr,1,eyLab); + MatrColumn(beta_matr,2,ezLab); + /* define the second vector (x) of the scattering plane (through prop and ez) + * if prop is along ez, we choose the original ex (that is - we generally assume that beam propagates from the left + * side and positive x-direction is in the right side (left and right is with respect to the possible rotation of + * ex,ey over the z-axis + */ + if (prop_0[2]==1) vCopy(incPolX_0,ex); + else if (prop_0[2]==-1) vMultScal(-1,incPolX_0,ex); + else { + ex[0]=prop_0[0]; + ex[1]=prop_0[1]; + ex[2]=0; + vNormalize(ex); + } + MatrVec(beta_matr,ex,exSP); + // if needed rotate beam center if (beam_asym) MatrVec(beta_matr,beam_center_0,beam_center); } @@ -582,7 +605,7 @@ double ExtCross(const double * restrict incPol) double sum; size_t i; - if (beamtype==B_PLANE) { + if (beamtype==B_PLANE && !surface) { CalcField (ebuff,prop); sum=crDotProd_Re(ebuff,incPol); // incPol is real, so no conjugate is needed MyInnerProduct(&sum,double_type,1,&Timing_ScatQuanComm); @@ -703,13 +726,60 @@ double AbsCross(void) //====================================================================================================================== +void SetScatPlane(const double ct,const double st,const double phi,double robs[static restrict 3], + double polPer[static restrict 3]) +/* Given theta (cos,sin) and phi, calculates scattering direction and unit vector perpendicular to the scattering plane. + * Currently, two alternative definitions are used: theta and phi are either angles in the laboratory reference frame, + * or in the beam reference frame. Generally, polPer = robs x prop (normalized to unit vector), but cases when the two + * latter vectors are aligned requires special treatment. + * + * For special cases with prop||ez we assume that the result is continuous for theta changing from 0 to pi. In other + * words th=0 and pi is effective replaced by th=+0 and pi-0 respectively. Then phi determines the scattering plane. + * + * For other special cases (prop not aligned with ez) the scattering plane is taken to include ez (and hence IncPolX). + * + * Overall, we spend a lot of effort to make consistent definition of polPer. However, its sign is irrelevant because + * change of sign (simultaneously for polPer and both incident and scattering polPar) doesn't change the amplitude + * (and Mueller) matrix. + */ +{ + double cphi,sphi; + + cphi=cos(phi); + sphi=sin(phi); + if (surface) { // definitions are in the laboratory reference frame + // robs = cos(theta)*ezLab + sin(theta)*[cos(phi)*exLab + sin(phi)*eyLab]; + LinComb(exLab,eyLab,cphi,sphi,robs); + LinComb(ezLab,robs,ct,st,robs); + } + else { // definitions are with respect to the incident beam + // robs = cos(theta)*prop + sin(theta)*[cos(phi)*incPolX + sin(phi)*incPolY]; + LinComb(incPolX,incPolY,cphi,sphi,robs); + LinComb(prop,robs,ct,st,robs); + } + // set unit vectors which defines Eper + // two special cases; testing ct is more accurate than st==0 + // 1) robs has +0 component along exLab (for phi=0): incPolper = (+-)(sin(phi)*exLab - cos(phi)*eyLab) + if (surface && propAlongZ && fabs(ct)==1) LinComb(exLab,eyLab,sphi*propAlongZ,-cphi*propAlongZ,polPer); + // 2) robs has +0 component along incPolX (for phi=0): incPolper = sin(phi)*incPolX - cos(phi)*incPolY + else if (!surface && fabs(ct)==1) LinComb(incPolX,incPolY,sphi,-cphi,polPer); + else { // general case: incPolPer = robserver x prop /||...|| + CrossProd(robs,prop,polPer); + // when robs || prop scattering plane contains ez, then polPer = prop x ez/||...|| = -incPolY + // currently this may only occur for surface, since otherwise it is handled by a special case above + if (DotProd(polPer,polPer)==0) vMultScal(-1,incPolY,polPer); + else vNormalize(polPer); + } +} +//====================================================================================================================== + void CalcAlldir(void) // calculate scattered field in many directions { int index,npoints,point; size_t i,j; TIME_TYPE tstart; - double robserver[3],incPolpar[3],incPolper[3],cthet,sthet,cphi,sphi,th,ph; + double robserver[3],incPolpar[3],incPolper[3],cthet,sthet,th,ph; doublecomplex ebuff[3]; // Calculate field @@ -722,11 +792,13 @@ void CalcAlldir(void) sthet=sin(th); for (j=0;jDipoleCoord[i3+2]) minZco=DipoleCoord[i3+2]; // crude way to find the minimum on the way } - + /* test that particle is wholly above the substrate; strictly speaking, we test dipole centers to be above the + * substrate - hsub+minZco>0, while the geometric boundary of the particle may still intersect with the substrate. + * However, the current test is sufficient to ensure that corresponding routines to calculate reflected Green's + * tensor do not fail. And accuracy of the DDA itself is anyway questionable when some of the dipoles are very close + * to the substrate (whether they cross it or not). + */ + if (surface && hsub<=-minZco) LogError(ALL_POS,"The particle must be entirely above the substrate. There exist a" + "dipole with z="GFORMDEF" (relative to the center), making specified height of the center ("GFORMDEF") too " + "small",minZco,hsub); // save geometry if (save_geom) #ifndef SPARSE diff --git a/src/param.c b/src/param.c index f71b08a8..3e7b466d 100644 --- a/src/param.c +++ b/src/param.c @@ -101,7 +101,6 @@ 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_mueller; // Calculate and write Mueller matrix 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 @@ -109,7 +108,6 @@ bool avg_inc_pol; // whether to average CC over incident polariz 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 prop_0[3]; // initial incident direction (in laboratory reference frame) 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 @@ -162,6 +160,10 @@ static const char *avg_parms; // name of file with orientation averaging param 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 +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 @@ -330,6 +332,7 @@ 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); @@ -351,6 +354,7 @@ PARSE_FUNC(save_geom); PARSE_FUNC(scat); PARSE_FUNC(scat_grid_inp); PARSE_FUNC(scat_matr); +PARSE_FUNC(scat_plane); #ifndef SPARSE PARSE_FUNC(sg_format); #endif @@ -364,6 +368,7 @@ 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; @@ -461,6 +466,14 @@ static struct opt_struct options[]={ "!!! All options except 'poi' incur a significant slowing down in sparse mode.\n" #endif "Default: poi",UNDEF,NULL}, + {PAR(int_surf),"{img|som}", + "Sets prescription to calculate the interaction term.\n" + "'img' - approximation based on a single image dipole (fast but inaccurate).\n" + "'som' - direct evaluation of Sommerfeld integrals.\n" + #ifdef SPARSE + "!!! In sparse mode 'som' is expected to be very slow.\n" + #endif + "Default: som",1,NULL}, {PAR(iter),"{bcgs2|bicg|bicgstab|cgnr|csym|qmr|qmr2}","Sets the iterative solver.\n" "Default: qmr",1,NULL}, /* TO ADD NEW ITERATIVE SOLVER @@ -546,6 +559,9 @@ static struct opt_struct options[]={ {PAR(scat_matr),"{muel|ampl|both|none}","Specifies which scattering matrices (from Mueller and amplitude) should " "be saved to file. Amplitude matrix is never integrated (in combination with '-orient avg' or '-phi_integr').\n" "Default: muel",1,NULL}, + {PAR(scat_plane),"","Explicitly enables calculation of the scattering in the plane through e_z (in laboratory " + "reference frame) and propagation direction of incident wave. For default incidence, this is the xz-plane. It " + "can also be implicitly enabled by other options.",0,NULL}, #ifndef SPARSE {PAR(sg_format),"{text|text_ext|ddscat6|ddscat7}","Specifies format for saving geometry files. First two are ADDA " "default formats for single- and multi-domain particles respectively. 'text' is automatically changed to " @@ -573,6 +589,10 @@ static struct opt_struct options[]={ #endif {PAR(store_int_field),"","Save internal fields to a file",0,NULL}, {PAR(store_scat_grid),"","Calculate Mueller matrix for a grid of scattering angles and save it to a file.",0,NULL}, + {PAR(surf)," ","Specifies that scatterer is located above the plane surface, parallel to the " + "xy-plane. First two arguments specify the refractive index of the substrate (below the surface), assuming that " + "the vacuum is above the surface. specifies the height of particle center above the surface (along the " + "z-axis, in um). Particle must be entirely above the substrate.",3,NULL}, {PAR(sym),"{auto|no|enf}","Automatically determine particle symmetries ('auto'), do not take them into account " "('no'), or enforce them ('enf').\n" "Default: auto",1,NULL}, @@ -1107,6 +1127,12 @@ PARSE_FUNC(int) else NotSupported("Interaction term prescription",argv[1]); if (Narg>1 && strcmp(argv[1],"igt")!=0) PrintErrorHelp("Additional arguments are allowed only for 'igt'"); } +PARSE_FUNC(int_surf) +{ + if (strcmp(argv[1],"img")==0) ReflRelation=GR_IMG; + else if (strcmp(argv[1],"som")==0) ReflRelation=GR_SOM; + else NotSupported("Interaction term prescription",argv[1]); +} PARSE_FUNC(iter) { if (strcmp(argv[1],"bcgs2")==0) IterMethod=IT_BCGS2; @@ -1189,6 +1215,10 @@ PARSE_FUNC(orient) ScanDoubleError(argv[2],&bet_deg); ScanDoubleError(argv[3],&gam_deg); } + /* this will catch even trivial cases like '-orient 0 0 0', but we still want to point the user to existent + * incompatibilities. This may help to prevent misunderstanding on the user's side. + */ + orient_used=true; } PARSE_FUNC(phi_integr) { @@ -1231,10 +1261,8 @@ PARSE_FUNC(prop) ScanDoubleError(argv[3],prop_0+2); tmp=DotProd(prop_0,prop_0); if (tmp==0) PrintErrorHelp("Given propagation vector is null"); - tmp=1/sqrt(tmp); - prop_0[0]*=tmp; - prop_0[1]*=tmp; - prop_0[2]*=tmp; + vMultScal(1/sqrt(tmp),prop_0,prop_0); + prop_used=true; } PARSE_FUNC(recalc_resid) { @@ -1274,6 +1302,11 @@ PARSE_FUNC(scat_matr) else if (strcmp(argv[1],"none")==0) store_mueller=store_ampl=false; else NotSupported("Scattering matrix specifier",argv[1]); } +PARSE_FUNC(scat_plane) +{ + scat_plane_used = true; + scat_plane = true; +} #ifndef SPARSE PARSE_FUNC(sg_format) { @@ -1357,6 +1390,18 @@ PARSE_FUNC(store_scat_grid) { store_scat_grid = true; } +PARSE_FUNC(surf) +{ + double mre,mim; + + ScanDoubleError(argv[1],&mre); + ScanDoubleError(argv[2],&mim); + msub = mre + I*mim; + ScanDoubleError(argv[3],&hsub); + TestPositive(hsub,"height above surface"); + surface = true; + symZ=false; +} PARSE_FUNC(sym) { if (strcmp(argv[1],"auto")==0) sym_type=SYM_AUTO; @@ -1548,6 +1593,7 @@ PARSE_FUNC(vec) } PARSE_FUNC(yz) { + yz_used = true; yzplane = true; } /* TO ADD NEW COMMAND LINE OPTION @@ -1647,9 +1693,8 @@ static void RemoveLockFile(FILEHANDLE fd ONLY_FOR_LOCK,const char * restrict fna void InitVariables(void) // some defaults are specified also in const.h { - prop_0[0]=0; // by default beam propagates along z-axis - prop_0[1]=0; - prop_0[2]=1; + prop_used=false; + orient_used=false; directory=""; lambda=TWO_PI; // initialize ref_index of scatterer @@ -1657,7 +1702,6 @@ void InitVariables(void) ref_index[0]=1.5; // initialize to null to determine further whether it is initialized logfile=NULL; - boxX=boxY=boxZ=UNDEF; sizeX=UNDEF; a_eq=UNDEF; @@ -1692,6 +1736,9 @@ void InitVariables(void) save_geom=false; save_geom_fname=""; yzplane=false; + yz_used=false; + scat_plane=false; + scat_plane_used=false; all_dir=false; scat_grid=false; phi_integr=false; @@ -1719,6 +1766,8 @@ void InitVariables(void) igt_eps=UNDEF; InitField=IF_AUTO; recalc_resid=false; + surface=false; + ReflRelation=GR_SOM; // sometimes the following two are left uninitialized beam_fnameX=NULL; infi_fnameX=NULL; @@ -1788,6 +1837,11 @@ void VariablesInterconnect(void) // initialize WaveNum ASAP WaveNum = TWO_PI/lambda; + // set default incident direction, which is +z for all configurations + if (!prop_used) { + prop_0[0]=prop_0[1]=0; + prop_0[2] = 1; + } // parameter interconnections /* very unlikely that calc_Cabs will ever be false, but strictly speaking dCabs should be calculated before Cext, * when SQ_FINDIP is used @@ -1795,14 +1849,18 @@ void VariablesInterconnect(void) if (ScatRelation==SQ_FINDIP && calc_Cext) calc_Cabs=true; if (IntRelation==G_SO) reduced_FFT=false; if (calc_Csca || calc_vec) all_dir = true; - // yzplane is true except when this case is true and when not forced by -yz option + // by default, one of the scattering options is activated if (store_scat_grid || phi_integr) scat_grid = true; - else yzplane = true; + else if (!scat_plane_used && !yz_used) { // default values when no scattering plane is explicitly specified + if (surface) scat_plane = true; + else yzplane = true; + } // if not initialized before, IGT precision is set to that of the iterative solver if (igt_eps==UNDEF) igt_eps=iter_eps; // parameter incompatibilities + if (scat_plane && yzplane) PrintError("Currently '-scat_plane' and '-yz' cannot be used together."); if (orient_avg) { - if (prop_0[2]!=1) PrintError("'-prop' and '-orient avg' can not be used together"); + if (prop_used) PrintError("'-prop' and '-orient avg' can not be used together"); if (store_int_field) PrintError("'-store_int_field' and '-orient avg' can not be used together"); if (store_dip_pol) PrintError("'-store_dip_pol' and '-orient avg' can not be used together"); if (store_beam) PrintError("'-store_beam' and '-orient avg' can not be used together"); @@ -1817,6 +1875,8 @@ void VariablesInterconnect(void) LogWarning(EC_WARN,ONE_POS,"Amplitude matrix can not be averaged over orientations. So switching off " "calculation of Mueller matrix results in no calculation of scattering matrices at all."); } + if (scat_plane && yzplane) PrintError("Only one scattering plane can be used during orientation averaging, so " + "a combination of '-scat_plane' and '-yz' cannot be used together with '-orient avg'."); } if (phi_integr && !store_mueller) PrintError("Integration over phi can only be performed for Mueller matrix. " "Hence, '-phi_integr' is incompatible with '-scat_matr {ampl|none}'"); @@ -1825,8 +1885,13 @@ void VariablesInterconnect(void) "'-ntheta' irrelevant. Hence, these two options are incompatible."); if (store_scat_grid) PrintError("'-scat_matr none' turns off calculation of angle-resolved quantities. Hence, " "it is incompatible with '-store_scat_grid'."); + if (yz_used) PrintError("'-scat_matr none' turns off calculation of angle-resolved quantities. Hence, it is " + "incompatible with '-yz'."); + if (scat_plane_used) PrintError("'-scat_matr none' turns off calculation of angle-resolved quantities. Hence, " + "it is incompatible with '-scat_plane'."); nTheta=0; yzplane=false; + scat_plane=false; scat_grid=false; } if (anisotropy) { @@ -1846,6 +1911,18 @@ void VariablesInterconnect(void) if (sizeX!=UNDEF && a_eq!=UNDEF) PrintError("'-size' and '-eq_rad' can not be used together"); if (calc_mat_force && beamtype!=B_PLANE) PrintError("Currently radiation forces can not be calculated for non-plane incident wave"); + if (InitField==IF_WKB) { + if (prop_used) PrintError("Currently '-init_field wkb' and '-prop' can not be used together"); + if (beamtype!=B_PLANE) PrintError("'-init_field wkb' is incompatible with non-plane incident wave"); + } + if (surface) { // currently a lot of limitations for particles near surface + if (orient_used) PrintError("Currently '-orient' and '-surf' can not be used together"); + if (calc_mat_force) PrintError("Currently calculation of radiation forces is incompatible with '-surf'"); + if (beamtype!=B_PLANE) PrintError("Currently non-plane incident wave is incompatible with '-surf'"); + if (InitField==IF_WKB) PrintError("'-init_field wkb' and '-surf' can not be used together"); + 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"); + } #ifdef SPARSE if (shape==SH_SPHERE) PrintError("Sparse mode requires shape to be read from file (-shape read ...)"); #endif @@ -1866,12 +1943,13 @@ void VariablesInterconnect(void) /* Determine two incident polarizations. Equivalent to rotation of X,Y,Z basis by angles theta and phi from (0,0,1) * to given propagation vector. */ - if (fabs(prop_0[2])>=1) { // can not be >1 except for machine precision - incPolX_0[0]=prop_0[2]; + if ((prop_0[0]==0 && prop_0[1]==0) || fabs(prop_0[2])==1) { // a robust test of propagation along the z-axis + propAlongZ=incPolX_0[0]=copysign(1.0,prop_0[2]); incPolY_0[1]=1; incPolX_0[1]=incPolX_0[2]=incPolY_0[0]=incPolY_0[2]=0.0; } else { + propAlongZ=0; temp=sqrt(1-prop_0[2]*prop_0[2]); incPolX_0[0]=prop_0[0]*prop_0[2]/temp; incPolX_0[1]=prop_0[1]*prop_0[2]/temp; @@ -1891,10 +1969,11 @@ void VariablesInterconnect(void) else { // else - initialize rotation stuff InitRotation(); - /* if not default incidence, break the symmetry completely. This can be improved to account for some special - * cases, however, then symmetry of Gaussian beam should be treated more thoroughly than now. + /* if not default incidence (generally, along z-axis), break the symmetry completely. This can be improved to + * account for some special cases, however, then symmetry of Gaussian beam should be treated more thoroughly + * than now. */ - if (prop[2]!=1 && sym_type==SYM_AUTO) sym_type=SYM_NO; + if (!propAlongZ && sym_type==SYM_AUTO) sym_type=SYM_NO; } ipr_required=(IterMethod==IT_BICGSTAB || IterMethod==IT_CGNR); /* TO ADD NEW ITERATIVE SOLVER @@ -2044,29 +2123,27 @@ void PrintInfo(void) fprintf(logfile,"box dimensions: %ix%ix%i\n",boxX,boxY,boxZ); if (anisotropy) { fprintf(logfile,"refractive index (diagonal elements of the tensor):\n"); - if (Nmat==1) fprintf(logfile," "CFORM3V"\n",creal(ref_index[0]),cimag(ref_index[0]),creal(ref_index[1]), - cimag(ref_index[1]),creal(ref_index[2]),cimag(ref_index[2])); + if (Nmat==1) fprintf(logfile," "CFORM3V"\n",REIM3V(ref_index)); else { for (i=0;i&2 fi diff --git a/tests/2exec/suite b/tests/2exec/suite index eea1fe1b..afcaf277 100644 --- a/tests/2exec/suite +++ b/tests/2exec/suite @@ -75,9 +75,6 @@ all -h chp_load # the following may cause errors in mpi_seq mode all -chp_dir chp_tmp -chp_load ;mgn; -all -h Cpr -all -Cpr ;sep; ;mgn; - all -h Cpr all -Cpr ;mgn; all -Cpr ;sep; ;mgn; @@ -161,6 +158,7 @@ all -opt mem ;mgn; all -h orient all -orient 10 20 30 ;se; ;mgn; +all -orient 10 20 30 ;se; ;mgn; ;p; -scat_matr both all -orient avg ;se; ;mg4n; all -orient avg ap.dat ;se; ;mg4n; diff --git a/tests/2exec/suite_sparse b/tests/2exec/suite_sparse index 58be2cd6..107821a4 100644 --- a/tests/2exec/suite_sparse +++ b/tests/2exec/suite_sparse @@ -143,6 +143,7 @@ all -opt mem ;mgn; all -h orient all -orient 10 20 30 ;se; ;mn; +all -orient 10 20 30 ;se; ;mn; ;p; -scat_matr both all -orient avg ;se; ;mn; all -orient avg ap.dat ;se; ;mn; diff --git a/tests/2exec/suite_surf b/tests/2exec/suite_surf new file mode 100644 index 00000000..391861f8 --- /dev/null +++ b/tests/2exec/suite_surf @@ -0,0 +1,311 @@ +#------------------------- Variable definitons ------------------------------------------------------------------------- +# All variable names should be enclosed in ';' and start the line - its value is assigned to the rest of the line. It +# can be defined through other variables defined below. Each variable must be defined only once. + +# default for most tests +;mgn; ;m; ;g; ;n; +# default for two-domain particles +;2mgn; ;2m; ;g; ;n; +# for large computations (e.g. orientation averaging) +;mg4n; ;m; -grid 4 ;n; +# default addition to make particle completely non-symmetric (can be replaced by ;sc; with two refractive indices below) +;sep; ;se; ;p; +# default for testing of different grids +;smn; -size 8 ;m; ;n; + +# Elementary variables +;m; -m 1.1 0.1 +;2m; -m 1.1 0.1 1.2 0.2 +;3m; -m 1.05 0.1 1.1 0.2 1.2 0.3 +;g; -grid 8 +;n; -ntheta 5 +;p; -prop 1 2 3 +;se; -shape ellipsoid 0.5 1.5 +;sc; -shape coated 0.3 0.2 0.2 0.2 + +#----------------------------- List of tests --------------------------------------------------------------------------- +# The format is the following: ' ' the first one is coma-separated list of files to +# compare or 'all' (which compares all produced files. is everything after the first space and it is passed +# directly to ADDA. + +all + +# testing of different grids - relevant for FFT methods +# to remove redundant warnings for ocl_seq, only (2,3,5) numbers are used +all -grid 2 ;smn; +all -grid 4 ;smn; +all -grid 6 ;smn; +all -grid 8 ;smn; +all -grid 10 ;smn; +all -grid 12 ;smn; +all -grid 16 ;smn; +all -grid 18 ;smn; +all -grid 20 ;smn; +all -grid 24 ;smn; +all -grid 30 ;smn; +all -grid 32 ;smn; + +all -h alldir_inp +all -alldir_inp adp.dat -Csca ;mgn; + +all -h anisotr +all -anisotr ;3m; ;g; ;n; + +all -h asym +# nontrivial propagation direction make the results different +#all -asym ;sep; ;mgn; +all -asym ;sc; ;2mgn; + +all -h beam +all -h beam plane +all -beam plane ;mgn; +all -h beam lminus +#all -beam lminus 2 1 2 3 ;mgn; +all -h beam davis3 +#all -beam davis3 2 1 2 3 ;mgn; +all -h beam barton5 +#all -beam barton5 2 1 2 3 ;mgn; +all -h beam read +#all -beam read IncBeam-Y IncBeam-X ;se; ;mgn; + +all -h chpoint +all -chpoint 1s -eps 3 ;mgn; +all -h chp_type +all -chp_type normal -chpoint 1s -eps 3 ;mgn; +all -chp_dir chp_tmp -chp_type regular -chpoint 1s -eps 3 ;mgn; +all -h chp_dir +all -chp_dir chp_tmp -chp_type always -eps 3 ;mgn; +all -h chp_load +# the following may cause errors in mpi_seq mode +all -chp_dir chp_tmp -chp_load ;mgn; + +all -h Cpr +# radiative forces are not yet supported with surf +#all -Cpr ;mgn; +#all -Cpr ;sep; ;mgn; + +all -h Csca +all -Csca ;se; ;mgn; + +all -h dir + +all -h dpl +all -dpl 20 ;mgn; + +all -h eps +all -eps 10 ;mgn; + +all -h eq_rad +all -eq_rad 1 ;mgn; + +# It is hard to make meaningful comparison of stdout and log for random placement of granules. However, optical +# properties are compared using rather large tolerances +all -h granul +CrossSec-Y,CrossSec-X,mueller -granul 0.2 0.5 2 -size 8 -shape coated 0.5 ;3m; ;n; +CrossSec-Y,CrossSec-X,mueller -granul 0.2 2 2 -size 8 -shape coated 0.5 ;3m; ;n; + +all -h grid +all -grid 4 6 8 ;m; ;n; + +all -h +all -h h + +all -h init_field +all -init_field auto ;mgn; +all -init_field inc ;mgn; +all -init_field read IncBeam-Y IncBeam-X ;se; ;mgn; +# WKB incident field is not yet suppoted with surf +#all -init_field wkb ;mgn; +all -init_field zero ;mgn; + +all -h int +all -int fcd ;mgn; +all -int fcd_st ;mgn; +all -int igt ;mgn; +all -int igt 3 ;mg4n; +all -int igt 3 0.01 ;mgn; +all -int igt_so ;mgn; +all -int poi ;mgn; +all -int so ;mgn; + +all -h iter +all -iter bcgs2 ;mgn; +all -iter bicg ;mgn; +all -iter bicgstab ;mgn; +all -iter cgnr ;mgn; +all -iter csym ;mgn; +all -iter qmr ;mgn; +all -iter qmr2 ;mgn; + +all -h jagged +all -jagged 2 ;mg4n; + +all -h lambda +all -lambda 1 ;mgn; + +all -h m +all -m 1.2 0.2 ;g; ;n; + +all -h maxiter +all -maxiter 5 ;mgn; + +all -h no_reduced_fft +all -no_reduced_fft ;mgn; + +all -h no_vol_cor +all -no_vol_cor -size 3 ;mgn; + +all -h ntheta +all -ntheta 10 ;m; ;g; + +all -h opt +all -opt speed ;mgn; +all -opt mem ;mgn; + +all -h orient +# changing particle orientation is not yet suppoted with surf +#all -orient 10 20 30 ;se; ;mgn; +#all -orient 10 20 30 ;se; ;mgn; ;p; -scat_matr both +#all -orient avg ;se; ;mg4n; +#all -orient avg ap.dat ;se; ;mg4n; + +all -h phi_integr +# nontrivial propagation direction make the results different +#all -phi_integr 31 ;sep; ;mgn; +all -phi_integr 31 ;sc; ;2mgn; + +all -h pol +all -pol cldr ;p; ;mgn; +all -pol cm ;mgn; +all -pol dgf ;mgn; +all -pol fcd ;mgn; +all -pol igt_so ;mgn; +all -pol lak ;mgn; +all -pol ldr ;p; ;mgn; +all -pol ldr avgpol ;p; ;mgn; +all -pol rrc ;mgn; +all -pol so ;p; ;mgn; + +all -h prognosis +all -prognosis + +all -h prop +all -prop 5 1 3 ;se; ;mgn; + +all -h recalc_resid +all -recalc_resid ;mgn; + +all -h save_geom +all -save_geom -shape coated 0.4 0.1 0.15 0.2 -prognosis +all -h sg_format +all -save_geom -shape ellipsoid 0.5 0.25 -prognosis -sg_format text +all -save_geom -shape ellipsoid 0.5 0.25 -prognosis -sg_format text_ext +all -save_geom -shape ellipsoid 0.5 0.25 -prognosis -sg_format ddscat6 +all -save_geom -shape ellipsoid 0.5 0.25 -prognosis -sg_format ddscat7 + +all -h scat +all -scat dr ;mgn; +all -scat fin ;mgn; +all -scat igt_so ;mgn; +# the results for so are indeed different, because surf uses formulae for non-plane beams +#all -scat so ;mgn; + +all -h scat_grid_inp +all -scat_grid_inp sp.dat ;mgn; + +all -h scat_matr +all -scat_matr muel ;mgn; +all -scat_matr ampl ;mgn; +all -scat_matr both ;mgn; +# for equivalent runs with surface we additionally add -yz, which is incompatible with the following +#all -scat_matr none ;m; ;g; + +all -h shape +all -h shape axisymmetric +all -shape axisymmetric 196.txt ;mgn; +all -h shape bicoated +all -shape bicoated 3 0.5 ;2mgn; +all -h shape biellipsoid +all -shape biellipsoid 0.5 1.5 0.75 0.5 1.5 ;2mgn; +all -h shape bisphere +all -shape bisphere 2 ;mgn; +all -h shape box +all -shape box ;mgn; +all -shape box 0.5 1.5 ;mgn; +all -h shape capsule +all -shape capsule 1.5 ;mgn; +all -h shape chebyshev +all -shape chebyshev 0.3 5 ;mgn; +all -h shape coated +all -shape coated 0.5 ;2mgn; +all -shape coated 0.4 0.1 0.15 0.2 ;2mgn; +all -h shape cylinder +all -shape cylinder 1 ;mgn; +all -h shape egg +all -shape egg 0.5 0.2 ;mgn; +all -h shape ellipsoid +all -shape ellipsoid 0.25 2 ;mgn; +all -h shape line +all -shape line -grid 16 ;m; ;n; +all -h shape plate +all -shape plate 0.5 ;mgn; +all -h shape prism +all -shape prism 5 1.5 ;mgn; +all -h shape rbc +all -shape rbc 0.3 0.1 0.3 ;mgn; +all -h shape read +all -shape read ellipsoid.geom ;m; ;n; +all -shape read coated.geom ;2m; ;n; +all -shape read ell_ddscat6.dat ;m; ;n; +all -shape read ell_ddscat7.dat ;m; ;n; +all -h shape sphere +all -shape sphere ;mgn; +all -h shape spherebox +all -shape spherebox 0.5 ;2mgn; + +all -h size +all -size 8 ;mgn; + +all -h store_beam +all -store_beam ;se; ;mgn; + +all -h store_dip_pol +all -store_dip_pol ;se; ;mgn; + +all -h store_force +# radiative forces are not yet supported with surf +#all -store_force ;sep; ;mgn; + +all -h store_grans +granules -store_grans -granul 0.2 1 -size 4 ;2mgn; + +all -h store_int_field +all -store_int_field ;se; ;mgn; + +all -h store_scat_grid +# nontrivial propagation direction make the results different +#all -store_scat_grid ;sep; ;mgn; +all -store_scat_grid ;sc; ;2mgn; + +all -h sym +all -sym auto ;mgn; +all -sym no ;mgn; +all -sym enf ;mgn; +all -sym auto ;sep; ;mgn; +all -sym no ;sep; ;mgn; +all -sym enf ;sep; ;mgn; + +all -h test +all -test ;mgn; + +all -h V +all -V + +all -h vec +# nontrivial propagation direction make the results different +#all -vec ;sep; ;mgn; +all -vec ;sc; ;2mgn; + +all -h yz +# for equivalent runs with surface we additionally add -yz, which is incompatible with the following +#all -yz -store_scat_grid ;mgn;