diff --git a/DiffractionGratings/grating3D_data_bisin.geo b/DiffractionGratings/grating3D_data_bisin.geo
new file mode 100644
index 0000000000000000000000000000000000000000..de9787fec277e70dfa0bb738e320c6f8fcd6ea1d
--- /dev/null
+++ b/DiffractionGratings/grating3D_data_bisin.geo
@@ -0,0 +1,90 @@
+nm = 1;
+pp1 = "1Incident Plane Wave";
+pp2 = "2Layers Thicknesses";
+pp3 = "3Scatterer Properties";
+pp4 = "4Layer Materials";
+pp5 = "5Computational Paramameters";
+pp6 = "6Output";
+DefineConstant[
+ lambda0 = {830 , Name StrCat[pp1,"/1lambda0 [nm]"]},
+ thetadeg = {0 , Name StrCat[pp1,"/2theta0 [deg]"]},
+ phideg = {0 , Name StrCat[pp1,"/3phi0 [deg]"]},
+ psideg = {0 , Name StrCat[pp1,"/4psi0 [deg]"]},
+ period_x = {1000 , Name StrCat[pp2,"/1X period [nm]"]},
+ period_y = {1000 , Name StrCat[pp2,"/2Y period [nm]"]},
+ thick_L_1 = {100 , Name StrCat[pp2,"/3thickness layer 1 [nm] (superstrate)"]},
+ thick_L_2 = {100 , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
+ thick_L_3 = {500, Name StrCat[pp2,"/5thickness layer 3 [nm]"]},
+ 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)"]},
+
+ tag_geom = { 6 , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal"}},
+ rx = {1.25*lambda0, Name StrCat[pp3,"/1rx"]},
+ ry = {1.25*lambda0, Name StrCat[pp3,"/2ry"]},
+ rz = { 200 , Name StrCat[pp3,"/3rz"]},
+ flag_mat_scat = { 0 , Name StrCat[pp3,"/4Scatterer permittivity model"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ eps_re_Scat = { 4 , Name StrCat[pp3,"/7eps_re_Scat"]},
+ eps_im_Scat = { 0 , Name StrCat[pp3,"/8eps_im_Scat"]},
+
+ flag_mat_1 = { 0 , Name StrCat[pp4,"/1Layer 1"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_2 = { 0 , Name StrCat[pp4,"/2Layer 2"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_3 = { 0 , Name StrCat[pp4,"/3Layer 3"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_4 = { 0 , Name StrCat[pp4,"/4Layer 4"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_5 = { 0 , Name StrCat[pp4,"/5Layer 5"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_6 = { 0 , Name StrCat[pp4,"/6Layer 6"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ eps_re_L_1 = {1 , Name StrCat[pp4,"/layer 1: real part of relative permittivity"]},
+ eps_im_L_1 = {0 , Name StrCat[pp4,"/layer 1: imag part of relative permittivity"]},
+ eps_re_L_2 = {1 , Name StrCat[pp4,"/layer 2: real part of relative permittivity"]},
+ eps_im_L_2 = {0 , Name StrCat[pp4,"/layer 2: imag part of relative permittivity"]},
+ eps_re_L_3 = {1 , Name StrCat[pp4,"/layer 3: real part of relative permittivity"]},
+ eps_im_L_3 = {0 , Name StrCat[pp4,"/layer 3: imag part of relative permittivity"]},
+ eps_re_L_4 = {4 , Name StrCat[pp4,"/layer 4: real part of relative permittivity"]},
+ eps_im_L_4 = {0 , Name StrCat[pp4,"/layer 4: imag part of relative permittivity"]},
+ eps_re_L_5 = {4 , Name StrCat[pp4,"/layer 5: real part of relative permittivity"]},
+ eps_im_L_5 = {0 , Name StrCat[pp4,"/layer 5: imag part of relative permittivity"]},
+ eps_re_L_6 = {4 , Name StrCat[pp4,"/layer 6: real part of relative permittivity"]},
+ eps_im_L_6 = {0 , Name StrCat[pp4,"/layer 6: imag part of relative permittivity"]},
+
+ paramaille = {7 , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
+ scat_lc = {lambda0/(2*paramaille) , Name StrCat[pp5,"/2Scatterer absolute mesh size [nm]"]},
+ PML_top = {lambda0, Name StrCat[pp5,"/4PML top thickness [nm]"]},
+ PML_bot = {lambda0*1.2, Name StrCat[pp5,"/5PML bot thickness [nm]"]},
+ Nmax = {2 , Name StrCat[pp5,"/6Number of non specular order to output [-]"]},
+ refine_mesh_L1= {1 , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]},
+ refine_mesh_L2= {1 , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
+ refine_mesh_L3= {1 , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]},
+ refine_mesh_L4= {1 , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
+ refine_mesh_L5= {1 , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
+ refine_mesh_L6= {1 , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+
+ InterpSampling = { 30 , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
+ Flag_interp_cubic = { 1 , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
+ FlagOutEtotCuts = { 1 , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
+ FlagOutEscaCuts = { 1 , Name StrCat[pp6,"/3Output Scattered Electric Field cuts?"] , Choices {0,1} },
+ FlagOutPoyCut = { 1 , Name StrCat[pp6,"/4Output Poynting cuts?"] , Choices {0,1} },
+ FlagOutEtotFull = { 0 , Name StrCat[pp6,"/5Total Electric Field Full Output?"] , Choices {0,1} },
+ FlagOutEscaFull = { 0 , Name StrCat[pp6,"/6Scattered Electric Field Full Output?"] , Choices {0,1} },
+ FlagOutPoyFull = { 0 , Name StrCat[pp6,"/7Poynting Full Output?"] , Choices {0,1} }
+];
+
+lambda_m = lambda0;
+
+hh_L_6 = -thick_L_6;
+For k In {1:5}
+ hh_L~{6-k} = hh_L~{7-k}+thick_L~{7-k};
+EndFor
+
+PML_bot_hh = hh_L_6-PML_bot;
+PML_top_hh = hh_L_1+thick_L_1;
+
+theta0 = thetadeg*Pi/180;
+phi0 = phideg*Pi/180;
+psi0 = psideg*Pi/180;
+
+DomainZsizeSca = PML_top_hh+PML_bot-(hh_L_6-PML_bot);
+DomainZsizeTot = PML_top_hh-hh_L_6;
+npts_interpX = period_x/InterpSampling;
+npts_interpY = period_y/InterpSampling;
+npts_interpZSca = DomainZsizeSca/InterpSampling;
+npts_interpZTot = DomainZsizeTot/InterpSampling;
diff --git a/DiffractionGratings/grating3D_data_checker.geo b/DiffractionGratings/grating3D_data_checker.geo
new file mode 100644
index 0000000000000000000000000000000000000000..8dddcb078d879e1542f9dc70d43890f0a0450b3a
--- /dev/null
+++ b/DiffractionGratings/grating3D_data_checker.geo
@@ -0,0 +1,90 @@
+nm = 1;
+pp1 = "1Incident Plane Wave";
+pp2 = "2Layers Thicknesses";
+pp3 = "3Scatterer Properties";
+pp4 = "4Layer Materials";
+pp5 = "5Computational Paramameters";
+pp6 = "6Output";
+DefineConstant[
+ lambda0 = {500 , Name StrCat[pp1,"/1lambda0 [nm]"]},
+ thetadeg = {0 , Name StrCat[pp1,"/2theta0 [deg]"]},
+ phideg = {0 , Name StrCat[pp1,"/3phi0 [deg]"]},
+ psideg = {45 , Name StrCat[pp1,"/4psi0 [deg]"]},
+ period_x = {1.25*Sqrt[2]*lambda0, Name StrCat[pp2,"/1X period [nm]"]},
+ period_y = {1.25*Sqrt[2]*lambda0, Name StrCat[pp2,"/2Y period [nm]"]},
+ thick_L_1 = {50 , Name StrCat[pp2,"/3thickness layer 1 [nm] (superstrate)"]},
+ thick_L_2 = {50 , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
+ thick_L_3 = {1.2*lambda0, Name StrCat[pp2,"/5thickness layer 3 [nm]"]},
+ 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)"]},
+
+ tag_geom = { 5 , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal"}},
+ rx = {1.25*lambda0, Name StrCat[pp3,"/1rx"]},
+ ry = {1.25*lambda0, Name StrCat[pp3,"/2ry"]},
+ rz = {lambda0 , Name StrCat[pp3,"/3rz"]},
+ flag_mat_scat = { 0 , Name StrCat[pp3,"/4Scatterer permittivity model"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ eps_re_Scat = { 2.25 , Name StrCat[pp3,"/7eps_re_Scat"]},
+ eps_im_Scat = { 0 , Name StrCat[pp3,"/8eps_im_Scat"]},
+
+ flag_mat_1 = { 0 , Name StrCat[pp4,"/1Layer 1"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_2 = { 0 , Name StrCat[pp4,"/2Layer 2"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_3 = { 0 , Name StrCat[pp4,"/3Layer 3"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_4 = { 0 , Name StrCat[pp4,"/4Layer 4"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_5 = { 0 , Name StrCat[pp4,"/5Layer 5"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_6 = { 0 , Name StrCat[pp4,"/6Layer 6"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ eps_re_L_1 = {2.25 , Name StrCat[pp4,"/layer 1: real part of relative permittivity"]},
+ eps_im_L_1 = {0 , Name StrCat[pp4,"/layer 1: imag part of relative permittivity"]},
+ eps_re_L_2 = {2.25 , Name StrCat[pp4,"/layer 2: real part of relative permittivity"]},
+ eps_im_L_2 = {0 , Name StrCat[pp4,"/layer 2: imag part of relative permittivity"]},
+ eps_re_L_3 = {1 , Name StrCat[pp4,"/layer 3: real part of relative permittivity"]},
+ eps_im_L_3 = {0 , Name StrCat[pp4,"/layer 3: imag part of relative permittivity"]},
+ eps_re_L_4 = {1 , Name StrCat[pp4,"/layer 4: real part of relative permittivity"]},
+ eps_im_L_4 = {0 , Name StrCat[pp4,"/layer 4: imag part of relative permittivity"]},
+ eps_re_L_5 = {1 , Name StrCat[pp4,"/layer 5: real part of relative permittivity"]},
+ eps_im_L_5 = {0 , Name StrCat[pp4,"/layer 5: imag part of relative permittivity"]},
+ eps_re_L_6 = {1 , Name StrCat[pp4,"/layer 6: real part of relative permittivity"]},
+ eps_im_L_6 = {0 , Name StrCat[pp4,"/layer 6: imag part of relative permittivity"]},
+
+ paramaille = {6 , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
+ scat_lc = {lambda0/(1.5*paramaille) , Name StrCat[pp5,"/2Scatterer absolute mesh size [nm]"]},
+ PML_top = {lambda0, Name StrCat[pp5,"/4PML top thickness [nm]"]},
+ PML_bot = {lambda0, Name StrCat[pp5,"/5PML bot thickness [nm]"]},
+ Nmax = {2 , Name StrCat[pp5,"/6Number of non specular order to output [-]"]},
+ refine_mesh_L1= {1.5 , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]},
+ refine_mesh_L2= {1.5 , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
+ refine_mesh_L3= {1 , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]},
+ refine_mesh_L4= {1 , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
+ refine_mesh_L5= {1 , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
+ refine_mesh_L6= {1 , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+
+ InterpSampling = { 30 , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
+ Flag_interp_cubic = { 1 , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
+ FlagOutEtotCuts = { 1 , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
+ FlagOutEscaCuts = { 1 , Name StrCat[pp6,"/3Output Scattered Electric Field cuts?"] , Choices {0,1} },
+ FlagOutPoyCut = { 1 , Name StrCat[pp6,"/4Output Poynting cuts?"] , Choices {0,1} },
+ FlagOutEtotFull = { 0 , Name StrCat[pp6,"/5Total Electric Field Full Output?"] , Choices {0,1} },
+ FlagOutEscaFull = { 0 , Name StrCat[pp6,"/6Scattered Electric Field Full Output?"] , Choices {0,1} },
+ FlagOutPoyFull = { 0 , Name StrCat[pp6,"/7Poynting Full Output?"] , Choices {0,1} }
+];
+
+lambda_m = lambda0;
+
+hh_L_6 = -thick_L_6;
+For k In {1:5}
+ hh_L~{6-k} = hh_L~{7-k}+thick_L~{7-k};
+EndFor
+
+PML_bot_hh = hh_L_6-PML_bot;
+PML_top_hh = hh_L_1+thick_L_1;
+
+theta0 = thetadeg*Pi/180;
+phi0 = phideg*Pi/180;
+psi0 = psideg*Pi/180;
+
+DomainZsizeSca = PML_top_hh+PML_bot-(hh_L_6-PML_bot);
+DomainZsizeTot = PML_top_hh-hh_L_6;
+npts_interpX = period_x/InterpSampling;
+npts_interpY = period_y/InterpSampling;
+npts_interpZSca = DomainZsizeSca/InterpSampling;
+npts_interpZTot = DomainZsizeTot/InterpSampling;
diff --git a/DiffractionGratings/grating3D_data_halfellipsoid.geo b/DiffractionGratings/grating3D_data_halfellipsoid.geo
new file mode 100644
index 0000000000000000000000000000000000000000..bc4951062fe672ba3bf63542c9fefed6724ae73b
--- /dev/null
+++ b/DiffractionGratings/grating3D_data_halfellipsoid.geo
@@ -0,0 +1,89 @@
+nm = 1;
+pp1 = "1Incident Plane Wave";
+pp2 = "2Layers Thicknesses";
+pp3 = "3Scatterer Properties";
+pp4 = "4Layer Materials";
+pp5 = "5Computational Paramameters";
+pp6 = "6Output";
+DefineConstant[
+ lambda0 = {495 , Name StrCat[pp1,"/1lambda0 [nm]"]},
+ thetadeg = {40 , Name StrCat[pp1,"/2theta0 [deg]"]},
+ phideg = {36 , Name StrCat[pp1,"/3phi0 [deg]"]},
+ psideg = {72 , Name StrCat[pp1,"/4psi0 [deg]"]},
+ period_x = {250 , Name StrCat[pp2,"/1X period [nm]"]},
+ period_y = {250 , Name StrCat[pp2,"/2Y period [nm]"]},
+ thick_L_1 = {50 , Name StrCat[pp2,"/3thickness layer 1 [nm] (superstrate)"]},
+ thick_L_2 = {50 , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
+ thick_L_3 = {100 , Name StrCat[pp2,"/5thickness layer 3 [nm]"]},
+ 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)"]},
+
+ tag_geom = { 4 , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal"}},
+ rx = {107 , Name StrCat[pp3,"/1rx"]},
+ ry = {47 , Name StrCat[pp3,"/2ry"]},
+ rz = {40 , Name StrCat[pp3,"/3rz"]},
+ flag_mat_scat = { 0 , Name StrCat[pp3,"/4Scatterer permittivity model"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ eps_re_Scat = {-2.23, Name StrCat[pp3,"/7eps_re_Scat"]},
+ eps_im_Scat = { 3.89, Name StrCat[pp3,"/8eps_im_Scat"]},
+
+ flag_mat_1 = { 0 , Name StrCat[pp4,"/1Layer 1"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_2 = { 0 , Name StrCat[pp4,"/2Layer 2"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_3 = { 0 , Name StrCat[pp4,"/3Layer 3"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_4 = { 0 , Name StrCat[pp4,"/4Layer 4"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_5 = { 0 , Name StrCat[pp4,"/5Layer 5"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_6 = { 0 , Name StrCat[pp4,"/6Layer 6"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ eps_re_L_1 = {1 , Name StrCat[pp4,"/layer 1: real part of relative permittivity"]},
+ eps_im_L_1 = {0 , Name StrCat[pp4,"/layer 1: imag part of relative permittivity"]},
+ eps_re_L_2 = {1 , Name StrCat[pp4,"/layer 2: real part of relative permittivity"]},
+ eps_im_L_2 = {0 , Name StrCat[pp4,"/layer 2: imag part of relative permittivity"]},
+ eps_re_L_3 = {1 , Name StrCat[pp4,"/layer 3: real part of relative permittivity"]},
+ eps_im_L_3 = {0 , Name StrCat[pp4,"/layer 3: imag part of relative permittivity"]},
+ eps_re_L_4 = {1 , Name StrCat[pp4,"/layer 4: real part of relative permittivity"]},
+ eps_im_L_4 = {0 , Name StrCat[pp4,"/layer 4: imag part of relative permittivity"]},
+ eps_re_L_5 = {1 , Name StrCat[pp4,"/layer 5: real part of relative permittivity"]},
+ eps_im_L_5 = {0 , Name StrCat[pp4,"/layer 5: imag part of relative permittivity"]},
+ eps_re_L_6 = {4 , Name StrCat[pp4,"/layer 6: real part of relative permittivity"]},
+ eps_im_L_6 = {0 , Name StrCat[pp4,"/layer 6: imag part of relative permittivity"]},
+
+ paramaille = {8 , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
+ scat_lc = {10 , Name StrCat[pp5,"/2metal mesh size [nm]"]},
+ PML_top = {lambda0 , Name StrCat[pp5,"/4PML top thickness [nm]"]},
+ PML_bot = {lambda0 , Name StrCat[pp5,"/5PML bot thickness [nm]"]},
+ Nmax = {1 , Name StrCat[pp5,"/6Number of non specular order to output [-]"]},
+ refine_mesh_L1= {1 , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]},
+ refine_mesh_L2= {1 , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
+ refine_mesh_L3= {1 , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]},
+ refine_mesh_L4= {1 , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
+ refine_mesh_L5= {1 , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
+ refine_mesh_L6= {1 , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+
+ InterpSampling = { 10 , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
+ Flag_interp_cubic = { 1 , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
+ FlagOutEtotCuts = { 1 , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
+ FlagOutEscaCuts = { 1 , Name StrCat[pp6,"/3Output Scattered Electric Field cuts?"] , Choices {0,1} },
+ FlagOutPoyCut = { 1 , Name StrCat[pp6,"/4Output Poynting cuts?"] , Choices {0,1} },
+ FlagOutEtotFull = { 0 , Name StrCat[pp6,"/5Total Electric Field Full Output?"] , Choices {0,1} },
+ FlagOutEscaFull = { 0 , Name StrCat[pp6,"/6Scattered Electric Field Full Output?"] , Choices {0,1} },
+ FlagOutPoyFull = { 0 , Name StrCat[pp6,"/7Poynting Full Output?"] , Choices {0,1} }
+];
+
+lambda_m = lambda0;
+
+hh_L_6 = -thick_L_6;
+For k In {1:5}
+ hh_L~{6-k} = hh_L~{7-k}+thick_L~{7-k};
+EndFor
+PML_bot_hh = hh_L_6-PML_bot;
+PML_top_hh = hh_L_1+thick_L_1;
+
+theta0 = thetadeg*Pi/180;
+phi0 = phideg*Pi/180;
+psi0 = psideg*Pi/180;
+
+DomainZsizeSca = PML_top_hh+PML_bot-(hh_L_6-PML_bot);
+DomainZsizeTot = PML_top_hh-hh_L_6;
+npts_interpX = period_x/InterpSampling;
+npts_interpY = period_y/InterpSampling;
+npts_interpZSca = DomainZsizeSca/InterpSampling;
+npts_interpZTot = DomainZsizeTot/InterpSampling;
diff --git a/DiffractionGratings/grating3D_data_hole.geo b/DiffractionGratings/grating3D_data_hole.geo
new file mode 100644
index 0000000000000000000000000000000000000000..1b91f97d9cb5afa82136c7032ef52a755161b09a
--- /dev/null
+++ b/DiffractionGratings/grating3D_data_hole.geo
@@ -0,0 +1,90 @@
+nm = 1;
+pp1 = "1Incident Plane Wave";
+pp2 = "2Layers Thicknesses";
+pp3 = "3Scatterer Properties";
+pp4 = "4Layer Materials";
+pp5 = "5Computational Paramameters";
+pp6 = "6Output";
+DefineConstant[
+ lambda0 = {500 , Name StrCat[pp1,"/1lambda0 [nm]"]},
+ thetadeg = {0 , Name StrCat[pp1,"/2theta0 [deg]"]},
+ phideg = {0 , Name StrCat[pp1,"/3phi0 [deg]"]},
+ psideg = {45 , Name StrCat[pp1,"/4psi0 [deg]"]},
+ period_x = {999.9 , Name StrCat[pp2,"/1X period [nm]"]},
+ period_y = {999.9 , Name StrCat[pp2,"/2Y period [nm]"]},
+ thick_L_1 = {100 , Name StrCat[pp2,"/3thickness layer 1 [nm] (superstrate)"]},
+ thick_L_2 = {100 , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
+ thick_L_3 = {50 , Name StrCat[pp2,"/5thickness layer 3 [nm]"]},
+ 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 = {100 , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
+
+ tag_geom = { 2 , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal"}},
+ rx = {250 , Name StrCat[pp3,"/1rx"]},
+ ry = {47 , Name StrCat[pp3,"/2ry"]},
+ rz = {500 , Name StrCat[pp3,"/3rz"]},
+ flag_mat_scat = { 0 , Name StrCat[pp3,"/4Scatterer permittivity model"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ eps_re_Scat = {1 , Name StrCat[pp3,"/7eps_re_Scat"]},
+ eps_im_Scat = {0 , Name StrCat[pp3,"/8eps_im_Scat"]},
+
+ flag_mat_1 = { 0 , Name StrCat[pp4,"/1Layer 1"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_2 = { 0 , Name StrCat[pp4,"/2Layer 2"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_3 = { 0 , Name StrCat[pp4,"/3Layer 3"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_4 = { 0 , Name StrCat[pp4,"/4Layer 4"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_5 = { 0 , Name StrCat[pp4,"/5Layer 5"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_6 = { 0 , Name StrCat[pp4,"/6Layer 6"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ eps_re_L_1 = {1 , Name StrCat[pp4,"/layer 1: real part of relative permittivity"]},
+ eps_im_L_1 = {0 , Name StrCat[pp4,"/layer 1: imag part of relative permittivity"]},
+ eps_re_L_2 = {1 , Name StrCat[pp4,"/layer 2: real part of relative permittivity"]},
+ eps_im_L_2 = {0 , Name StrCat[pp4,"/layer 2: imag part of relative permittivity"]},
+ eps_re_L_3 = {0.8125, Name StrCat[pp4,"/layer 3: real part of relative permittivity"]},
+ eps_im_L_3 = {5.25 , Name StrCat[pp4,"/layer 3: imag part of relative permittivity"]},
+ eps_re_L_4 = {2.25 , Name StrCat[pp4,"/layer 4: real part of relative permittivity"]},
+ eps_im_L_4 = {0 , Name StrCat[pp4,"/layer 4: imag part of relative permittivity"]},
+ eps_re_L_5 = {2.25 , Name StrCat[pp4,"/layer 5: real part of relative permittivity"]},
+ eps_im_L_5 = {0 , Name StrCat[pp4,"/layer 5: imag part of relative permittivity"]},
+ eps_re_L_6 = {2.25 , Name StrCat[pp4,"/layer 6: real part of relative permittivity"]},
+ eps_im_L_6 = {0 , Name StrCat[pp4,"/layer 6: imag part of relative permittivity"]},
+
+ paramaille = {5 , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
+ scat_lc = {lambda0/(3*paramaille) , Name StrCat[pp5,"/2Scatterer absolute mesh size [nm]"]},
+ PML_top = {lambda0, Name StrCat[pp5,"/4PML top thickness [nm]"]},
+ PML_bot = {lambda0, Name StrCat[pp5,"/5PML bot thickness [nm]"]},
+ Nmax = {2 , Name StrCat[pp5,"/6Number of non specular order to output [-]"]},
+ refine_mesh_L1= {1 , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]},
+ refine_mesh_L2= {1 , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
+ refine_mesh_L3= {5 , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]},
+ refine_mesh_L4= {1 , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
+ refine_mesh_L5= {1 , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
+ refine_mesh_L6= {1 , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+
+ InterpSampling = { 30 , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
+ Flag_interp_cubic = { 0 , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
+ FlagOutEtotCuts = { 1 , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
+ FlagOutEscaCuts = { 1 , Name StrCat[pp6,"/3Output Scattered Electric Field cuts?"] , Choices {0,1} },
+ FlagOutPoyCut = { 1 , Name StrCat[pp6,"/4Output Poynting cuts?"] , Choices {0,1} },
+ FlagOutEtotFull = { 0 , Name StrCat[pp6,"/5Total Electric Field Full Output?"] , Choices {0,1} },
+ FlagOutEscaFull = { 0 , Name StrCat[pp6,"/6Scattered Electric Field Full Output?"] , Choices {0,1} },
+ FlagOutPoyFull = { 0 , Name StrCat[pp6,"/7Poynting Full Output?"] , Choices {0,1} }
+];
+
+lambda_m = lambda0;
+
+hh_L_6 = -thick_L_6;
+For k In {1:5}
+ hh_L~{6-k} = hh_L~{7-k}+thick_L~{7-k};
+EndFor
+
+PML_bot_hh = hh_L_6-PML_bot;
+PML_top_hh = hh_L_1+thick_L_1;
+
+theta0 = thetadeg*Pi/180;
+phi0 = phideg*Pi/180;
+psi0 = psideg*Pi/180;
+
+DomainZsizeSca = PML_top_hh+PML_bot-(hh_L_6-PML_bot);
+DomainZsizeTot = PML_top_hh-hh_L_6;
+npts_interpX = period_x/InterpSampling;
+npts_interpY = period_y/InterpSampling;
+npts_interpZSca = DomainZsizeSca/InterpSampling;
+npts_interpZTot = DomainZsizeTot/InterpSampling;
diff --git a/DiffractionGratings/grating3D_data_pyramid.geo b/DiffractionGratings/grating3D_data_pyramid.geo
new file mode 100644
index 0000000000000000000000000000000000000000..05c3754acf72a5ff5ea62a6fb98c5d157e20099d
--- /dev/null
+++ b/DiffractionGratings/grating3D_data_pyramid.geo
@@ -0,0 +1,89 @@
+nm = 1;
+pp1 = "1Incident Plane Wave";
+pp2 = "2Layers Thicknesses";
+pp3 = "3Scatterer Properties";
+pp4 = "4Layer Materials";
+pp5 = "5Computational Paramameters";
+pp6 = "6Output";
+DefineConstant[
+ lambda0 = {1533 , Name StrCat[pp1,"/1lambda0 [nm]"]},
+ thetadeg = {30 , Name StrCat[pp1,"/2theta0 [deg]"]},
+ phideg = {45 , Name StrCat[pp1,"/3phi0 [deg]"]},
+ psideg = { 0 , Name StrCat[pp1,"/4psi0 [deg]"]},
+ period_x = {1500 , Name StrCat[pp2,"/1X period [nm]"]},
+ period_y = {1000 , Name StrCat[pp2,"/2Y period [nm]"]},
+ thick_L_1 = {50 , Name StrCat[pp2,"/3thickness layer 1 [nm] (superstrate)"]},
+ thick_L_2 = {50 , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
+ thick_L_3 = {300 , Name StrCat[pp2,"/5thickness layer 3 [nm]"]},
+ 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)"]},
+
+ tag_geom = { 1 , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal"}},
+ rx = {107 , Name StrCat[pp3,"/1rx"]},
+ ry = {250 , Name StrCat[pp3,"/2ry"]},
+ rz = {250 , Name StrCat[pp3,"/3rz"]},
+ flag_mat_scat = { 0 , Name StrCat[pp3,"/4Scatterer permittivity model"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ eps_re_Scat = {2.25 , Name StrCat[pp3,"/5Custom real part of relative permittivity"]},
+ eps_im_Scat = { 0 , Name StrCat[pp3,"/6Custom real part of relative permittivity"]},
+
+ flag_mat_1 = { 0 , Name StrCat[pp4,"/1Layer 1"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_2 = { 0 , Name StrCat[pp4,"/2Layer 2"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_3 = { 0 , Name StrCat[pp4,"/3Layer 3"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_4 = { 0 , Name StrCat[pp4,"/4Layer 4"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_5 = { 0 , Name StrCat[pp4,"/5Layer 5"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ flag_mat_6 = { 0 , Name StrCat[pp4,"/6Layer 6"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+ eps_re_L_1 = {1 , Name StrCat[pp4,"/layer 1: real part of relative permittivity"]},
+ eps_im_L_1 = {0 , Name StrCat[pp4,"/layer 1: imag part of relative permittivity"]},
+ eps_re_L_2 = {1 , Name StrCat[pp4,"/layer 2: real part of relative permittivity"]},
+ eps_im_L_2 = {0 , Name StrCat[pp4,"/layer 2: imag part of relative permittivity"]},
+ eps_re_L_3 = {1 , Name StrCat[pp4,"/layer 3: real part of relative permittivity"]},
+ eps_im_L_3 = {0 , Name StrCat[pp4,"/layer 3: imag part of relative permittivity"]},
+ eps_re_L_4 = {2.25 , Name StrCat[pp4,"/layer 4: real part of relative permittivity"]},
+ eps_im_L_4 = {0 , Name StrCat[pp4,"/layer 4: imag part of relative permittivity"]},
+ eps_re_L_5 = {2.25 , Name StrCat[pp4,"/layer 5: real part of relative permittivity"]},
+ eps_im_L_5 = {0 , Name StrCat[pp4,"/layer 5: imag part of relative permittivity"]},
+ eps_re_L_6 = {2.25 , Name StrCat[pp4,"/layer 6: real part of relative permittivity"]},
+ eps_im_L_6 = {0 , Name StrCat[pp4,"/layer 6: imag part of relative permittivity"]},
+
+ paramaille = {8 , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
+ scat_lc = {lambda0/(1.5*paramaille) , Name StrCat[pp5,"/2Scatterer absolute mesh size [nm]"]},
+ PML_top = {lambda0 , Name StrCat[pp5,"/4PML top thickness [nm]"]},
+ PML_bot = {lambda0 , Name StrCat[pp5,"/5PML bot thickness [nm]"]},
+ Nmax = {1 , Name StrCat[pp5,"/6Number of non specular order to output [-]"]},
+ refine_mesh_L1= {1 , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]},
+ refine_mesh_L2= {1 , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
+ refine_mesh_L3= {1 , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]},
+ refine_mesh_L4= {1 , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
+ refine_mesh_L5= {1 , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
+ refine_mesh_L6= {1 , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+
+ InterpSampling = { 50 , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
+ Flag_interp_cubic = { 1 , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
+ FlagOutEtotCuts = { 1 , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
+ FlagOutEscaCuts = { 1 , Name StrCat[pp6,"/3Output Scattered Electric Field cuts?"] , Choices {0,1} },
+ FlagOutPoyCut = { 1 , Name StrCat[pp6,"/4Output Poynting cuts?"] , Choices {0,1} },
+ FlagOutEtotFull = { 0 , Name StrCat[pp6,"/5Total Electric Field Full Output?"] , Choices {0,1} },
+ FlagOutEscaFull = { 0 , Name StrCat[pp6,"/6Scattered Electric Field Full Output?"] , Choices {0,1} },
+ FlagOutPoyFull = { 0 , Name StrCat[pp6,"/7Poynting Full Output?"] , Choices {0,1} }
+];
+
+lambda_m = lambda0;
+
+hh_L_6 = -thick_L_6;
+For k In {1:5}
+ hh_L~{6-k} = hh_L~{7-k}+thick_L~{7-k};
+EndFor
+PML_bot_hh = hh_L_6-PML_bot;
+PML_top_hh = hh_L_1+thick_L_1;
+
+theta0 = thetadeg*Pi/180;
+phi0 = phideg*Pi/180;
+psi0 = psideg*Pi/180;
+
+DomainZsizeSca = PML_top_hh+PML_bot-(hh_L_6-PML_bot);
+DomainZsizeTot = PML_top_hh-hh_L_6;
+npts_interpX = period_x/InterpSampling;
+npts_interpY = period_y/InterpSampling;
+npts_interpZSca = DomainZsizeSca/InterpSampling;
+npts_interpZTot = DomainZsizeTot/InterpSampling;