Skip to content

Commit

Permalink
New equiv tests + fixed TML bb
Browse files Browse the repository at this point in the history
  • Loading branch information
stefaniagl committed Feb 4, 2022
1 parent 63c59c8 commit 62be280
Show file tree
Hide file tree
Showing 4 changed files with 184 additions and 10 deletions.
3 changes: 1 addition & 2 deletions examples/bessel/bb_examples.py
Original file line number Diff line number Diff line change
Expand Up @@ -160,7 +160,6 @@ def plotField(xd,yd,ed,z0,mode):
# data save
os.makedirs('saved', exist_ok=True)

plt.savefig('saved/option_1.pdf', bbox_inches='tight')
plt.savefig('saved/option_2.pdf', bbox_inches='tight')
plt.savefig('saved/results.pdf', bbox_inches='tight')

plt.show()
16 changes: 8 additions & 8 deletions src/GenerateB.c
Original file line number Diff line number Diff line change
Expand Up @@ -72,10 +72,10 @@ static doublecomplex besM[4]; // Matrix M defining the generalized Bessel be
* afterwards. If you need local, intermediate variables, put them into the beginning of the corresponding function.
* Add descriptive comments, use 'static'.
*/

// EXTERNAL FUNCTIONS

#ifndef NO_FORTRAN
#ifndef NO_FORTRAN
void bjndd_(const int *n,const double *x,double *bj,double *dj,double *fj);
#endif

