diff --git a/DiffractionGratings/grating2D_clean_all.sh b/DiffractionGratings/grating2D_clean_all.sh
index d0513df214000bd16a631cf0ac235efb5caaef19..9d322ec29f6677c2a06e191a3e0adfd7855cf563 100644
--- a/DiffractionGratings/grating2D_clean_all.sh
+++ b/DiffractionGratings/grating2D_clean_all.sh
@@ -1,4 +1,5 @@
 rm -rf ./run_results
+rm -rf ./run_results3D
 rm *.msh
 rm *.pre
 rm *.db
diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index 3d7013cd897e504df2c956ec6bb51f3573d6c0f6..44d7b27b96db6ff66661e049011fbf190a57ccf9 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -1,3 +1,4 @@
+myDir = "run_results3D/";
 Include "grating3D_data_torus.geo";
 // Include "grating3D_data_hole.geo";
 // Include "grating3D_data_pyramid.geo";
@@ -55,22 +56,50 @@ Function{
 	om0         = 2*Pi*cel/lambda0;
 	k0          = 2.*Pi/lambda0;
 	Ae          = 1;
-	Ah          = Ae*Sqrt[ep0/mu0];
 	Pinc        =  0.5*Ae^2*Sqrt[eps_re_L_1*ep0/mu0] * Cos[theta0];
 
+  // permittivities
+  For i In {1:6}
+    If (flag_mat~{i}==0)
+      epsr[L~{i}]  = Complex[eps_re_L~{i} , eps_im_L~{i}] * TensorDiag[1,1,1];
+      If (i==1)
+        epsr1[] = Complex[eps_re_L~{i} , eps_im_L~{i}];
+      EndIf
+      If (i==6)
+        epsr2[] = Complex[eps_re_L~{i} , eps_im_L~{i}];
+      EndIf
+    Else
+      For j In {2:nb_available_materials}
+        If(flag_mat~{i}==j-1)
+          epsr[L~{i}]  = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]] * TensorDiag[1,1,1];
+          If (i==1)
+            epsr1[] = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]];
+          EndIf
+          If (i==6)
+            epsr2[] = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]];
+          EndIf
+        EndIf
+      EndFor
+    EndIf
+  EndFor
+
+  If (flag_mat_scat==0)
+    epsr[Scat] = Complex[eps_re_Scat , eps_im_Scat] * TensorDiag[1,1,1];
+  Else
+    For j In {flag_mat_scat:flag_mat_scat}
+      epsr[Scat] = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]] * TensorDiag[1,1,1];
+    EndFor
+  EndIf
+
 	For i In {1:6}
-		epsr[L~{i}]  = Complex[eps_re_L~{i} , eps_im_L~{i}] * TensorDiag[1,1,1];
 		mur[L~{i}]   = TensorDiag[1,1,1];
 	EndFor
-  epsr[Scat] = Complex[eps_re_Scat , eps_im_Scat] * TensorDiag[1,1,1];
 	mur[Scat]  = TensorDiag[1,1,1];
 
   ////// PMLS
 	a_pml           = 1.;
 	b_pml           = 1.;
   // bermu
-  epsr1[] = Complex[eps_re_L_1 , eps_im_L_1];
- 	epsr2[] = Complex[eps_re_L_6 , eps_im_L_6];
   n1[]    = Sqrt[epsr1[]];
   n2[]    = Sqrt[epsr2[]];
   k1norm[]= k0*n1[];
