diff --git a/DiffractionGratings/grating3D.geo b/DiffractionGratings/grating3D.geo
index 249fcebcf9e25eaa5f182ebb512b470b32c2fd47..782394f99fad2e395c7675ef881645f99251869c 100644
--- a/DiffractionGratings/grating3D.geo
+++ b/DiffractionGratings/grating3D.geo
@@ -28,10 +28,10 @@ Macro SetPBCs
   //   Printf("slaveY surf %g",slaveY(k));
   // EndFor
   If (tag_geom==6) // bi-sin : BoundingBox does not catch BSpline Surfaces?
-    masterX()+={48,52};
-    slaveX()+={50,54};
-    masterY()+={51,55};
-    slaveY()+={49,53};
+    masterX()+={17,12};
+    slaveX()+={19,14};
+    masterY()+={15,20};
+    slaveY()+={13,18};
   EndIf
   Periodic Surface{slaveX()} = {masterX()} Translate{period_x,        0,        0};
   Periodic Surface{slaveY()} = {masterY()} Translate{    dys, dyc,        0};
@@ -40,6 +40,7 @@ Return
 If (tag_geom==6)
   // Bi-sinusoidal profile by Spline surface filling
   // not the most general way to implement a freesurface height=f(x,y) and relies on multiple Coherence calls
+  // several hacks are used, one of them is to do all this at the very begining
   nsin = 20;
   For i In {0:nsin}
       xx1~{i} = -period_x/2+i*period_x/2/(nsin-1);
@@ -107,22 +108,26 @@ If (tag_geom==6)
   Curve Loop(30) = {31, -23, 32, 4, 10};
   Surface(18) = {30};
   Surface Loop(3) = {5, 16, 18, 12, 11, 9, 8, 7, 10};
-  Volume(9) = {3};    
+  Volume(1) = {3};    
   Surface Loop(2) = {6, 15, 14, 10, 7, 8, 9, 13, 17};
-  Volume(10) = {2};
+  Volume(2) = {2};
 EndIf
 
 For k In {7:0:-1}
-  If (tag_geom==6 && k!=3)
-    p=newp; l=newl; ll=newll; s=news;
-    Point(p)   = {0.5*(-period_x-dys), -dyc/2, hh_L~{k}};
-    Point(p+1) = {0.5*( period_x-dys), -dyc/2, hh_L~{k}};
-    Point(p+2) = {0.5*( period_x+dys),  dyc/2, hh_L~{k}};
-    Point(p+3) = {0.5*(-period_x+dys),  dyc/2, hh_L~{k}};
-    Line(l)={p,p+1};Line(l+1)={p+1,p+2};Line(l+2)={p+2,p+3};Line(l+3)={p+3,p};
-    Curve Loop(ll) = {l,l+1,l+2,l+3};
-    Surface(s) = {ll};
-    Extrude {0, 0, thick_L~{k}} {Surface{s};}
+  If (tag_geom==6)
+    If (k!=3)
+      p=newp; l=newl; ll=newll; s=news;
+      Point(p)   = {0.5*(-period_x-dys), -dyc/2, hh_L~{k}};
+      Point(p+1) = {0.5*( period_x-dys), -dyc/2, hh_L~{k}};
+      Point(p+2) = {0.5*( period_x+dys),  dyc/2, hh_L~{k}};
+      Point(p+3) = {0.5*(-period_x+dys),  dyc/2, hh_L~{k}};
+      Line(l)={p,p+1};Line(l+1)={p+1,p+2};Line(l+2)={p+2,p+3};Line(l+3)={p+3,p};
+      Curve Loop(ll) = {l,l+1,l+2,l+3};
+      Surface(s) = {ll};
+      Extrude {0, 0, thick_L~{k}} {Surface{s};}
+    Else
+      v=newv;
+    EndIf
   Else
     p=newp; l=newl; ll=newll; s=news;
     Point(p)   = {0.5*(-period_x-dys), -dyc/2, hh_L~{k}};
@@ -192,20 +197,28 @@ EndIf
 
 
 Coherence;
-// // // Physical Surface("ENDPML",770) = {69,37};
-
 Call SetPBCs;
