new loop over theta

DiffractionGratings/grating3D.geo

@@ -213,6 +213,11 @@
     Box(9) = {-period_x/2,-ry/2, hh_L_3, period_x, ry, rz};
   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;
   Call SetPBCs;

DiffractionGratings/grating3D_data_pillar.geo

@@ -4,7 +4,7 @@
 // solve wall time pastix : 966.046s - R00=0.2432916553175064
 // beware : very near a Rayleigh anomaly
 
-nm  = 1000;
+nm  = 1;
 pp1 = "1Incident Plane Wave";
 pp2 = "2Layers Thicknesses";
 pp3 = "3Scatterer Properties";
@@ -27,7 +27,7 @@ DefineConstant[
     thick_L_6     = {50   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
     xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
-    tag_geom      = { 2    , 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"}},
+    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"]},
@@ -60,10 +60,10 @@ DefineConstant[
     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          = {2      , Name StrCat[pp5,"/6Number of non specular order to output [-]"]},
+    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= {2      , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
-    refine_mesh_L_3= {6      , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]},
+    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 [-]"]},

DiffractionGratings/grating3D_parallel_Mmatrix.sh

@@ -2,50 +2,61 @@
 export OPENBLAS_NUM_THREADS=1
 export OMP_NUM_THREADS=1
 export NPROC=96
-export nb_lam=50
-export nb_phi=50
-# export nb_lam=10
-# export nb_phi=10
+
+export flag_angle_study="theta"
+export loop_angle_max=50
+export fixed_angle=90
+export GRATING_CASE="pillar"
+
+# export flag_angle_study="phi"
+# export loop_angle_max=360
+# export fixed_angle=50
+# export GRATING_CASE="half_ellipsoid"
+
+export nb_angle=2
+export nb_lam=2
 export lambda_min=400
 export lambda_max=1200
 export FLAG_TOTAL=0
-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_"$flag_angle_study$nb_angle"_total"$FLAG_TOTAL
 
 rm -r $myDir
 mkdir $myDir
 gmsh grating3D.geo -setstring test_case $GRATING_CASE -3 -o grating3D.msh
 
 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
     cp grating3D.msh grating3D.pro grating3D_data_$GRATING_CASE.geo grating3D_data.geo grating3D_materials.pro $mysubDir
-    if (( $nb_lam == 1 ))
-    then
-        local lam=$(echo "scale=10;$lambda_min" | bc )
+    local lam=$(echo "scale=10;$lambda_min+($lambda_max-$lambda_min)/($nb_lam-1)*$1" | bc )
+    if [ $flag_angle_study = "phi" ]; then
+        local phi=$(echo "scale=10;$loop_angle_max/($nb_angle-1)*$2" | bc )
+        local theta=$(echo "scale=10;$fixed_angle" | bc )
     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
-    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 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
-    if [ $3 -eq 1 ]; then cp res3D/rp.txt ../r_sin_pout_lam$1_phi$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/tp.txt ../t_pin_pout_lam$1_phi$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/tp.txt ../t_sin_pout_lam$1_phi$2.out ; fi
-    cp res3D/eff_t1.txt ../eff_t1_lam$1_phi$2_psi$3.out
-    cp res3D/eff_r1.txt ../eff_r1_lam$1_phi$2_psi$3.out
-    cp res3D/eff_t2.txt ../eff_t2_lam$1_phi$2_psi$3.out
-    cp res3D/eff_r2.txt ../eff_r2_lam$1_phi$2_psi$3.out
-    cp res3D/temp-Q_scat.txt ../Q_scat_lam$1_phi$2_psi$3.out
+    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_$flag_angle_study$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_$flag_angle_study$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_$flag_angle_study$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_$flag_angle_study$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_$flag_angle_study$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_$flag_angle_study$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_$flag_angle_study$2_psi$3.out
     cd ../..
     rm -r $mysubDir
 }
 
 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

@@ -2,6 +2,29 @@ import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib as mpl
 import scipy.constants as scc
+import argparse
+
+parser = argparse.ArgumentParser()
+parser.add_argument("-flag_angle_study" , help="flag_angle_study" , type=str   , default='phi')
+parser.add_argument("-nb_angle"         , help="nb_angle"         , type=int   , default=2)
+parser.add_argument("-loop_angle_max"   , help="loop_angle_max"   , type=float , default=360)
+parser.add_argument("-fixed_angle"      , help="fixed_angle"      , type=float , default=50)
+parser.add_argument("-nb_lam"           , help="nb_lam"           , type=int   , default=2)
+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
+
+
 plt.rcParams.update({"text.usetex": True, "font.family": "serif"})
 
 def add_colorbar(mappable):
@@ -23,14 +46,7 @@ def load_getdp_integral(fname):
     else:
         return temp[:,1]+1j*temp[:,2]
 
-nb_lam = 50
-nb_angle = 51
-FLAG_TOT = 0
-
-flag_study = 'theta'
-# flag_study = 'phi'
-
-respath  = 'res_Matrix_nb_lam%g_nb_%s%g_total0/'%(nb_lam,flag_study,nb_angle)
+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)
@@ -47,13 +63,8 @@ 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))
 
