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

M matrix plots

parent d4765a95
No related branches found
No related tags found
No related merge requests found
Pipeline #10105 passed
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib as mpl import matplotlib as mpl
plt.rcParams.update({"text.usetex": True})
def add_colorbar(mappable): def add_colorbar(mappable):
from mpl_toolkits.axes_grid1 import make_axes_locatable from mpl_toolkits.axes_grid1 import make_axes_locatable
...@@ -22,9 +23,10 @@ def load_getdp_integral(fname): ...@@ -22,9 +23,10 @@ def load_getdp_integral(fname):
return temp[:,1]+1j*temp[:,2] return temp[:,1]+1j*temp[:,2]
nb_lam = 100 nb_lam = 100
nb_phi = 50 nb_phi = 51
FLAG_TOT = 0 FLAG_TOT = 0
respath = 'res_Matrix_nb_phi50_total%g/'%(FLAG_TOT) # respath = 'res_Matrix_nb_phi50_total%g/'%(FLAG_TOT)
respath = 'res_Matrix_nb_phi50_total0_bermu1/'
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)
...@@ -55,7 +57,8 @@ M = np.zeros((4,4,len(tab_lam),len(tab_phi)),dtype=complex) ...@@ -55,7 +57,8 @@ M = np.zeros((4,4,len(tab_lam),len(tab_phi)),dtype=complex)
# 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)
for i in range(nb_lam): for i in range(nb_lam):
for j in range(nb_phi): 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)) 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)) 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)) rps[i,j] = load_getdp_integral(respath+'r_sin_pout_lam%g_phi%g.out'%(i,j))
...@@ -71,6 +74,12 @@ for i in range(nb_lam): ...@@ -71,6 +74,12 @@ for i in range(nb_lam):
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_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))) 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]
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)
R1_sin = np.sum(effr1_sin,axis=0) R1_sin = np.sum(effr1_sin,axis=0)
...@@ -101,8 +110,6 @@ M[4-1,2-1,:,:] =-np.imag(rpp*np.conj(rsp)-rps*np.conj(rss)) ...@@ -101,8 +110,6 @@ 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(2, 2, figsize=(12,12))
# L,P = np.meshgrid(tab_lam,tab_phi,indexing='ij') # L,P = np.meshgrid(tab_lam,tab_phi,indexing='ij')
...@@ -126,7 +133,7 @@ M[4-1,4-1,:,:] = np.real(rpp*np.conj(rss)-rps*np.conj(rsp)) ...@@ -126,7 +133,7 @@ M[4-1,4-1,:,:] = np.real(rpp*np.conj(rss)-rps*np.conj(rsp))
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$ (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):
...@@ -137,18 +144,20 @@ for i in range(4): ...@@ -137,18 +144,20 @@ for i in range(4):
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([])
plt.colorbar(sm, ax=ax) cbar = plt.colorbar(sm, ax=ax, fraction=0.046, pad=0.04)
cbar.ax.locator_params(nbins=5)
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_phi*np.pi/180,which_r,M[i,j]/M[0,0]-1,cmap=plt.cm.bwr,vmin=-1,vmax=1)
label_position=ax.get_rlabel_position() label_position=ax.get_rlabel_position()
ax.text(np.radians(label_position+30),ax.get_rmax()/1.3,rlabel, ax.text(np.radians(label_position+30),ax.get_rmax()/1.3,rlabel,
rotation=label_position,ha='center',va='center',fontsize=16) rotation=label_position,ha='center',va='center',fontsize=16)
ax.text(np.radians(-15),tab_lam.max()*1.01,r"$\varphi$",fontsize=16) ax.text(np.radians(-15),tab_lam.max()*1.01,r"$\varphi_0$",fontsize=16)
norm = mpl.colors.Normalize(vmin=-1,vmax=1) norm = mpl.colors.Normalize(vmin=-1,vmax=1)
sm = mpl.cm.ScalarMappable(cmap=plt.cm.bwr, norm=norm) sm = mpl.cm.ScalarMappable(cmap=plt.cm.bwr, norm=norm)
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.jpg')
plt.savefig('Mmatrix.pdf', bbox_inches='tight', pad_inches=0)
plt.show() plt.show()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment