diff --git a/DiffractionGratings/grating3D_postplot_Mmatrix.py b/DiffractionGratings/grating3D_postplot_Mmatrix.py
new file mode 100644
index 0000000000000000000000000000000000000000..f209b638a953679cc6911fe862db91e61d1e190a
--- /dev/null
+++ b/DiffractionGratings/grating3D_postplot_Mmatrix.py
@@ -0,0 +1,155 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+
+def add_colorbar(mappable):
+ from mpl_toolkits.axes_grid1 import make_axes_locatable
+ import matplotlib.pyplot as plt
+ last_axes = plt.gca()
+ ax = mappable.axes
+ fig = ax.figure
+ divider = make_axes_locatable(ax)
+ cax = divider.append_axes("right", size="5%", pad=0.05)
+ cbar = fig.colorbar(mappable, cax=cax)
+ plt.sca(last_axes)
+ return cbar
+
+def load_getdp_integral(fname):
+ temp = np.loadtxt(fname)
+ if len(temp.shape)==1:
+ return temp[1]+1j*temp[2]
+ else:
+ return temp[:,1]+1j*temp[:,2]
+
+nb_lam = 100
+nb_phi = 50
+FLAG_TOT = 0
+respath = 'res_Matrix_nb_phi50_total%g/'%(FLAG_TOT)
+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(400,800,nb_lam)
+tab_phi = np.linspace(0,360,nb_phi)
+M = np.zeros((4,4,len(tab_lam),len(tab_phi)),dtype=complex)
+
+
+### Convention used e.g. in https://doi.org/10.1364/JOSAB.36.000E78
+# [Ep] [ rpp rps ] [Ep]
+# [ ] = [ ] [ ]
+# [Es]ref [ rsp rss ] [Es]inc
+# where Ep (resp. Es) are the plane wave electric field components
+# that are parallel (resp. perpendicular) to the plane of incidence
+# rps : send Es_inc (s-in), project the reflected field along p (p-out)
+
+for i in range(nb_lam):
+ for j in range(nb_phi):
+ 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))
+ rps[i,j] = load_getdp_integral(respath+'r_sin_pout_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))
+ 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)))
+
+T1_sin = np.sum(efft1_sin,axis=0)
+T2_sin = np.sum(efft2_sin,axis=0)
+R1_sin = np.sum(effr1_sin,axis=0)
+R2_sin = np.sum(effr2_sin,axis=0)
+T1_pin = np.sum(efft1_pin,axis=0)
+T2_pin = np.sum(efft2_pin,axis=0)
+R1_pin = np.sum(effr1_pin,axis=0)
+R2_pin = np.sum(effr2_pin,axis=0)
+TOT1_pin = T1_pin+R1_pin+Qscat_pin
+TOT2_pin = T2_pin+R2_pin+Qscat_pin
+TOT1_sin = T1_sin+R1_sin+Qscat_sin
+TOT2_sin = T2_sin+R2_sin+Qscat_sin
+
+M[1-1,1-1,:,:] = 0.5*(np.abs(rpp)**2+np.abs(rsp)**2+np.abs(rps)**2+np.abs(rss)**2)
+M[1-1,2-1,:,:] = 0.5*(np.abs(rpp)**2+np.abs(rsp)**2-np.abs(rps)**2-np.abs(rss)**2)
+M[1-1,3-1,:,:] = np.real(rpp*np.conj(rps)+rsp*np.conj(rss)) #*-1
+M[1-1,4-1,:,:] = np.imag(rpp*np.conj(rps)+rsp*np.conj(rss))
+M[2-1,1-1,:,:] = 0.5*(np.abs(rpp)**2-np.abs(rsp)**2+np.abs(rps)**2-np.abs(rss)**2)
+M[2-1,2-1,:,:] = 0.5*(np.abs(rpp)**2-np.abs(rsp)**2-np.abs(rps)**2+np.abs(rss)**2)
+M[2-1,3-1,:,:] = np.real(rpp*np.conj(rps)-rsp*np.conj(rss)) #*-1
+M[2-1,4-1,:,:] = np.imag(rpp*np.conj(rps)-rsp*np.conj(rss))
+M[3-1,1-1,:,:] = np.real(rpp*np.conj(rsp)+rps*np.conj(rss)) #*-1
+M[3-1,2-1,:,:] = np.real(rpp*np.conj(rsp)-rps*np.conj(rss)) #*-1
+M[3-1,3-1,:,:] = np.real(rpp*np.conj(rss)+rps*np.conj(rsp))
+M[3-1,4-1,:,:] = np.imag(rpp*np.conj(rss)-rps*np.conj(rsp)) #*-1
+M[4-1,1-1,:,:] =-np.imag(rpp*np.conj(rsp)+rps*np.conj(rss))
+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)) #*-1
+M[4-1,4-1,:,:] = np.real(rpp*np.conj(rss)-rps*np.conj(rsp))
+
+
+
+# ### Energy balance
+# fig, axes = plt.subplots(2, 2, figsize=(12,12))
+# L,P = np.meshgrid(tab_lam,tab_phi,indexing='ij')
+# for form in range(2):
+# for pol in range(2):
+# if form==0 and pol==0:
+# data=TOT1_pin
+# title = 'TOT1_pin'
+# if form==0 and pol==1:
+# data=TOT1_sin
+# title = 'TOT1_sin'
+# if form==1 and pol==0:
+# data=TOT2_pin
+# title = 'TOT2_pin'
+# if form==1 and pol==1:
+# data=TOT2_sin
+# title = 'TOT2_sin'
+# zplot = axes[form,pol].pcolormesh(L,P,data.real,vmin=0.99,vmax=1.01)
+# add_colorbar(zplot)
+# plt.savefig('BALANCE_FLAG_TOT%g.jpg'%FLAG_TOT)
+
+fig, axes = plt.subplots(4, 4, subplot_kw=dict(projection='polar') ,figsize=(12,12))
+flag_lam = True
+rlabel = r"$\lambda$ (nm)"
+which_r = 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)
+ ax.text(0. , 1. , r"$M_{%d%d}$"%(i+1,j+1) , fontsize=16 , transform=ax.transAxes)
+ ax.set_xticks([])
+ ax.set_yticks([])
+ plt.colorbar(sm, ax=ax)
+ 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)
+ label_position=ax.get_rlabel_position()
+ ax.text(np.radians(label_position+30),ax.get_rmax()/1.3,rlabel,
+ rotation=label_position,ha='center',va='center',fontsize=16)
+ ax.text(np.radians(-15),tab_lam.max()*1.01,r"$\varphi$",fontsize=16)
+ norm = mpl.colors.Normalize(vmin=-1,vmax=1)
+ sm = mpl.cm.ScalarMappable(cmap=plt.cm.bwr, norm=norm)
+ sm.set_array([])
+ # plt.colorbar(sm, ticks=np.linspace(-1,1,3))
+ 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.jpg')
+plt.show()
+