-
-lambda_min = 400
-lambda_max = 1200
-angledeg_max = 50
-
 tab_lam   = np.linspace(lambda_min,lambda_max,nb_lam)
-tab_angle = np.linspace(0,angledeg_max,nb_angle)
+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)
 
@@ -71,26 +82,20 @@ str_reftrans = 'r'
 for i in range(nb_lam):
     print(i)
     for j in range(nb_angle):
-        rpp[i,j] = load_getdp_integral(respath+str_reftrans+'_pin_pout_lam%g_%s%g.out'%(i,flag_study,j))
-        rss[i,j] = load_getdp_integral(respath+str_reftrans+'_sin_sout_lam%g_%s%g.out'%(i,flag_study,j))
-        rps[i,j] = load_getdp_integral(respath+str_reftrans+'_sin_pout_lam%g_%s%g.out'%(i,flag_study,j))
-        rsp[i,j] = load_getdp_integral(respath+str_reftrans+'_pin_sout_lam%g_%s%g.out'%(i,flag_study,j))
-        efft1_pin[:,i,j] = load_getdp_integral(respath+'eff_t1_lam%g_%s%g_psi0.out'%(i,flag_study,j))
-        efft2_pin[:,i,j] = load_getdp_integral(respath+'eff_t2_lam%g_%s%g_psi0.out'%(i,flag_study,j))
-        effr1_pin[:,i,j] = load_getdp_integral(respath+'eff_r1_lam%g_%s%g_psi0.out'%(i,flag_study,j))
-        effr2_pin[:,i,j] = load_getdp_integral(respath+'eff_r2_lam%g_%s%g_psi0.out'%(i,flag_study,j))
-        Qscat_pin[i,j]   = np.real(load_getdp_integral(respath+'Q_scat_lam%g_%s%g_psi0.out'%(i,flag_study,j)))
-        efft1_sin[:,i,j] = load_getdp_integral(respath+'eff_t1_lam%g_%s%g_psi1.out'%(i,flag_study,j))
-        efft2_sin[:,i,j] = load_getdp_integral(respath+'eff_t2_lam%g_%s%g_psi1.out'%(i,flag_study,j))
-        effr1_sin[:,i,j] = load_getdp_integral(respath+'eff_r1_lam%g_%s%g_psi1.out'%(i,flag_study,j))
-        effr2_sin[:,i,j] = load_getdp_integral(respath+'eff_r2_lam%g_%s%g_psi1.out'%(i,flag_study,j))
-        Qscat_sin[i,j]   = np.real(load_getdp_integral(respath+'Q_scat_lam%g_%s%g_psi1.out'%(i,flag_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]
+        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_%s%g.out'%(i,flag_angle_study,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_%s%g.out'%(i,flag_angle_study,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_%s%g_psi0.out'%(i,flag_angle_study,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_%s%g_psi0.out'%(i,flag_angle_study,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_%s%g_psi1.out'%(i,flag_angle_study,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_%s%g_psi1.out'%(i,flag_angle_study,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_%s%g_psi1.out'%(i,flag_angle_study,j)))
 
 T1_sin = np.sum(efft1_sin,axis=0)
 T2_sin = np.sum(efft2_sin,axis=0)
@@ -167,7 +172,7 @@ plt.savefig('BALANCE_new_code.jpg')
 fig, axes = plt.subplots(4, 4, subplot_kw=dict(projection='polar') ,figsize=(12,12))
 flag_lam = False
 rlabel = r"$\hbar \nu$" if not flag_lam else r"$\lambda_0$"
-anglelabel = r"$\varphi_0$" if flag_study=='phi' else r"$\theta_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_orig = 0  #if not flag_lam else 350
 for i in range(4):
@@ -198,7 +203,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_%s_%s.jpg'%(flag_study,str_reftrans))
+plt.savefig('Mmatrix_%s_%s.jpg'%(flag_angle_study,str_reftrans))
 # plt.savefig('Mmatrix.pdf',bbox_inches='tight',pad_inches=0)
-plt.show()
+# plt.show()