Skip to content

Commit

Permalink
Added rigorous calculations of Sommerfeld integrals. Imported somnec.…
Browse files Browse the repository at this point in the history 
…c from NEC2C, which computes the integrals themselves. The remaining code is in interaction.c - it creates a table of integrals (for different rho and z) and further indexes it to calculate the reflected Green's tensor. There is a lot that can be improved, but in sparse mode it has the main functionality associated with substrate. Moreover, the '-int_surf som' is only marginally (10-20%) slower than the '-int_surf img'.

Small comparison of simulation results of a sphere on the surface was performed to compare with (Loke&Menguc,2010) - Figs.10,11. Qualitatively the results are very similar, but there are certain quantitative differences. We will do a more careful benchmarking with the FFT version

Options '-surf <mre> <mim> <h>' was changed to '-surf <h> {<mre> <mim>|inf}'. So the order of parameters has been changed and option to specify perfectly reflecting substrate is now available. The results of '-surf ... inf' were checked to correspond to the limit of increasing <mre> to very large values. For perfectly reflecting substrate the default (and only) method to calculate reflected Green's tensor ('-int surf ...') is the image-dipole one.

Version incremented to 1.3a2. Fixed minor bug in sparse_ops.h, introduced recently. Added function sizetVector to memory.c/h. Call to imExp in CalcReflTerm_img was replaced by accImExp.
  • Loading branch information
@myurkin
myurkin committed Sep 13, 2013
1 parent d852cde commit c870ace
Show file tree
Hide file tree
Showing 13 changed files with 1,274 additions and 69 deletions.
43 changes: 27 additions & 16 deletions src/GenerateB.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,7 @@ void InitBeam(void)
if (surface) {
// 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
ki=msub*prop_0[2];
kt=cSqrtCut(1 - msub*msub*(prop_0[0]*prop_0[0]+prop_0[1]*prop_0[1]));
// determine propagation direction and full wavevector of wave transmitted into substrate
Expand All @@ -101,17 +102,21 @@ void InitBeam(void)
else if (prop_0[2]<0) { // beam comes from above the substrate
vRefl(prop_0,prIncRefl);
ki=-prop_0[2];
kt=cSqrtCut(msub*msub - (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];
ktVec[2]=-kt;
if (!msubInf) {
kt=cSqrtCut(msub*msub - (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];
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);
if (!msubInf) {
vReal(ktVec,prIncTran);
vNormalize(prIncTran);
}
}
return;
case B_LMINUS:
Expand Down Expand Up @@ -229,7 +234,7 @@ void GenerateB (const enum incpol which, // x - or y polarized incident light
*/
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
// 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);
Expand All @@ -256,20 +261,26 @@ 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);
cvBuildRe(ex,eIncTran);
rc=FresnelRS(ki,kt);
tc=FresnelTS(ki,kt);
if (msubInf) rc=-1;
else {
cvBuildRe(ex,eIncTran);
rc=FresnelRS(ki,kt);
tc=FresnelTS(ki,kt);
}
}
else { // p-polarized
vInvRefl_cr(ex,eIncRefl);
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);
if (msubInf) rc=1;
else {
crCrossProd(ey,ktVec,eIncTran);
cvMultScal_cmplx(1/msub,eIncTran,eIncTran); // normalize eIncTran by ||ktVec||=msub
rc=FresnelRP(ki,kt,msub);
tc=FresnelTP(ki,kt,msub);
}
}
// phase shift due to the origin at height hsub
cvMultScal_cmplx(rc*imExp(2*WaveNum*ki*hsub),eIncRefl,eIncRefl);
cvMultScal_cmplx(tc*cexp(I*WaveNum*(ki-kt)*hsub),eIncTran,eIncTran);
if (!msubInf) cvMultScal_cmplx(tc*cexp(I*WaveNum*(ki-kt)*hsub),eIncTran,eIncTran);
// main part
for (i=0;i<local_nvoid_Ndip;i++) {
j=3*i;
Expand Down
4 changes: 2 additions & 2 deletions src/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,8 +171,8 @@ override EXTRA_FLAGS +=

# C files are located in source folder (src/), other files may be added below
CSOURCE := ADDAmain.c CalculateE.c calculator.c chebyshev.c comm.c crosssec.c GenerateB.c interaction.c io.c \
iterative.c linalg.c make_particle.c matvec.c memory.c mt19937ar.c param.c Romberg.c sinint.c timing.c \
vars.c
iterative.c linalg.c make_particle.c matvec.c memory.c mt19937ar.c param.c Romberg.c sinint.c somnec.c \
timing.c vars.c
# Fortran files are located in src/fort folder, other files may be added below
FSOURCE := d07hre.f d09hre.f d113re.f d132re.f dadhre.f dchhre.f dcuhre.f dfshre.f dinhre.f drlhre.f dtrhre.f \
propaesplibreintadda.f
Expand Down
2 changes: 1 addition & 1 deletion src/const.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
#define __const_h

// version number (string)
#define ADDA_VERSION "1.3a1"
#define ADDA_VERSION "1.3a2"

/* 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
Expand Down
51 changes: 31 additions & 20 deletions src/crosssec.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -633,6 +633,10 @@ static void CalcFieldSurf(doublecomplex ebuff[static restrict 3], // where to wr
nN[2]*=-1;
}
else { // transmission
if (msubInf) { // no transmission for perfectly reflecting substrate => zero result
cvInit(ebuff);
return;
}
doublecomplex eps=msub*msub;
// effective eps and (real part of) refractive index
doublecomplex epsEff=(creal(eps) + I*cimag(eps)/nF[2]);
Expand Down Expand Up @@ -728,21 +732,28 @@ static void CalcFieldSurf(doublecomplex ebuff[static restrict 3], // where to wr
* quantities, like Csca (if sufficient very number of integration points is chosen)
*/
ki=-nN[2];
if (cabs(msub-1)<ROUND_ERR && fabs(ki)<ROUND_ERR) kt=ki; // special case to avoid randomness due to round-off errors
else kt=cSqrtCut(msub*msub - (nN[0]*nN[0]+nN[1]*nN[1]));
if (above) {
/* here we test only the exact zero, since for other cases (when msub=1, and very small values of ki,kt) the
* above assignment kt=ki guarantees correct results through standard functions
*/
if (ki==0 && kt==0) cs=cp=0;
else {
cs=FresnelRS(ki,kt);
cp=FresnelRP(ki,kt,msub);
}
if (msubInf) {
cs=-1;
cp=1;
kt=NAN; // redundant to remove warnings below
}
else { // below surface
cs=FresnelTS(ki,kt);
cp=FresnelTP(ki,kt,msub);
else {
if (cabs(msub-1)<ROUND_ERR && fabs(ki)<ROUND_ERR) kt=ki; // special case to avoid randomness due to round-off errors
else kt=cSqrtCut(msub*msub - (nN[0]*nN[0]+nN[1]*nN[1]));
if (above) {
/* here we test only the exact zero, since for other cases (when msub=1, and very small values of ki,kt) the
* above assignment kt=ki guarantees correct results through standard functions
*/
if (ki==0 && kt==0) cs=cp=0;
else {
cs=FresnelRS(ki,kt);
cp=FresnelRP(ki,kt,msub);
}
}
else { // below surface
cs=FresnelTS(ki,kt);
cp=FresnelTP(ki,kt,msub);
}
}
// set unit vectors for s- and p-polarizations; es is the same for all cases
if (vAlongZ(nF)) { // special case - es=ey
Expand All @@ -756,7 +767,7 @@ static void CalcFieldSurf(doublecomplex ebuff[static restrict 3], // where to wr
}
CrossProd(es,nN,epN); // epN = es x nN
// epF is different for reflected and transmitted (similar to the code in GenerateB.c)
if (above || cimag(msub)==0) { // epF = es x nF; also includes simple case with real ktVec
if (above || cimag(msub)==0) { // epF = es x nF; also includes simple cases with real ktVec and msubInf
double epFRe[3];
CrossProd(es,nF,epFRe);
cvBuildRe(epFRe,epF);
Expand Down Expand Up @@ -1093,7 +1104,7 @@ static double CscaIntegrand(const int theta,const int phi,double * restrict res)
* still unclear (especially for strong absorption). So we use Re[msub]*||E||^2 in all cases, which should be
* exact for real msub, and a good approximation in the case of (weakly) absorbing medium.
*/
if (surface && cos(Deg2Rad(theta_int.val[theta]))<=-ROUND_ERR) res[0]*=creal(msub);
if (surface && !msubInf && cos(Deg2Rad(theta_int.val[theta]))<=-ROUND_ERR) res[0]*=creal(msub);
return 0;
}

Expand Down Expand Up @@ -1134,7 +1145,7 @@ static double gIntegrand(const int theta,const int phi,double * restrict res)
*
* Moreover, after such modification g = (gCsca)/Csca has very little sense, in particular it is not <cos(th)>
*/
if (surface && cos(th)<=-ROUND_ERR) E_square*=creal(msub)*creal(msub);
if (surface && !msubInf && cos(th)<=-ROUND_ERR) E_square*=creal(msub)*creal(msub);
res[0] = E_square*sin(th)*cos(ph);
res[1] = E_square*sin(th)*sin(ph);
res[2] = E_square*cos(th);
Expand Down Expand Up @@ -1173,7 +1184,7 @@ static double gxIntegrand(const int theta,const int phi,double * restrict res)
else {
res[0]=E2_alldir[AlldirIndex(theta,phi)]*sin(Deg2Rad(th))*cos(Deg2Rad(phi_int.val[phi]));
// the following correction is explained in gIntegrand()
if (surface && cos(Deg2Rad(th))<=-ROUND_ERR) res[0]*=creal(msub)*creal(msub);
if (surface && !msubInf && cos(Deg2Rad(th))<=-ROUND_ERR) res[0]*=creal(msub)*creal(msub);
}
return 0;
}
Expand Down Expand Up @@ -1208,7 +1219,7 @@ static double gyIntegrand(const int theta,const int phi,double * restrict res)
else {
res[0]=E2_alldir[AlldirIndex(theta,phi)]*sin(Deg2Rad(th))*sin(Deg2Rad(phi_int.val[phi]));
// the following correction is explained in gIntegrand()
if (surface && cos(Deg2Rad(th))<=-ROUND_ERR) res[0]*=creal(msub)*creal(msub);
if (surface && !msubInf && cos(Deg2Rad(th))<=-ROUND_ERR) res[0]*=creal(msub)*creal(msub);
}
return 0;
}
Expand Down Expand Up @@ -1237,7 +1248,7 @@ static double gzIntegrand(const int theta,const int phi,double * restrict res)
double th=Deg2Rad(theta_int.val[theta]);
res[0]=E2_alldir[AlldirIndex(theta,phi)]*cos(th);
// the following correction is explained in gIntegrand()
if (surface && cos(th)<=-ROUND_ERR) res[0]*=creal(msub)*creal(msub);
if (surface && !msubInf && cos(th)<=-ROUND_ERR) res[0]*=creal(msub)*creal(msub);
return 0;
}

Expand Down
Loading

0 comments on commit c870ace

Please sign in to comment.