diff --git a/DiffractionGratings/grating2D.pro b/DiffractionGratings/grating2D.pro
index 1f73550b2791d676c33d90d08832e5bc39486074..e5e9b8dec4ed049e3b3aee758222117d9ce398d0 100644
--- a/DiffractionGratings/grating2D.pro
+++ b/DiffractionGratings/grating2D.pro
@@ -118,13 +118,13 @@ Function{
   k0         = 2.*Pi/lambda0;
   epsr_sup[] = Complex[epsr_sup_re[],epsr_sup_im[]];
   epsr_sub[] = Complex[epsr_sub_re[],epsr_sub_im[]];
-  n_sup[]    = (epsr_sup[])^(0.5);
-  n_sub[]    = (epsr_sub[])^(0.5);
+  n_sup[]    = Sqrt[epsr_sup[]];
+  n_sub[]    = Sqrt[epsr_sub[]];
   k_sup[]    = k0*n_sup[];
   k_sub[]    = k0*n_sub[];
   alpha[]    =  k_sup[]*Sin[theta_deg*deg2rad];
   beta_sup[] =  k_sup[]*Cos[theta_deg*deg2rad];
-  beta_sub[] = (k0^2*epsr_sub[]-alpha[]^2)^(0.5);
+  beta_sub[] = Sqrt[k0^2*epsr_sub[]-alpha[]^2];
   
   If (flag_Hparallel==1)
     beta_pol_sup[] = beta_sup[]/epsr_sup[];
@@ -143,8 +143,8 @@ Function{
   For i In {0:2*nb_orders}
   alpha_orders~{i}[]  = Complex[alpha[]  + 2.*Pi/d*(i-nb_orders),0];
   expialpha_orders~{i}[]= Complex[ Cos[Re[alpha_orders~{i}[]]*X[]] , Sin[Re[alpha_orders~{i}[]]*X[]] ];
-  betat_sup~{i}[] = (k_sup[]^2-alpha_orders~{i}[]^2)^(0.5);
-  betat_sub~{i}[] = (k_sub[]^2-alpha_orders~{i}[]^2)^(0.5);
+  betat_sup~{i}[] = Sqrt[k_sup[]^2-alpha_orders~{i}[]^2];
+  betat_sub~{i}[] = Sqrt[k_sub[]^2-alpha_orders~{i}[]^2];
   EndFor
   
   a_pml        = 1.;
@@ -371,112 +371,112 @@ PostProcessing {
   { Name postpro_energy; NameOfFormulation helmoltz_scalar;
     Quantity {
       If (flag_Hparallel==1)
-  { Name debr         ; Value { Local { [ r[]                 ]; In Omega; Jacobian JVol; } } }
-      { Name debt         ; Value { Local { [ t[]                 ]; In Omega; Jacobian JVol; } } }
-      { Name testtm       ; Value { Local { [ {u2d}               ]; In Omega; Jacobian JVol; } } }
-      { Name deb_beta_sub ; Value { Local { [ beta_sub[]          ]; In Omega; Jacobian JVol; } } }
-      { Name epsr         ; Value { Local { [ CompZZ[epsilonr[]]  ]; In Omega; Jacobian JVol; } } }
-      { Name Hz_diff      ; Value { Local { [ {u2d}+u1d[]         ]; In Omega; Jacobian JVol; } } }
-      { Name Hz_tot       ; Value { Local { [ {u2d}+u1[]          ]; In Omega; Jacobian JVol; } } }
-      { Name NormHz_tot   ; Value { Local { [ Norm[{u2d}+u1[]]          ]; In Omega; Jacobian JVol; } } }
-      { Name Hz_totp1     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 1*-alpha[]*d],Sin[ 1*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name Hz_totp2     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 2*-alpha[]*d],Sin[ 2*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name Hz_totp3     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 3*-alpha[]*d],Sin[ 3*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name Hz_totp4     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 4*-alpha[]*d],Sin[ 4*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name Hz_totm1     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-1*-alpha[]*d],Sin[-1*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name Hz_totm2     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-2*-alpha[]*d],Sin[-2*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name Hz_totm3     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-3*-alpha[]*d],Sin[-3*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name Hz_totm4     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-4*-alpha[]*d],Sin[-4*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name boundary     ; Value { Local { [ bndCol[] ] ; In Plot_bnd ; Jacobian JVol ; } } }
+        { Name debr         ; Value { Local { [ r[]                 ]; In Omega; Jacobian JVol; } } }
+        { Name debt         ; Value { Local { [ t[]                 ]; In Omega; Jacobian JVol; } } }
+        { Name testtm       ; Value { Local { [ {u2d}               ]; In Omega; Jacobian JVol; } } }
+        { Name deb_beta_sub ; Value { Local { [ beta_sub[]          ]; In Omega; Jacobian JVol; } } }
+        { Name epsr         ; Value { Local { [ CompZZ[epsilonr[]]  ]; In Omega; Jacobian JVol; } } }
+        { Name Hz_diff      ; Value { Local { [ {u2d}+u1d[]         ]; In Omega; Jacobian JVol; } } }
+        { Name Hz_tot       ; Value { Local { [ {u2d}+u1[]          ]; In Omega; Jacobian JVol; } } }
+        { Name NormHz_tot   ; Value { Local { [ Norm[{u2d}+u1[]]          ]; In Omega; Jacobian JVol; } } }
+        { Name Hz_totp1     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 1*-alpha[]*d],Sin[ 1*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+        { Name Hz_totp2     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 2*-alpha[]*d],Sin[ 2*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+        { Name Hz_totp3     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 3*-alpha[]*d],Sin[ 3*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+        { Name Hz_totp4     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 4*-alpha[]*d],Sin[ 4*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+        { Name Hz_totm1     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-1*-alpha[]*d],Sin[-1*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+        { Name Hz_totm2     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-2*-alpha[]*d],Sin[-2*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+        { Name Hz_totm3     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-3*-alpha[]*d],Sin[-3*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+        { Name Hz_totm4     ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-4*-alpha[]*d],Sin[-4*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+        { Name boundary     ; Value { Local { [ bndCol[] ] ; In Plot_bnd ; Jacobian JVol ; } } }
         
-      { Name Ex_diff_re; Value { Local { [ Re[Exi[]]         ]; In Omega; Jacobian JVol; } } }
-      { Name Ex_diff_im; Value { Local { [ Im[Exi[]]         ]; In Omega; Jacobian JVol; } } }
-      { Name Ey_diff_re; Value { Local { [ Re[Eyi[]]         ]; In Omega; Jacobian JVol; } } }
-      { Name Ey_diff_im; Value { Local { [ Im[Eyi[]]         ]; In Omega; Jacobian JVol; } } }
-      { Name En_tot    ; Value { Local { [ Log10[Sqrt[SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Exr[]/CompXX[epsilonr[]]+Exi[]/CompXX[epsilonr[]] + Ext[] ] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Eyr[]/CompXX[epsilonr[]]+Eyi[]/CompXX[epsilonr[]] + Eyt[]] ] ] ]; In Omega; Jacobian JVol; } } }
+        { Name Ex_diff_re; Value { Local { [ Re[Exi[]]         ]; In Omega; Jacobian JVol; } } }
+        { Name Ex_diff_im; Value { Local { [ Im[Exi[]]         ]; In Omega; Jacobian JVol; } } }
+        { Name Ey_diff_re; Value { Local { [ Re[Eyi[]]         ]; In Omega; Jacobian JVol; } } }
+        { Name Ey_diff_im; Value { Local { [ Im[Eyi[]]         ]; In Omega; Jacobian JVol; } } }
+        { Name En_tot    ; Value { Local { [ Log10[Sqrt[SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Exr[]/CompXX[epsilonr[]]+Exi[]/CompXX[epsilonr[]] + Ext[] ] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Eyr[]/CompXX[epsilonr[]]+Eyi[]/CompXX[epsilonr[]] + Eyt[]] ] ] ]; In Omega; Jacobian JVol; } } }
 
-      For i In {0:2*nb_orders}
-      { Name s_r~{i} ; Value { 
-    Integral{ [ expialpha_orders~{i}[] * ({u2d}+u1d[])/d ] ;
-      In SurfCutSuper1 ; Jacobian JSur ; Integration Int_1 ; } } }
-      { Name s_t~{i} ; Value { 
-    Integral{ [ expialpha_orders~{i}[] * ({u2d}+u1d[])/d ] ;
-      In SurfCutSubs1  ; Jacobian JSur ; Integration Int_1 ; } } }
-      { Name order_t_angle~{i} ; Value { 
-    Local{ [-Atan2[Re[alpha_orders~{i}[]],Re[betat_sub~{i}[]]]/deg2rad ] ;
-      In Omega; Jacobian JVol; } } }
-      { Name order_r_angle~{i} ; Value { 
-    Local{ [ Atan2[Re[alpha_orders~{i}[]],Re[betat_sup~{i}[]]]/deg2rad ] ;
-      In Omega; Jacobian JVol; } } }
-      EndFor
         For i In {0:2*nb_orders}
-      { Name eff_r~{i} ; Value {
-    Term{ Type Global; [ SquNorm[#i]*betat_sup~{i}[]/beta_sup[] ] ;
-      In SurfCutSuper1 ; } } }
-      { Name eff_t~{i} ; Value {
-    Term{ Type Global; [ SquNorm[#(2*nb_orders+1+i)]*(betat_sub~{i}[]/beta_sup[])*(epsr_sup[]/epsr_sub[]) ] ;
-      In SurfCutSubs1 ; } } }
-      EndFor
-  For i In {0:N_rods-1:1}
-      { Name Q_rod~{i}   ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_rods_im[]]        * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]]] ) / (Pinc[]*d) ] ; In rod~{i}   ; Integration Int_1 ; Jacobian JVol ; } } }
-      EndFor
-  { Name Q_sub         ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_sub_im[]      ]     * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]] ] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In sub         ; Integration Int_1 ; Jacobian JVol ; } } }
-      { Name Q_rod_out     ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_rod_out_im[]]     * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]] ] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In rod_out   ; Integration Int_1 ; Jacobian JVol ; } } }
-      { Name Q_layer_dep   ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_layer_dep_im[]]     * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]] ] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In layer_dep   ; Integration Int_1 ; Jacobian JVol ; } } }
-      { Name Q_layer_cov   ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_layer_cov_im[]]     * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]] ] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In layer_cov   ; Integration Int_1 ; Jacobian JVol ; } } }
-      { Name Q_tot         ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[Im[CompZZ[epsilonr[]]]]  * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]] ] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In Plot_domain ; Integration Int_1 ; Jacobian JVol ; } } }
-      { Name lambda_step   ; Value { Local { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
+          { Name s_r~{i} ; Value { 
+            Integral{ [ expialpha_orders~{i}[] * ({u2d}+u1d[])/d ] ;
+            In SurfCutSuper1 ; Jacobian JSur ; Integration Int_1 ; } } }
+          { Name s_t~{i} ; Value { 
+            Integral{ [ expialpha_orders~{i}[] * ({u2d}+u1d[])/d ] ;
+            In SurfCutSubs1  ; Jacobian JSur ; Integration Int_1 ; } } }
+          { Name order_t_angle~{i} ; Value { 
+            Local{ [-Atan2[Re[alpha_orders~{i}[]],Re[betat_sub~{i}[]]]/deg2rad ] ;
+            In Omega; Jacobian JVol; } } }
+          { Name order_r_angle~{i} ; Value { 
+            Local{ [ Atan2[Re[alpha_orders~{i}[]],Re[betat_sup~{i}[]]]/deg2rad ] ;
+            In Omega; Jacobian JVol; } } }
+        EndFor
+        For i In {0:2*nb_orders}
+          { Name eff_r~{i} ; Value {
+            Term{ Type Global; [ SquNorm[#i]*betat_sup~{i}[]/beta_sup[] ] ;
+            In SurfCutSuper1 ; } } }
+          { Name eff_t~{i} ; Value {
+            Term{ Type Global; [ SquNorm[#(2*nb_orders+1+i)]*(betat_sub~{i}[]/beta_sup[])*(epsr_sup[]/epsr_sub[]) ] ;
+            In SurfCutSubs1 ; } } }
+        EndFor
+        For i In {0:N_rods-1:1}
+          { Name Q_rod~{i}   ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_rods_im[]]         * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]]] ) / (Pinc[]*d) ] ; In rod~{i}   ; Integration Int_1 ; Jacobian JVol ; } } }
+        EndFor
+        { Name Q_sub         ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_sub_im[]]           * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]] ] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In sub         ; Integration Int_1 ; Jacobian JVol ; } } }
+        { Name Q_rod_out     ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_rod_out_im[]]       * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]] ] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In rod_out   ; Integration Int_1 ; Jacobian JVol ; } } }
+        { Name Q_layer_dep   ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_layer_dep_im[]]     * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]] ] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In layer_dep   ; Integration Int_1 ; Jacobian JVol ; } } }
+        { Name Q_layer_cov   ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_layer_cov_im[]]     * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]] ] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In layer_cov   ; Integration Int_1 ; Jacobian JVol ; } } }
+        { Name Q_tot         ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[Im[CompZZ[epsilonr[]]]]  * ( SquNorm[-CompY[{Grad u2d}]*I[]/(omega0*epsilon0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]] ] + SquNorm[CompX[{Grad u2d}]*I[]/(omega0*epsilon0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In Plot_domain ; Integration Int_1 ; Jacobian JVol ; } } }
+        { Name lambda_step   ; Value { Local { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
       EndIf
   If (flag_Hparallel==0)
-      { Name testte    ; Value { Local { [ {u2d} ]; In Omega; Jacobian JVol; } } }
-      { Name epsr      ; Value { Local { [ CompZZ[epsilonr[]] ]; In Omega; Jacobian JVol; } } }
-      { Name Ez_diff   ; Value { Local { [ {u2d}+u1d[]        ]; In Omega; Jacobian JVol; } } }
-      { Name Ez_tot    ; Value { Local { [ {u2d}+u1[]          ]; In Omega; Jacobian JVol; } } }
-      { Name Ez_totp1  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 1*-alpha[]*d],Sin[ 1*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name Ez_totp2  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 2*-alpha[]*d],Sin[ 2*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name Ez_totp3  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 3*-alpha[]*d],Sin[ 3*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name Ez_totp4  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 4*-alpha[]*d],Sin[ 4*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name Ez_totm1  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-1*-alpha[]*d],Sin[-1*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name Ez_totm2  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-2*-alpha[]*d],Sin[-2*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name Ez_totm3  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-3*-alpha[]*d],Sin[-3*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name Ez_totm4  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-4*-alpha[]*d],Sin[-4*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
-      { Name boundary  ; Value { Local { [ bndCol[] ] ; In Plot_bnd ; Jacobian JVol ; } } }
-      { Name u         ; Value { Local { [ {u2d}    ]; In Omega; Jacobian JVol; } } }
-        
-      // modif effic
-      For i In {0:2*nb_orders}
+    { Name testte    ; Value { Local { [ {u2d} ]; In Omega; Jacobian JVol; } } }
+    { Name epsr      ; Value { Local { [ CompZZ[epsilonr[]] ]; In Omega; Jacobian JVol; } } }
+    { Name Ez_diff   ; Value { Local { [ {u2d}+u1d[]        ]; In Omega; Jacobian JVol; } } }
+    { Name Ez_tot    ; Value { Local { [ {u2d}+u1[]          ]; In Omega; Jacobian JVol; } } }
+    { Name Ez_totp1  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 1*-alpha[]*d],Sin[ 1*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+    { Name Ez_totp2  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 2*-alpha[]*d],Sin[ 2*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+    { Name Ez_totp3  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 3*-alpha[]*d],Sin[ 3*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+    { Name Ez_totp4  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[ 4*-alpha[]*d],Sin[ 4*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+    { Name Ez_totm1  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-1*-alpha[]*d],Sin[-1*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+    { Name Ez_totm2  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-2*-alpha[]*d],Sin[-2*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+    { Name Ez_totm3  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-3*-alpha[]*d],Sin[-3*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+    { Name Ez_totm4  ; Value { Local { [ ({u2d}+u1[])*Complex[Cos[-4*-alpha[]*d],Sin[-4*-alpha[]*d]] ]; In Omega; Jacobian JVol; } } }
+    { Name boundary  ; Value { Local { [ bndCol[] ] ; In Plot_bnd ; Jacobian JVol ; } } }
+    { Name u         ; Value { Local { [ {u2d}    ]; In Omega; Jacobian JVol; } } }
+      
+    // modif effic
+    For i In {0:2*nb_orders}
       { Name s_r~{i} ; Value {
-    Integral{ [ expialpha_orders~{i}[] * ({u2d}+u1d[])/d ] ;
-      In SurfCutSuper1 ; Jacobian JSur ; Integration Int_1 ; } } }
+        Integral{ [ expialpha_orders~{i}[] * ({u2d}+u1d[])/d ] ;
+        In SurfCutSuper1 ; Jacobian JSur ; Integration Int_1 ; } } }
       { Name s_t~{i} ; Value { 
-    Integral{ [ expialpha_orders~{i}[] * ({u2d}+u1d[])/d ] ;
-      In SurfCutSubs1  ; Jacobian JSur ; Integration Int_1 ; } } }
+        Integral{ [ expialpha_orders~{i}[] * ({u2d}+u1d[])/d ] ;
+        In SurfCutSubs1  ; Jacobian JSur ; Integration Int_1 ; } } }
       { Name order_t_angle~{i} ; Value { 
-    Local{ [-Atan2[Re[alpha_orders~{i}[]],Re[betat_sub~{i}[]]]/deg2rad ] ;
-      In Omega; Jacobian JVol; } } }
+        Local{ [-Atan2[Re[alpha_orders~{i}[]],Re[betat_sub~{i}[]]]/deg2rad ] ;
+        In Omega; Jacobian JVol; } } }
       { Name order_r_angle~{i} ; Value { 
-    Local{ [ Atan2[Re[alpha_orders~{i}[]],Re[betat_sup~{i}[]]]/deg2rad ] ;
-      In Omega; Jacobian JVol; } } }
-      EndFor        
-        For i In {0:2*nb_orders}
+        Local{ [ Atan2[Re[alpha_orders~{i}[]],Re[betat_sup~{i}[]]]/deg2rad ] ;
+        In Omega; Jacobian JVol; } } }
+    EndFor        
+    For i In {0:2*nb_orders}
       { Name eff_r~{i} ; Value {
-    Term{ Type Global; [ SquNorm[#i]*betat_sup~{i}[]/beta_sup[] ] ;
-      In SurfCutSuper1 ; } } }
+        Term{ Type Global; [ SquNorm[#i]*betat_sup~{i}[]/beta_sup[] ] ;
+        In SurfCutSuper1 ; } } }
       { Name eff_t~{i} ; Value {
-    Term{ Type Global; [ SquNorm[#(2*nb_orders+1+i)]*(betat_sub~{i}[]/beta_sup[])] ;
-      In SurfCutSubs1 ; } } }
-      EndFor
-  For i In {0:N_rods-1:1}
+        Term{ Type Global; [ SquNorm[#(2*nb_orders+1+i)]*(betat_sub~{i}[]/beta_sup[])] ;
+        In SurfCutSubs1 ; } } }
+    EndFor
+    For i In {0:N_rods-1:1}
       { Name Q_rod~{i}  ; Value { Integral { [0.5 * epsilon0*omega0*Fabs[epsr_rods_im[]]            * (SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In rod~{i}    ; Integration Int_1 ; Jacobian JVol ; } } }
-      EndFor
-        { Name Q_sub        ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_sub_im[]]           * (SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In sub         ; Integration Int_1 ; Jacobian JVol ; } } }
-      { Name Q_rod_out    ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_rod_out_im[]]       * (SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In rod_out   ; Integration Int_1 ; Jacobian JVol ; } } }
-      { Name Q_layer_dep  ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_layer_dep_im[]]     * (SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In layer_dep   ; Integration Int_1 ; Jacobian JVol ; } } }
-      { Name Q_layer_cov  ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_layer_cov_im[]]     * (SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In layer_cov   ; Integration Int_1 ; Jacobian JVol ; } } }
-      { Name Q_tot        ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[Im[CompXX[epsilonr[]]]]  * (SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In Plot_domain ; Integration Int_1 ; Jacobian JVol ; } } }
-      { Name lambda_step  ; Value { Local    { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
-      EndIf
-  }
+    EndFor
+    { Name Q_sub        ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_sub_im[]]           * (SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In sub         ; Integration Int_1 ; Jacobian JVol ; } } }
+    { Name Q_rod_out    ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_rod_out_im[]]       * (SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In rod_out   ; Integration Int_1 ; Jacobian JVol ; } } }
+    { Name Q_layer_dep  ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_layer_dep_im[]]     * (SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In layer_dep   ; Integration Int_1 ; Jacobian JVol ; } } }
+    { Name Q_layer_cov  ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[epsr_layer_cov_im[]]     * (SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In layer_cov   ; Integration Int_1 ; Jacobian JVol ; } } }
+    { Name Q_tot        ; Value { Integral { [  0.5 * epsilon0*omega0*Fabs[Im[CompXX[epsilonr[]]]]  * (SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In Plot_domain ; Integration Int_1 ; Jacobian JVol ; } } }
+    { Name lambda_step  ; Value { Local    { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
+    EndIf
+    }
   }
 }
 
@@ -499,9 +499,9 @@ PostOperation {
       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) ]];           
