diff --git a/DiffractionGratings/grating2D.geo b/DiffractionGratings/grating2D.geo
index e56fe9e8c7d5cfb583a88af589a981356b272a90..db87e8e1ed5751096e727cf398dfd0d0b9baf626 100644
--- a/DiffractionGratings/grating2D.geo
+++ b/DiffractionGratings/grating2D.geo
@@ -270,3 +270,4 @@ EndIf
 Solver.AutoMesh=1;
 Geometry.Points = 0;
 Mesh.SurfaceEdges = 0;
+// Mesh.MshFileVersion = 2.0;
\ No newline at end of file
diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index d474b5e4022220ae020cf1ef02a613a3c1aeb5e5..1a4741b4e42e1c0c1628d35b95791015aad4a1e9 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -200,8 +200,8 @@ Function{
   H1d[Omega_super] = Hr[];
   H1d[Omega_subs]  = Ht[];
 
-  source[] = (om0/cel)^2*(epsr[]-epsr_annex[])*E1[];
-
+  source_vol_scat[] = (om0/cel)^2*(epsr[]-epsr_annex[])*E1[];
+  source_surf_tot[] = -2*I[]*om0*mu0*Vector[0,0,1] /\ Hi[];
   // Bloch phase shifts
   skx1[] =  k1x[];
   // sky1[] = -k0*n1[]*Sin[theta0]*Sin[phi0+xsi];
@@ -324,7 +324,11 @@ Formulation {
 		Equation{
       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; }
+      If (FLAG_TOTAL==0)
+        Galerkin { [source_vol_scat[] ,      {u}]; In Omega_source; Jacobian JVol; Integration I1; }
+      Else
+        Galerkin { [source_surf_tot[] ,      {u}]; In SurfIntTop; Jacobian JVol; Integration I1; }
+      EndIf
     }
   }
 }
