diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index b9396c364e3f77acc2edd4528094d9a1131fe6c9..8bebf2071c850ce8c26627a7efb4336febe2ebe2 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -125,27 +125,37 @@ Function{
 	a_pml           = 1.;
 	b_pml           = 1.;
   // bermu
-  n1[]    = Sqrt[epsr1[]];
-  n2[]    = Sqrt[epsr2[]];
-  k1norm[]= k0*n1[];
-  k2norm[]= k0*n2[];
-  Zmax = PML_top_hh;
-  Zmin = hh_L_6;
-  Damp_pml_top[] = 1/Abs[Zmax + PML_top - Z[]];
-  Damp_pml_bot[] = 1/Abs[Zmin - PML_bot - Z[]];
-  // DampingProfileX[] = 1/(Xmax + SizePMLX - Fabs[X[]]) - 1/(SizePMLX);
-
-  sx            = 1.;
-	sz[]          = Complex[a_pml,b_pml];
-	sz_bermutop[] = Complex[1,Damp_pml_top[]/k1norm[]];
-	sz_bermubot[] = Complex[1,Damp_pml_bot[]/k2norm[]];
-	sy            = 1.;
+  n1[]     = Sqrt[epsr1[]];
+  n2[]     = Sqrt[epsr2[]];
+  k1norm[] = k0*n1[];
+  k2norm[] = k0*n2[];
+
+  Zmax     = PML_top_hh;
+  Zmin     = hh_L_6;
+  Damp_pml_top[] = 1/(Zmax + PML_top - Fabs[Z[]]) - 1/(PML_top);
+  Damp_pml_bot[] = 1/(Zmin + PML_top - Fabs[Z[]]) - 1/(PML_bot);
+  Sigma_top[] = 0.5*(Damp_pml_top[] + Fabs[Damp_pml_top[]]);
+  Sigma_bot[] = 0.5*(Damp_pml_bot[] + Fabs[Damp_pml_bot[]]);
+
+	If (PML_TYPE==0)
+    sz[]          = Complex[a_pml,b_pml];
+	Else
+    sz[PMLtop] = Complex[1,Damp_pml_top[]/k1norm[]];
+	  sz[PMLbot] = Complex[1,Damp_pml_bot[]/k2norm[]];
+  EndIf
+  sx = 1.;
+	sy = 1.;
 
 	epsr[PMLtop]  = Re[epsr1[]]*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
 	mur[PMLtop]   =             TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
   epsr[PMLbot]  = Re[epsr2[]]*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
   mur[PMLbot]   =             TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
 
+  // epsr[PMLtop]  = Re[epsr1[]]*TensorDiag[sz_bermutop[]*sy/sx,sx*sz_bermutop[]/sy,sx*sy/sz_bermutop[]];
+	// mur[PMLtop]   =             TensorDiag[sz_bermutop[]*sy/sx,sx*sz_bermutop[]/sy,sx*sy/sz_bermutop[]];
+  // epsr[PMLbot]  = Re[epsr2[]]*TensorDiag[sz_bermubot[]*sy/sx,sx*sz_bermubot[]/sy,sx*sy/sz_bermubot[]];
+  // mur[PMLbot]   =             TensorDiag[sz_bermubot[]*sy/sx,sx*sz_bermubot[]/sy,sx*sy/sz_bermubot[]];
+
   // epsr[PMLbot] = epsr2[];
   // mur[PMLbot]  = TensorDiag[1,1,1];
 
diff --git a/DiffractionGratings/grating3D_data_2Dlamellar.geo b/DiffractionGratings/grating3D_data_2Dlamellar.geo
index 27dc260a162f77b443dd8c429349f7c94488473d..ed3c4a729bdb24d5b302ed45622ce8c951700c7b 100644
--- a/DiffractionGratings/grating3D_data_2Dlamellar.geo
+++ b/DiffractionGratings/grating3D_data_2Dlamellar.geo
@@ -62,6 +62,7 @@ DefineConstant[
     refine_mesh_L_5= {0.3       , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {0.3       , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
     
     InterpSampling     = { 30   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 0    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_bisin.geo b/DiffractionGratings/grating3D_data_bisin.geo
index 65840b288e8c965f05dce0f04aa60128b6d30285..216ef0b2710c139bc04ae748e9b7f8ba6b3596d8 100644
--- a/DiffractionGratings/grating3D_data_bisin.geo
+++ b/DiffractionGratings/grating3D_data_bisin.geo
@@ -62,7 +62,8 @@ DefineConstant[
     refine_mesh_L_5= {1      , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {1      , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
-    
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
+
     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} },
diff --git a/DiffractionGratings/grating3D_data_checker.geo b/DiffractionGratings/grating3D_data_checker.geo
index 169cb1fa72d1f288ab9aa4f83512ab90711131ed..2e43fc200373707d35adaf8da31d31146fedf512 100644
--- a/DiffractionGratings/grating3D_data_checker.geo
+++ b/DiffractionGratings/grating3D_data_checker.geo
@@ -62,7 +62,8 @@ DefineConstant[
     refine_mesh_L_5= {1      , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {1      , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
-    
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
+
     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} },
diff --git a/DiffractionGratings/grating3D_data_conv.geo b/DiffractionGratings/grating3D_data_conv.geo
index 1e227c44c2e34dd160d2d31b9716d5a4d69477f6..ebe61a71689da2fdf3d373f51877b446207c1d0c 100644
--- a/DiffractionGratings/grating3D_data_conv.geo
+++ b/DiffractionGratings/grating3D_data_conv.geo
@@ -63,7 +63,8 @@ DefineConstant[
     refine_mesh_L_5 = {1       , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6 = {1       , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
     FlagLinkFacets  = {0       , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
-        
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
+
     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} },
diff --git a/DiffractionGratings/grating3D_data_halfellipsoid.geo b/DiffractionGratings/grating3D_data_halfellipsoid.geo
index 04cd83c5bc2f3803dcaa6d18aa256834beebdd58..e27f2d1b66a6c790473680cae94cd59cbcdcc90c 100644
--- a/DiffractionGratings/grating3D_data_halfellipsoid.geo
+++ b/DiffractionGratings/grating3D_data_halfellipsoid.geo
@@ -7,8 +7,7 @@ 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]"]},
-    lambda0       = {400  , Min 400, Max 800, Step 200, Name StrCat[pp1, "0wavelength [nm]"] , Loop 1},
+    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]"]},
@@ -63,6 +62,7 @@ DefineConstant[
     refine_mesh_L_5= {1     , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {1     , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
         
     InterpSampling     = { 10   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 0    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_hole.geo b/DiffractionGratings/grating3D_data_hole.geo
index d05712415eb148937abb803f6dd63dac47964632..32bf71279b1423c4b2368d8e9975b12b9db39d3b 100644
--- a/DiffractionGratings/grating3D_data_hole.geo
+++ b/DiffractionGratings/grating3D_data_hole.geo
@@ -68,7 +68,8 @@ DefineConstant[
     refine_mesh_L_5= {1      , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {1      , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
-    
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
+
     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} },
diff --git a/DiffractionGratings/grating3D_data_pyramid.geo b/DiffractionGratings/grating3D_data_pyramid.geo
index fb97f520459ccaf20660eddfc4f59d934d42dc56..579361b100742a5e1cce49e1295d8b12c95b4721 100644
--- a/DiffractionGratings/grating3D_data_pyramid.geo
+++ b/DiffractionGratings/grating3D_data_pyramid.geo
@@ -62,6 +62,7 @@ DefineConstant[
     refine_mesh_L_5= {1     , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {1     , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
 
     InterpSampling     = { 50   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_skew.geo b/DiffractionGratings/grating3D_data_skew.geo
index 7ceee288e05637a3136cc0fb06e57dd90d619d9a..2bf61c73e048e313d0e685bd570ef1235ca1f4c9 100644
--- a/DiffractionGratings/grating3D_data_skew.geo
+++ b/DiffractionGratings/grating3D_data_skew.geo
@@ -63,7 +63,8 @@ DefineConstant[
     refine_mesh_L_5= {1     , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {1     , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
-        
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
+
     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} },
diff --git a/DiffractionGratings/grating3D_data_solarcell.geo b/DiffractionGratings/grating3D_data_solarcell.geo
index 7c26bb74ea80a0617bb6e56552e23966e741c1f0..ccfd3c968465d3930db4930155f5fa1efa384ec5 100644
--- a/DiffractionGratings/grating3D_data_solarcell.geo
+++ b/DiffractionGratings/grating3D_data_solarcell.geo
@@ -63,7 +63,7 @@ DefineConstant[
     refine_mesh_L_5= {2      , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {2      , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
     FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
-
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
     
     InterpSampling     = { 20   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_torus.geo b/DiffractionGratings/grating3D_data_torus.geo
index 643534525c88cb4a8fc36887b0c4178c4907e101..bc27227d377d6b0b47cb37a3576cf86583b3a451 100644
--- a/DiffractionGratings/grating3D_data_torus.geo
+++ b/DiffractionGratings/grating3D_data_torus.geo
@@ -62,7 +62,8 @@ DefineConstant[
     refine_mesh_L_5= {1       , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {1       , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
     FlagLinkFacets = {1      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0,1}, Visible 0},
-    
+    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
+
     InterpSampling     = {15   , 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} },