-
-Physical Volume("PMLBOT"        ,1) = {1};
-Physical Volume("LAYER_L6_SUBS" ,2) = {2};
-Physical Volume("LAYER_L5"      ,3) = {3};
-Physical Volume("LAYER_L4"      ,4) = {4};
-Physical Volume("LAYER_L3"      ,5) = {10};
-Physical Volume("LAYER_L2"      ,6) = {6};
-Physical Volume("LAYER_L1_SUPER",7) = {7};
-Physical Volume("PMLTOP"        ,8) = {8};
-Physical Volume("SCAT"          ,9) = {9};
-
+If (tag_geom==6)
+  Physical Volume("PMLBOT"        ,1) = {3};
+  Physical Volume("LAYER_L6_SUBS" ,2) = {4};
+  Physical Volume("LAYER_L5"      ,3) = {5};
+  Physical Volume("LAYER_L4"      ,4) = {6};
+  Physical Volume("LAYER_L3"      ,5) = {2};
+  Physical Volume("LAYER_L2"      ,6) = {7};
+  Physical Volume("LAYER_L1_SUPER",7) = {8};
+  Physical Volume("PMLTOP"        ,8) = {9};
+  Physical Volume("SCAT"          ,9) = {1};
+Else
+  Physical Volume("PMLBOT"        ,1) = {1};
+  Physical Volume("LAYER_L6_SUBS" ,2) = {2};
+  Physical Volume("LAYER_L5"      ,3) = {3};
+  Physical Volume("LAYER_L4"      ,4) = {4};
+  Physical Volume("LAYER_L3"      ,5) = {10};
+  Physical Volume("LAYER_L2"      ,6) = {6};
+  Physical Volume("LAYER_L1_SUPER",7) = {7};
+  Physical Volume("PMLTOP"        ,8) = {8};
+  Physical Volume("SCAT"          ,9) = {9};
+EndIf
 Physical Surface("BXM" ,101 ) = masterX();
 Physical Surface("BXP" ,102 ) = slaveX();
 Physical Surface("BYM" ,201 ) = masterY();