@@ -87,10 +116,10 @@ Function{
 	sz_bermubot[] = Complex[1,Damp_pml_bot[]/k2norm[]];
 	sy            = 1.;
   
-	epsr[PMLtop]  = eps_re_L_1*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
-	mur[PMLtop]   =            TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
-	epsr[PMLbot]  = eps_re_L_6*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
-	mur[PMLbot]   =            TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
+	epsr[PMLtop]  = Re[epsr1[]]*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
+	mur[PMLtop]   =             TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
+	epsr[PMLbot]  = Re[epsr2[]]*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
+	mur[PMLbot]   =             TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
 	// epsr[PMLtop]     = TensorDiag[sz_bermutop[]*sy/sx,sx*sz_bermutop[]/sy,sx*sy/sz_bermutop[]];
 	// mur[PMLtop]      = TensorDiag[sz_bermutop[]*sy/sx,sx*sz_bermutop[]/sy,sx*sy/sz_bermutop[]];
 	// epsr[PMLbot]     = TensorDiag[sz_bermubot[]*sy/sx,sx*sz_bermubot[]/sy,sx*sy/sz_bermubot[]];
@@ -98,7 +127,7 @@ Function{
 	
 	epsr_annex[PMLbot]       = epsr[];
 	epsr_annex[PMLtop]       = epsr[];
-	epsr_annex[Omega_source] = Complex[eps_re_L_1 , eps_im_L_1] * TensorDiag[1,1,1];
+	epsr_annex[Omega_source] = epsr1[] * TensorDiag[1,1,1];
 	epsr_annex[L_1]          = epsr[];
 	epsr_annex[L_6]          = epsr[];
 
@@ -126,12 +155,12 @@ Function{
   AmplEr[] = rp[]*Cos[psi0]*k1rp[] + rs[]*Sin[psi0]*k1rs[];
   AmplEt[] = tp[]*Cos[psi0]* k2p[] + ts[]*Sin[psi0]* k2s[];
   
-  Ei[]     = AmplEi[] * Exp[I[]*k1[]*XYZ[]];
-  Et[]     = AmplEt[] * Exp[I[]*k2[]*XYZ[]];
-  Er[]     = AmplEr[] * Exp[I[]*k1r[]*XYZ[]];
-  Hi[]     = 1/(om0*mu0*mur[]) * k1[]  /\ Ei[];
-  Hr[]     = 1/(om0*mu0*mur[]) * k1r[] /\ Er[];
-  Ht[]     = 1/(om0*mu0*mur[]) * k2[]  /\ Et[];
+  Ei[]     = Ae*AmplEi[] * Exp[I[]*k1[]*XYZ[]];
+  Et[]     = Ae*AmplEt[] * Exp[I[]*k2[]*XYZ[]];
+  Er[]     = Ae*AmplEr[] * Exp[I[]*k1r[]*XYZ[]];
+  Hi[]     =  1/(om0*mu0*mur[]) * k1[]  /\ Ei[];
+  Hr[]     =  1/(om0*mu0*mur[]) * k1r[] /\ Er[];
+  Ht[]     =  1/(om0*mu0*mur[]) * k2[]  /\ Et[];
 
   E1[Omega_super]  = Ei[]+Er[];
   E1[Omega_subs]   = Et[];
@@ -206,7 +235,7 @@ Jacobian {
 }
 
 Integration {
-  { Name Int_1 ;
+  { Name I1 ;
     Case { 
       { Type Gauss ;
         Case { 
@@ -257,9 +286,9 @@ Formulation {
       { Name u; Type Local; NameOfSpace Hcurl;}
     }
 		Equation{
-      Galerkin { [-1/mur[]    * Dof{Curl u} , {Curl u}]; In Omega       ; Jacobian JVol; Integration Int_1; }
-      Galerkin { [k0^2*epsr[] * Dof{u}      ,      {u}]; In Omega       ; Jacobian JVol; Integration Int_1; }
-      Galerkin { [source[]                  ,      {u}]; In Omega_source; Jacobian JVol; Integration Int_1; }
+      Galerkin { [-1/mur[]    * Dof{Curl u} , {Curl u}]; In Omega       ; Jacobian JVol; Integration I1; }
+      Galerkin { [k0^2*epsr[] * Dof{u}      ,      {u}]; In Omega       ; Jacobian JVol; Integration I1; }
+      Galerkin { [source[]                  ,      {u}]; In Omega_source; Jacobian JVol; Integration I1; }
     }
   }
 }
@@ -270,7 +299,9 @@ Resolution {
       { Name M; NameOfFormulation helmholtz_vector; Type ComplexValue; }
     }
     Operation { 
-      Generate[M]; Solve[M]; //SaveSolution[M];  
+      CreateDir[Str[myDir]];
+      Generate[M]; 
+      Solve[M]; //SaveSolution[M];  
     }
   }
 }
@@ -289,16 +320,16 @@ PostProcessing {
       { Name Poy_tot; Value { Local { [  0.5*Re[Cross[{u}+E1[] , Conj[ H1[]-I[]/(mur[]*mu0*om0)*{Curl u}]]]  ]; In Omega; Jacobian JVol; } } }
       
       For k In {2:6}
-        { Name Abs_L~{k} ; Value { Integral { [ ep0*om0 * 0.5*Fabs[eps_im_L~{k}]*(SquNorm[{u}+E1[]]) / (Pinc*period_x*period_y) ] ; In L~{k} ; Integration Int_1 ; Jacobian JVol ; } } }
+        { Name Abs_L~{k} ; Value { Integral { [ ep0*om0 * 0.5*Fabs[eps_im_L~{k}]*(SquNorm[{u}+E1[]]) / (Pinc*period_x*period_y) ] ; In L~{k} ; Integration I1 ; Jacobian JVol ; } } }
       EndFor
-      { Name Abs_scat ; Value { Integral { [ ep0*om0 * 0.5*Fabs[eps_im_Scat]*(SquNorm[{u}+E1[]]) / (Pinc*period_x*period_y) ] ; In Scat ; Integration Int_1 ; Jacobian JVol ; } } }
+      { Name Abs_scat ; Value { Integral { [ ep0*om0 * 0.5*Fabs[eps_im_Scat]*(SquNorm[{u}+E1[]]) / (Pinc*period_x*period_y) ] ; In Scat ; Integration I1 ; Jacobian JVol ; } } }
 
       For i In {0:Nb_ordre-1}
         For j In {0:Nb_ordre-1}
-          { Name int_x_t~{i}~{j} ; Value { Integral { [ CompX[{u}+E1[] ]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration Int_1 ; Jacobian JSur ; } } }
-          { Name int_y_t~{i}~{j} ; Value { Integral { [ CompY[{u}+E1[] ]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration Int_1 ; Jacobian JSur ; } } }
-          { Name int_x_r~{i}~{j} ; Value { Integral { [ CompX[{u}+E1d[]]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration Int_1 ; Jacobian JSur ; } } }
-          { Name int_y_r~{i}~{j} ; Value { Integral { [ CompY[{u}+E1d[]]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration Int_1 ; Jacobian JSur ; } } }
+          { Name int_x_t~{i}~{j} ; Value { Integral { [ CompX[{u}+E1[] ]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
+          { Name int_y_t~{i}~{j} ; Value { Integral { [ CompY[{u}+E1[] ]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
+          { Name int_x_r~{i}~{j} ; Value { Integral { [ CompX[{u}+E1d[]]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
+          { Name int_y_r~{i}~{j} ; Value { Integral { [ CompY[{u}+E1d[]]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
           { Name eff_t~{i}~{j}   ; Value { Term{Type Global; [
                  1/(gammat~{i}~{j}[]*-CompZ[k1[]]) *((gammat~{i}~{j}[]^2+alpha~{i}[]^2)*SquNorm[$int_x_t~{i}~{j}]+
                                                        (gammat~{i}~{j}[]^2+ beta~{j}[]^2)*SquNorm[$int_y_t~{i}~{j}]+
@@ -307,8 +338,7 @@ PostProcessing {
                  1/(gammar~{i}~{j}[]*-CompZ[k1[]])*((gammar~{i}~{j}[]^2+alpha~{i}[]^2)*SquNorm[$int_x_r~{i}~{j}]+
                                                       (gammar~{i}~{j}[]^2+ beta~{j}[]^2)*SquNorm[$int_y_r~{i}~{j}]+
                                                       2*alpha~{i}[]*beta~{j}[]*Re[$int_x_r~{i}~{j}*Conj[$int_y_r~{i}~{j}]]) ] ; In SurfIntTop ; } } }
-          { Name gammat~{i}~{j}   ; Value { Term{Type Global; [gammat~{i}~{j}[]] ; In SurfIntBot ; } } }
-          { Name gammar~{i}~{j}   ; Value { Term{Type Global; [gammar~{i}~{j}[]] ; In SurfIntTop ; } } }          
+          { Name numbering_ij~{i}~{j}   ; Value { Term{Type Global; [Vector[i-Nb_ordre,j-Nb_ordre,0]] ; In SurfIntBot ; } } }
         EndFor
       EndFor
 
@@ -329,51 +359,51 @@ PostOperation {
   { Name postop_helmholtz_vector; NameOfPostProcessing postpro_helmholtz_vector ;
     Operation {
       If (FlagOutEscaCuts==1)
-        Print [ u , OnPlane { {-period_x/2,0,hh_L_6-PML_bot} {period_x/2,0,hh_L_6-PML_bot} {-period_x/2,0,hh_L_1+thick_L_1+PML_top} } {npts_interpX,npts_interpZSca} , File "u_cut_Y=0.pos", Name "u_cut_Y=0"];
-        Print [ u , OnPlane { {0,-period_y/2,hh_L_6-PML_bot} {0,period_y/2,hh_L_6-PML_bot} {0,-period_y/2,hh_L_1+thick_L_1+PML_top} } {npts_interpY,npts_interpZSca} , File "u_cut_X=0.pos", Name "u_cut_X=0"];
+        Print [ u , OnPlane { {-period_x/2,0,hh_L_6-PML_bot} {period_x/2,0,hh_L_6-PML_bot} {-period_x/2,0,hh_L_1+thick_L_1+PML_top} } {npts_interpX,npts_interpZSca} , File StrCat[myDir,"u_cut_Y=0.pos"], Name "u_cut_Y=0"];
+        Print [ u , OnPlane { {0,-period_y/2,hh_L_6-PML_bot} {0,period_y/2,hh_L_6-PML_bot} {0,-period_y/2,hh_L_1+thick_L_1+PML_top} } {npts_interpY,npts_interpZSca} , File StrCat[myDir,"u_cut_X=0.pos"], Name "u_cut_X=0"];
       EndIf
       If (FlagOutEtotCuts==1)
-        Print [ Etot , OnPlane { {-period_x/2,0,hh_L_6} {period_x/2,0,hh_L_6} {-period_x/2,0,hh_L_1+thick_L_1} } {npts_interpX,npts_interpZTot} , File "Etot_cut_Y=0.pos", Name "Etot_cut_Y=0"];
-        Print [ Etot , OnPlane { {0,-period_y/2,hh_L_6} {0,period_y/2,hh_L_6} {0,-period_y/2,hh_L_1+thick_L_1} } {npts_interpY,npts_interpZTot} , File "Etot_cut_X=0.pos", Name "Etot_cut_X=0"];
+        Print [ Etot , OnPlane { {-period_x/2,0,hh_L_6} {period_x/2,0,hh_L_6} {-period_x/2,0,hh_L_1+thick_L_1} } {npts_interpX,npts_interpZTot} , File StrCat[myDir,"Etot_cut_Y=0.pos"], Name "Etot_cut_Y=0"];
+        Print [ Etot , OnPlane { {0,-period_y/2,hh_L_6} {0,period_y/2,hh_L_6} {0,-period_y/2,hh_L_1+thick_L_1} } {npts_interpY,npts_interpZTot} , File StrCat[myDir,"Etot_cut_X=0.pos"], Name "Etot_cut_X=0"];
       EndIf
       If (FlagOutPoyCut==1)
-        Print [ Poy_tot , OnPlane { {-period_x/2,0,hh_L_6} {period_x/2,0,hh_L_6} {-period_x/2,0,hh_L_1+thick_L_1} } {npts_interpX,npts_interpZSca} , File "Poy_tot_cut_Y=0.pos", Name "Poy_tot_cut_Y=0"];
-        Print [ Poy_tot , OnPlane { {0,-period_y/2,hh_L_6} {0,period_y/2,hh_L_6} {0,-period_y/2,hh_L_1+thick_L_1} } {npts_interpY,npts_interpZSca} , File "Poy_tot_cut_X=0.pos", Name "Poy_tot_cut_X=0"];
+        Print [ Poy_tot , OnPlane { {-period_x/2,0,hh_L_6} {period_x/2,0,hh_L_6} {-period_x/2,0,hh_L_1+thick_L_1} } {npts_interpX,npts_interpZSca} , File StrCat[myDir,"Poy_tot_cut_Y=0.pos"], Name "Poy_tot_cut_Y=0"];
+        Print [ Poy_tot , OnPlane { {0,-period_y/2,hh_L_6} {0,period_y/2,hh_L_6} {0,-period_y/2,hh_L_1+thick_L_1} } {npts_interpY,npts_interpZSca} , File StrCat[myDir,"Poy_tot_cut_X=0.pos"], Name "Poy_tot_cut_X=0"];
       EndIf
       If (FlagOutEscaFull==1)
         If (Flag_interp_cubic==1)
-          Print [ u , OnBox { {-period_x/2,-period_y/2,hh_L_6-PML_bot} {period_x/2,-period_y/2,hh_L_6-PML_bot} {-period_x/2,period_y/2,hh_L_6-PML_bot} {-period_x/2,-period_y/2,hh_L_1+thick_L_1+PML_top} } {npts_interpX,npts_interpY,npts_interpZSca} , File "u_grid.pos", Name "u_grid"];
+          Print [ u , OnBox { {-period_x/2,-period_y/2,hh_L_6-PML_bot} {period_x/2,-period_y/2,hh_L_6-PML_bot} {-period_x/2,period_y/2,hh_L_6-PML_bot} {-period_x/2,-period_y/2,hh_L_1+thick_L_1+PML_top} } {npts_interpX,npts_interpY,npts_interpZSca} , File StrCat[myDir,"u_grid.pos"], Name "u_grid"];
         Else
-          Print [ u , OnElementsOf Omega, File "Etot.pos"];
+          Print [ u , OnElementsOf Omega, File StrCat[myDir,"Etot.pos"]];
         EndIf
       EndIf
       If (FlagOutEtotFull==1)
         If (Flag_interp_cubic==1)
-          Print [ Etot , OnBox { {-period_x/2,-period_y/2,hh_L_6} {period_x/2,-period_y/2,hh_L_6} {-period_x/2,period_y/2,hh_L_6} {-period_x/2,-period_y/2,hh_L_1+thick_L_1} } {npts_interpX,npts_interpY,npts_interpZTot} , File "Etot_grid.pos", Name "Etot_grid"];
+          Print [ Etot , OnBox { {-period_x/2,-period_y/2,hh_L_6} {period_x/2,-period_y/2,hh_L_6} {-period_x/2,period_y/2,hh_L_6} {-period_x/2,-period_y/2,hh_L_1+thick_L_1} } {npts_interpX,npts_interpY,npts_interpZTot} , File StrCat[myDir,"Etot_grid.pos"], Name "Etot_grid"];
         Else
-          Print [ Etot , OnElementsOf Omega_plot, File "Etot.pos"];
+          Print [ Etot , OnElementsOf Omega_plot, File StrCat[myDir,"Etot.pos"]];
         EndIf
       EndIf
       If (FlagOutPoyFull==1)
         If (Flag_interp_cubic==1)
-          Print [ Poy_tot , OnBox { {-period_x/2,-period_y/2,hh_L_6} {period_x/2,-period_y/2,hh_L_6} {-period_x/2,period_y/2,hh_L_6} {-period_x/2,-period_y/2,hh_L_1+thick_L_1} } {npts_interpX,npts_interpY,npts_interpZTot} , File "Poy_tot_grid.pos", Name "Poy_tot_grid"];
+          Print [ Poy_tot , OnBox { {-period_x/2,-period_y/2,hh_L_6} {period_x/2,-period_y/2,hh_L_6} {-period_x/2,period_y/2,hh_L_6} {-period_x/2,-period_y/2,hh_L_1+thick_L_1} } {npts_interpX,npts_interpY,npts_interpZTot} , File StrCat[myDir,"Poy_tot_grid.pos", Name "Poy_tot_grid"];
         Else
-          Print [ Poy_tot , OnElementsOf Omega_plot, File "Poy_tot.pos"];
+          Print [ Poy_tot , OnElementsOf Omega_plot, File StrCat[myDir,"Poy_tot.pos"]];
         EndIf
       EndIf
       // // For DEBUG
-      // Print [ epsr_xx , OnElementsOf Omega, File "epsr_xx.pos"];
-      // Print [ Ecm , OnElementsOf Omega_plot, File "Ecm.pos"];
-      // Print [ E1, OnElementsOf Omega_plot, File "E1.pos"];
-      // Print [ Damp_pml_top, OnElementsOf PMLtop, File "Damp_pml_top.pos"];
+      // Print [ epsr_xx , OnElementsOf Omega, File StrCat[myDir,"epsr_xx.pos"];
+      // Print [ Ecm , OnElementsOf Omega_plot, File StrCat[myDir,"Ecm.pos"];
+      // Print [ E1, OnElementsOf Omega_plot, File StrCat[myDir,"E1.pos"];
+      // Print [ Damp_pml_top, OnElementsOf PMLtop, File StrCat[myDir,"Damp_pml_top.pos"];
       // For i In {0:nb_slice-1}
-      //   Print [ Etot , OnPlane { {-period_x/2+small_delta,-period_y/2+small_delta,zcut_sub~{i}} {-period_x/2+small_delta,period_y/2-small_delta,zcut_sub~{i}} {period_x/2-small_delta,-period_y/2+small_delta,zcut_sub~{i}} } {ninterv_integ,ninterv_integ} , File > "./Views/Etot_XYcut.out",Format Table];
-      //   Print [ Edif , OnPlane { {-period_x/2+small_delta,-period_y/2+small_delta,zcut_sup~{i}} {-period_x/2+small_delta,period_y/2-small_delta,zcut_sup~{i}} {period_x/2-small_delta,-period_y/2+small_delta,zcut_sup~{i}} } {ninterv_integ,ninterv_integ} , File > "./Views/Edif_XYcut.out",Format Table];
+      //   Print [ Etot , OnPlane { {-period_x/2+small_delta,-period_y/2+small_delta,zcut_sub~{i}} {-period_x/2+small_delta,period_y/2-small_delta,zcut_sub~{i}} {period_x/2-small_delta,-period_y/2+small_delta,zcut_sub~{i}} } {ninterv_integ,ninterv_integ} , File StrCat[myDir,> "./Views/Etot_XYcut.out",Format Table];
+      //   Print [ Edif , OnPlane { {-period_x/2+small_delta,-period_y/2+small_delta,zcut_sup~{i}} {-period_x/2+small_delta,period_y/2-small_delta,zcut_sup~{i}} {period_x/2-small_delta,-period_y/2+small_delta,zcut_sup~{i}} } {ninterv_integ,ninterv_integ} , File StrCat[myDir,> "./Views/Edif_XYcut.out",Format Table];
       // EndFor
       For k In {2:6}
-        Print[ Abs_L~{k}[L~{k}], OnGlobal, File Sprintf["temp-Q_L_%g.txt",k], Format Table ];
+        Print[ Abs_L~{k}[L~{k}], OnGlobal, File > StrCat[myDir,Sprintf["temp-Q_L_%g.txt",k]], Format Table ];
       EndFor
-      Print[ Abs_scat[Scat]  , OnGlobal, File "temp-Q_scat.txt", Format Table ];
+      Print[ Abs_scat[Scat]  , OnGlobal, File > StrCat[myDir,"temp-Q_scat.txt"], Format Table ];
 
       For i In {0:Nb_ordre-1}
         For j In {0:Nb_ordre-1}
@@ -386,14 +416,13 @@ PostOperation {
 
       For i In {0:Nb_ordre-1}
         For j In {0:Nb_ordre-1}
-          Print[ eff_t~{i}~{j}[SurfIntBot], OnRegion SurfIntBot, File > "eff_t.txt", Format Table ];
-          Print[ eff_r~{i}~{j}[SurfIntTop], OnRegion SurfIntTop, File > "eff_r.txt", Format Table ];
+          Print[ eff_t~{i}~{j}[SurfIntBot], OnRegion SurfIntBot, File > StrCat[myDir, "eff_t.txt"], Format Table ];
+          Print[ eff_r~{i}~{j}[SurfIntTop], OnRegion SurfIntTop, File > StrCat[myDir, "eff_r.txt"], Format Table ];
         EndFor
       EndFor
       For i In {0:Nb_ordre-1}
         For j In {0:Nb_ordre-1}
-          Print[ gammat~{i}~{j}[SurfIntBot], OnRegion SurfIntBot, File > "gammat.txt", Format Table ];
-          Print[ gammar~{i}~{j}[SurfIntBot], OnRegion SurfIntTop, File > "gammar.txt", Format Table ];
+          Print[ numbering_ij~{i}~{j}[SurfIntBot], OnRegion SurfIntBot, File StrCat[myDir,"numbering_ij.txt"], Format Table ];
         EndFor
       EndFor
     }
diff --git a/DiffractionGratings/gratings_common_materials.pro b/DiffractionGratings/gratings_common_materials.pro
index 032217703c8f581cb59a17576100237443153de9..352842b94dcb5da19100a4ffef8c0015b06c8ebc 100644
--- a/DiffractionGratings/gratings_common_materials.pro
+++ b/DiffractionGratings/gratings_common_materials.pro
@@ -68,3 +68,30 @@ epsr_mat_im_12 = {-9.93483e-01,-8.48454e-01,-7.26134e-01,-6.23625e-01,-5.38097e-
 lambdamat_13   = {1.90700e-07,2.06600e-07,2.13800e-07,2.11400e-07,2.29600e-07,2.38400e-07,2.48000e-07,2.58300e-07,2.69500e-07,2.81800e-07,2.95200e-07,3.09900e-07,3.26300e-07,3.44400e-07,3.64600e-07,3.87400e-07,4.13300e-07,4.42800e-07,4.76800e-07,5.16600e-07,5.39000e-07,5.63500e-07,5.90400e-07,6.19900e-07,6.52500e-07,6.70200e-07,6.88000e-07,7.08400e-07,7.29300e-07,8.26500e-07,1.24000e-06,1.26500e-06,1.29100e-06,1.31900e-06,1.34800e-06,1.37800e-06,1.40900e-06,1.44200e-06,1.47600e-06,1.51200e-06,1.55000e-06,1.58900e-06,1.63100e-06,1.67500e-06,1.72200e-06,1.77100e-06,1.82300e-06,1.87800e-06,1.93700e-06,2.00000e-06,2.06600e-06,2.13800e-06,2.21400e-06,2.29600e-06,2.38400e-06,2.48000e-06,2.58300e-06,2.69500e-06,2.81800e-06,2.95200e-06,3.10000e-06,3.26300e-06,3.44400e-06,3.64600e-06};
 epsr_mat_re_13 = {-9.59136e-01,-1.44650e+00,-1.57890e+00,-1.63520e+00,-1.53000e+00,-1.33560e+00,-1.00750e+00,-5.83200e-01,-4.78500e-01,-4.69500e-01,-6.73200e-01,-1.16280e+00,-1.48050e+00,-1.78080e+00,-2.18960e+00,-3.49170e+00,-3.49170e+00,-4.20070e+00,-4.92750e+00,-5.50560e+00,-5.62650e+00,-6.07772e+00,-7.67708e+00,-1.04236e+01,-1.34231e+01,-1.48534e+01,-1.63571e+01,-1.79318e+01,-1.95752e+01,-2.76000e+01,-7.13841e+01,-4.57224e+01,-4.76314e+01,-4.95784e+01,-5.17085e+01,-5.38823e+01,-5.64017e+01,-5.89775e+01,-6.17672e+01,-6.46202e+01,-6.78604e+01,-7.13478e+01,-7.50949e+01,-7.91148e+01,-8.34198e+01,-8.80254e+01,-9.31369e+01,-9.87894e+01,-1.05426e+02,-1.11638e+02,-1.20208e+02,-1.26820e+02,-1.38280e+02,-1.47779e+02,-1.60102e+02,-1.72917e+02,-1.86202e+02,-2.05696e+02,-2.23123e+02,-2.44329e+02,-2.69722e+02,-2.99767e+02,-3.34950e+02,-3.79792e+02};
 epsr_mat_im_13 = {-2.62492e+00,-3.30720e+00,-3.67400e+00,-4.10640e+00,-4.55680e+00,-4.96800e+00,-5.23320e+00,-5.23260e+00,-5.07680e+00,-4.88720e+00,-4.65760e+00,-4.60960e+00,-4.85080e+00,-4.89940e+00,-4.95300e+00,-5.21560e+00,-5.21560e+00,-5.52240e+00,-5.75000e+00,-5.82400e+00,-5.38720e+00,-4.29520e+00,-2.63016e+00,-1.76256e+00,-1.57076e+00,-1.65980e+00,-1.72530e+00,-1.81472e+00,-1.97578e+00,-2.73520e+00,-7.32636e+00,-6.72576e+00,-6.98920e+00,-7.27180e+00,-7.57050e+00,-7.88992e+00,-8.23782e+00,-8.60860e+00,-9.01472e+00,-9.44632e+00,-1.00111e+01,-1.06214e+01,-1.12796e+01,-1.19885e+01,-1.27690e+01,-1.36069e+01,-1.45587e+01,-1.55931e+01,-1.67890e+01,-1.80200e+01,-1.95800e+01,-2.10858e+01,-2.31280e+01,-2.51320e+01,-2.76860e+01,-3.03600e+01,-3.34280e+01,-3.71520e+01,-4.11000e+01,-4.61580e+01,-5.24700e+01,-6.02040e+01,-6.99200e+01,-8.19280e+01};
+
+Function{
+  epsr_re_interp_mat_2[]  = InterpolationLinear[$1]{ListAlt[lambdamat_2 ,epsr_mat_re_2 ]};
+  epsr_im_interp_mat_2[]  = InterpolationLinear[$1]{ListAlt[lambdamat_2 ,epsr_mat_im_2 ]};
+  epsr_re_interp_mat_3[]  = InterpolationLinear[$1]{ListAlt[lambdamat_3 ,epsr_mat_re_3 ]};
+  epsr_im_interp_mat_3[]  = InterpolationLinear[$1]{ListAlt[lambdamat_3 ,epsr_mat_im_3 ]};
+  epsr_re_interp_mat_4[]  = InterpolationLinear[$1]{ListAlt[lambdamat_4 ,epsr_mat_re_4 ]};
+  epsr_im_interp_mat_4[]  = InterpolationLinear[$1]{ListAlt[lambdamat_4 ,epsr_mat_im_4 ]};
+  epsr_re_interp_mat_5[]  = InterpolationLinear[$1]{ListAlt[lambdamat_5 ,epsr_mat_re_5 ]};
+  epsr_im_interp_mat_5[]  = InterpolationLinear[$1]{ListAlt[lambdamat_5 ,epsr_mat_im_5 ]};
+  epsr_re_interp_mat_6[]  = InterpolationLinear[$1]{ListAlt[lambdamat_6 ,epsr_mat_re_6 ]};
+  epsr_im_interp_mat_6[]  = InterpolationLinear[$1]{ListAlt[lambdamat_6 ,epsr_mat_im_6 ]};
+  epsr_re_interp_mat_7[]  = InterpolationLinear[$1]{ListAlt[lambdamat_7 ,epsr_mat_re_7 ]};
+  epsr_im_interp_mat_7[]  = InterpolationLinear[$1]{ListAlt[lambdamat_7 ,epsr_mat_im_7 ]};
+  epsr_re_interp_mat_8[]  = InterpolationLinear[$1]{ListAlt[lambdamat_8 ,epsr_mat_re_8 ]};
+  epsr_im_interp_mat_8[]  = InterpolationLinear[$1]{ListAlt[lambdamat_8 ,epsr_mat_im_8 ]};
+  epsr_re_interp_mat_9[]  = InterpolationLinear[$1]{ListAlt[lambdamat_9 ,epsr_mat_re_9 ]};
+  epsr_im_interp_mat_9[]  = InterpolationLinear[$1]{ListAlt[lambdamat_9 ,epsr_mat_im_9 ]};
+  epsr_re_interp_mat_10[] = InterpolationLinear[$1]{ListAlt[lambdamat_10,epsr_mat_re_10]};
+  epsr_im_interp_mat_10[] = InterpolationLinear[$1]{ListAlt[lambdamat_10,epsr_mat_im_10]};
+  epsr_re_interp_mat_11[] = InterpolationLinear[$1]{ListAlt[lambdamat_11,epsr_mat_re_11]};
+  epsr_im_interp_mat_11[] = InterpolationLinear[$1]{ListAlt[lambdamat_11,epsr_mat_im_11]};
+  epsr_re_interp_mat_12[] = InterpolationLinear[$1]{ListAlt[lambdamat_12,epsr_mat_re_12]};
+  epsr_im_interp_mat_12[] = InterpolationLinear[$1]{ListAlt[lambdamat_12,epsr_mat_im_12]};
+  epsr_re_interp_mat_13[] = InterpolationLinear[$1]{ListAlt[lambdamat_13,epsr_mat_re_13]};
+  epsr_im_interp_mat_13[] = InterpolationLinear[$1]{ListAlt[lambdamat_13,epsr_mat_im_13]};
+}
\ No newline at end of file