// Copyright (C) 2018 Guillaume Demésy
//
// This file is part of the model NonLinearEVP.pro.
// 
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License
// along with This program. If not, see <https://www.gnu.org/licenses/>.

pp0        = "1Geometry/0";
pp1        = "2Polarization-Bloch/0";
pp3        = "3Eigenvalue problem parameters/0";
pp2        = "4Drude Model parameters/0";
pp4        = "5Discretization/0";
pp5        = "6Formulations/0";
pp6        = "7Output/0";
close_menu = 0;
colorpp0   = "MintCream";
colorpp1   = "Ivory";
colorpp2   = "PaleGoldenRod";
colorpp3   = "Ivory";
colorpp4   = "LightSalmon";
colorpp5   = "Ivory";

DefineConstant[ a_lat = {50     , Name StrCat[pp0  , "1grating period d [nm]"]          , Highlight Str[colorpp0]  , Closed close_menu} ];


DefineConstant[
  d_sq           = {0.806 , Name StrCat[pp0 , "2square side size (fraction of period) [-]"] , Highlight Str[colorpp0]  , Min 0.01, Max 0.99} , 
  space2pml      = {1     , Name StrCat[pp0 , "3PML distance to object [d]"]  , Highlight Str[colorpp0], Min 0.1, Max 2} , 
  pmlsize        = {3     , Name StrCat[pp0 , "4PML thickness [d]"] , Highlight Str[colorpp0], Min 0.5, Max 4.} , 
  
  flag_Hparallel = {1     , Name StrCat[pp1 , "1polarization case"] , Choices {0="E //",1="H //"} },
  kx             = {0.75  , Name StrCat[pp1 , "2kx [Pi\a]"] , Highlight Str[colorpp1] , Min 0., Max 1.} , 

  eps_oo_1       = {1     , Name StrCat[pp2 , "0epsilon_inf_Drude [ - ]"] , Highlight Str[colorpp2]  , Closed close_menu} , 
  om_d_1         = {1.1   , Name StrCat[pp2 , "1Omega_Drude [2Pic\a]"]    , Highlight Str[colorpp2]  , Closed close_menu} , 
  gam_1          = {0.05  , Name StrCat[pp2 , "2Gamma_Drude [2Pic\a]"]    , Highlight Str[colorpp2]  , Closed close_menu} , 

  neig          =  {1     , Name StrCat[pp3 , "0Number of eigenvalues [int]"]  , Highlight Str[colorpp3]  , Min 1, Max 200} ,
  eig_target_re  = {0.0077, Name StrCat[pp3 , "1EV real part target [2Pic\a]"] , Highlight Str[colorpp3]  , Closed close_menu} , 
  eig_target_im  = {0.2598, Name StrCat[pp3 , "2EV imag part target [2Pic\a]"] , Highlight Str[colorpp3]  , Closed close_menu} , 
  eig_min_re     = {0.    , Name StrCat[pp3 , "3EV real min [2Pic\a]"] , Highlight Str[colorpp3]  , Closed close_menu} , 
  eig_max_re     = {0.1   , Name StrCat[pp3 , "4EV real max [2Pic\a]"] , Highlight Str[colorpp3]  , Closed close_menu} , 
  eig_min_im     = {0.0005, Name StrCat[pp3 , "5EV imag min [2Pic\a]"] , Highlight Str[colorpp3]  , Closed close_menu} , 
  eig_max_im     = {4.    , Name StrCat[pp3 , "6EV imag max [2Pic\a]"] , Highlight Str[colorpp3]  , Closed close_menu} , 

  paramaille     = {4     , Name StrCat[pp4 , "0number of mesh elements per period []"]  , Highlight Str[colorpp4], Min 2, Max 10} , 
  flag_Tmesh     = {0     , Name StrCat[pp4 , "2locally structured mesh?"] , Choices {0="unstruct",1="struct"} },
  flag_o2g       = {1     , Name StrCat[pp4 , "3Geometry order"] , Choices {0="o1 (linear)",1="o2 (curved)"} },
  flag_o2i       = {1     , Name StrCat[pp4 , "4Interpolation order"] , Choices {0="o1",1="o2"} },
  flag_rounding  = {1     , Name StrCat[pp4 , "5Corner rounding/0Do it!"], Choices{0,1}, ServerAction "ResetDatabase"},
  corner_rad_frac= {0.1   , Name StrCat[pp4 , "5Corner rounding/1corner radius (fraction of square side)"] , Highlight Str[colorpp2]  , Min 0.01, Max 0.49},
  flag_outEigvec = {1     , Name StrCat[pp4 , "7output eigenvector?"], Choices{0,1}},
  
  flag_res       = {4     , Name StrCat[pp5  , "0resolution type"], 
                        // Choices {0="Aux_E" ,1="PEP_E" ,2="NEP_E" ,3="Lag_E" ,4="PEP_h", 5="all"},ServerAction "ResetDatabase"},
                        Choices {0="Aux_E" ,1="PEP_E" ,2="NEP_E" ,3="Lag_E" ,4="PEP_h"},ServerAction "ResetDatabase"}
];

// Normalized units so that 2*pi*c/a=1
cel      = a_lat/(2*Pi);
epsf     = 8.854187817e-3;
muf      = 400.*Pi;
nm       = 2*Pi/(a_lat*Sqrt[epsf*muf]);
epsilon0 = epsf*nm;
mu0      = muf*nm;
norm     = a_lat/(2.*Pi*cel);

// Normalize
d_sq           = d_sq          * a_lat;
space2pml      = space2pml     * a_lat;
pmlsize        = pmlsize       * a_lat;
kx             = kx            * Pi/a_lat;
eig_target_re  = eig_target_re / norm;
eig_target_im  = eig_target_im / norm;
eig_min_re     = eig_min_re    / norm;
eig_max_re     = eig_max_re    / norm;
eig_min_im     = eig_min_im    / norm;
eig_max_im     = eig_max_im    / norm;
om_d_1         = om_d_1        / norm;
gam_1          = gam_1         / norm;
corner_rad     = d_sq*corner_rad_frac;

eps_oo_2       = 1;
om_d_2         = 0;
gam_2          = 0;

slepc_options_rg = StrCat(" -rg_interval_endpoints ",
                     Sprintf("%.8lf,",eig_min_re),
                     Sprintf("%.8lf,",eig_max_re),
                     Sprintf("%.8lf,",eig_min_im),
                     Sprintf("%.8lf",eig_max_im));