diff --git a/DiffractionGratings/grating3D_data_2Dlamellar.geo b/DiffractionGratings/grating3D_data_2Dlamellar.geo
new file mode 100644
index 0000000000000000000000000000000000000000..84257371523303dd55d32efa8ea042ebb48268b2
--- /dev/null
+++ b/DiffractionGratings/grating3D_data_2Dlamellar.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 = {1000 , Name StrCat[pp1,"/1lambda0 [nm]"]},
+ thetadeg = {30 , 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 = {20 , Name StrCat[pp2,"/2Y period [nm]"]},
+ thick_L_1 = {200 , Name StrCat[pp2,"/3thickness layer 1 [nm] (superstrate)"]},
+ thick_L_2 = {200 , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
+ thick_L_3 = {1050 , Name StrCat[pp2,"/5thickness layer 3 [nm]"]},
+ thick_L_4 = {250 , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
+ thick_L_5 = {250 , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
+ thick_L_6 = {500 , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
+
+ 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 = {500 , Name StrCat[pp3,"/1rx"]},
+ ry = {0 , Name StrCat[pp3,"/2ry"]},
+ rz = {1000 , 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 = {-44.9757 , Name StrCat[pp3,"/7eps_re_Scat"]},
+ eps_im_Scat = {2.9524 , 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 = {-44.9757 , Name StrCat[pp4,"/layer 4: real part of relative permittivity"]},
+ eps_im_L_4 = {2.9524 , Name StrCat[pp4,"/layer 4: imag part of relative permittivity"]},
+ eps_re_L_5 = {-44.9757 , Name StrCat[pp4,"/layer 5: real part of relative permittivity"]},
+ eps_im_L_5 = {2.9524 , Name StrCat[pp4,"/layer 5: imag part of relative permittivity"]},
+ eps_re_L_6 = {-44.9757 , Name StrCat[pp4,"/layer 6: real part of relative permittivity"]},
+ eps_im_L_6 = {2.9524 , 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,"/2Scatterer absolute mesh size [nm]"]},
+ PML_top = {lambda0 , Name StrCat[pp5,"/4PML top thickness [nm]"]},
+ PML_bot = {lambda0/10, 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= {2 , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
+ refine_mesh_L5= {2 , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
+ refine_mesh_L6= {2 , 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 = { 0 , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
+ FlagOutEscaCuts = { 0 , Name StrCat[pp6,"/3Output Scattered Electric Field cuts?"] , Choices {0,1} },
+ FlagOutPoyCut = { 0 , Name StrCat[pp6,"/4Output Poynting cuts?"] , Choices {0,1} },
+ FlagOutEtotFull = { 0 , Name StrCat[pp6,"/5Total Electric Field Full Output?"] , Choices {0,1} },
+ FlagOutEscaFull = { 1 , 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;