diff --git a/DiffractionGratings/grating3D_data_2Dlamellar.geo b/DiffractionGratings/grating3D_data_2Dlamellar.geo
index 86daafe6256df2c5c96c5a3a1bc1473d05d043ef..bdc6a2d6cc25636904658fddd9047b19e990fb1b 100644
--- a/DiffractionGratings/grating3D_data_2Dlamellar.geo
+++ b/DiffractionGratings/grating3D_data_2Dlamellar.geo
@@ -18,7 +18,7 @@ DefineConstant[
     thick_L_4     = {150     , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
     thick_L_5     = {150     , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
     thick_L_6     = {100    , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
-    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 0},
+    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
     tag_geom      = {  7      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = { 0     , Name StrCat[pp3,"/1rx"]},
diff --git a/DiffractionGratings/grating3D_data_bisin.geo b/DiffractionGratings/grating3D_data_bisin.geo
index d3d18c0673dbadb03e541c8db70e7aebb6f9adab..4b2296a4bc62ba35927bfbdf200ad1f3f7bf786e 100644
--- a/DiffractionGratings/grating3D_data_bisin.geo
+++ b/DiffractionGratings/grating3D_data_bisin.geo
@@ -18,7 +18,7 @@ DefineConstant[
     thick_L_4     = {100   , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
     thick_L_5     = {100   , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
     thick_L_6     = {200   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
-    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 0},
+    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
     tag_geom      = {  6      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = {1.25*lambda0, Name StrCat[pp3,"/1rx"]},
diff --git a/DiffractionGratings/grating3D_data_checker.geo b/DiffractionGratings/grating3D_data_checker.geo
index bfab429e67f0aa6bbbccf0e7a4e2dc4acb751afd..e85f3c2832039e12f64f87b29e5ce449916cca82 100644
--- a/DiffractionGratings/grating3D_data_checker.geo
+++ b/DiffractionGratings/grating3D_data_checker.geo
@@ -18,7 +18,7 @@ DefineConstant[
     thick_L_4     = {50    , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
     thick_L_5     = {50    , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
     thick_L_6     = {50    , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
-    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 0},
+    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
     tag_geom      = {  5      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = {1.25*lambda0, Name StrCat[pp3,"/1rx"]},
diff --git a/DiffractionGratings/grating3D_data_conv.geo b/DiffractionGratings/grating3D_data_conv.geo
index 9dfef1a302e711dc2c2cec75d80a743fccda6eaf..3c7d5019245c0ea228fd228e98ec7a62832c8553 100644
--- a/DiffractionGratings/grating3D_data_conv.geo
+++ b/DiffractionGratings/grating3D_data_conv.geo
@@ -18,7 +18,7 @@ DefineConstant[
     thick_L_4     = {50   , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
     thick_L_5     = {50   , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
     thick_L_6     = {50   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
-    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 0},
+    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
     tag_geom      = {  4      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = {107  , Name StrCat[pp3,"/1rx"]},
diff --git a/DiffractionGratings/grating3D_data_halfellipsoid.geo b/DiffractionGratings/grating3D_data_halfellipsoid.geo
index f7735245b13692a92c596938f1d97b5b19648931..ac0efcd99dc4eb684293c71ffaf14cde4716a6bb 100644
--- a/DiffractionGratings/grating3D_data_halfellipsoid.geo
+++ b/DiffractionGratings/grating3D_data_halfellipsoid.geo
@@ -18,7 +18,7 @@ DefineConstant[
     thick_L_4     = {50   , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
     thick_L_5     = {50   , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
     thick_L_6     = {50   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
-    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 0},
+    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
     tag_geom      = {  4      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = {107  , Name StrCat[pp3,"/1rx"]},
diff --git a/DiffractionGratings/grating3D_data_hole.geo b/DiffractionGratings/grating3D_data_hole.geo
index bdb23d20adcacdf8b8a93d5d1ce2e4dab51c6386..3a8c822d822428e8e64af0c4ffcc2d5989378b9a 100644
--- a/DiffractionGratings/grating3D_data_hole.geo
+++ b/DiffractionGratings/grating3D_data_hole.geo
@@ -18,7 +18,7 @@ DefineConstant[
     thick_L_4     = {50   , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
     thick_L_5     = {100   , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
     thick_L_6     = {100   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
-    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 0},
+    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
     tag_geom      = { 2      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = {250   , Name StrCat[pp3,"/1rx"]},
diff --git a/DiffractionGratings/grating3D_data_pyramid.geo b/DiffractionGratings/grating3D_data_pyramid.geo
index d6af7609de8feeda472cb60a3980eedb2fe54464..e9f0338751d1f6acc9ddbc56daa6cdca7719de9b 100644
--- a/DiffractionGratings/grating3D_data_pyramid.geo
+++ b/DiffractionGratings/grating3D_data_pyramid.geo
@@ -18,7 +18,7 @@ DefineConstant[
     thick_L_4     = {50    , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
     thick_L_5     = {50    , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
     thick_L_6     = {50    , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
-    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 0},
+    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
     tag_geom      = {  1      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = {107   , Name StrCat[pp3,"/1rx"]},
diff --git a/DiffractionGratings/grating3D_data_skew.geo b/DiffractionGratings/grating3D_data_skew.geo
index 44e2aa992c646ceab2d6940402e016ea9d39e6ac..31edf8340f6fe7ed00cd52fda137f789229d3447 100644
--- a/DiffractionGratings/grating3D_data_skew.geo
+++ b/DiffractionGratings/grating3D_data_skew.geo
@@ -18,7 +18,7 @@ DefineConstant[
     thick_L_4     = {50   , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
     thick_L_5     = {50   , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
     thick_L_6     = {50   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
-    xsideg        = {30   , Name StrCat[pp2,"/9skew angle [deg]"],Visible 0},
+    xsideg        = {30   , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
     tag_geom      = {  4      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = {100 , Name StrCat[pp3,"/1rx"]},
diff --git a/DiffractionGratings/grating3D_data_solarcell.geo b/DiffractionGratings/grating3D_data_solarcell.geo
index 173efae138c115daca4e81e26d9d93c22adf8776..3501b60712480d5373babc3f4c54d66016cdc07a 100644
--- a/DiffractionGratings/grating3D_data_solarcell.geo
+++ b/DiffractionGratings/grating3D_data_solarcell.geo
@@ -20,7 +20,7 @@ DefineConstant[
     thick_L_4     = {25   , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
     thick_L_5     = {25   , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
     thick_L_6     = {25   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
-    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 0},
+    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
     tag_geom      = { 2      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = { 50   , Name StrCat[pp3,"/1rx"]},
diff --git a/DiffractionGratings/grating3D_data_torus.geo b/DiffractionGratings/grating3D_data_torus.geo
index 01aab56699944e6ec5cb365bc8472425300f96f4..220bd6477f09798549f8f866314e43f4346b8491 100644
--- a/DiffractionGratings/grating3D_data_torus.geo
+++ b/DiffractionGratings/grating3D_data_torus.geo
@@ -18,7 +18,7 @@ DefineConstant[
     thick_L_4     = {100   , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
     thick_L_5     = {50    , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
     thick_L_6     = {50    , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
-    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 0},
+    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
     tag_geom      = {  3      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = {150/2 , Name StrCat[pp3,"/1rx"]},
diff --git a/DiffractionGratings/grating3D_postplot.py b/DiffractionGratings/grating3D_postplot.py
index 67d1f9d22027c22bd5d4b34c96f89e237f120902..83a8e154623c36ff973681cf9ff48b4bca35f497 100644
--- a/DiffractionGratings/grating3D_postplot.py
+++ b/DiffractionGratings/grating3D_postplot.py
@@ -11,10 +11,9 @@ def dtrap_poy(fname_in,nx,ny):
 
 myDir = sys.argv[1]
 
-intpoyz_tot =  dtrap_poy(myDir+'/Poy_tot_gd.pos',50,50) #*np.cos(float(sys.argv[2])*np.pi/180)
-intpoyz_ref =  dtrap_poy(myDir+'/Poy_ref_gd.pos',50,50) #*np.cos(float(sys.argv[2])*np.pi/180)
-intpoyz_inc =  dtrap_poy(myDir+'/Poy_inc_gd.pos',50,50) #*np.cos(float(sys.argv[2])*np.pi/180)
-Ascat2      = np.loadtxt(myDir+'/temp-Q_scat2.txt')[1]
+intpoyz_tot = abs(dtrap_poy(myDir+'/Poy_tot_gd.pos',50,50))
+intpoyz_ref = abs(dtrap_poy(myDir+'/Poy_ref_gd.pos',50,50))
+intpoyz_inc = abs(dtrap_poy(myDir+'/Poy_inc_gd.pos',50,50))
 
 Rnm = np.loadtxt(myDir+'/eff_r.txt',ndmin=2)[:,1] + 1j*np.loadtxt(myDir+'/eff_r.txt',ndmin=2)[:,2]
 Tnm = np.loadtxt(myDir+'/eff_t.txt',ndmin=2)[:,1] + 1j*np.loadtxt(myDir+'/eff_t.txt',ndmin=2)[:,2]
@@ -24,6 +23,7 @@ Q=np.array(Q)
 TOT = Rnm.real.sum()+Tnm.real.sum()+Q.sum()
 
 if myDir[6:]=='solarcell':
+    print('cf pdf')
     Nmax=2
     tab_lambdas=np.loadtxt(myDir+'/temp_lambda_step.txt',ndmin=2)[:,8]
     nb_lambdas=len(tab_lambdas)
@@ -44,13 +44,21 @@ if myDir[6:]=='solarcell':
     pl.legend()
     pl.xlabel(r'$\lambda$ [nm]')
     pl.ylabel('fraction of incident energy')
-    pl.savefig('solar_balance.pdf')
+    pl.savefig('fig_solar_balance.pdf')
+elif myDir[6:]=='conv':
+    print('cf pdf')
+    pl.figure()
+    pl.plot(np.loadtxt(myDir+'/temp-Q_scat.txt',ndmin=2)[:,1])
+    pl.savefig('fig_convergence_absorption.pdf')
 else:
+    print('===> Computed from diffraction efficiencies')
     print('Rtot ',Rnm.real.sum())
-    print('Rtot2',-intpoyz_ref/intpoyz_inc)
     print('Ttot ',Tnm.real.sum())
+    print('Atot ',Q.sum())
+    print('TOT  ',TOT)
+    print('===> Computed from trapz on Poynting vector')
+    print('(expected to be less precise)')
+    print('Rtot2',intpoyz_ref/intpoyz_inc)
     print('Ttot2',intpoyz_tot/intpoyz_inc)
     print('Atot ',Q.sum())
-    print('Atot2 ',Ascat2/-intpoyz_inc)
-    print('TOT   ',TOT)
-    print('TOT2  ',(-Ascat2+intpoyz_tot-intpoyz_ref)/intpoyz_inc)
+    print('TOT2 ',Q.sum()+(intpoyz_tot+intpoyz_ref)/intpoyz_inc)
diff --git a/DiffractionGratings/grating3D_runall.sh b/DiffractionGratings/grating3D_runall.sh
index 5b2562e00e5eb12cca93ece640c90cbe6c1079bc..368e0a387f64d2fb326893b40a10ca4f974d01a3 100644
--- a/DiffractionGratings/grating3D_runall.sh
+++ b/DiffractionGratings/grating3D_runall.sh
@@ -8,6 +8,6 @@ done
 for t in bisin checker halfellipsoid hole pyramid torus 2Dlamellar solarcell conv skew
 do
     $GMSHDIR/gmsh grating3D.geo -setstring test_case $t res3D_$t/*.pos &
-    echo "\n"$t
+    echo "____________\n"$t
     python grating3D_postplot.py res3D_$t
 done
\ No newline at end of file