Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Multi-surface mode] #312

Draft
wants to merge 12 commits into
base: master
Choose a base branch
from
2 changes: 1 addition & 1 deletion src/CalculateE.c
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,7 @@ static void ComputeMuellerMatrix(double matrix[4][4], const doublecomplex s1,con

if (surface) {
double scale=inc_scale;
if (TestBelowDeg(theta)) scale*=creal(msub);
if (TestBelowDeg(theta)) scale*=creal(sub.m[sub.N-1]);
if (fabs(scale-1)>ROUND_ERR) for (int i=0;i<4;i++) for (int j=0;j<4;j++) matrix[i][j]*=scale;
}
}
Expand Down
76 changes: 38 additions & 38 deletions src/GenerateB.c
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ double C0dipole,C0dipole_refl; // inherent cross sections of exciting dipole (in
// 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
* The transmitted wave can be inhomogeneous wave (when reflactive index of the substrate 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
Expand Down Expand Up @@ -95,31 +95,31 @@ void InitBeam(void)
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 "
if (sub.mInf && 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
// here sub.m[sub.N-1] should always be defined
inc_scale=1/creal(sub.m[sub.N-1]);
ki=sub.m[sub.N-1]*prop_0[2];
/* Special case for sub.m[sub.N-1] near 1 to remove discontinuities for near-grazing incidence. The details
* are discussed in CalcFieldSurf() in crosssec.c.
*/
if (cabs(msub-1)<ROUND_ERR && cabs(ki)<SQRT_RND_ERR) kt=ki;
else kt=cSqrtCut(1 - msub*msub*(prop_0[0]*prop_0[0]+prop_0[1]*prop_0[1]));
if (cabs(sub.m[sub.N-1]-1)<ROUND_ERR && cabs(ki)<SQRT_RND_ERR) kt=ki;
else kt=cSqrtCut(1 - sub.m[sub.N-1]*sub.m[sub.N-1]*(prop_0[0]*prop_0[0]+prop_0[1]*prop_0[1]));
// determine propagation direction and full wavevector of wave transmitted into substrate
ktVec[0]=msub*prop_0[0];
ktVec[1]=msub*prop_0[1];
ktVec[0]=sub.m[sub.N-1]*prop_0[0];
ktVec[1]=sub.m[sub.N-1]*prop_0[1];
ktVec[2]=kt;
}
else if (prop_0[2]<0) { // beam comes from above the substrate
inc_scale=1;
vRefl(prop_0,prIncRefl);
ki=-prop_0[2]; // always real
if (!msubInf) {
if (!sub.mInf) {
// same special case as above
if (cabs(msub-1)<ROUND_ERR && cabs(ki)<SQRT_RND_ERR) kt=ki;
else kt=cSqrtCut(msub*msub - (prop_0[0]*prop_0[0]+prop_0[1]*prop_0[1]));
if (cabs(sub.m[sub.N-1]-1)<ROUND_ERR && cabs(ki)<SQRT_RND_ERR) kt=ki;
else kt=cSqrtCut(sub.m[sub.N-1]*sub.m[sub.N-1] - (prop_0[0]*prop_0[0]+prop_0[1]*prop_0[1]));
// determine propagation direction of wave transmitted into substrate
ktVec[0]=prop_0[0];
ktVec[1]=prop_0[1];
Expand All @@ -129,7 +129,7 @@ void InitBeam(void)
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);
if (!msubInf) {
if (!sub.mInf) {
vReal(ktVec,prIncTran);
vNormalize(prIncTran);
}
Expand All @@ -138,7 +138,7 @@ void InitBeam(void)
case B_DIPOLE:
vCopy(beam_pars,beam_center_0);
if (surface) {
if (beam_center_0[2]<=-hsub)
if (beam_center_0[2]<=-sub.hP)
PrintErrorHelp("External dipole should be placed strictly above the surface");
inc_scale=1; // but scaling of Mueller matrix is weird anyway
}
Expand Down Expand Up @@ -257,14 +257,14 @@ void GenerateB (const enum incpol which, // x - or y polarized incident light
if (surface) {
/* With respect to normalization we use here the same assumption as in the free space - the origin is in
* the particle center, and amplitude of incoming plane wave is equal to 1. Then irradiance of the beam
* coming from below is c*Re(msub)/(8pi), different from that coming from above.
* coming from below is c*Re(sub.m[sub.N-1])/(8pi), different from that coming from above.
* Original incident (incoming) beam propagating from the vacuum (above) is Exp(i*k*r.a), while - from
* the substrate (below) is Exp(i*k*msub*r.a). We assume that the incoming beam is homogeneous in its
* the substrate (below) is Exp(i*k*sub.m[sub.N-1]*r.a). We assume that the incoming beam is homogeneous in its
* original medium.
*/
doublecomplex rc,tc; // reflection and transmission coefficients
if (prop[2]>0) { // beam comes from the substrate (below)
// determine amplitude of the reflected and transmitted waves; here msub is always defined
// determine amplitude of the reflected and transmitted waves; here sub.m[sub.N-1] is always defined
if (which==INCPOL_Y) { // s-polarized
cvBuildRe(ex,eIncRefl);
cvBuildRe(ex,eIncTran);
Expand All @@ -274,12 +274,12 @@ void GenerateB (const enum incpol which, // x - or y polarized incident light
else { // p-polarized
vInvRefl_cr(ex,eIncRefl);
crCrossProd(ey,ktVec,eIncTran);
rc=FresnelRP(ki,kt,1/msub);
tc=FresnelTP(ki,kt,1/msub);
rc=FresnelRP(ki,kt,1/sub.m[sub.N-1]);
tc=FresnelTP(ki,kt,1/sub.m[sub.N-1]);
}
// 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);
// phase shift due to the origin at height sub.hP
cvMultScal_cmplx(rc*cexp(-2*I*WaveNum*ki*sub.hP),eIncRefl,eIncRefl);
cvMultScal_cmplx(tc*cexp(I*WaveNum*(kt-ki)*sub.hP),eIncTran,eIncTran);
// main part
for (i=0;i<local_nvoid_Ndip;i++) {
j=3*i;
Expand All @@ -291,7 +291,7 @@ void GenerateB (const enum incpol which, // x - or y polarized incident light
// determine amplitude of the reflected and transmitted waves
if (which==INCPOL_Y) { // s-polarized
cvBuildRe(ex,eIncRefl);
if (msubInf) {
if (sub.mInf) {
rc=-1;
tc=0; // to remove compiler warnings
}
Expand All @@ -303,20 +303,20 @@ void GenerateB (const enum incpol which, // x - or y polarized incident light
}
else { // p-polarized
vInvRefl_cr(ex,eIncRefl);
if (msubInf) {
if (sub.mInf) {
rc=1;
tc=0; // to remove compiler warnings
}
else {
crCrossProd(ey,ktVec,eIncTran);
cvMultScal_cmplx(1/msub,eIncTran,eIncTran); // normalize eIncTran by ||ktVec||=msub
rc=FresnelRP(ki,kt,msub);
tc=FresnelTP(ki,kt,msub);
cvMultScal_cmplx(1/sub.m[sub.N-1],eIncTran,eIncTran); // normalize eIncTran by ||ktVec||=sub.m[sub.N-1]
rc=FresnelRP(ki,kt,sub.m[sub.N-1]);
tc=FresnelTP(ki,kt,sub.m[sub.N-1]);
}
}
// phase shift due to the origin at height hsub
cvMultScal_cmplx(rc*imExp(2*WaveNum*creal(ki)*hsub),eIncRefl,eIncRefl); // assumes real ki
if (!msubInf) cvMultScal_cmplx(tc*cexp(I*WaveNum*(ki-kt)*hsub),eIncTran,eIncTran);
// phase shift due to the origin at height sub.hP
cvMultScal_cmplx(rc*imExp(2*WaveNum*creal(ki)*sub.hP),eIncRefl,eIncRefl); // assumes real ki
if (!sub.mInf) cvMultScal_cmplx(tc*cexp(I*WaveNum*(ki-kt)*sub.hP),eIncTran,eIncTran);
// main part
for (i=0;i<local_nvoid_Ndip;i++) {
j=3*i;
Expand All @@ -341,7 +341,7 @@ void GenerateB (const enum incpol which, // x - or y polarized incident light
(*InterTerm_real)(r1,gt);
cSymMatrVecReal(gt,dip_p,b+j);
if (surface) { // add reflected field
r1[2]=DipoleCoord[j+2]+beam_center[2]+2*hsub;
r1[2]=DipoleCoord[j+2]+beam_center[2]+2*sub.hP;
(*ReflTerm_real)(r1,gt);
cReflMatrVecReal(gt,dip_p,v1);
cvAdd(v1,b+j,b+j);
Expand All @@ -358,7 +358,7 @@ void GenerateB (const enum incpol which, // x - or y polarized incident light
C0dipole=2*FOUR_PI_OVER_THREE*temp*temp;
if (surface) {
r1[0]=r1[1]=0;
r1[2]=2*(beam_center[2]+hsub);
r1[2]=2*(beam_center[2]+sub.hP);
(*ReflTerm_real)(r1,gt);
double tmp;
/* the following expression uses that dip_p is real and a specific (anti-)symmetry of the gt
Expand Down Expand Up @@ -449,11 +449,11 @@ void GenerateB (const enum incpol which, // x - or y polarized incident light
* add a case above. Identifier ('B_...') should be defined inside 'enum beam' in const.h. This case should set
* complex vector 'b', describing the incident field in the particle reference frame. It is set inside the cycle for
* each dipole of the particle and is calculated using
* 1) 'DipoleCoord' � array of dipole coordinates;
* 2) 'prop' � propagation direction of the incident field;
* 3) 'ex' � direction of incident polarization;
* 4) 'ey' � complementary unity vector of polarization (orthogonal to both 'prop' and 'ex');
* 5) 'beam_center' � beam center in the particle reference frame (automatically calculated from 'beam_center_0'
* 1) 'DipoleCoord' � array of dipole coordinates;
* 2) 'prop' � propagation direction of the incident field;
* 3) 'ex' � direction of incident polarization;
* 4) 'ey' � complementary unity vector of polarization (orthogonal to both 'prop' and 'ex');
* 5) 'beam_center' � beam center in the particle reference frame (automatically calculated from 'beam_center_0'
* defined in InitBeam).
* If the new beam type is compatible with '-surf', include here the corresponding code. For that you will need
* the variables, related to surface - see vars.c after "// related to a nearby surface".
Expand Down
3 changes: 2 additions & 1 deletion src/calculator.c
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,8 @@ double * restrict muel_alpha; // mueller matrix for different values of alpha

// used in crosssec.c
doublecomplex * restrict E_ad; // complex field E, calculated for alldir
double * restrict E2_alldir; // square of E (scaled with msub, so ~ Poynting vector or dC/dOmega), calculated for alldir
double * restrict E2_alldir; /* square of E (scaled with the refractive index of the last layer, so ~ Poynting vector or
dC/dOmega), calculated for alldir */
palatni marked this conversation as resolved.
Show resolved Hide resolved
doublecomplex cc[MAX_NMAT][3]; // couple constants
#ifndef SPARSE
doublecomplex * restrict expsX,* restrict expsY,* restrict expsZ; // arrays of exponents along 3 axes (for calc_field)
Expand Down
1 change: 1 addition & 0 deletions src/const.h
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,7 @@ the compilation may fail or produce wrong results. If you still want to try, ena
#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
#define MAX_N_LAYERS 10 // maximum number of layers in multi-surface mode
myurkin marked this conversation as resolved.
Show resolved Hide resolved

// 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
Expand Down
42 changes: 21 additions & 21 deletions src/crosssec.c
Original file line number Diff line number Diff line change
Expand Up @@ -652,12 +652,12 @@ static void CalcFieldSurf(doublecomplex ebuff[static restrict 3], // where to wr
cvInit(sumN);
if (above) cvInit(sumF); //additional storage for directly propagated scattering

/* There is an inherent discontinuity for msub approaching 1 and scattering angle 90 degrees (nF[2]=0). The problem
/* There is an inherent discontinuity for sub.m[sub.N-1] approaching 1 and scattering angle 90 degrees (nF[2]=0). The problem
* is that for m=1+-0, but |m-1|>>(nF[2])^2, ki<<kt<<1 => rs=rp=-1
* while for m=1 (exactly) the limit of nF[2]->0 results in kt=ki => rs=rp=0
* Therefore, below is a certain logic, which behaves in an intuitively expected way, for common special cases.
* However, it is still not expected to be continuous for fine-changing parameters (like msub approaching 1).
* In particular, the jump occurs when msub crosses 1+-ROUND_ERR boundary.
* However, it is still not expected to be continuous for fine-changing parameters (like sub.m[sub.N-1] approaching 1).
* In particular, the jump occurs when sub.m[sub.N-1] crosses 1+-ROUND_ERR boundary.
* Still, the discontinuity should apply only to scattering at exactly 90 degrees, but not to, e.g., integral
* quantities, like Csca (if sufficient large number of integration points is chosen).
*/
Expand All @@ -668,44 +668,44 @@ static void CalcFieldSurf(doublecomplex ebuff[static restrict 3], // where to wr
* a particle at a planar interface," J. Opt. Soc. Am. A 30, 2519-2525 (2013) for theoretical discussion of
* this fact.
*/
if (fabs(nF[2])<ROUND_ERR && cabs(msub-1)>ROUND_ERR) {
if (fabs(nF[2])<ROUND_ERR && cabs(sub.m[sub.N-1]-1)>ROUND_ERR) {
cvInit(ebuff);
return;
}
cvBuildRe(nF,nN);
nN[2]*=-1;
ki=nF[2]; // real
if (msubInf) {
if (sub.mInf) {
cs=-1;
cp=1;
}
// since kt is not further needed, we directly calculate cs and cp (equivalent to kt=ki)
else if (cabs(msub-1)<ROUND_ERR && cabs(ki)<SQRT_RND_ERR) cs=cp=0;
else if (cabs(sub.m[sub.N-1]-1)<ROUND_ERR && cabs(ki)<SQRT_RND_ERR) cs=cp=0;
else { // no special treatment here, since other cases, including 90deg-scattering, are taken care above.
kt=cSqrtCut(msub*msub - (nN[0]*nN[0]+nN[1]*nN[1]));
kt=cSqrtCut(sub.m[sub.N-1]*sub.m[sub.N-1] - (nN[0]*nN[0]+nN[1]*nN[1]));
cs=FresnelRS(ki,kt);
cp=FresnelRP(ki,kt,msub);
cp=FresnelRP(ki,kt,sub.m[sub.N-1]);
}
phSh=imExp(2*WaveNum*hsub*creal(ki)); // assumes real ki
phSh=imExp(2*WaveNum*sub.hP*creal(ki)); // assumes real ki
}
else { // transmission; here nF[2] is negative
// formulae correspond to plane wave incoming from below, but with change ki<->kt
if (msubInf) { // no transmission for perfectly reflecting substrate => zero result
if (sub.mInf) { // no transmission for perfectly reflecting substrate => zero result
cvInit(ebuff);
return;
}
kt=-msub*nF[2];
if (cabs(msub-1)<ROUND_ERR && cabs(kt)<SQRT_RND_ERR) ki=kt;
else ki=cSqrtCut(1 - msub*msub*(nF[0]*nF[0]+nF[1]*nF[1]));
kt=-sub.m[sub.N-1]*nF[2];
if (cabs(sub.m[sub.N-1]-1)<ROUND_ERR && cabs(kt)<SQRT_RND_ERR) ki=kt;
else ki=cSqrtCut(1 - sub.m[sub.N-1]*sub.m[sub.N-1]*(nF[0]*nF[0]+nF[1]*nF[1]));
// here nN may be complex, but normalized to n.n=1
nN[0]=msub*nF[0];
nN[1]=msub*nF[1];
nN[0]=sub.m[sub.N-1]*nF[0];
nN[1]=sub.m[sub.N-1]*nF[1];
nN[2]=-ki;
// these formulae works fine for ki=kt (even very small), and ki=kt=0 is impossible here
cs=FresnelTS(kt,ki);
cp=FresnelTP(kt,ki,1/msub);
cp=FresnelTP(kt,ki,1/sub.m[sub.N-1]);
// coefficient comes from k0->k in definition of F(n) (in denominator)
phSh=msub*cexp(I*WaveNum*hsub*(ki-kt));
phSh=sub.m[sub.N-1]*cexp(I*WaveNum*sub.hP*(ki-kt));
}
#ifndef SPARSE
// prepare values of exponents, along each of the coordinates
Expand Down Expand Up @@ -1050,13 +1050,13 @@ void CalcAlldir(void)
Accumulate(E_ad,cmplx_type,2*npoints,&Timing_EFieldADComm);
// calculate square of the field
for (point=0;point<npoints;point++) E2_alldir[point] = cAbs2(E_ad[2*point]) + cAbs2(E_ad[2*point+1]);
/* when below surface we scale E2 by Re(1/msub) in accordance with formula for the Poynting vector (and factor of
/* when below surface we scale E2 by Re(1/sub.m[sub.N-1]) in accordance with formula for the Poynting vector (and factor of
* k_sca^2). After that Csca (and g) computed using the standard formula should correctly describe the energy and
* momentum balance for any (even complex) msub (since the scattered wave is homogeneous at far-field), but doesn't
* momentum balance for any (even complex) sub.m[sub.N-1] (since the scattered wave is homogeneous at far-field), but doesn't
* include energy or momentum absorbed (obtained) by the medium.
*/
if (surface && !msubInf) {
double scale=creal(1/msub); // == Re(msub)/|msub|^2
if (surface && !sub.mInf) {
double scale=creal(1/sub.m[sub.N-1]); // == Re(sub.m[sub.N-1])/|sub.m[sub.N-1]|^2
for (i=0;i<theta_int.N;++i) if (TestBelowDeg(theta_int.val[i]))
for (j=0,point=phi_int.N*i;j<phi_int.N;j++,point++) E2_alldir[point]*=scale;
}
Expand Down
4 changes: 2 additions & 2 deletions src/interaction.c
Original file line number Diff line number Diff line change
Expand Up @@ -1423,7 +1423,7 @@ void InitInteraction(void)
case GR_IMG: SET_FUNC_POINTERS(ReflTerm,img); break;
case GR_SOM:
SET_FUNC_POINTERS(ReflTerm,som);
if (!prognosis) som_init(msub*msub);
if (!prognosis) som_init(sub.m[sub.N-1]*sub.m[sub.N-1]);
CalcSomTable();
break;
/* TO ADD NEW REFLECTION FORMULATION
Expand All @@ -1438,7 +1438,7 @@ void InitInteraction(void)
default: LogError(ONE_POS, "Invalid reflection term calculation method: %d",(int)ReflRelation);
// no break
}
surfRCn=msubInf ? -1 : ((1-msub*msub)/(1+msub*msub));
surfRCn=sub.mInf ? -1 : ((1-sub.m[sub.N-1]*sub.m[sub.N-1])/(1+sub.m[sub.N-1]*sub.m[sub.N-1]));
}

#ifdef USE_SSE3
Expand Down
10 changes: 5 additions & 5 deletions src/make_particle.c
Original file line number Diff line number Diff line change
Expand Up @@ -2437,14 +2437,14 @@ void MakeParticle(void)
if (minZco>DipoleCoord[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.
* substrate - sub.hP+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 "
if (surface && sub.hP<=-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);
"small",minZco,sub.hP);
// save geometry
if (save_geom)
#ifndef SPARSE
Expand All @@ -2469,10 +2469,10 @@ void MakeParticle(void)
box_origin_unif[1]=-dsY*cY;
#ifndef SPARSE
box_origin_unif[2]=dsZ*(local_z0_unif-cZ);
if (surface) ZsumShift=2*(hsub+(local_z0-cZ)*dsZ);
if (surface) ZsumShift=2*(sub.hP+(local_z0-cZ)*dsZ);
#else
box_origin_unif[2]=-dsZ*cZ;
if (surface) ZsumShift=2*(hsub-cZ*dsZ);
if (surface) ZsumShift=2*(sub.hP-cZ*dsZ);
# ifdef PARALLEL
AllGather(NULL,position_full,int3_type,NULL);
# endif
Expand Down
Loading