Skip to content
Snippets Groups Projects
Commit 4005a31a authored by Guillaume Demesy's avatar Guillaume Demesy
Browse files

up

parent f463f63e
No related branches found
No related tags found
No related merge requests found
Pipeline #10490 passed
Stage: test
Hide whitespace changes
Inline Side-by-side
DiffractionGratings/grating3D.pro
+ 6
− 5
View file @ 4005a31a
...@@ -182,7 +182,8 @@ Function{ ...@@ -182,7 +182,8 @@ Function{
tp[] = (2.*k1z[]*epsr2[])/(k1z[]*epsr2[]+k2z[]*epsr1[]); tp[] = (2.*k1z[]*epsr2[])/(k1z[]*epsr2[]+k2z[]*epsr1[]);
spol[] = Vector[Sin[phi0],-Cos[phi0],0]; 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[]; AmplEis[] = spol[];
AmplErs[] = rs[]*spol[]; AmplErs[] = rs[]*spol[];
...@@ -444,10 +445,10 @@ PostProcessing { ...@@ -444,10 +445,10 @@ PostProcessing {
$int_z_t~{ispecular}~{jspecular}] ] ; In SurfIntBot ; } } } $int_z_t~{ispecular}~{jspecular}] ] ; In SurfIntBot ; } } }
// Project er_specular on the (s,p) basis // Project er_specular on the (s,p) basis
{ Name rp ; Value { Term{ Type Global; [ ppol[] * $er_specular] ; In SurfIntTop ; } } } { Name rp ; Value { Term{ Type Global; [ ppol_r[] * $er_specular] ; In SurfIntTop ; } } }
{ Name rs ; Value { Term{ Type Global; [ spol[] * $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 tp ; Value { Term{ Type Global; [ ppol_t[] * $et_specular] ; In SurfIntBot ; } } }
{ Name ts ; Value { Term{ Type Global; [ spol[] * $et_specular] ; In SurfIntBot ; } } } { Name ts ; Value { Term{ Type Global; [ spol[] * $et_specular] ; In SurfIntBot ; } } }
} }
} }
......
DiffractionGratings/grating3D_parallel_Mmatrix.sh
+ 3
− 3
View file @ 4005a31a
...@@ -6,8 +6,8 @@ export NPROC=96 ...@@ -6,8 +6,8 @@ export NPROC=96
# export nb_phi=1 # export nb_phi=1
export nb_lam=100 export nb_lam=100
export nb_phi=59 export nb_phi=59
export lambda_min=350 export lambda_min=400
export lambda_max=800 export lambda_max=1600
export FLAG_TOTAL=0 export FLAG_TOTAL=0
export GRATING_CASE="half_ellipsoid" export GRATING_CASE="half_ellipsoid"
export myDir="res_Matrix_nb_lam"$nb_lam"_nb_phi"$nb_phi"_total"$FLAG_TOTAL export myDir="res_Matrix_nb_lam"$nb_lam"_nb_phi"$nb_phi"_total"$FLAG_TOTAL
...@@ -24,7 +24,7 @@ myfunc() { ...@@ -24,7 +24,7 @@ myfunc() {
local phi=$(echo "scale=10;360/($nb_phi)*$2" | bc ) local phi=$(echo "scale=10;360/($nb_phi)*$2" | bc )
local psi=$(echo "scale=10;90*$3" | bc ) local psi=$(echo "scale=10;90*$3" | bc )
cd $mysubDir 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/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 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 if [ $3 -eq 1 ]; then cp res3D/rs.txt ../r_sin_sout_lam$1_phi$2.out ; fi
......
DiffractionGratings/grating3D_postplot_Mmatrix.py
+ 38
− 31
View file @ 4005a31a
...@@ -22,10 +22,10 @@ def load_getdp_integral(fname): ...@@ -22,10 +22,10 @@ def load_getdp_integral(fname):
else: else:
return temp[:,1]+1j*temp[:,2] return temp[:,1]+1j*temp[:,2]
nb_lam = 100 nb_lam = 200
nb_phi = 60 nb_phi = 60
FLAG_TOT = 0 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) rpp = np.zeros((nb_lam,nb_phi),dtype=complex)
rps = 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) rsp = np.zeros((nb_lam,nb_phi),dtype=complex)
...@@ -42,7 +42,7 @@ effr2_sin = np.zeros((9,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)) 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) tab_phi = np.linspace(0,360,nb_phi)
M = np.zeros((4,4,len(tab_lam),len(tab_phi)),dtype=complex) 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) ...@@ -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 # that are parallel (resp. perpendicular) to the plane of incidence
# rps : send Es_inc (s-in), project the reflected field along p (p-out) # rps : send Es_inc (s-in), project the reflected field along p (p-out)
str_reftrans = 'r'
for i in range(nb_lam): for i in range(nb_lam):
print(i)
for j in range(nb_phi-1): for j in range(nb_phi-1):
print(i,j) rpp[i,j] = load_getdp_integral(respath+str_reftrans+'_pin_pout_lam%g_phi%g.out'%(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+str_reftrans+'_sin_sout_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+str_reftrans+'_sin_pout_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+str_reftrans+'_pin_sout_lam%g_phi%g.out'%(i,j))
rsp[i,j] = load_getdp_integral(respath+'r_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)) 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)) 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)) 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)) ...@@ -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,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)) M[4-1,4-1,:,:] = np.real(rpp*np.conj(rss)-rps*np.conj(rsp))
# ### Energy balance ### Energy balance
# fig, axes = plt.subplots(2, 2, figsize=(12,12)) fig, axes = plt.subplots(3, 2, figsize=(12,12))
# L,P = np.meshgrid(tab_lam,tab_phi,indexing='ij') L,P = np.meshgrid(tab_lam,tab_phi,indexing='ij')
# for form in range(2): for form in range(2):
# for pol in range(2): for pol in range(2):
# if form==0 and pol==0: if form==0 and pol==0:
# data=TOT1_pin data=TOT1_pin
# title = 'TOT1_pin' title = 'TOT1 pin'
# if form==0 and pol==1: if form==0 and pol==1:
# data=TOT1_sin data=TOT1_sin
# title = 'TOT1_sin' title = 'TOT1 sin'
# if form==1 and pol==0: if form==1 and pol==0:
# data=TOT2_pin data=TOT2_pin
# title = 'TOT2_pin' title = 'TOT2 pin'
# if form==1 and pol==1: if form==1 and pol==1:
# data=TOT2_sin data=TOT2_sin
# title = 'TOT2_sin' title = 'TOT2 sin'
# zplot = axes[form,pol].pcolormesh(L,P,data.real,vmin=0.99,vmax=1.01) zplot = axes[form,pol].pcolormesh(L,P,data.real,vmin=0.99,vmax=1.01)
# add_colorbar(zplot) axes[form,pol].title.set_text(title)
# plt.savefig('BALANCE_FLAG_TOT%g.jpg'%FLAG_TOT) 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)) fig, axes = plt.subplots(4, 4, subplot_kw=dict(projection='polar') ,figsize=(12,12))
flag_lam = True flag_lam = True
rlabel = r"$\lambda_0$ (nm)" rlabel = r"$\lambda_0$ (nm)"
which_r = tab_lam 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 i in range(4):
for j in range(4): for j in range(4):
ax=axes[i,j] ax=axes[i,j]
...@@ -162,7 +169,7 @@ for i in range(4): ...@@ -162,7 +169,7 @@ for i in range(4):
sm.set_array([]) sm.set_array([])
ax.set_rorigin(which_orig) 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.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_%s.jpg'%str_reftrans)
plt.savefig('Mmatrix.pdf',bbox_inches='tight',pad_inches=0) # plt.savefig('Mmatrix.pdf',bbox_inches='tight',pad_inches=0)
plt.show() plt.show()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment