diff --git a/DiffractionGratings/grating2D_conical.pro b/DiffractionGratings/grating2D_conical.pro
index 4e7a9f8d79f886f966f4f32d83df44165ee64847..b7251d35833b527252147ba312b8fe6a4b57489a 100644
--- a/DiffractionGratings/grating2D_conical.pro
+++ b/DiffractionGratings/grating2D_conical.pro
@@ -439,6 +439,7 @@ Resolution {
 PostProcessing {
   { Name postpro_energy; NameOfFormulation helmoltz_conical;
     Quantity {
+    { Name lambda_step   ; Value { Local { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
     { Name E2d    ; Value { Local { [ {Et}+{El}      ]; In Omega; Jacobian JVol; } } }
     { Name Etot   ; Value { Local { [ {Et}+{El}+E1[] ]; In Omega; Jacobian JVol; } } }
     { Name H2d    ; Value { Local { [ (-I[]/(mur[]*mu0*om0))*({Curl Et}+{Curl El}) ]; In Omega; Jacobian JVol; } } }
@@ -497,6 +498,7 @@ PostProcessing {
 PostOperation {
   { Name postop_energy; NameOfPostProcessing postpro_energy ;
     Operation {
+      Print[ lambda_step, OnPoint{0,0,0}, Format Table, File > StrCat[myDir, "temp_lambda_step.txt"], SendToServer "GetDP/Lambda_step" ] ;
       Print [ Etot , OnElementsOf Omega, File "Etot.pos" ];
       For i In {0:2*nb_orders}
         Print[ int_x_t~{i}[SurfCutSubs1] , OnGlobal, StoreInVariable $int_x_t~{i}, File > StrCat[myDir, "int_x_t.txt"], Format Table];
@@ -526,10 +528,6 @@ DefineConstant[
   P_ = {"postop_energy", Name "GetDP/2PostOperationChoices", Visible 1}
 ];
 
-DefineConstant[
-  M_ = {"grating2D.msh", Name "Gmsh/MshFileName", Visible 1}
-];
-
 If(plotRTgraphs)
   DefineConstant[
     refl_  = {0, Name "GetDP/R0", ReadOnly 1, Graph "02000000", Visible 1},
diff --git a/DiffractionGratings/grating2D_data.geo b/DiffractionGratings/grating2D_data.geo
index f2f238e179bc827359ea076427cfd3cfbc171fac..8f4b1417f230f5a8696f2cbb31c02819a6e8a763 100644
--- a/DiffractionGratings/grating2D_data.geo
+++ b/DiffractionGratings/grating2D_data.geo
@@ -8,4 +8,6 @@ DefineConstant[
   }
 ];
 
+nb_plot_periods = 3;
+
 Include StrCat["grating2D_data_",test_case,".geo"];
diff --git a/DiffractionGratings/grating2D_scalar.pro b/DiffractionGratings/grating2D_scalar.pro
index cccf5dd541e18b7d7a043a7afd013fa15f087cb9..7e3e68f98a974bff99067ece2e1e326706369405 100644
--- a/DiffractionGratings/grating2D_scalar.pro
+++ b/DiffractionGratings/grating2D_scalar.pro
@@ -372,6 +372,9 @@ Resolution {
 PostProcessing {
   { Name postpro_energy; NameOfFormulation helmoltz_scalar;
     Quantity {
+      For i In {-nb_plot_periods:nb_plot_periods}
+        { Name u_tot~{i}     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*i*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
+      EndFor
       If (flag_polar==1)
         { Name debr         ; Value { Local { [ r[]                 ]; In Omega; Jacobian JVol; } } }
         { Name debt         ; Value { Local { [ t[]                 ]; In Omega; Jacobian JVol; } } }
@@ -382,14 +385,6 @@ PostProcessing {
         { Name NormHz_tot   ; Value { Local { [ Norm[{u2d}+u1[]]    ]; In Omega; Jacobian JVol; } } }
         { Name E1           ; Value { Local { [ E1[]                ]; In Omega; Jacobian JVol; } } }
         { Name u1           ; Value { Local { [ u1[]                ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totp1     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 1*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totp2     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 2*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totp3     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 3*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totp4     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 4*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totm1     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-1*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totm2     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-2*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totm3     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-3*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totm4     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-4*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
         { Name boundary     ; Value { Local { [ bndCol[] ] ; In Plot_bnd ; Jacobian JVol ; } } }
 
         For i In {0:2*nb_orders}
@@ -454,59 +449,60 @@ PostProcessing {
 PostOperation {
   { Name postop_energy; NameOfPostProcessing postpro_energy ;
     Operation {
+      Print[ lambda_step, OnPoint{0,0,0}, Format Table, File > StrCat[myDir, "temp_lambda_step.txt"], SendToServer "GetDP/Lambda_step" ] ;
+      For i In {0:2*nb_orders}
+        Print[ s_r~{i}[SurfCutSuper1], OnGlobal, Store i                , Format Table , File > StrCat[myDir, Sprintf("temp_s_r_%g.txt", i-nb_orders)]];
+        Print[ s_t~{i}[SurfCutSubs1] , OnGlobal, Store (2*nb_orders+1+i), Format Table , File > StrCat[myDir, Sprintf("temp_s_t_%g.txt", i-nb_orders)]];
+      EndFor
+      For i In {0:2*nb_orders}
+        Print[ eff_r~{i}[SurfCutSuper1], OnRegion SurfCutSuper1,Store (4*nb_orders+1+i), Format FrequencyTable, File > StrCat[myDir, Sprintf("efficiency_r_%g.txt", i-nb_orders)]];
+        Print[ eff_t~{i}[SurfCutSubs1] , OnRegion SurfCutSubs1 ,Store (6*nb_orders+1+i), Format FrequencyTable, File > StrCat[myDir, Sprintf("efficiency_t_%g.txt", i-nb_orders)]];
+        Print[ order_r_angle~{i}     , OnPoint{0,0,0}, Format Table , File > StrCat[myDir, Sprintf("order_r_angle_%g.txt", i-nb_orders)]];
+        Print[ order_t_angle~{i}     , OnPoint{0,0,0}, Format Table , File > StrCat[myDir, Sprintf("order_t_angle_%g.txt", i-nb_orders)]];          
+      EndFor
+      Print[ eff_r~{nb_orders}[SurfCutSuper1], OnRegion SurfCutSuper1, Format Table, SendToServer "GetDP/R0", File StrCat[myDir, "temp_R0.txt"]];
+      Print[ eff_t~{nb_orders}[SurfCutSubs1] , OnRegion SurfCutSubs1 , Format Table, SendToServer "GetDP/T0", File StrCat[myDir, "temp_T0.txt"]];
+      Print[ Q_tot[Plot_domain]             , OnGlobal  , Format FrequencyTable ,SendToServer "GetDP/total absorption", File > StrCat[myDir, "absorption-Q_tot.txt"]];
+      For i In {0:N_rods-1:1}
+        Print[ Q_rod~{i}[rod~{i}] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, Sprintf("absorption-Q_rod_%g.txt", i+1) ]];           
+      EndFor  
+      Print[ Q_tot[Plot_domain]     , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_tot.txt"]];
+      Print[ Q_subs[sub]            , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_subs.txt"]];
+      Print[ Q_rod_out[rod_out]     , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_rod_out.txt"]];
+      Print[ Q_layer_dep[layer_dep] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_dep.txt"]];
+      Print[ Q_layer_cov[layer_cov] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_cov.txt"]];
       If (flag_polar==1)
-        Print[ lambda_step, OnPoint{0,0,0}, Format Table, File > StrCat[myDir, "temp_lambda_step.txt"], SendToServer "GetDP/Lambda_step" ] ;
-        For i In {0:2*nb_orders}
-          Print[ s_r~{i}[SurfCutSuper1], OnGlobal, Store i                , Format Table , File > StrCat[myDir, Sprintf("temp_s_r_%g.txt", i-nb_orders)]];
-          Print[ s_t~{i}[SurfCutSubs1] , OnGlobal, Store (2*nb_orders+1+i), Format Table , File > StrCat[myDir, Sprintf("temp_s_t_%g.txt", i-nb_orders)]];
-        EndFor
-        For i In {0:2*nb_orders}
-          Print[ eff_r~{i}[SurfCutSuper1], OnRegion SurfCutSuper1,Store (4*nb_orders+1+i), Format FrequencyTable, File > StrCat[myDir, Sprintf("efficiency_r_%g.txt", i-nb_orders)]];
-          Print[ eff_t~{i}[SurfCutSubs1] , OnRegion SurfCutSubs1 ,Store (6*nb_orders+1+i), Format FrequencyTable, File > StrCat[myDir, Sprintf("efficiency_t_%g.txt", i-nb_orders)]];
-          Print[ order_r_angle~{i}     , OnPoint{0,0,0}, Format Table , File > StrCat[myDir, Sprintf("order_r_angle_%g.txt", i-nb_orders)]];
-          Print[ order_t_angle~{i}     , OnPoint{0,0,0}, Format Table , File > StrCat[myDir, Sprintf("order_t_angle_%g.txt", i-nb_orders)]];          
-        EndFor
-        Print[ eff_r~{nb_orders}[SurfCutSuper1], OnRegion SurfCutSuper1, Format Table, SendToServer "GetDP/R0", File StrCat[myDir, "temp_R0.txt"]];
-        Print[ eff_t~{nb_orders}[SurfCutSubs1] , OnRegion SurfCutSubs1 , Format Table, SendToServer "GetDP/T0", File StrCat[myDir, "temp_T0.txt"]];
-        Print[ Q_tot[Plot_domain]             , OnGlobal  , Format FrequencyTable ,SendToServer "GetDP/total absorption", File > StrCat[myDir, "absorption-Q_tot.txt"]];
-        For i In {0:N_rods-1:1}
-          Print[ Q_rod~{i}[rod~{i}] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, Sprintf("absorption-Q_rod_%g.txt", i+1) ]];           
-        EndFor  
-        Print[ Q_tot[Plot_domain]     , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_tot.txt"]];
-        Print[ Q_subs[sub]            , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_subs.txt"]];
-        Print[ Q_rod_out[rod_out]     , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_rod_out.txt"]];
-        Print[ Q_layer_dep[layer_dep] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_dep.txt"]];
-        Print[ Q_layer_cov[layer_cov] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_cov.txt"]];
         Print[ Hz_tot    , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm)];
-        Print[ u, OnElementsOf Omega, File StrCat[myDir, Sprintf("u2d_%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("u2d_%.2fnm.pos", lambda0/nm)];
-
+        // Print[ u, OnElementsOf Omega, File StrCat[myDir, Sprintf("u2d_%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("u2d_%.2fnm.pos", lambda0/nm)];
       EndIf
       If (flag_polar==0)
-        Print[ lambda_step, OnPoint{0,0,0}, Format Table, File StrCat[myDir, "temp_lambda_step.txt"], SendToServer "GetDP/Lambda_step" ] ;
-        For i In {0:2*nb_orders}
-          Print[ s_r~{i}[SurfCutSuper1], OnGlobal, Store i                , Format Table , File > StrCat[myDir, Sprintf("temp_s_r_%g.txt", i-nb_orders)]];
-          Print[ s_t~{i}[SurfCutSubs1] , OnGlobal, Store (2*nb_orders+1+i), Format Table , File > StrCat[myDir, Sprintf("temp_s_t_%g.txt", i-nb_orders)]];
-        EndFor
-        For i In {0:2*nb_orders}
-          Print[ eff_r~{i}[SurfCutSuper1], OnRegion SurfCutSuper1,Store (4*nb_orders+1+i), Format FrequencyTable, File > StrCat[myDir, Sprintf("efficiency_r_%g.txt", i-nb_orders)]];
-          Print[ eff_t~{i}[SurfCutSubs1] , OnRegion SurfCutSubs1 ,Store (6*nb_orders+1+i), Format FrequencyTable, File > StrCat[myDir, Sprintf("efficiency_t_%g.txt", i-nb_orders)]];
-          Print[ order_r_angle~{i}       , OnPoint{0,0,0}, Format Table , File > StrCat[myDir, Sprintf("order_r_angle_%g.txt", i-nb_orders)]];
-          Print[ order_t_angle~{i}       , OnPoint{0,0,0}, Format Table , File > StrCat[myDir, Sprintf("order_t_angle_%g.txt", i-nb_orders)]];          
-        EndFor
-        Print[ eff_r~{nb_orders}[SurfCutSuper1], OnRegion SurfCutSuper1, Format Table, SendToServer "GetDP/R0", File StrCat[myDir, "temp_R0.txt"]];
-        Print[ eff_t~{nb_orders}[SurfCutSubs1] , OnRegion SurfCutSubs1 , Format Table, SendToServer "GetDP/T0", File StrCat[myDir, "temp_T0.txt"]];
-        Print[ Q_tot[Plot_domain]              , OnGlobal  , Format FrequencyTable ,SendToServer "GetDP/total absorption", File > StrCat[myDir, "absorption-Q_tot.txt"]];
-        For i In {0:N_rods-1:1}
-          Print[ Q_rod~{i}[rod~{i}] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, Sprintf("absorption-Q_rod_%g.txt", i+1) ]];           
-        EndFor
-        Print[ Q_subs[sub]             , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_subs.txt"        ] ];
-        Print[ Q_rod_out[rod_out] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_rod_out.txt"] ];
-        Print[ Q_layer_dep[layer_dep] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_dep.txt"] ];
-        Print[ Q_layer_cov[layer_cov] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_cov.txt"] ];
         Print[ Ez_tot  , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm)];
-
-        Print[ u, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("u2d_%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("u2d_%.2fnm.pos", lambda0/nm)];
+        // Print[ u, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("u2d_%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("u2d_%.2fnm.pos", lambda0/nm)];
+      EndIf
+      If(multiplot)
+        Echo[ Str["For i In {PostProcessing.NbViews-1:0:-1}",
+                    "  If(!StrCmp(View[i].Name, 'boundary') || !StrCmp(View[i].Name, 'boundary_Combine'))",
+                    "    Delete View[i];",
+                    "  EndIf",
+                    "EndFor"], File StrCat[myDir,"tmp1.geo"]] ;
+        For i In {-nb_plot_periods:nb_plot_periods}
+          If (i!=0)
+            Print [ u_tot~{i}, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("u_tot_lambda%.2fnm_%g.pos", lambda0/nm,i+nb_plot_periods)], ChangeOfCoordinates {$X+i*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
+          EndIf
+          Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,Sprintf("boundary_%g.pos",i+nb_plot_periods)], ChangeOfCoordinates {$X+i*d,$Y,$Z} ] ;
+        EndFor
+        Echo[ Str["Combine ElementsByViewName;",
+                  "Hide {",
+                  "Point{1,2,7,8,9,10,20,22};",
+                  "Line{1,7,8,9,10,30,32,34,2,3,4,5,6,12,16,20,24,28};",
+                  "Surface{36,48};}",
+                  "Geometry.Color.Lines = {0,0,0};",
+                  "l=PostProcessing.NbViews-1; View[l].ColorTable={Black}; ",
+                  "View[l-1].Visible=1; View[l-1].ShowScale=0;",
+                  "View[l].ShowScale=0; View[l].LineWidth=1.5; View[l].LineType=1;Geometry.LineWidth=0;"],
+                  File StrCat[myDir,"tmp3.geo" ]] ;
       EndIf
+
     }
   }
 }