@@ -346,40 +350,54 @@ PostProcessing {
   { Name postpro_helmholtz_vector; NameOfFormulation helmholtz_vector; NameOfSystem M;
     Quantity {
       { Name u      ; Value { Local { [ {u}       ]; In Omega; Jacobian JVol; } } }
-      { Name Etot   ; Value { Local { [ {u}+E1[]  ]; In Omega; Jacobian JVol; } } }
-      { Name Htot   ; Value { Local { [ H1[]-I[]/(mur[]*mu0*om0)*{Curl u}]; In Omega; Jacobian JVol; } } }
-      { Name Edif   ; Value { Local { [ {u}+E1d[] ]; In Omega; Jacobian JVol; } } }
-      { Name E1     ; Value { Local { [     E1[]  ]; In Omega; Jacobian JVol; } } }
-      { Name Em     ; Value { Local { [     E1d[] ]; In Omega; Jacobian JVol; } } }
-      { Name source ; Value { Local { [  source[] ]; In Omega; Jacobian JVol; } } }
-      { Name epsr_xx; Value { Local { [  CompXX[epsr[]] ]; In Omega; Jacobian JVol; } } }
       { Name Damp_pml_top; Value { Local { [Damp_pml_top[]  ]; In Omega; Jacobian JVol; } } }
-      { Name Poy_tot; Value { Local { [ 0.5*Re[Cross[{u}+E1[] , Conj[ H1[]-I[]/(mur[]*mu0*om0)*{Curl u}]]] ]; In Omega; Jacobian JVol; } } }
-      { Name Poy_ref; Value { Local { [ 0.5*Re[Cross[{u}+E1d[], Conj[H1d[]-I[]/(mur[]*mu0*om0)*{Curl u}]]] ]; In Omega; Jacobian JVol; } } }
-      { Name Poy_inc; Value { Local { [ 0.5*Re[Cross[    Ei[] , Conj[Hi[]]]                              ] ]; In Omega; Jacobian JVol; } } }
+      { Name epsr_xx; Value { Local { [  CompXX[epsr[]] ]; In Omega; Jacobian JVol; } } }
+      { Name Poy_inc; Value { Local { [ 0.5*Re[Cross[    Ei[] , Conj[Hi[]]]] ]; In Omega; Jacobian JVol; } } }
+      { Name E1     ; Value { Local { [     E1[]  ]; In Omega; Jacobian JVol; } } }
       { Name lambda_step   ; Value { Local { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
-
-      For k In {2:6}
-        { Name Abs_L~{k} ; Value { Integral { [ ep0*om0 * 0.5*Im[CompXX[epsr[]]]*(SquNorm[{u}+E1[]]) / (Pinc*period_x*dyc) ] ; In L~{k} ; Integration I1 ; Jacobian JVol ; } } }
-      EndFor
-      { Name Abs_scat ; Value { Integral { [ ep0*om0 * 0.5*Im[CompXX[epsr[]]]*(SquNorm[{u}+E1[]]) / (Pinc*period_x*dyc) ] ; In Scat ; Integration I1 ; Jacobian JVol ; } } }
-      { Name Abs_scat2; Value { Integral { [ ep0*om0 * 0.5*Im[CompXX[epsr[]]]*(SquNorm[{u}+E1[]]) ] ; In Scat ; Integration I1 ; Jacobian JVol ; } } }
+      
+      If (FLAG_TOTAL==0)
+        { Name Etot   ; Value { Local { [ {u}+E1[]  ]; In Omega; Jacobian JVol; } } }
+        { Name Htot   ; Value { Local { [ H1[]-I[]/(mur[]*mu0*om0)*{Curl u}]; In Omega; Jacobian JVol; } } }
+        { Name Poy_tot; Value { Local { [ 0.5*Re[Cross[{u}+E1[] , Conj[ H1[]-I[]/(mur[]*mu0*om0)*{Curl u}]]] ]; In Omega; Jacobian JVol; } } }
+        { Name Poy_ref; Value { Local { [ 0.5*Re[Cross[{u}+E1d[], Conj[H1d[]-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*Im[CompXX[epsr[]]]*(SquNorm[{u}+E1[]]) / (Pinc*period_x*dyc) ] ; In L~{k} ; Integration I1 ; Jacobian JVol ; } } }
+        EndFor
+        { Name Abs_scat ; Value { Integral { [ ep0*om0 * 0.5*Im[CompXX[epsr[]]]*(SquNorm[{u}+E1[]]) / (Pinc*period_x*dyc) ] ; In Scat ; Integration I1 ; Jacobian JVol ; } } }
+      Else
+        { Name Etot   ; Value { Local { [ {u}  ]; In Omega; Jacobian JVol; } } }
+        { Name Htot   ; Value { Local { [ -I[]/(mur[]*mu0*om0)*{Curl u}]; In Omega; Jacobian JVol; } } }
+        { Name Poy_tot; Value { Local { [ 0.5*Re[Cross[{u} , Conj[ -I[]/(mur[]*mu0*om0)*{Curl u}]]] ]; In Omega; Jacobian JVol; } } }
+        { Name Poy_ref; Value { Local { [ 0.5*Re[Cross[{u}-Ei[], Conj[-Hi[]-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*Im[CompXX[epsr[]]]*(SquNorm[{u}]) / (Pinc*period_x*dyc) ] ; In L~{k} ; Integration I1 ; Jacobian JVol ; } } }
+        EndFor
+        { Name Abs_scat ; Value { Integral { [ ep0*om0 * 0.5*Im[CompXX[epsr[]]]*(SquNorm[{u}]) / (Pinc*period_x*dyc) ] ; In Scat ; Integration I1 ; Jacobian JVol ; } } }
+      EndIf
 
       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[] ]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
-          { Name int_y_t~{i}~{j} ; Value { Integral { [ CompY[{u}+E1[] ]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
-          { Name int_x_r~{i}~{j} ; Value { Integral { [ CompX[{u}+E1d[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
-          { Name int_y_r~{i}~{j} ; Value { Integral { [ CompY[{u}+E1d[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
-          { Name eff_t~{i}~{j}   ; Value { Term{Type Global; [
-                 1/(gammat~{i}~{j}[]*-k1z[]*Cos[xsi]^2) * ( (gammat~{i}~{j}[]^2+alpha~{i}~{j}[]^2)*SquNorm[$int_x_t~{i}~{j}]+
+          If (FLAG_TOTAL==0)
+            { Name int_x_t~{i}~{j} ; Value { Integral { [ CompX[{u}+E1[] ]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
+            { Name int_y_t~{i}~{j} ; Value { Integral { [ CompY[{u}+E1[] ]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
+            { Name int_x_r~{i}~{j} ; Value { Integral { [ CompX[{u}+E1d[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
+            { Name int_y_r~{i}~{j} ; Value { Integral { [ CompY[{u}+E1d[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
+          Else
+            { Name int_x_t~{i}~{j} ; Value { Integral { [ CompX[{u}     ]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
+            { Name int_y_t~{i}~{j} ; Value { Integral { [ CompY[{u}     ]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
+            { Name int_x_r~{i}~{j} ; Value { Integral { [ CompX[{u}-Ei[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
+            { Name int_y_r~{i}~{j} ; Value { Integral { [ CompY[{u}-Ei[]]*expialphaxy~{i}~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
+          EndIf
+            { Name eff_t~{i}~{j}   ; Value { Term{Type Global; [
+                  1/(gammat~{i}~{j}[]*-k1z[]*Cos[xsi]^2) * ((gammat~{i}~{j}[]^2+alpha~{i}~{j}[]^2)*SquNorm[$int_x_t~{i}~{j}]+
                                                             (gammat~{i}~{j}[]^2+ beta~{j}~{j}[]^2)*SquNorm[$int_y_t~{i}~{j}]+
-                                                             2*alpha~{i}~{j}[]*beta~{i}~{j}[]*Re[$int_x_t~{i}~{j}*Conj[$int_y_t~{i}~{j}]] ) ] ; In SurfIntBot ; } } }
-          { Name eff_r~{i}~{j}   ; Value { Term{Type Global; [
-                 1/(gammar~{i}~{j}[]*-k1z[]*Cos[xsi]^2) * ((gammar~{i}~{j}[]^2+alpha~{i}~{j}[]^2)*SquNorm[$int_x_r~{i}~{j}]+
-                                                           (gammar~{i}~{j}[]^2+ beta~{i}~{j}[]^2)*SquNorm[$int_y_r~{i}~{j}]+
-                                                            2*alpha~{i}~{j}[]*beta~{i}~{j}[]*Re[$int_x_r~{i}~{j}*Conj[$int_y_r~{i}~{j}]]) ] ; In SurfIntTop ; } } }
-          { Name numbering_ij~{i}~{j}   ; Value { Term{Type Global; [Vector[i-Nmax,j-Nmax,0]] ; In SurfIntBot ; } } }
+                                                              2*alpha~{i}~{j}[]*beta~{i}~{j}[]*Re[$int_x_t~{i}~{j}*Conj[$int_y_t~{i}~{j}]] ) ] ; In SurfIntBot ; } } }
+            { Name eff_r~{i}~{j}   ; Value { Term{Type Global; [
+                  1/(gammar~{i}~{j}[]*-k1z[]*Cos[xsi]^2) * ((gammar~{i}~{j}[]^2+alpha~{i}~{j}[]^2)*SquNorm[$int_x_r~{i}~{j}]+
+                                                            (gammar~{i}~{j}[]^2+ beta~{i}~{j}[]^2)*SquNorm[$int_y_r~{i}~{j}]+
+                                                              2*alpha~{i}~{j}[]*beta~{i}~{j}[]*Re[$int_x_r~{i}~{j}*Conj[$int_y_r~{i}~{j}]]) ] ; In SurfIntTop ; } } }
+            { Name numbering_ij~{i}~{j}   ; Value { Term{Type Global; [Vector[i-Nmax,j-Nmax,0]] ; In SurfIntBot ; } } }
         EndFor
       EndFor
     }
@@ -426,34 +444,6 @@ PostOperation {
         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_interpZTot} , 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_interpZTot} , File StrCat[myDir,"Poy_tot_cut_X=0.pos"], Name "Poy_tot_cut_X=0"];
       EndIf
-      // // // // For DEBUG
-      // Print [ epsr_xx , OnElementsOf Omega, File StrCat[myDir,"epsr_xx.pos"]];
-      // Print [ H1x , OnElementsOf Omega, File StrCat[myDir,"H1x.pos"]];
-      // Print [ Etot , OnElementsOf Omega, File StrCat[myDir,"Etot.pos"]];
-      // Print [ u , OnElementsOf Omega, File StrCat[myDir,"u.pos"]];
-      // Print [ E1 , OnElementsOf Omega, File StrCat[myDir,"E1.pos"]];
-      // Print [ Htotx , OnElementsOf Omega, File StrCat[myDir,"Htotx.pos"]];
-      // Print [ uper   , OnElementsOf Omega, File StrCat[myDir,"uper.pos"]];
-      // Print [ uperx  , OnElementsOf Omega, File StrCat[myDir,"uperx.pos"]];
-      // Print [ upery  , OnElementsOf Omega, File StrCat[myDir,"upery.pos"]];
-      // Print [ uperz  , OnElementsOf Omega, File StrCat[myDir,"uperz.pos"]];
-
-      // Print [ E1perx , OnElementsOf Omega, File StrCat[myDir,"E1perx.pos"]];
-      // Print [ E1pery , OnElementsOf Omega, File StrCat[myDir,"E1pery.pos"]];
-      // Print [ Etotpery , OnElementsOf Omega, File StrCat[myDir,"Etotpery.pos"]];
-      // Print [ Poy_tot , OnPlane { {0.5*(-period_x-dys), -dyc/2,(hh_L_6+hh_L_5)/2}
-      //                             {0.5*( period_x-dys), -dyc/2,(hh_L_6+hh_L_5)/2}
-      //                             {0.5*(-period_x+dys),  dyc/2,(hh_L_6+hh_L_5)/2} }
-      //                             {npts_checkpoyX-1,npts_checkpoyY-1} , File StrCat[myDir,"Poy_tot.pos"]];
-      // Print [ Poy_ref , OnPlane { {0.5*(-period_x-dys), -dyc/2, hh_L_1+thick_L_1/2}
-      //                             {0.5*( period_x-dys), -dyc/2, hh_L_1+thick_L_1/2}
-      //                             {0.5*(-period_x+dys),  dyc/2, hh_L_1+thick_L_1/2} }
-      //                             {npts_checkpoyX-1,npts_checkpoyY-1} , File StrCat[myDir,"Poy_ref.pos"]];
-      // Print [ Poy_inc , OnPlane { {0.5*(-period_x-dys), -dyc/2, hh_L_1+thick_L_1/2}
-      //                             {0.5*( period_x-dys), -dyc/2, hh_L_1+thick_L_1/2}
-      //                             {0.5*(-period_x+dys),  dyc/2, hh_L_1+thick_L_1/2} }
-      //                             {npts_checkpoyX-1,npts_checkpoyY-1} , File StrCat[myDir,"Poy_inc.pos"]];
-      //////////// END DEBUG  /////////
 
       Print [ Poy_tot , OnPlane { {0.5*(-period_x-dys), -dyc/2,(hh_L_6+hh_L_5)/2}
                                   {0.5*( period_x-dys), -dyc/2,(hh_L_6+hh_L_5)/2}
diff --git a/DiffractionGratings/grating3D_data_2Dlamellar.geo b/DiffractionGratings/grating3D_data_2Dlamellar.geo
index b6d19e07c254b7e797047d111fd9931f9ca24768..27dc260a162f77b443dd8c429349f7c94488473d 100644
--- a/DiffractionGratings/grating3D_data_2Dlamellar.geo
+++ b/DiffractionGratings/grating3D_data_2Dlamellar.geo
@@ -6,6 +6,7 @@ pp4 = "4Layer Materials";
 pp5 = "5Computational Parameters";
 pp6 = "6Output";
 DefineConstant[
+    FLAG_TOTAL    = {0     , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
     lambda0       = {1000  , Name StrCat[pp1,"/1lambda0 [nm]"]},
     thetadeg      = {30    , Name StrCat[pp1,"/2theta0 [deg]"]},
     phideg        = {90     , Name StrCat[pp1,"/3phi0 [deg]"]},
diff --git a/DiffractionGratings/grating3D_data_bisin.geo b/DiffractionGratings/grating3D_data_bisin.geo
index f46e34327bd0c44ca4178eab847d9cec4318a275..65840b288e8c965f05dce0f04aa60128b6d30285 100644
--- a/DiffractionGratings/grating3D_data_bisin.geo
+++ b/DiffractionGratings/grating3D_data_bisin.geo
@@ -6,6 +6,7 @@ pp4 = "4Layer Materials";
 pp5 = "5Computational Parameters";
 pp6 = "6Output";
 DefineConstant[
+    FLAG_TOTAL    = {0     , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
     lambda0       = {830   , Name StrCat[pp1,"/1lambda0 [nm]"]},
     thetadeg      = {0     , Name StrCat[pp1,"/2theta0 [deg]"]},
     phideg        = {0     , Name StrCat[pp1,"/3phi0 [deg]"]},
diff --git a/DiffractionGratings/grating3D_data_checker.geo b/DiffractionGratings/grating3D_data_checker.geo
index f10ec4d5143cc96cd90fe0869ebbd5b4337312db..169cb1fa72d1f288ab9aa4f83512ab90711131ed 100644
--- a/DiffractionGratings/grating3D_data_checker.geo
+++ b/DiffractionGratings/grating3D_data_checker.geo
@@ -6,6 +6,7 @@ pp4 = "4Layer Materials";
 pp5 = "5Computational Parameters";
 pp6 = "6Output";
 DefineConstant[
+    FLAG_TOTAL    = {0     , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
     lambda0       = {500   , Name StrCat[pp1,"/1lambda0 [nm]"]},
     thetadeg      = {0     , Name StrCat[pp1,"/2theta0 [deg]"]},
     phideg        = {0     , Name StrCat[pp1,"/3phi0 [deg]"]},
diff --git a/DiffractionGratings/grating3D_data_conv.geo b/DiffractionGratings/grating3D_data_conv.geo
index 334fa7ba6bd249a2bcb4771247a7b0f6e95eb2dc..1e227c44c2e34dd160d2d31b9716d5a4d69477f6 100644
--- a/DiffractionGratings/grating3D_data_conv.geo
+++ b/DiffractionGratings/grating3D_data_conv.geo
@@ -6,6 +6,7 @@ pp4 = "4Layer Materials";
 pp5 = "5Computational Paramameters";
 pp6 = "6Output";
 DefineConstant[
+    FLAG_TOTAL    = {0     , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
     lambda0       = {495  , Name StrCat[pp1,"/1lambda0 [nm]"]},
     thetadeg      = {40   , Name StrCat[pp1,"/2theta0 [deg]"]},
     phideg        = {36   , Name StrCat[pp1,"/3phi0 [deg]"]},
diff --git a/DiffractionGratings/grating3D_data_halfellipsoid.geo b/DiffractionGratings/grating3D_data_halfellipsoid.geo
index 0dd7f2ff1cc27ccfd92a774029c4fd2e3c137f13..9075c042b0ed9f507a20ab60cfc3ee88e2e5427e 100644
--- a/DiffractionGratings/grating3D_data_halfellipsoid.geo
+++ b/DiffractionGratings/grating3D_data_halfellipsoid.geo
@@ -6,6 +6,7 @@ pp4 = "4Layer Materials";
 pp5 = "5Computational Parameters";
 pp6 = "6Output";
 DefineConstant[
+    FLAG_TOTAL    = {0     , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
     lambda0       = {495  , Name StrCat[pp1,"/1lambda0 [nm]"]},
     thetadeg      = {40   , Name StrCat[pp1,"/2theta0 [deg]"]},
     phideg        = {36   , Name StrCat[pp1,"/3phi0 [deg]"]},
diff --git a/DiffractionGratings/grating3D_data_hole.geo b/DiffractionGratings/grating3D_data_hole.geo
index 593f69ec14f92aed5ba833aea3db759a2f255f77..d05712415eb148937abb803f6dd63dac47964632 100644
--- a/DiffractionGratings/grating3D_data_hole.geo
+++ b/DiffractionGratings/grating3D_data_hole.geo
@@ -12,6 +12,7 @@ pp4 = "4Layer Materials";
 pp5 = "5Computational Parameters";
 pp6 = "6Output";
 DefineConstant[
+    FLAG_TOTAL    = {0     , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
     lambda0       = {500   , Name StrCat[pp1,"/1lambda0 [nm]"]},
     thetadeg      = {0     , Name StrCat[pp1,"/2theta0 [deg]"]},
     phideg        = {0     , Name StrCat[pp1,"/3phi0 [deg]"]},
diff --git a/DiffractionGratings/grating3D_data_pyramid.geo b/DiffractionGratings/grating3D_data_pyramid.geo
index 1382421140200dd841e4403069d25973dec10451..fb97f520459ccaf20660eddfc4f59d934d42dc56 100644
--- a/DiffractionGratings/grating3D_data_pyramid.geo
+++ b/DiffractionGratings/grating3D_data_pyramid.geo
@@ -6,6 +6,7 @@ pp4 = "4Layer Materials";
 pp5 = "5Computational Parameters";
 pp6 = "6Output";
 DefineConstant[
+    FLAG_TOTAL    = {0     , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
     lambda0       = {1533  , Name StrCat[pp1,"/1lambda0 [nm]"]},
     thetadeg      = {30    , Name StrCat[pp1,"/2theta0 [deg]"]},
     phideg        = {45    , Name StrCat[pp1,"/3phi0 [deg]"]},
diff --git a/DiffractionGratings/grating3D_data_skew.geo b/DiffractionGratings/grating3D_data_skew.geo
index 064c68a21f70499c59f55635fc9da3e332d854bd..7ceee288e05637a3136cc0fb06e57dd90d619d9a 100644
--- a/DiffractionGratings/grating3D_data_skew.geo
+++ b/DiffractionGratings/grating3D_data_skew.geo
@@ -6,6 +6,8 @@ pp4 = "4Layer Materials";
 pp5 = "5Computational Parameters";
 pp6 = "6Output";
 DefineConstant[
+    FLAG_TOTAL    = {0     , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
+    FLAG_TOTAL    = {0     , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
     lambda0       = {495  , Name StrCat[pp1,"/1lambda0 [nm]"]},
     thetadeg      = {40   , Name StrCat[pp1,"/2theta0 [deg]"]},
     phideg        = {36   , Name StrCat[pp1,"/3phi0 [deg]"]},
diff --git a/DiffractionGratings/grating3D_data_solarcell.geo b/DiffractionGratings/grating3D_data_solarcell.geo
index 74a0c161091262cbb3a403e39385f8e94a122b53..7c26bb74ea80a0617bb6e56552e23966e741c1f0 100644
--- a/DiffractionGratings/grating3D_data_solarcell.geo
+++ b/DiffractionGratings/grating3D_data_solarcell.geo
@@ -7,8 +7,8 @@ pp5 = "5Computational Parameters";
 pp6 = "6Output";
 
 DefineConstant[
-    // lambda0       = {500   , Name StrCat[pp1,"/1lambda0 [nm]"]},
-    lambda0 = {500 , Min 500, Max 700, Step 50, Name StrCat[pp1, "0wavelength [nm]"] , Loop 1},
+    FLAG_TOTAL    = {0     , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
+    lambda0       = {500 , Min 500, Max 700, Step 50, Name StrCat[pp1, "0wavelength [nm]"] , Loop 1},
     thetadeg      = {0     , Name StrCat[pp1,"/2theta0 [deg]"]},
     phideg        = {0     , Name StrCat[pp1,"/3phi0 [deg]"]},
     psideg        = {45    , Name StrCat[pp1,"/4psi0 [deg]"]},
diff --git a/DiffractionGratings/grating3D_data_torus.geo b/DiffractionGratings/grating3D_data_torus.geo
index a22f891c05b2749d54f051ed6498b708bbb58510..643534525c88cb4a8fc36887b0c4178c4907e101 100644
--- a/DiffractionGratings/grating3D_data_torus.geo
+++ b/DiffractionGratings/grating3D_data_torus.geo
@@ -6,6 +6,7 @@ pp4 = "4Layer Materials";
 pp5 = "5Computational Parameters";
 pp6 = "6Output";
 DefineConstant[
+    FLAG_TOTAL    = {0     , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
     lambda0       = {1000  , Name StrCat[pp1,"/1lambda0 [nm]"]},
     thetadeg      = {0     , Name StrCat[pp1,"/2theta0 [deg]"]},
     phideg        = {0     , Name StrCat[pp1,"/3phi0 [deg]"]},
diff --git a/DiffractionGratings/grating3D_runall.sh b/DiffractionGratings/grating3D_runall.sh
index 368e0a387f64d2fb326893b40a10ca4f974d01a3..b296018bc832d559220deaf1949515422b001486 100644
--- a/DiffractionGratings/grating3D_runall.sh
+++ b/DiffractionGratings/grating3D_runall.sh
@@ -1,13 +1,22 @@
-GMSHDIR="/home/demesy/programs/gmsh-4.4.1-Linux64/bin"
-for t in bisin checker halfellipsoid hole pyramid torus 2Dlamellar solarcell conv skew
+rm -r res3D*
+
+# for t in bisin checker halfellipsoid hole pyramid torus 2Dlamellar solarcell conv skew
+for t in halfellipsoid torus  
 do
-    $GMSHDIR/gmsh grating3D.pro -setstring test_case $t -
-    mv res3D res3D_$t
+    for f in 0 1
+    do
+        gmsh grating3D.pro -setstring test_case $t -setnumber FLAG_TOTAL $f -
+        mv res3D res3D_$t"_FLAG_TOTAL_"$f
+    done
 done
 
-for t in bisin checker halfellipsoid hole pyramid torus 2Dlamellar solarcell conv skew
+# for t in bisin checker halfellipsoid hole pyramid torus 2Dlamellar solarcell conv skew
+for t in halfellipsoid torus  
 do
-    $GMSHDIR/gmsh grating3D.geo -setstring test_case $t res3D_$t/*.pos &
-    echo "____________\n"$t
-    python grating3D_postplot.py res3D_$t
+    for f in 0 1
+    do
+        # gmsh grating3D.geo -setstring test_case $t res3D_$t"_FLAG_TOTAL_"$f"/*.pos" &
+        echo "____________\ntest case "$t" - FLAG_TOT="$f
+        python grating3D_postplot.py res3D_$t"_FLAG_TOTAL_"$f
+    done
 done
\ No newline at end of file
diff --git a/ElectromagneticScattering/scattering.pro b/ElectromagneticScattering/scattering.pro
index b4dde40106c78d868bf25f980451d950f0d1da2c..789f6bcd253278f7cd422e20c736b11b33b9201a 100644
--- a/ElectromagneticScattering/scattering.pro
+++ b/ElectromagneticScattering/scattering.pro
@@ -369,11 +369,11 @@ Resolution {
       }
       Operation {
         CreateDir[Str[myDir]];
-        Evaluate[Python[]{"scattering_init.py"}];
+        // Evaluate[Python[]{"scattering_init.py"}];
         Generate[P]; 
         Solve[P];
         PostOperation[PW_postop];
-        Evaluate[Python[]{"scattering_post.py"}];
+        // Evaluate[Python[]{"scattering_post.py"}];
       }
     }
   EndIf
@@ -384,7 +384,7 @@ Resolution {
       }
       Operation {
         CreateDir[Str[myDir]];
-        Evaluate[Python[]{"scattering_init.py"}];
+        // Evaluate[Python[]{"scattering_init.py"}];
         
         For pe In {1:p_max}
           Evaluate[ $isN = 0 ];
@@ -403,7 +403,7 @@ Resolution {
           SolveAgain[T];
           PostOperation[VPWN_postop~{pe}];
         EndFor
-        Evaluate[Python[]{"scattering_post.py"}];
+        // Evaluate[Python[]{"scattering_post.py"}];
       }
     }
 
@@ -417,13 +417,13 @@ Resolution {
       }
       Operation {
         CreateDir[Str[myDir]];
-        Evaluate[Python[]{"scattering_init.py"}];
+        // Evaluate[Python[]{"scattering_init.py"}];
         For ncomp In {0:2}
           Generate[M~{ncomp}];
           Solve[M~{ncomp}];
           PostOperation[GreenAll_postop~{ncomp}];
         EndFor
-        Evaluate[Python[]{"scattering_post.py"}];
+        // Evaluate[Python[]{"scattering_post.py"}];
       }
     }    
   EndIf
@@ -434,11 +434,11 @@ Resolution {
       }
       Operation{
         CreateDir[Str[myDir]];
-        Evaluate[Python[]{"scattering_init.py"}];
+        // Evaluate[Python[]{"scattering_init.py"}];
         GenerateSeparate[E];
         EigenSolve[E,neig,eigtarget,0,EigFilter[]];
         PostOperation[QNM_postop];
-        Evaluate[Python[]{"scattering_post.py"}];
+        // Evaluate[Python[]{"scattering_post.py"}];
       }
     }
   EndIf