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

Merge branch 'pillar' into 'master' 

pillar + check Mmatrix again

See merge request !7
parents e50893d1 93ff16c1
No related branches found
No related tags found
1 merge request!7pillar + check Mmatrix again
Pipeline #10540 passed
Stage: test
Hide whitespace changes
Inline Side-by-side
DiffractionGratings/grating3D.geo
+ 5
− 0
View file @ 4d43092a
...@@ -213,6 +213,11 @@ ...@@ -213,6 +213,11 @@
Box(9) = {-period_x/2,-ry/2, hh_L_3, period_x, ry, rz}; Box(9) = {-period_x/2,-ry/2, hh_L_3, period_x, ry, rz};
EndIf EndIf
If (tag_geom==8)
Cylinder(9) = {0,0,hh_L_3 , 0,0,rz , rx};
Dilate { { 0,0,hh_L_3 }, { 1, ry/rx, 1 } } { Volume{9}; }
EndIf
Coherence; Coherence;
Call SetPBCs; Call SetPBCs;
......
DiffractionGratings/grating3D_data_pillar.geo 0 → 100644
+ 131
− 0
View file @ 4d43092a
// note : invcreasing paramaille to 7 leads to:
// Ndofs 1116912
// solve wall time mumps : 634.318s - R00=0.2432916553175037
// solve wall time pastix : 966.046s - R00=0.2432916553175064
// beware : very near a Rayleigh anomaly
nm = 1;
pp1 = "1Incident Plane Wave";
pp2 = "2Layers Thicknesses";
pp3 = "3Scatterer Properties";
pp4 = "4Layer Materials";
pp5 = "5Computational Parameters";
pp6 = "6Output";
DefineConstant[
FLAG_TOTAL = {0 , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
lambda0 = {500 , Name StrCat[pp1,"/1lambda0 [nm]"]},
thetadeg = {0 , Name StrCat[pp1,"/2theta0 [deg]"]},
phideg = {0 , Name StrCat[pp1,"/3phi0 [deg]"]},
psideg = {45 , Name StrCat[pp1,"/4psi0 [deg]"]},
period_x = {390 , Name StrCat[pp2,"/1X period [nm]"]},
period_y = {390 , Name StrCat[pp2,"/2Y period [nm]"]},
thick_L_1 = {50 , Name StrCat[pp2,"/3thickness layer 1 [nm] (superstrate)"]},
thick_L_2 = {50 , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
thick_L_3 = {200 , Name StrCat[pp2,"/5thickness layer 3 [nm]"]},
thick_L_4 = {50 , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
thick_L_5 = {50 , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
thick_L_6 = {50 , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
xsideg = {0 , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
tag_geom = { 8 , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="U",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar",8="cylindrical pillar"}},
rx = {60 , Name StrCat[pp3,"/1rx"]},
ry = {60 , Name StrCat[pp3,"/2ry"]},
rz = {130 , Name StrCat[pp3,"/3rz"]},
flag_mat_scat = {10 , Name StrCat[pp3,"/4Scatterer permittivity model"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
eps_re_Scat = {1 , Name StrCat[pp3,"/7eps_re_Scat"]},
eps_im_Scat = {0 , Name StrCat[pp3,"/8eps_im_Scat"]},
flag_mat_1 = { 0 , Name StrCat[pp4,"/1Layer 1"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
flag_mat_2 = { 0 , Name StrCat[pp4,"/2Layer 2"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
flag_mat_3 = { 0 , Name StrCat[pp4,"/3Layer 3"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
flag_mat_4 = { 0 , Name StrCat[pp4,"/4Layer 4"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
flag_mat_5 = { 0 , Name StrCat[pp4,"/5Layer 5"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
flag_mat_6 = { 0 , Name StrCat[pp4,"/6Layer 6"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
eps_re_L_1 = {1 , Name StrCat[pp4,"/layer 1: real part of relative permittivity"]},
eps_im_L_1 = {0 , Name StrCat[pp4,"/layer 1: imag part of relative permittivity"]},
eps_re_L_2 = {1 , Name StrCat[pp4,"/layer 2: real part of relative permittivity"]},
eps_im_L_2 = {0 , Name StrCat[pp4,"/layer 2: imag part of relative permittivity"]},
eps_re_L_3 = {1 , Name StrCat[pp4,"/layer 3: real part of relative permittivity"]},
eps_im_L_3 = {0 , Name StrCat[pp4,"/layer 3: imag part of relative permittivity"]},
eps_re_L_4 = {2.25 , Name StrCat[pp4,"/layer 4: real part of relative permittivity"]},
eps_im_L_4 = {0 , Name StrCat[pp4,"/layer 4: imag part of relative permittivity"]},
eps_re_L_5 = {2.25 , Name StrCat[pp4,"/layer 5: real part of relative permittivity"]},
eps_im_L_5 = {0 , Name StrCat[pp4,"/layer 5: imag part of relative permittivity"]},
eps_re_L_6 = {2.25 , Name StrCat[pp4,"/layer 6: real part of relative permittivity"]},
eps_im_L_6 = {0 , Name StrCat[pp4,"/layer 6: imag part of relative permittivity"]},
og = {0 , Name StrCat[pp5,"/0geometrical order [-]"] , Choices {0="1",1="2"} },
oi = {1 , Name StrCat[pp5,"/0interpolation order [-]"], Choices {0="1",1="2"} },
paramaille = {5. , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
lc_scat = {lambda0/(5*paramaille) , Name StrCat[pp5,"/2Scatterer absolute mesh size [nm]"]},
PML_top = {lambda0, Name StrCat[pp5,"/4PML top thickness [nm]"]},
PML_bot = {lambda0, Name StrCat[pp5,"/5PML bot thickness [nm]"]},
Nmax = {1 , Name StrCat[pp5,"/6Number of non specular order to output [-]"]},
refine_mesh_L_1= {1 , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]},
refine_mesh_L_2= {1 , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
refine_mesh_L_3= {1 , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]},
refine_mesh_L_4= {1 , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
refine_mesh_L_5= {1 , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
refine_mesh_L_6= {1 , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
FlagLinkFacets = {0 , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
PML_TYPE = {0 , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
InterpSampling = { 30 , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
Flag_interp_cubic = { 1 , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
FlagOutEtotCuts = { 1 , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
FlagOutHtotCuts = { 0 , Name StrCat[pp6,"/3Output Total Magnetic Field cuts?"] , Choices {0,1} },
FlagOutEscaCuts = { 1 , Name StrCat[pp6,"/4Output Scattered Electric Field cuts?"] , Choices {0,1} },
FlagOutPoyCut = { 1 , Name StrCat[pp6,"/5Output Poynting cuts?"] , Choices {0,1} },
FlagOutEtotFull = { 0 , Name StrCat[pp6,"/6Total Electric Field Full Output?"] , Choices {0,1} },
FlagOutEscaFull = { 0 , Name StrCat[pp6,"/7Scattered Electric Field Full Output?"] , Choices {0,1} },
FlagOutPoyFull = { 0 , Name StrCat[pp6,"/8Poynting Full Output?"] , Choices {0,1} }
];
lambda0 = nm * lambda0;
period_x = nm * period_x;
period_y = nm * period_y;
thick_L_1 = nm * thick_L_1;
thick_L_2 = nm * thick_L_2;
thick_L_3 = nm * thick_L_3;
thick_L_4 = nm * thick_L_4;
thick_L_5 = nm * thick_L_5;
thick_L_6 = nm * thick_L_6;
rx = nm * rx;
ry = nm * ry;
rz = nm * rz;
lc_scat = nm * lc_scat;
PML_top = nm * PML_top;
PML_bot = nm * PML_bot;
InterpSampling= nm * InterpSampling;
lambda_m = lambda0;
og+=1;
oi+=1;
hh_L_6 = -thick_L_6;
For k In {1:5}
hh_L~{6-k} = hh_L~{7-k}+thick_L~{7-k};
EndFor
PML_bot_hh = hh_L_6-PML_bot;
PML_top_hh = hh_L_1+thick_L_1;
hh_L_7 = PML_bot_hh;
hh_L_0 = PML_top_hh;
thick_L_7 = PML_bot;
thick_L_0 = PML_top;
theta0 = thetadeg*Pi/180;
phi0 = phideg*Pi/180;
psi0 = psideg*Pi/180;
xsi = xsideg*Pi/180;
dyc = period_y*Cos[xsi];
dys = period_y*Sin[xsi];
DomainZsizeSca = PML_top_hh+PML_bot-(hh_L_6-PML_bot);
DomainZsizeTot = PML_top_hh-hh_L_6;
npts_interpX = period_x/InterpSampling;
npts_interpY = period_y/InterpSampling;
npts_checkpoyX = 50;
npts_checkpoyY = 50;
npts_interpZSca = DomainZsizeSca/InterpSampling;
npts_interpZTot = DomainZsizeTot/InterpSampling;
DiffractionGratings/grating3D_parallel_Mmatrix.sh
+ 42
− 26
View file @ 4d43092a
...@@ -2,50 +2,66 @@ ...@@ -2,50 +2,66 @@
export OPENBLAS_NUM_THREADS=1 export OPENBLAS_NUM_THREADS=1
export OMP_NUM_THREADS=1 export OMP_NUM_THREADS=1
export NPROC=96 export NPROC=96
### choose configuration
### theta loop study for the "U" case
### config 1
export flag_angle_study="theta"
export loop_angle_max=50
export fixed_angle=90
export GRATING_CASE="U"
### phi loop study for the "U" case
### config 2
# export flag_angle_study="phi"
# export loop_angle_max=360
# export fixed_angle=50
# export GRATING_CASE="U"
export nb_angle=40
export nb_lam=50 export nb_lam=50
export nb_phi=50
# export nb_lam=10
# export nb_phi=10
export lambda_min=400 export lambda_min=400
export lambda_max=1200 export lambda_max=1200
export FLAG_TOTAL=0 export FLAG_TOTAL=0
export GRATING_CASE="half_ellipsoid" export myDir="res_Matrix_nb_lam"$nb_lam"_nb_"$flag_angle_study$nb_angle"_total"$FLAG_TOTAL
export myDir="res_Matrix_nb_lam"$nb_lam"_nb_phi"$nb_phi"_total"$FLAG_TOTAL
rm -r $myDir rm -r $myDir
mkdir $myDir mkdir $myDir
gmsh grating3D.geo -setstring test_case $GRATING_CASE -3 -o grating3D.msh gmsh grating3D.geo -setstring test_case $GRATING_CASE -3 -o grating3D.msh
myfunc() { myfunc() {
local mysubDir=$myDir/run_lam$1_phi$2_psi$3 local mysubDir=$myDir/run_lam$1_$flag_angle_study$2_psi$3
mkdir $mysubDir mkdir $mysubDir
cp grating3D.msh grating3D.pro grating3D_data_$GRATING_CASE.geo grating3D_data.geo grating3D_materials.pro $mysubDir cp grating3D.msh grating3D.pro grating3D_data_$GRATING_CASE.geo grating3D_data.geo grating3D_materials.pro $mysubDir
if (( $nb_lam == 1 )) local lam=$(echo "scale=10;$lambda_min+($lambda_max-$lambda_min)/($nb_lam-1)*$1" | bc )
then if [ $flag_angle_study = "phi" ]; then
local lam=$(echo "scale=10;$lambda_min" | bc ) local phi=$(echo "scale=10;$loop_angle_max/($nb_angle-1)*$2" | bc )
local theta=$(echo "scale=10;$fixed_angle" | bc )
else else
local lam=$(echo "scale=10;$lambda_min+($lambda_max-$lambda_min)/($nb_lam-1)*$1" | bc ) local theta=$(echo "scale=10;$loop_angle_max/($nb_angle-1)*$2" | bc )
local phi=$(echo "scale=10;$fixed_angle" | bc )
fi fi
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 50 -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 $theta -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_$flag_angle_study$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_$flag_angle_study$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_$flag_angle_study$2.out ; fi
if [ $3 -eq 1 ]; then cp res3D/rp.txt ../r_sin_pout_lam$1_phi$2.out ; fi if [ $3 -eq 1 ]; then cp res3D/rp.txt ../r_sin_pout_lam$1_$flag_angle_study$2.out ; fi
if [ $3 -eq 0 ]; then cp res3D/ts.txt ../t_pin_sout_lam$1_phi$2.out ; fi if [ $3 -eq 0 ]; then cp res3D/ts.txt ../t_pin_sout_lam$1_$flag_angle_study$2.out ; fi
if [ $3 -eq 0 ]; then cp res3D/tp.txt ../t_pin_pout_lam$1_phi$2.out ; fi if [ $3 -eq 0 ]; then cp res3D/tp.txt ../t_pin_pout_lam$1_$flag_angle_study$2.out ; fi
if [ $3 -eq 1 ]; then cp res3D/ts.txt ../t_sin_sout_lam$1_phi$2.out ; fi if [ $3 -eq 1 ]; then cp res3D/ts.txt ../t_sin_sout_lam$1_$flag_angle_study$2.out ; fi
if [ $3 -eq 1 ]; then cp res3D/tp.txt ../t_sin_pout_lam$1_phi$2.out ; fi if [ $3 -eq 1 ]; then cp res3D/tp.txt ../t_sin_pout_lam$1_$flag_angle_study$2.out ; fi
cp res3D/eff_t1.txt ../eff_t1_lam$1_phi$2_psi$3.out cp res3D/eff_t1.txt ../eff_t1_lam$1_$flag_angle_study$2_psi$3.out
cp res3D/eff_r1.txt ../eff_r1_lam$1_phi$2_psi$3.out cp res3D/eff_r1.txt ../eff_r1_lam$1_$flag_angle_study$2_psi$3.out
cp res3D/eff_t2.txt ../eff_t2_lam$1_phi$2_psi$3.out cp res3D/eff_t2.txt ../eff_t2_lam$1_$flag_angle_study$2_psi$3.out
cp res3D/eff_r2.txt ../eff_r2_lam$1_phi$2_psi$3.out cp res3D/eff_r2.txt ../eff_r2_lam$1_$flag_angle_study$2_psi$3.out
cp res3D/temp-Q_scat.txt ../Q_scat_lam$1_phi$2_psi$3.out cp res3D/temp-Q_scat.txt ../Q_scat_lam$1_$flag_angle_study$2_psi$3.out
cd ../.. cd ../..
rm -r $mysubDir rm -r $mysubDir
} }
export -f myfunc export -f myfunc
parallel -j $NPROC myfunc ::: $(seq 0 $(($nb_lam-1))) ::: $(seq 0 $(($nb_phi-1))) ::: 0 1 parallel -j $NPROC myfunc ::: $(seq 0 $(($nb_lam-1))) ::: $(seq 0 $(($nb_angle-1))) ::: 0 1
python grating3D_postplot_Mmatrix.py -flag_angle_study $flag_angle_study -nb_angle $nb_angle -loop_angle_max $loop_angle_max -fixed_angle $fixed_angle -nb_lam $nb_lam -lambda_min $lambda_min -lambda_max $lambda_max
DiffractionGratings/grating3D_postplot_Mmatrix.py
+ 66
− 57
View file @ 4d43092a
...@@ -2,6 +2,7 @@ import numpy as np ...@@ -2,6 +2,7 @@ import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib as mpl import matplotlib as mpl
import scipy.constants as scc import scipy.constants as scc
import argparse
plt.rcParams.update({"text.usetex": True, "font.family": "serif"}) plt.rcParams.update({"text.usetex": True, "font.family": "serif"})
def add_colorbar(mappable): def add_colorbar(mappable):
...@@ -23,31 +24,48 @@ def load_getdp_integral(fname): ...@@ -23,31 +24,48 @@ def load_getdp_integral(fname):
else: else:
return temp[:,1]+1j*temp[:,2] return temp[:,1]+1j*temp[:,2]
nb_lam = 50 parser = argparse.ArgumentParser()
nb_phi = 51 parser.add_argument("-flag_angle_study" , help="flag_angle_study" , type=str , default='theta')
parser.add_argument("-nb_angle" , help="nb_angle" , type=int , default=20)
parser.add_argument("-loop_angle_max" , help="loop_angle_max" , type=float , default=50)
parser.add_argument("-fixed_angle" , help="fixed_angle" , type=float , default=90)
parser.add_argument("-nb_lam" , help="nb_lam" , type=int , default=100)
parser.add_argument("-lambda_min" , help="lambda_min" , type=float , default=400)
parser.add_argument("-lambda_max" , help="lambda_max" , type=float , default=1200)
args = parser.parse_args()
flag_angle_study = args.flag_angle_study
nb_angle = args.nb_angle
loop_angle_max = args.loop_angle_max
fixed_angle = args.fixed_angle
nb_lam = args.nb_lam
lambda_min = args.lambda_min
lambda_max = args.lambda_max
FLAG_TOT = 0 FLAG_TOT = 0
respath = 'res_Matrix_nb_lam%g_nb_phi%g_total0/'%(nb_lam,nb_phi-1) str_reftrans = 'r'
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)
rss = np.zeros((nb_lam,nb_phi),dtype=complex)
efft1_pin = np.zeros((9,nb_lam,nb_phi),dtype=complex)
efft2_pin = np.zeros((9,nb_lam,nb_phi),dtype=complex)
effr1_pin = np.zeros((9,nb_lam,nb_phi),dtype=complex)
effr2_pin = np.zeros((9,nb_lam,nb_phi),dtype=complex)
Qscat_pin = np.zeros((nb_lam,nb_phi))
efft1_sin = np.zeros((9,nb_lam,nb_phi),dtype=complex)
efft2_sin = np.zeros((9,nb_lam,nb_phi),dtype=complex)
effr1_sin = np.zeros((9,nb_lam,nb_phi),dtype=complex)
effr2_sin = np.zeros((9,nb_lam,nb_phi),dtype=complex)
Qscat_sin = np.zeros((nb_lam,nb_phi))
tab_lam = np.linspace(495,1200,nb_lam)
tab_phi = np.linspace(0,360,nb_phi)
tab_hnu = 2*scc.pi*scc.c/(tab_lam*1e-9/scc.hbar*scc.eV)
M = np.zeros((4,4,len(tab_lam),len(tab_phi)),dtype=complex)
respath = 'res_Matrix_nb_lam%g_nb_%s%g_total0/'%(nb_lam,flag_angle_study,nb_angle)
rpp = np.zeros((nb_lam,nb_angle),dtype=complex)
rps = np.zeros((nb_lam,nb_angle),dtype=complex)
rsp = np.zeros((nb_lam,nb_angle),dtype=complex)
rss = np.zeros((nb_lam,nb_angle),dtype=complex)
efft1_pin = np.zeros((9,nb_lam,nb_angle),dtype=complex)
efft2_pin = np.zeros((9,nb_lam,nb_angle),dtype=complex)
effr1_pin = np.zeros((9,nb_lam,nb_angle),dtype=complex)
effr2_pin = np.zeros((9,nb_lam,nb_angle),dtype=complex)
Qscat_pin = np.zeros((nb_lam,nb_angle))
efft1_sin = np.zeros((9,nb_lam,nb_angle),dtype=complex)
efft2_sin = np.zeros((9,nb_lam,nb_angle),dtype=complex)
effr1_sin = np.zeros((9,nb_lam,nb_angle),dtype=complex)
effr2_sin = np.zeros((9,nb_lam,nb_angle),dtype=complex)
Qscat_sin = np.zeros((nb_lam,nb_angle))
tab_lam = np.linspace(lambda_min,lambda_max,nb_lam)
tab_angle = np.linspace(0,loop_angle_max,nb_angle)
tab_hnu = 2*scc.pi*scc.c/(tab_lam*1e-9/scc.hbar*scc.eV)
M = np.zeros((4,4,nb_lam,nb_angle),dtype=complex)
### Convention used e.g. in https://doi.org/10.1364/JOSAB.36.000E78 ### Convention used e.g. in https://doi.org/10.1364/JOSAB.36.000E78
# [Ep] [ rpp rps ] [Ep] # [Ep] [ rpp rps ] [Ep]
...@@ -57,31 +75,23 @@ M = np.zeros((4,4,len(tab_lam),len(tab_phi)),dtype=complex) ...@@ -57,31 +75,23 @@ 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) print(i)
for j in range(nb_phi-1): for j in range(nb_angle):
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+str_reftrans+'_pin_pout_lam%g_%s%g.out'%(i,flag_angle_study,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+str_reftrans+'_sin_sout_lam%g_%s%g.out'%(i,flag_angle_study,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+str_reftrans+'_sin_pout_lam%g_%s%g.out'%(i,flag_angle_study,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+str_reftrans+'_pin_sout_lam%g_%s%g.out'%(i,flag_angle_study,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_%s%g_psi0.out'%(i,flag_angle_study,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_%s%g_psi0.out'%(i,flag_angle_study,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_%s%g_psi0.out'%(i,flag_angle_study,j))
effr2_pin[:,i,j] = load_getdp_integral(respath+'eff_r2_lam%g_phi%g_psi0.out'%(i,j)) effr2_pin[:,i,j] = load_getdp_integral(respath+'eff_r2_lam%g_%s%g_psi0.out'%(i,flag_angle_study,j))
Qscat_pin[i,j] = np.real(load_getdp_integral(respath+'Q_scat_lam%g_phi%g_psi0.out'%(i,j))) Qscat_pin[i,j] = np.real(load_getdp_integral(respath+'Q_scat_lam%g_%s%g_psi0.out'%(i,flag_angle_study,j)))
efft1_sin[:,i,j] = load_getdp_integral(respath+'eff_t1_lam%g_phi%g_psi1.out'%(i,j)) efft1_sin[:,i,j] = load_getdp_integral(respath+'eff_t1_lam%g_%s%g_psi1.out'%(i,flag_angle_study,j))
efft2_sin[:,i,j] = load_getdp_integral(respath+'eff_t2_lam%g_phi%g_psi1.out'%(i,j)) efft2_sin[:,i,j] = load_getdp_integral(respath+'eff_t2_lam%g_%s%g_psi1.out'%(i,flag_angle_study,j))
effr1_sin[:,i,j] = load_getdp_integral(respath+'eff_r1_lam%g_phi%g_psi1.out'%(i,j)) effr1_sin[:,i,j] = load_getdp_integral(respath+'eff_r1_lam%g_%s%g_psi1.out'%(i,flag_angle_study,j))
effr2_sin[:,i,j] = load_getdp_integral(respath+'eff_r2_lam%g_phi%g_psi1.out'%(i,j)) effr2_sin[:,i,j] = load_getdp_integral(respath+'eff_r2_lam%g_%s%g_psi1.out'%(i,flag_angle_study,j))
Qscat_sin[i,j] = np.real(load_getdp_integral(respath+'Q_scat_lam%g_phi%g_psi1.out'%(i,j))) Qscat_sin[i,j] = np.real(load_getdp_integral(respath+'Q_scat_lam%g_%s%g_psi1.out'%(i,flag_angle_study,j)))
## replicate data for phi=2\pi from phi=0
rpp[:,-1] = rpp[:,0]
rss[:,-1] = rss[:,0]
rps[:,-1] = rps[:,0]
rsp[:,-1] = rsp[:,0]
T1_sin = np.sum(efft1_sin,axis=0) T1_sin = np.sum(efft1_sin,axis=0)
T2_sin = np.sum(efft2_sin,axis=0) T2_sin = np.sum(efft2_sin,axis=0)
...@@ -119,7 +129,7 @@ M[4-1,4-1,:,:] = np.real(rpp*np.conj(rss)-rps*np.conj(rsp)) ...@@ -119,7 +129,7 @@ M[4-1,4-1,:,:] = np.real(rpp*np.conj(rss)-rps*np.conj(rsp))
### Energy balance ### Energy balance
fig, axes = plt.subplots(3, 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_angle,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:
...@@ -144,7 +154,7 @@ for form in range(2): ...@@ -144,7 +154,7 @@ for form in range(2):
plt.savefig('BALANCE_FLAG_TOT%g.jpg'%FLAG_TOT) plt.savefig('BALANCE_FLAG_TOT%g.jpg'%FLAG_TOT)
fig, axes = plt.subplots(3, 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_angle,indexing='ij')
ax=axes[0,0];zplot = ax.pcolormesh(L,P,T00_sin.real );add_colorbar(zplot);ax.title.set_text('T00 sin') ax=axes[0,0];zplot = ax.pcolormesh(L,P,T00_sin.real );add_colorbar(zplot);ax.title.set_text('T00 sin')
ax=axes[1,0];zplot = ax.pcolormesh(L,P,R00_sin.real );add_colorbar(zplot);ax.title.set_text('R00 sin') ax=axes[1,0];zplot = ax.pcolormesh(L,P,R00_sin.real );add_colorbar(zplot);ax.title.set_text('R00 sin')
ax=axes[2,0];zplot = ax.pcolormesh(L,P,Qscat_sin.real);add_colorbar(zplot);ax.title.set_text('Abs sin') ax=axes[2,0];zplot = ax.pcolormesh(L,P,Qscat_sin.real);add_colorbar(zplot);ax.title.set_text('Abs sin')
...@@ -153,19 +163,18 @@ ax=axes[1,1];zplot = ax.pcolormesh(L,P,R00_pin.real );add_colorbar(zplot);ax.ti ...@@ -153,19 +163,18 @@ ax=axes[1,1];zplot = ax.pcolormesh(L,P,R00_pin.real );add_colorbar(zplot);ax.ti
ax=axes[2,1];zplot = ax.pcolormesh(L,P,Qscat_pin.real);add_colorbar(zplot);ax.title.set_text('Abs pin') ax=axes[2,1];zplot = ax.pcolormesh(L,P,Qscat_pin.real);add_colorbar(zplot);ax.title.set_text('Abs pin')
plt.savefig('BALANCE_new_code.jpg') plt.savefig('BALANCE_new_code.jpg')
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 = False flag_lam = True
rlabel = r"$\hbar \nu$" if not flag_lam else r"$\lambda_0$" rlabel = r"$\hbar \nu$" if not flag_lam else r"$\lambda_0$"
anglelabel = r"$\varphi_0$" if flag_angle_study=='phi' else r"$\theta_0$"
which_r = tab_hnu if not flag_lam else tab_lam which_r = tab_hnu if not flag_lam else 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]
if i!=0 or j!=0: if i!=0 or j!=0:
sm=ax.contourf(tab_phi*np.pi/180,which_r,M[i,j]/M[0,0],cmap=plt.cm.bwr,levels=30,vmin=-1,vmax=1) # sm=ax.contourf(tab_angle*np.pi/180,which_r,M[i,j]/M[0,0],cmap=plt.cm.bwr,levels=30,vmin=-1,vmax=1)
# sm=ax.contourf(tab_phi*np.pi/180,which_r,M[i,j]/M[0,0],cmap=plt.cm.bwr,levels=30) sm=ax.contourf(tab_angle*np.pi/180,which_r,M[i,j]/M[0,0],cmap=plt.cm.bwr,levels=30)
ax.text(0. , 1. , r"$M_{%d%d}$"%(i+1,j+1) , fontsize=16 , transform=ax.transAxes) ax.text(0. , 1. , r"$M_{%d%d}$"%(i+1,j+1) , fontsize=16 , transform=ax.transAxes)
ax.set_xticks([]) ax.set_xticks([])
ax.set_yticks([]) ax.set_yticks([])
...@@ -173,14 +182,14 @@ for i in range(4): ...@@ -173,14 +182,14 @@ for i in range(4):
# cbar.ax.locator_params(nbins=5) # cbar.ax.locator_params(nbins=5)
# cbar.ax.tick_params(labelsize=14) # cbar.ax.tick_params(labelsize=14)
else: else:
p00 = ax.contourf(tab_phi*np.pi/180,which_r,M[i,j]/M[0,0]-1,cmap=plt.cm.bwr,vmin=-1,vmax=1) p00 = ax.contourf(tab_angle*np.pi/180,which_r,M[i,j]/M[0,0]-1,cmap=plt.cm.bwr,vmin=-1,vmax=1)
ax.xaxis.label.set_color('C0') #setting up X-axis label color to yellow ax.xaxis.label.set_color('C0') #setting up X-axis label color to yellow
ax.yaxis.label.set_color('C3') #setting up Y-axis label color to blue ax.yaxis.label.set_color('C3') #setting up Y-axis label color to blue
ax.tick_params(axis='x', colors='C0') #setting up X-axis tick color to red ax.tick_params(axis='x', colors='C0') #setting up X-axis tick color to red
ax.tick_params(axis='y', colors='C3') #setting up Y-axis tick color to black ax.tick_params(axis='y', colors='C3') #setting up Y-axis tick color to black
ax.set_rlabel_position(70) ax.set_rlabel_position(70)
ax.tick_params(axis='both', which='major', labelsize=14) ax.tick_params(axis='both', which='major', labelsize=14)
ax.text(np.radians(22.5),tab_lam.max()*1.03,r"$\varphi_0$",fontsize=16,color='C0') ax.text(np.radians(22.5),tab_lam.max()*1.03,anglelabel,fontsize=16,color='C0')
ax.text(np.radians(90),ax.get_rmax()/1.3,rlabel, ax.text(np.radians(90),ax.get_rmax()/1.3,rlabel,
rotation=70,ha='center',va='center',fontsize=16,color='C3') rotation=70,ha='center',va='center',fontsize=16,color='C3')
norm = mpl.colors.Normalize(vmin=-1,vmax=1) norm = mpl.colors.Normalize(vmin=-1,vmax=1)
...@@ -188,7 +197,7 @@ for i in range(4): ...@@ -188,7 +197,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_%s.jpg'%str_reftrans) plt.savefig('Mmatrix_%s_%s.jpg'%(flag_angle_study,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()
DiffractionGratings/gratings_clean_all.sh
+ 10
− 2
View file @ 4d43092a
rm -rf res2D* rm -r res*
rm -rf res3D* rm -r res2D*
rm -r res3D*
rm *.msh rm *.msh
rm *.pre rm *.pre
rm *.db rm *.db
rm *.pdf rm *.pdf
rm *.jpg
rm *.npz rm *.npz
rm *.out
rm *.tgz
rm *.*~
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment