diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index eb066f2b83fcb59fca1ac16316912f3a78663503..0348ee7954691006089154c2aae08569146b65d8 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -182,7 +182,8 @@ Function{
tp[] = (2.*k1z[]*epsr2[])/(k1z[]*epsr2[]+k2z[]*epsr1[]);
spol[] = Vector[Sin[phi0],-Cos[phi0],0];
- ppol[] = (k1r[]/Norm[k1r[]]) /\ spol[];
+ ppol_r[] = (k1r[]/Norm[k1r[]]) /\ spol[];
+ ppol_t[] = (k2[] /Norm[k2[]] ) /\ spol[];
AmplEis[] = spol[];
AmplErs[] = rs[]*spol[];
@@ -444,10 +445,10 @@ PostProcessing {
$int_z_t~{ispecular}~{jspecular}] ] ; In SurfIntBot ; } } }
// Project er_specular on the (s,p) basis
- { Name rp ; Value { Term{ Type Global; [ ppol[] * $er_specular] ; In SurfIntTop ; } } }
- { Name rs ; Value { Term{ Type Global; [ spol[] * $er_specular] ; In SurfIntTop ; } } }
- { Name tp ; Value { Term{ Type Global; [ ppol[] * $et_specular] ; In SurfIntBot ; } } }
- { Name ts ; Value { Term{ Type Global; [ spol[] * $et_specular] ; In SurfIntBot ; } } }
+ { Name rp ; Value { Term{ Type Global; [ ppol_r[] * $er_specular] ; In SurfIntTop ; } } }
+ { Name rs ; Value { Term{ Type Global; [ spol[] * $er_specular] ; In SurfIntTop ; } } }
+ { Name tp ; Value { Term{ Type Global; [ ppol_t[] * $et_specular] ; In SurfIntBot ; } } }
+ { Name ts ; Value { Term{ Type Global; [ spol[] * $et_specular] ; In SurfIntBot ; } } }
}
}
diff --git a/DiffractionGratings/grating3D_parallel_Mmatrix.sh b/DiffractionGratings/grating3D_parallel_Mmatrix.sh
index 4c5f7cfad09490f850b33bf6c413a1fbe52f5150..887d864165e8fa66fed6838ae31e15c1725e05ba 100644
--- a/DiffractionGratings/grating3D_parallel_Mmatrix.sh
+++ b/DiffractionGratings/grating3D_parallel_Mmatrix.sh
@@ -6,8 +6,8 @@ export NPROC=96
# export nb_phi=1
export nb_lam=100
export nb_phi=59
-export lambda_min=350
-export lambda_max=800
+export lambda_min=400
+export lambda_max=1600
export FLAG_TOTAL=0
export GRATING_CASE="half_ellipsoid"
export myDir="res_Matrix_nb_lam"$nb_lam"_nb_phi"$nb_phi"_total"$FLAG_TOTAL
@@ -24,7 +24,7 @@ myfunc() {
local phi=$(echo "scale=10;360/($nb_phi)*$2" | bc )
local psi=$(echo "scale=10;90*$3" | bc )
cd $mysubDir
- getdp grating3D.pro -pre helmholtz_vector -msh grating3D.msh -cal -pos postop_helmholtz_vector -petsc_prealloc 200 -setstring test_case $GRATING_CASE -setnumber lambda0 $lam -setnumber thetadeg 15 -setnumber phideg $phi -setnumber psideg $psi -setnumber FLAG_TOTAL $FLAG_TOTAL
+ getdp grating3D.pro -pre helmholtz_vector -msh grating3D.msh -cal -pos postop_helmholtz_vector -petsc_prealloc 200 -setstring test_case $GRATING_CASE -setnumber lambda0 $lam -setnumber thetadeg 50 -setnumber phideg $phi -setnumber psideg $psi -setnumber FLAG_TOTAL $FLAG_TOTAL
if [ $3 -eq 0 ]; then cp res3D/rs.txt ../r_pin_sout_lam$1_phi$2.out ; fi
if [ $3 -eq 0 ]; then cp res3D/rp.txt ../r_pin_pout_lam$1_phi$2.out ; fi
if [ $3 -eq 1 ]; then cp res3D/rs.txt ../r_sin_sout_lam$1_phi$2.out ; fi
diff --git a/DiffractionGratings/grating3D_postplot_Mmatrix.py b/DiffractionGratings/grating3D_postplot_Mmatrix.py
index 5499df969898e1c4051e114376bafce1b87d5c0e..d3a0975feec718638789c90e8893e50e80e73ebd 100644
--- a/DiffractionGratings/grating3D_postplot_Mmatrix.py
+++ b/DiffractionGratings/grating3D_postplot_Mmatrix.py
@@ -22,10 +22,10 @@ def load_getdp_integral(fname):
else:
return temp[:,1]+1j*temp[:,2]
-nb_lam = 100
+nb_lam = 200
nb_phi = 60
FLAG_TOT = 0
-respath = 'res_Matrix_nb_lam100_nb_phi59_total0/'
+respath = 'res_Matrix_nb_lam200_nb_phi59_total0/'
rpp = np.zeros((nb_lam,nb_phi),dtype=complex)
rps = np.zeros((nb_lam,nb_phi),dtype=complex)
rsp = np.zeros((nb_lam,nb_phi),dtype=complex)
@@ -42,7 +42,7 @@ effr2_sin = np.zeros((9,nb_lam,nb_phi),dtype=complex)
Qscat_sin = np.zeros((nb_lam,nb_phi))
-tab_lam = np.linspace(350,800,nb_lam)
+tab_lam = np.linspace(350,1600,nb_lam)
tab_phi = np.linspace(0,360,nb_phi)
M = np.zeros((4,4,len(tab_lam),len(tab_phi)),dtype=complex)
@@ -55,13 +55,15 @@ M = np.zeros((4,4,len(tab_lam),len(tab_phi)),dtype=complex)
# that are parallel (resp. perpendicular) to the plane of incidence
# rps : send Es_inc (s-in), project the reflected field along p (p-out)
+str_reftrans = 'r'
+
for i in range(nb_lam):
+ print(i)
for j in range(nb_phi-1):
- print(i,j)
- rpp[i,j] = load_getdp_integral(respath+'r_pin_pout_lam%g_phi%g.out'%(i,j))
- rss[i,j] = load_getdp_integral(respath+'r_sin_sout_lam%g_phi%g.out'%(i,j))
- rps[i,j] = load_getdp_integral(respath+'r_sin_pout_lam%g_phi%g.out'%(i,j))
- rsp[i,j] = load_getdp_integral(respath+'r_pin_sout_lam%g_phi%g.out'%(i,j))
+ rpp[i,j] = load_getdp_integral(respath+str_reftrans+'_pin_pout_lam%g_phi%g.out'%(i,j))
+ rss[i,j] = load_getdp_integral(respath+str_reftrans+'_sin_sout_lam%g_phi%g.out'%(i,j))
+ rps[i,j] = load_getdp_integral(respath+str_reftrans+'_sin_pout_lam%g_phi%g.out'%(i,j))
+ rsp[i,j] = load_getdp_integral(respath+str_reftrans+'_pin_sout_lam%g_phi%g.out'%(i,j))
efft1_pin[:,i,j] = load_getdp_integral(respath+'eff_t1_lam%g_phi%g_psi0.out'%(i,j))
efft2_pin[:,i,j] = load_getdp_integral(respath+'eff_t2_lam%g_phi%g_psi0.out'%(i,j))
effr1_pin[:,i,j] = load_getdp_integral(respath+'eff_r1_lam%g_phi%g_psi0.out'%(i,j))
@@ -109,32 +111,37 @@ M[4-1,2-1,:,:] =-np.imag(rpp*np.conj(rsp)-rps*np.conj(rss))
M[4-1,3-1,:,:] =-np.imag(rpp*np.conj(rss)+rps*np.conj(rsp))
M[4-1,4-1,:,:] = np.real(rpp*np.conj(rss)-rps*np.conj(rsp))
-# ### Energy balance
-# fig, axes = plt.subplots(2, 2, figsize=(12,12))
-# L,P = np.meshgrid(tab_lam,tab_phi,indexing='ij')
-# for form in range(2):
-# for pol in range(2):
-# if form==0 and pol==0:
-# data=TOT1_pin
-# title = 'TOT1_pin'
-# if form==0 and pol==1:
-# data=TOT1_sin
-# title = 'TOT1_sin'
-# if form==1 and pol==0:
-# data=TOT2_pin
-# title = 'TOT2_pin'
-# if form==1 and pol==1:
-# data=TOT2_sin
-# title = 'TOT2_sin'
-# zplot = axes[form,pol].pcolormesh(L,P,data.real,vmin=0.99,vmax=1.01)
-# add_colorbar(zplot)
-# plt.savefig('BALANCE_FLAG_TOT%g.jpg'%FLAG_TOT)
+### Energy balance
+fig, axes = plt.subplots(3, 2, figsize=(12,12))
+L,P = np.meshgrid(tab_lam,tab_phi,indexing='ij')
+for form in range(2):
+ for pol in range(2):
+ if form==0 and pol==0:
+ data=TOT1_pin
+ title = 'TOT1 pin'
+ if form==0 and pol==1:
+ data=TOT1_sin
+ title = 'TOT1 sin'
+ if form==1 and pol==0:
+ data=TOT2_pin
+ title = 'TOT2 pin'
+ if form==1 and pol==1:
+ data=TOT2_sin
+ title = 'TOT2 sin'
+ zplot = axes[form,pol].pcolormesh(L,P,data.real,vmin=0.99,vmax=1.01)
+ axes[form,pol].title.set_text(title)
+ add_colorbar(zplot)
+ zplot = axes[2,0].pcolormesh(L,P,T1_sin.real);add_colorbar(zplot)
+ axes[2,0].title.set_text('efft1 pin')
+ zplot = axes[2,1].pcolormesh(L,P,R1_sin.real);add_colorbar(zplot)
+ axes[2,1].title.set_text('effr1 pin')
+plt.savefig('BALANCE_FLAG_TOT%g.jpg'%FLAG_TOT)
fig, axes = plt.subplots(4, 4, subplot_kw=dict(projection='polar') ,figsize=(12,12))
flag_lam = True
rlabel = r"$\lambda_0$ (nm)"
which_r = tab_lam
-which_orig = 0 if not flag_lam else 350
+which_orig = 0 #if not flag_lam else 350
for i in range(4):
for j in range(4):
ax=axes[i,j]
@@ -162,7 +169,7 @@ for i in range(4):
sm.set_array([])
ax.set_rorigin(which_orig)
plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=0.95, hspace=0.25,wspace=0.35)
-# plt.savefig('Mmatrix.jpg')
-plt.savefig('Mmatrix.pdf',bbox_inches='tight',pad_inches=0)
+plt.savefig('Mmatrix_%s.jpg'%str_reftrans)
+# plt.savefig('Mmatrix.pdf',bbox_inches='tight',pad_inches=0)
plt.show()