+        Print[ Q_rod~{i}[rod~{i}] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, Sprintf("absorption-Q_rod_%g.txt", i+1) ]];           
       EndFor  
-  Print[ Q_sub[sub]         , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_sub.txt"      ]];
+      Print[ Q_sub[sub]         , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_sub.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"]];
@@ -511,106 +511,106 @@ PostOperation {
       //                   "EndFor"], File StrCat[myDir,"tmp0.geo"]] ;
       // // View[i].RangeType = 2;View[i].CustomMin=0.;View[i].CustomMax=1.;
       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"]] ;
-      Print [ Hz_totp1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_2.pos", lambda0/nm)], ChangeOfCoordinates {$X+1*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Print [ Hz_totm1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_3.pos", lambda0/nm)], ChangeOfCoordinates {$X-1*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Print [ Hz_totp2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_4.pos", lambda0/nm)], ChangeOfCoordinates {$X+2*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Print [ Hz_totm2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_5.pos", lambda0/nm)], ChangeOfCoordinates {$X-2*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Print [ Hz_totp3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_6.pos", lambda0/nm)], ChangeOfCoordinates {$X+3*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Print [ Hz_totm3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_7.pos", lambda0/nm)], ChangeOfCoordinates {$X-3*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Print [ Hz_totp4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_8.pos", lambda0/nm)], ChangeOfCoordinates {$X+4*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Print [ Hz_totm4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_9.pos", lambda0/nm)], ChangeOfCoordinates {$X-4*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo"]] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary1.pos"], ChangeOfCoordinates {$X+0*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary2.pos"], ChangeOfCoordinates {$X+1*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary3.pos"], ChangeOfCoordinates {$X-1*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary4.pos"], ChangeOfCoordinates {$X+2*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary5.pos"], ChangeOfCoordinates {$X-2*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary6.pos"], ChangeOfCoordinates {$X+3*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary7.pos"], ChangeOfCoordinates {$X-3*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary8.pos"], ChangeOfCoordinates {$X+4*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary9.pos"], ChangeOfCoordinates {$X-4*d,$Y,$Z} ];
-      Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo" ]] ;
-      Echo[ Str["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" ]] ;
+        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"]] ;
+        Print [ Hz_totp1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_2.pos", lambda0/nm)], ChangeOfCoordinates {$X+1*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Print [ Hz_totm1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_3.pos", lambda0/nm)], ChangeOfCoordinates {$X-1*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Print [ Hz_totp2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_4.pos", lambda0/nm)], ChangeOfCoordinates {$X+2*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Print [ Hz_totm2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_5.pos", lambda0/nm)], ChangeOfCoordinates {$X-2*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Print [ Hz_totp3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_6.pos", lambda0/nm)], ChangeOfCoordinates {$X+3*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Print [ Hz_totm3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_7.pos", lambda0/nm)], ChangeOfCoordinates {$X-3*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Print [ Hz_totp4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_8.pos", lambda0/nm)], ChangeOfCoordinates {$X+4*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Print [ Hz_totm4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_9.pos", lambda0/nm)], ChangeOfCoordinates {$X-4*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo"]] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary1.pos"], ChangeOfCoordinates {$X+0*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary2.pos"], ChangeOfCoordinates {$X+1*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary3.pos"], ChangeOfCoordinates {$X-1*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary4.pos"], ChangeOfCoordinates {$X+2*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary5.pos"], ChangeOfCoordinates {$X-2*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary6.pos"], ChangeOfCoordinates {$X+3*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary7.pos"], ChangeOfCoordinates {$X-3*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary8.pos"], ChangeOfCoordinates {$X+4*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary9.pos"], ChangeOfCoordinates {$X-4*d,$Y,$Z} ];
+        Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo" ]] ;
+        Echo[ Str["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
   EndIf
   If (flag_Hparallel==0)
-        // Print [ testte  , OnElementsOf Omega, File "testte.pos"];
-        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 {1:2*nb_orders-1}
-        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 [ testte  , OnElementsOf Omega, File "testte.pos"];
+    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 {1:2*nb_orders-1}
+      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_sub[sub]             , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_sub.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)];
-      //         Echo[ Str["For i In {0:2}",
-      //          " View[i].LineWidth = 4;View[i].ColormapNumber = 15;View[0].AxesFormatX = '%.1f';",
-      //                   "EndFor"], File StrCat[myDir,"tmp0.geo"]] ;
-      // // View[i].RangeType = 2;View[i].CustomMin =0.;View[i].CustomMax =1.;
-      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"]] ;
-      Print [ Ez_totp1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_2.pos", lambda0/nm)], ChangeOfCoordinates {$X+1*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Print [ Ez_totm1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_3.pos", lambda0/nm)], ChangeOfCoordinates {$X-1*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Print [ Ez_totp2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_4.pos", lambda0/nm)], ChangeOfCoordinates {$X+2*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Print [ Ez_totm2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_5.pos", lambda0/nm)], ChangeOfCoordinates {$X-2*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Print [ Ez_totp3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_6.pos", lambda0/nm)], ChangeOfCoordinates {$X+3*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Print [ Ez_totm3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_7.pos", lambda0/nm)], ChangeOfCoordinates {$X-3*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Print [ Ez_totp4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_8.pos", lambda0/nm)], ChangeOfCoordinates {$X+4*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Print [ Ez_totm4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_9.pos", lambda0/nm)], ChangeOfCoordinates {$X-4*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-      Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo"]] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary1.pos"], ChangeOfCoordinates {$X+0*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary2.pos"], ChangeOfCoordinates {$X+1*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary3.pos"], ChangeOfCoordinates {$X-1*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary4.pos"], ChangeOfCoordinates {$X+2*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary5.pos"], ChangeOfCoordinates {$X-2*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary6.pos"], ChangeOfCoordinates {$X+3*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary7.pos"], ChangeOfCoordinates {$X-3*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary8.pos"], ChangeOfCoordinates {$X+4*d,$Y,$Z} ] ;
-      Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary9.pos"], ChangeOfCoordinates {$X-4*d,$Y,$Z} ];
-      Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo" ]] ;
-      Echo[ Str["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" ]] ;
+    EndFor
+    Print[ Q_sub[sub]             , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_sub.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)];
+    //         Echo[ Str["For i In {0:2}",
+    //          " View[i].LineWidth = 4;View[i].ColormapNumber = 15;View[0].AxesFormatX = '%.1f';",
+    //                   "EndFor"], File StrCat[myDir,"tmp0.geo"]] ;
+    // // View[i].RangeType = 2;View[i].CustomMin =0.;View[i].CustomMax =1.;
+    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"]] ;
+        Print [ Ez_totp1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_2.pos", lambda0/nm)], ChangeOfCoordinates {$X+1*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Print [ Ez_totm1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_3.pos", lambda0/nm)], ChangeOfCoordinates {$X-1*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Print [ Ez_totp2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_4.pos", lambda0/nm)], ChangeOfCoordinates {$X+2*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Print [ Ez_totm2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_5.pos", lambda0/nm)], ChangeOfCoordinates {$X-2*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Print [ Ez_totp3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_6.pos", lambda0/nm)], ChangeOfCoordinates {$X+3*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Print [ Ez_totm3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_7.pos", lambda0/nm)], ChangeOfCoordinates {$X-3*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Print [ Ez_totp4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_8.pos", lambda0/nm)], ChangeOfCoordinates {$X+4*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Print [ Ez_totm4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_9.pos", lambda0/nm)], ChangeOfCoordinates {$X-4*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
+        Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo"]] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary1.pos"], ChangeOfCoordinates {$X+0*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary2.pos"], ChangeOfCoordinates {$X+1*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary3.pos"], ChangeOfCoordinates {$X-1*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary4.pos"], ChangeOfCoordinates {$X+2*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary5.pos"], ChangeOfCoordinates {$X-2*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary6.pos"], ChangeOfCoordinates {$X+3*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary7.pos"], ChangeOfCoordinates {$X-3*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary8.pos"], ChangeOfCoordinates {$X+4*d,$Y,$Z} ] ;
+        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary9.pos"], ChangeOfCoordinates {$X-4*d,$Y,$Z} ];
+        Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo" ]] ;
+        Echo[ Str["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        
-      EndIf
-  }
+    EndIf
+    }
   }
 }
 DefineConstant[