Expand All @@ -89,7 +89,7 @@ void InitBeam(void)
/* 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
Expand All @@ -100,9 +100,9 @@ void InitBeam(void)
case B_PLANE:
if (IFROOT) beam_descr="plane wave";
if (surface) {
/* here we assume that prop_0 will not change further (e.g., by rotation of particle),
/* here we assume that prop_0 will not change further (e.g., by rotation of particle),
* i.e. prop=prop_0 in GenerateBeam() below
*/
*/
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 "
Expand Down Expand Up @@ -255,7 +255,7 @@ void InitBeam(void)
case B_BES_TML:
TestRangeNN(beam_pars[1],"half-cone angle for TML type",0,90);
besM[0]=0; besM[1]=besKz/besKt;
besM[2]=WaveNum/besKt; besM[2]=0;
besM[2]=WaveNum/besKt; besM[3]=0;
if (IFROOT) tmp_str="linear component of the TM";
break;
default: LogError(ONE_POS,"Incompatibility error in GenerateB");
Expand Down Expand Up @@ -359,7 +359,7 @@ void GenerateB (const enum incpol which, // x - or y polarized incident light
* the substrate (below) is Exp(i*k*msub*r.a). We assume that the incoming beam is homogeneous in its
* original medium.
*/
double hbeam=hsub+beam_center[2]; // height of beam center above the surface
double hbeam=hsub+beam_center[2]; // height of beam center above the surface
if (prop[2]>0) { // beam comes from the substrate (below)
doublecomplex tc; // transmission coefficients
// determine amplitude of the transmitted wave; here msub is always defined
Expand Down Expand Up @@ -635,7 +635,7 @@ void GenerateB (const enum incpol which, // x - or y polarized incident light
cvMultScal_cmplx(ctemp,v1,b+j);
}
return;
#endif // !NO_FORTRAN
#endif // !NO_FORTRAN
case B_READ:
if (which==INCPOL_Y) fname=beam_fnameY;
else fname=beam_fnameX; // which==INCPOL_X
Expand Down
5 changes: 5 additions & 0 deletions tests/equiv/bessel/README.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
The comparison of equivalent cases (command lines)
for the scattering of Bessel beams calculated in ADDA.
Differences are shown when |diff| > 1e-06

See results in results.txt
170 changes: 170 additions & 0 deletions tests/equiv/bessel/bb_equiv.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,170 @@
'''
# This code represents the comparison of equivalent cases (command lines)
# for the scattering of Bessel beams calculated in ADDA.
'''

import os,shutil,re,math


# path to adda executable
#adda_exec = "../../win64/adda.exe"
adda_exec = os.path.abspath(__file__ + "/../../../../src/seq/adda")

dirname = 'out'
fdiff = 1e-6 # fixed difference

pi = 3.14159265359


def extractCrosSeq(mode,pol):
dname = dirname + str(mode)
files = os.listdir(dname)
with open(dname + '/' + str(files[files.index('CrossSec-' + pol)])) as f:
file = f.read()
f.close()
dat = re.split('[\n \t]',file)
Cext = float(dat[2])
Qext = float(dat[5])
Cabs = float(dat[8])
Qabs = float(dat[11])

return Cext,Qext,Cabs,Qabs

# data generation (run of ADDA code)
def adda_run(mode,option):
dname = dirname + str(mode)
#os.makedirs(dname, exist_ok=True)
cmdline = adda_exec + ' -dir ' + dname + option
os.system(cmdline)

return 0

# print difference
def printdiff(val,x,y):
cdiff = math.fabs(x-y) # calculated difference
if (cdiff > fdiff):
print('\n\t\t'+val + ':\n\t\tcase 1:\t'+str(x)+'\n\t\tcase 2:\t'+str(y)+'\n\t\tdiff:\t'+str(cdiff)+'\n')
#print('\ncase 1:\t'+str(x)+'\ncase 2:\t'+str(y)+'\nvariation:\t'+str(cdiff)+' %')
fr.write('\n\t\t'+val + ':\n\t\tcase 1:\t'+str(x)+'\n\t\tcase 2:\t'+str(y)+'\n\t\tdiff:\t'+str(cdiff)+'\n')

return 0

def printsearch(pol):
print('\n\tCrosSec-'+pol+': search diff...')
fr.write('\n\n\tCrosSec-'+pol+': search diff')
cext1,qext1,cabs1,qabs1 = extractCrosSeq(1,pol)
cext2,qext2,cabs2,qabs2 = extractCrosSeq(2,pol)
printdiff('Cext',cext1,cext2)
printdiff('Qext',qext1,qext2)
printdiff('Cabs',cabs1,cabs2)
printdiff('Qabs',qabs1,qabs2)
print('\tCrosSec-'+pol+': done')
fr.write('\tCrosSec-'+pol+': done')

def compare(line1,line2):

adda_run(1,line1)
adda_run(2,line2)

print('\nCompare:\n'+line1+'\n'+line2)
fr.write('\n\nCompare:\n'+line1+'\n'+line2)
printsearch('X')
printsearch('Y')
print('\nDone\n___')
fr.write('\n\nDone\n___')

shutil.rmtree('out1')
shutil.rmtree('out2')

return 0

fr = open("results.txt", "w")
fr.write("The comparison of equivalent cases (command lines) \nfor the scattering of Bessel beams calculated in ADDA.")
fr.write("\nDifferences are shown when |diff| > " + str(fdiff) + "\n\n___")
fr.close()

fr = open("results.txt", "a")


print('\n\nPlane wave case of LE Bessel beam')
fr.write('\n\nPlane wave case of LE Bessel beam')
common = ' -sym no'
opt1 = common
opt2 = common + ' -beam besselLE 0 0'
compare(opt1,opt2)

print('\n\nPlane wave case of LM Bessel beam')
fr.write('\n\nPlane wave case of LM Bessel beam')
common = ' -sym no'
opt1 = common
opt2 = common + ' -beam besselLM 0 0'
compare(opt1,opt2)

print('\n\nPlane wave case of CS Bessel beam')
fr.write('\n\nPlane wave case of CS Bessel beam')
common = ' -sym no'
opt1 = common
opt2 = common + ' -beam besselCS 0 0'
compare(opt1,opt2)

print('\n\nGeneralized and LE Bessel beams')
fr.write('\n\nGeneralized and LE Bessel beams')
opt1 = ' -beam besselM 2 5 0 0 0 1'
opt2 = ' -beam besselLE 2 5'
compare(opt1,opt2)

opt1 = ' -beam besselM 2 85 0 0 0 1'
opt2 = ' -beam besselLE 2 85'
compare(opt1,opt2)

print('\n\nGeneralized and LM Bessel beams')
fr.write('\n\nGeneralized and LM Bessel beams')
opt1 = ' -beam besselM 2 5 0 1 0 0'
opt2 = ' -beam besselLM 2 5'
compare(opt1,opt2)

opt1 = ' -beam besselM 2 85 0 1 0 0'
opt2 = ' -beam besselLM 2 85'
compare(opt1,opt2)

print('\n\nGeneralized and CS Bessel beams')
fr.write('\n\nGeneralized and CS Bessel beams')
opt1 = ' -beam besselM 2 5 0.5 0 0 0.5'
opt2 = ' -beam besselCS 2 5'
compare(opt1,opt2)

opt1 = ' -beam besselM 2 85 0.5 0 0 0.5'
opt2 = ' -beam besselCS 2 85'
compare(opt1,opt2)

print("\n\nGeneralized and CS' Bessel beams")
fr.write("\n\nGeneralized and CS' Bessel beams")
opt1 = ' -beam besselM 2 5 0.5 0 0 -0.5'
opt2 = ' -beam besselCSp 2 5'
compare(opt1,opt2)

opt1 = ' -beam besselM 2 85 0.5 0 0 -0.5'
opt2 = ' -beam besselCSp 2 85'
compare(opt1,opt2)

print('\n\nGeneralized and TEL Bessel beams')
fr.write('\n\nGeneralized and TEL Bessel beams')
opt1 = ' -beam besselM 2 5 '+str(-1/math.sin(5*pi/180))+' 0 0 '+str(1/math.tan(5*pi/180))
opt2 = ' -beam besselTEL 2 5'
compare(opt1,opt2)

opt1 = ' -beam besselM 2 85 '+str(-1/math.sin(85*pi/180))+' 0 0 '+str(1/math.tan(85*pi/180))
opt2 = ' -beam besselTEL 2 85'
compare(opt1,opt2)

print('\n\nGeneralized and TML Bessel beams')
fr.write('\n\nGeneralized and TML Bessel beams')
opt1 = ' -beam besselM 2 5 0 '+str(1/math.tan(5*pi/180))+' '+str(1/math.sin(5*pi/180))+' 0'
opt2 = ' -beam besselTML 2 5'
compare(opt1,opt2)

opt1 = ' -beam besselM 2 85 0 '+str(1/math.tan(85*pi/180))+' '+str(1/math.sin(85*pi/180))+' 0'
opt2 = ' -beam besselTML 2 85'
compare(opt1,opt2)

fr.close()

0 comments on commit 62be280

Please sign in to comment.