From 05835d22c88f57e78ba4aa16b0293bf082816d1a Mon Sep 17 00:00:00 2001
From: Guillaume Demesy <guillaume.demesy@fresnel.fr>
Date: Tue, 14 Feb 2023 17:42:16 +0100
Subject: [PATCH] choose theta or phi plots
---
.../grating3D_postplot_Mmatrix.py | 108 ++++++++++--------
1 file changed, 59 insertions(+), 49 deletions(-)
diff --git a/DiffractionGratings/grating3D_postplot_Mmatrix.py b/DiffractionGratings/grating3D_postplot_Mmatrix.py
index 96f1697..5b7b783 100644
--- a/DiffractionGratings/grating3D_postplot_Mmatrix.py
+++ b/DiffractionGratings/grating3D_postplot_Mmatrix.py
@@ -24,29 +24,38 @@ def load_getdp_integral(fname):
return temp[:,1]+1j*temp[:,2]
nb_lam = 50
-nb_phi = 51
+nb_angle = 51
FLAG_TOT = 0
-respath = 'res_Matrix_nb_lam%g_nb_phi%g_total0/'%(nb_lam,nb_phi-1)
-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)
+
+flag_study = 'theta'
+# flag_study = 'phi'
+
+respath = 'res_Matrix_nb_lam%g_nb_%s%g_total0/'%(nb_lam,flag_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))
+
+
+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_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)
+M = np.zeros((4,4,nb_lam,nb_angle),dtype=complex)
### Convention used e.g. in https://doi.org/10.1364/JOSAB.36.000E78
@@ -61,27 +70,27 @@ str_reftrans = 'r'
for i in range(nb_lam):
print(i)
- for j in range(nb_phi-1):
- 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))
- effr2_pin[:,i,j] = load_getdp_integral(respath+'eff_r2_lam%g_phi%g_psi0.out'%(i,j))
- Qscat_pin[i,j] = np.real(load_getdp_integral(respath+'Q_scat_lam%g_phi%g_psi0.out'%(i,j)))
- efft1_sin[:,i,j] = load_getdp_integral(respath+'eff_t1_lam%g_phi%g_psi1.out'%(i,j))
- efft2_sin[:,i,j] = load_getdp_integral(respath+'eff_t2_lam%g_phi%g_psi1.out'%(i,j))
- effr1_sin[:,i,j] = load_getdp_integral(respath+'eff_r1_lam%g_phi%g_psi1.out'%(i,j))
- effr2_sin[:,i,j] = load_getdp_integral(respath+'eff_r2_lam%g_phi%g_psi1.out'%(i,j))
- Qscat_sin[i,j] = np.real(load_getdp_integral(respath+'Q_scat_lam%g_phi%g_psi1.out'%(i,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]
+ 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]
T1_sin = np.sum(efft1_sin,axis=0)
T2_sin = np.sum(efft2_sin,axis=0)
@@ -119,7 +128,7 @@ M[4-1,4-1,:,:] = np.real(rpp*np.conj(rss)-rps*np.conj(rsp))
### Energy balance
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 pol in range(2):
if form==0 and pol==0:
@@ -144,7 +153,7 @@ for form in range(2):
plt.savefig('BALANCE_FLAG_TOT%g.jpg'%FLAG_TOT)
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[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')
@@ -158,14 +167,15 @@ 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$"
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):
for j in range(4):
ax=axes[i,j]
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_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,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)
ax.text(0. , 1. , r"$M_{%d%d}$"%(i+1,j+1) , fontsize=16 , transform=ax.transAxes)
ax.set_xticks([])
ax.set_yticks([])
@@ -173,14 +183,14 @@ for i in range(4):
# cbar.ax.locator_params(nbins=5)
# cbar.ax.tick_params(labelsize=14)
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.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='y', colors='C3') #setting up Y-axis tick color to black
ax.set_rlabel_position(70)
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,
rotation=70,ha='center',va='center',fontsize=16,color='C3')
norm = mpl.colors.Normalize(vmin=-1,vmax=1)
@@ -188,7 +198,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.jpg'%str_reftrans)
+plt.savefig('Mmatrix_%s_%s.jpg'%(flag_study,str_reftrans))
# plt.savefig('Mmatrix.pdf',bbox_inches='tight',pad_inches=0)
plt.show()
--
GitLab