From 52c4ded2ea8f7b260c5283291b377d3b2fd4f9d8 Mon Sep 17 00:00:00 2001
From: Guillaume Demesy <guillaume.demesy@fresnel.fr>
Date: Fri, 30 Oct 2020 11:46:13 +0100
Subject: [PATCH] extra interface layer for grating2D

---
 DiffractionGratings/grating2D.pro | 595 +-----------------------------
 1 file changed, 5 insertions(+), 590 deletions(-)

diff --git a/DiffractionGratings/grating2D.pro b/DiffractionGratings/grating2D.pro
index bbfaf79..c2811b6 100644
--- a/DiffractionGratings/grating2D.pro
+++ b/DiffractionGratings/grating2D.pro
@@ -16,595 +16,10 @@
 // along with This program. If not, see <https://www.gnu.org/licenses/>.
 
 
-/** Select the config. (See also grating2D.geo) **/
-// Include "grating2D_data_AnisotropicGrating.geo";
-// Include "grating2D_data_LamellarGrating.geo";
-// Include "grating2D_data_PhotonicCrystalSlab.geo";
-// Include "grating2D_data_ResonantGrating.geo";
-Include "grating2D_data_plasmonics.geo";
-
+Include "grating2D_data.geo";
 Include "grating2D_materials.pro";
-
-myDir = "res2D/";
-DefineConstant[
-  lambda0 = {lambda_min , Min lambda_min, Max lambda_max, Step (lambda_max-lambda_min)/(nb_lambdas-1), Name StrCat[pp2, "0wavelength [nm]"] , Loop 1, Highlight Str[colorpp2],Graph "200000200020", ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"}
-];
-lambda0    = lambda0    * nm;
-lambda_min = lambda_min * nm;
-lambda_max = lambda_max * nm;
-A          = 1.;
-nb_orders  = nb_orders+1;
-
-Group {
-  // Boundaries
-  SurfBlochLeft   = Region[SURF_BLOCH_X_LEFT];
-  SurfBlochRight  = Region[SURF_BLOCH_X_RIGHT];
-  SurfCutSubs1    = Region[SURF_INTEG_SUB1];
-  SurfCutSuper1   = Region[SURF_INTEG_SUP1];
-  SurfCutSubs2    = Region[SURF_INTEG_SUB2];
-  SurfCutSuper2   = Region[SURF_INTEG_SUP2];
-
-  // Domains
-  pmlbot          = Region[PMLBOT];
-  sub             = Region[{SUB,SURF_INTEG_SUB1,SURF_INTEG_SUB2}];
-  layer_dep       = Region[LAYERDEP];
-  rod_out         = Region[RODOUT];
-  For i In {0:N_rods-1:1}
-    rod~{i} = Region[{ROD~{i}}];
-    rods   += Region[{ROD~{i}}];
-  EndFor
-  layer_cov       = Region[LAYERCOV];
-  sup             = Region[{SUP,SURF_INTEG_SUP1,SURF_INTEG_SUP2}];
-  pmltop          = Region[PMLTOP];
-  Omega_source    = Region[{layer_dep,rod_out,rods,layer_cov}];
-  Omega_nosource  = Region[{pmltop,pmlbot,sup,sub}];
-  Omega           = Region[{Omega_source,Omega_nosource}];
-  Omega_top       = Region[{layer_dep,rod_out,rods,layer_cov,sup}];
-  Omega_bot       = Region[{sub}];
-  Omega_pml       = Region[{pmltop,pmlbot}];
-
-  Plot_domain     = Region[{Omega,-pmltop,-pmlbot}];
-  Plot_bnd        = Region[SURF_PLOT];
-  Print_point     = Region[PRINT_POINT];
-}
-
-Function{
-  I[] = Complex[0.0,1.0];
-  bndCol[Plot_bnd] = Complex[1,1];
-
-  epsr_re_interp_mat_1[]  = InterpolationLinear[$1]{ListAlt[lambdamat_1 ,epsr_mat_re_1 ]};
-  epsr_im_interp_mat_1[]  = InterpolationLinear[$1]{ListAlt[lambdamat_1 ,epsr_mat_im_1 ]};
-  epsr_re_interp_mat_2[]  = InterpolationLinear[$1]{ListAlt[lambdamat_2 ,epsr_mat_re_2 ]};
-  epsr_im_interp_mat_2[]  = InterpolationLinear[$1]{ListAlt[lambdamat_2 ,epsr_mat_im_2 ]};
-  epsr_re_interp_mat_3[]  = InterpolationLinear[$1]{ListAlt[lambdamat_3 ,epsr_mat_re_3 ]};
-  epsr_im_interp_mat_3[]  = InterpolationLinear[$1]{ListAlt[lambdamat_3 ,epsr_mat_im_3 ]};
-  epsr_re_interp_mat_4[]  = InterpolationLinear[$1]{ListAlt[lambdamat_4 ,epsr_mat_re_4 ]};
-  epsr_im_interp_mat_4[]  = InterpolationLinear[$1]{ListAlt[lambdamat_4 ,epsr_mat_im_4 ]};
-  epsr_re_interp_mat_5[]  = InterpolationLinear[$1]{ListAlt[lambdamat_5 ,epsr_mat_re_5 ]};
-  epsr_im_interp_mat_5[]  = InterpolationLinear[$1]{ListAlt[lambdamat_5 ,epsr_mat_im_5 ]};
-  epsr_re_interp_mat_6[]  = InterpolationLinear[$1]{ListAlt[lambdamat_6 ,epsr_mat_re_6 ]};
-  epsr_im_interp_mat_6[]  = InterpolationLinear[$1]{ListAlt[lambdamat_6 ,epsr_mat_im_6 ]};
-  epsr_re_interp_mat_7[]  = InterpolationLinear[$1]{ListAlt[lambdamat_7 ,epsr_mat_re_7 ]};
-  epsr_im_interp_mat_7[]  = InterpolationLinear[$1]{ListAlt[lambdamat_7 ,epsr_mat_im_7 ]};
-  epsr_re_interp_mat_8[]  = InterpolationLinear[$1]{ListAlt[lambdamat_8 ,epsr_mat_re_8 ]};
-  epsr_im_interp_mat_8[]  = InterpolationLinear[$1]{ListAlt[lambdamat_8 ,epsr_mat_im_8 ]};
-  epsr_re_interp_mat_9[]  = InterpolationLinear[$1]{ListAlt[lambdamat_9 ,epsr_mat_re_9 ]};
-  epsr_im_interp_mat_9[]  = InterpolationLinear[$1]{ListAlt[lambdamat_9 ,epsr_mat_im_9 ]};
-  epsr_re_interp_mat_10[] = InterpolationLinear[$1]{ListAlt[lambdamat_10,epsr_mat_re_10]};
-  epsr_im_interp_mat_10[] = InterpolationLinear[$1]{ListAlt[lambdamat_10,epsr_mat_im_10]};
-  epsr_re_interp_mat_11[] = InterpolationLinear[$1]{ListAlt[lambdamat_11,epsr_mat_re_11]};
-  epsr_im_interp_mat_11[] = InterpolationLinear[$1]{ListAlt[lambdamat_11,epsr_mat_im_11]};
-  epsr_re_interp_mat_12[] = InterpolationLinear[$1]{ListAlt[lambdamat_12,epsr_mat_re_12]};
-  epsr_im_interp_mat_12[] = InterpolationLinear[$1]{ListAlt[lambdamat_12,epsr_mat_im_12]};
-  epsr_re_interp_mat_13[] = InterpolationLinear[$1]{ListAlt[lambdamat_13,epsr_mat_re_13]};
-  epsr_im_interp_mat_13[] = InterpolationLinear[$1]{ListAlt[lambdamat_13,epsr_mat_im_13]};
-  epsr_re_interp_mat_14[] = InterpolationLinear[$1]{ListAlt[lambdamat_14,epsr_mat_re_14]};
-  epsr_im_interp_mat_14[] = InterpolationLinear[$1]{ListAlt[lambdamat_14,epsr_mat_im_14]};
-  epsr_re_interp_mat_15[] = InterpolationLinear[$1]{ListAlt[lambdamat_15,epsr_mat_re_15]};
-  epsr_im_interp_mat_15[] = InterpolationLinear[$1]{ListAlt[lambdamat_15,epsr_mat_im_15]};
-  epsr_re_interp_mat_16[] = InterpolationLinear[$1]{ListAlt[lambdamat_16,epsr_mat_re_16]};
-  epsr_im_interp_mat_16[] = InterpolationLinear[$1]{ListAlt[lambdamat_16,epsr_mat_im_16]};
-
-  For i In {1:5}
-    For j In {1:nb_available_materials}
-      If(flag_mat~{i}==j-1)
-        epsr_re_dom~{i}[] = epsr_re_interp_mat~{j}[lambda0/nm*1e-9];
-        epsr_im_dom~{i}[] = epsr_im_interp_mat~{j}[lambda0/nm*1e-9];
-      EndIf
-    EndFor
-  EndFor
-  For k In {0:nb_available_lossless_materials-1}
-    If(flag_mat_6==lossless_material_list(k)-1)
-      epsr_re_dom_6[] = epsr_re_interp_mat~{lossless_material_list(k)}[lambda0/nm*1e-9];
-      epsr_im_dom_6[] = epsr_im_interp_mat~{lossless_material_list(k)}[lambda0/nm*1e-9];
-    EndIf
-  EndFor
-
-  epsr2_re[]       = epsr_re_dom_1[];
-  epsr2_im[]       = epsr_im_dom_1[];
-  epsr_layer_dep_re[] = epsr_re_dom_2[];
-  epsr_layer_dep_im[] = epsr_im_dom_2[];
-  epsr_rod_out_re[]   = epsr_re_dom_3[];
-  epsr_rod_out_im[]   = epsr_im_dom_3[];
-  epsr_rods_re[]      = epsr_re_dom_4[];
-  epsr_rods_im[]      = epsr_im_dom_4[];
-  epsr_layer_cov_re[] = epsr_re_dom_5[];
-  epsr_layer_cov_im[] = epsr_im_dom_5[];
-  epsr1_re[]       = epsr_re_dom_6[];
-  epsr1_im[]       = epsr_im_dom_6[];
-  
-  om0     = 2.*Pi*cel/lambda0;
-  k0      = 2.*Pi/lambda0;
-  epsr1[] = Complex[epsr1_re[],epsr1_im[]];
-  epsr2[] = Complex[epsr2_re[],epsr2_im[]];
-  n1[]    = Sqrt[epsr1[]];
-  n2[]    = Sqrt[epsr2[]];
-  k1norm[]= k0*n1[];
-  k2norm[]= k0*n2[];
-  alpha[] =  k0*n1[]*Sin[theta_deg*deg2rad];
-  beta1[] = -k0*n1[]*Cos[theta_deg*deg2rad];
-  beta2[] = -Sqrt[k0^2*epsr2[]-alpha[]^2];
-  
-  k1[] = Vector[alpha[], beta1[],0];
-  k1r[]= Vector[alpha[],-beta1[],0];
-  k2[] = Vector[alpha[], beta2[],0];
-
-  If (flag_Hparallel==1)
-    r[]    = (CompY[k1[]]*epsr2[]-CompY[k2[]]*epsr1[])/(CompY[k1[]]*epsr2[]+CompY[k2[]]*epsr1[]);
-    t[]    =                  (2.*CompY[k1[]]*epsr2[])/(CompY[k1[]]*epsr2[]+CompY[k2[]]*epsr1[]);
-    Pinc[] = 0.5*A^2*Sqrt[mu0/(ep0*epsr1_re[])] * Cos[theta_deg*deg2rad];
-  EndIf
-  If (flag_Hparallel==0)
-    r[]    = (CompY[k1[]]-CompY[k2[]])/(CompY[k1[]]+CompY[k2[]]);
-    t[]    =           2.*CompY[k1[]] /(CompY[k1[]]+CompY[k2[]]);
-    Pinc[] =  0.5*A^2*Sqrt[ep0*epsr1_re[]/mu0] * Cos[theta_deg*deg2rad]; 
-  EndIf
-
-  BlochX_phase_shift[] = Exp[I[]*alpha[]*d];
-
-  For i In {0:2*nb_orders}
-    alpha~{i}[]      = CompX[k1[]]  + 2.*Pi/d*(i-nb_orders);
-    expmialpha~{i}[] = Exp[-I[]*alpha~{i}[]*X[]];
-    beta_super~{i}[] = -Sqrt[k0^2*epsr1[]-alpha~{i}[]^2];
-    beta_subs~{i}[]  = -Sqrt[k0^2*epsr2[]-alpha~{i}[]^2];
-  EndFor
-  
-  a_pml        = 1.;
-  b_pml        = 1.;
-  sx           = 1.;
-  sy[]         = Complex[a_pml,b_pml]; 
-  sz           = 1.;
-  PML_Tensor[] = TensorDiag[sz*sy[]/sx,sx*sz/sy[],sx*sy[]/sz];
-
-  epsilonr[sup]           = epsr1[] * TensorDiag[1,1,1];
-  epsilonr[sub]           = epsr2[] * TensorDiag[1,1,1];
-  epsilonr[layer_dep]     = Complex[epsr_layer_dep_re[],epsr_layer_dep_im[]] * TensorDiag[1,1,1];
-  epsilonr[rod_out]       = Complex[epsr_rod_out_re[],epsr_rod_out_im[]    ] * TensorDiag[1,1,1];
-  If (Flag_aniso==0)
-    epsilonr[rods]        = Complex[epsr_rods_re[],epsr_rods_im[]] * TensorDiag[1,1,1];
-  Else
-    epsilonr[rods]        = Tensor[Complex[epsr_custom_anisoXX_re, epsr_custom_anisoXX_im],
-           Complex[epsr_custom_anisoXY_re, epsr_custom_anisoXY_im],
-           0,
-           Complex[epsr_custom_anisoXY_re,-epsr_custom_anisoXX_im],
-           Complex[epsr_custom_anisoYY_re,epsr_custom_anisoYY_im],
-           0,
-           0,
-           0,
-           Complex[epsr_custom_anisoZZ_re,epsr_custom_anisoZZ_im]
-           ];
-  EndIf
-  epsilonr[layer_cov]       = Complex[epsr_layer_cov_re[],epsr_layer_cov_im[]] * TensorDiag[1,1,1];
-  epsilonr[pmltop]          = epsr1_re[]*PML_Tensor[];
-  epsilonr[pmlbot]          = epsr2_re[]*PML_Tensor[]; 
-  
-  epsilonr_annex[sub]       = epsr2[] * TensorDiag[1,1,1];
-  epsilonr_annex[sup]       = epsr1[] * TensorDiag[1,1,1];
-  epsilonr_annex[layer_dep] = epsr1[] * TensorDiag[1,1,1];
-  epsilonr_annex[rod_out]   = epsr1[] * TensorDiag[1,1,1];
-  epsilonr_annex[rods]      = epsr1[] * TensorDiag[1,1,1];
-  epsilonr_annex[layer_cov] = epsr1[] * TensorDiag[1,1,1];
-  epsilonr_annex[pmltop]    = epsr1_re[]*PML_Tensor[];
-  epsilonr_annex[pmlbot]    = epsr2_re[]*PML_Tensor[];
-
-  mur[pmltop]     = PML_Tensor[];
-  mur[pmlbot]     = PML_Tensor[];
-  mur[sub]        = TensorDiag[1,1,1];
-  mur[sup]        = TensorDiag[1,1,1];
-  mur[layer_dep]  = TensorDiag[1,1,1];
-  mur[rods]       = TensorDiag[1,1,1];
-  mur[rod_out]    = TensorDiag[1,1,1];
-  mur[layer_cov]  = TensorDiag[1,1,1];
-
-  ui[pmltop]    = 0.;
-  ui[pmlbot]    = 0.;
-  ui[sup]       = A*Exp[I[]*k1[]*XYZ[]];
-  ui[layer_cov] = A*Exp[I[]*k1[]*XYZ[]];
-  ui[rod_out]   = A*Exp[I[]*k1[]*XYZ[]];
-  ui[rods]      = A*Exp[I[]*k1[]*XYZ[]];
-  ui[layer_dep] = A*Exp[I[]*k1[]*XYZ[]];
-  ui[sub]       = 0.;
-  
-  ur[pmltop]    = 0.;
-  ur[pmlbot]    = 0.;
-  ur[sup]       = r[]*Exp[I[]*k1r[]*XYZ[]];
-  ur[layer_dep] = r[]*Exp[I[]*k1r[]*XYZ[]];
-  ur[rod_out]   = r[]*Exp[I[]*k1r[]*XYZ[]];
-  ur[rods]      = r[]*Exp[I[]*k1r[]*XYZ[]];
-  ur[layer_cov] = r[]*Exp[I[]*k1r[]*XYZ[]];
-  ur[sub]       = 0;
-
-  ut[pmltop]    = 0.;
-  ut[pmlbot]    = 0.;
-  ut[sup]       = 0.;
-  ut[layer_dep] = 0.;
-  ut[rod_out]   = 0.;
-  ut[rods]      = 0.;
-  ut[layer_cov] = 0.;
-  ut[sub]       = t[]*Exp[I[]*k2[]*XYZ[]];
-  
-  u1[]          = ui[]+ur[]+ut[];
-  u1d[]         = ur[]+ut[];
-  
-  If (flag_Hparallel==1)
-    E1i[]  = -1/(om0*ep0*epsilonr_annex[]) * k1[]  /\ Vector[0,0,ui[]];
-    E1r[]  = -1/(om0*ep0*epsilonr_annex[]) * k1r[] /\ Vector[0,0,ur[]];
-    E1t[]  = -1/(om0*ep0*epsilonr_annex[]) * k2[]  /\ Vector[0,0,ut[]];
-    E1[]   = E1i[]+E1r[]+E1t[];
-    Ex1[]  = CompX[E1[]];
-    Ey1[]  = CompY[E1[]];
-
-    detepst[]       = CompXX[epsilonr[]]*CompYY[epsilonr[]]-CompXY[epsilonr[]]*Conj[CompYX[epsilonr[]]];
-    detepst_annex[] = CompXX[epsilonr_annex[]]*CompYY[epsilonr_annex[]]-CompXY[epsilonr_annex[]]*Conj[CompYX[epsilonr_annex[]]];
-    xsi[]           = Transpose[epsilonr[]]/detepst[];
-    xsi_annex[]     = Transpose[epsilonr_annex[]]/detepst_annex[];
-    chi[]           = CompZZ[mur[]];
-    source_xsi_r[]      = (xsi[]-xsi_annex[]) * -k1[]  * I[] * ui[];
-    source_xsi_i[]      = (xsi[]-xsi_annex[]) * -k1r[] * I[] * ur[];
-    source_chi_r[]      = 0;
-    source_chi_i[]      = 0;
-
-  Else
-    detmut[]    = CompXX[mur[]]*CompYY[mur[]]-CompXY[mur[]]*Conj[CompYX[mur[]]];
-    xsi[]       = Transpose[mur[]]/detmut[];
-    chi[]       = CompZZ[epsilonr[]];
-    chi_annex[] = CompZZ[epsilonr_annex[]];
-    source_xsi_r[] = Vector[0.,0.,0.];
-    source_xsi_i[] = Vector[0.,0.,0.];
-    source_chi_r[] = k0^2 * (chi[]-chi_annex[]) * ur[];
-    source_chi_i[] = k0^2 * (chi[]-chi_annex[]) * ui[];
-  EndIf
-}
-
-Constraint {
-  {Name Bloch;
-    Case {
-      { Region SurfBlochRight; Type LinkCplx ; RegionRef SurfBlochLeft;
-               Coefficient BlochX_phase_shift[]; Function Vector[$X-d,$Y,$Z] ;
-      }
-    }
-  }
-}
-
-Jacobian {
-  { Name JVol ;
-    Case { 
-      { Region All ; Jacobian Vol ; }
-    }
-  }
-  { Name JSur ;
-    Case {
-      { Region All ; Jacobian Sur ; }
-    }
-  }
-}
-
-Integration {
-  If (Flag_o2_geom==0)
-    { Name Int_1 ;
-      Case { 
-        { Type Gauss ;
-          Case { 
-            { GeoElement Point    ; NumberOfPoints  1 ; }
-            { GeoElement Line     ; NumberOfPoints  4 ; }
-            { GeoElement Triangle ; NumberOfPoints  12 ; }
-          }
-        }
-      }
-    }
-  EndIf
-  If (Flag_o2_geom==1)
-    { Name Int_1 ;
-      Case { 
-        { Type Gauss ;
-          Case { 
-            { GeoElement Point     ; NumberOfPoints  1 ; }
-            { GeoElement Line2     ; NumberOfPoints  4 ; }
-            { GeoElement Triangle2 ; NumberOfPoints  12 ; }
-          }
-        }
-      }
-    }
-  EndIf
-}
-
-FunctionSpace {
-  { Name Hgrad; Type Form0;
-    BasisFunction {
-      { Name sn;  NameOfCoef un;  Function BF_Node;    Support Region[Omega]; Entity NodesOf[Omega]; }
-      { Name sn2; NameOfCoef un2; Function BF_Node_2E; Support Region[Omega]; Entity EdgesOf[Omega]; }
-      If (Flag_o2_interp==1)
-        { Name un3; NameOfCoef un3; Function BF_Node_2F; Support Region[Omega]; Entity FacetsOf[Omega]; }
-        { Name un4; NameOfCoef un4; Function BF_Node_3E; Support Region[Omega]; Entity EdgesOf[Omega]; }
-        { Name un5; NameOfCoef un5; Function BF_Node_3F; Support Region[Omega]; Entity FacetsOf[Omega]; }
-      EndIf
-      
-    }
-    Constraint {
-      { NameOfCoef un;  EntityType NodesOf ; NameOfConstraint Bloch; }
-      { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Bloch; }
-      If (Flag_o2_interp==1)
-        { NameOfCoef un3; EntityType EdgesOf ; NameOfConstraint Bloch; }
-      EndIf
-    }
-  }
-}
-
-Formulation {
-  {Name helmoltz_scalar; Type FemEquation;
-    Quantity {
-      { Name u2d; Type Local; NameOfSpace Hgrad;}}
-    Equation {
-      Galerkin { [-xsi[]*Dof{Grad u2d} , {Grad u2d}]; In Omega; Jacobian JVol; Integration Int_1; }
-      Galerkin { [k0^2*chi[]*Dof{u2d}  ,      {u2d}]; In Omega; Jacobian JVol; Integration Int_1; }
-      Galerkin { [source_xsi_i[]       , {Grad u2d}]; In Omega; Jacobian JVol; Integration Int_1; }
-      Galerkin { [source_xsi_r[]       , {Grad u2d}]; In Omega; Jacobian JVol; Integration Int_1; }        
-      Galerkin { [source_chi_r[]       ,      {u2d}]; In Omega; Jacobian JVol; Integration Int_1; }
-      Galerkin { [source_chi_i[]       ,      {u2d}]; In Omega; Jacobian JVol; Integration Int_1; }
-    }
-  }
-}
-
-Resolution {
-  { Name helmoltz_scalar;
-    System {
-      { Name S; NameOfFormulation helmoltz_scalar; Type ComplexValue; }
-    }
-    Operation {
-      CreateDir[Str[myDir]];
-      Printf["%f",lambda0/nm];
-      Generate[S] ;
-      Solve[S] ;
-    }
-  }
-}
-
-PostProcessing {
-  { Name postpro_energy; NameOfFormulation helmoltz_scalar;
-    Quantity {
-      If (flag_Hparallel==1)
-        { Name debr         ; Value { Local { [ r[]                 ]; In Omega; Jacobian JVol; } } }
-        { Name debt         ; Value { Local { [ t[]                 ]; In Omega; Jacobian JVol; } } }
-        { Name u            ; Value { Local { [ {u2d}               ]; In Omega; Jacobian JVol; } } }
-        { Name epsr         ; Value { Local { [ CompZZ[epsilonr[]]  ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_diff      ; Value { Local { [ {u2d}+u1d[]         ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_tot       ; Value { Local { [ {u2d}+u1[]          ]; In Omega; Jacobian JVol; } } }
-        { Name NormHz_tot   ; Value { Local { [ Norm[{u2d}+u1[]]    ]; In Omega; Jacobian JVol; } } }
-        { Name E1           ; Value { Local { [ E1[]                ]; In Omega; Jacobian JVol; } } }
-        { Name u1           ; Value { Local { [ u1[]                ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totp1     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 1*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totp2     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 2*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totp3     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 3*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totp4     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 4*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totm1     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-1*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totm2     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-2*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totm3     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-3*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Hz_totm4     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-4*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name boundary     ; Value { Local { [ bndCol[] ] ; In Plot_bnd ; Jacobian JVol ; } } }
-
-        For i In {0:2*nb_orders}
-          { Name s_r~{i} ; Value { Integral{ [ expmialpha~{i}[] * ({u2d}+u1d[])/d ] ; In SurfCutSuper1 ; Jacobian JSur ; Integration Int_1 ; } } }
-          { Name s_t~{i} ; Value { Integral{ [ expmialpha~{i}[] * ({u2d}+u1d[])/d ] ; In SurfCutSubs1  ; Jacobian JSur ; Integration Int_1 ; } } }
-          { Name order_t_angle~{i} ; Value { Local{ [-Atan2[Re[alpha~{i}[]],Re[beta_subs~{i}[]]]/deg2rad ] ; In Omega; Jacobian JVol; } } }
-          { Name order_r_angle~{i} ; Value { Local{ [ Atan2[Re[alpha~{i}[]],Re[beta_super~{i}[]]]/deg2rad ] ; In Omega; Jacobian JVol; } } }
-        EndFor
-        For i In {0:2*nb_orders}
-          { Name eff_r~{i} ; Value { Term{ Type Global; [ SquNorm[#i]*beta_super~{i}[]/beta1[] ] ; In SurfCutSuper1 ; } } }
-          { Name eff_t~{i} ; Value { Term{ Type Global; [ SquNorm[#(2*nb_orders+1+i)]*(beta_subs~{i}[]/beta1[])*(epsr1[]/epsr2[]) ] ; In SurfCutSubs1 ; } } }
-        EndFor
-        For i In {0:N_rods-1:1}
-          { Name Q_rod~{i}   ; Value { Integral { [  0.5 * ep0*om0*Fabs[epsr_rods_im[]]          * ( SquNorm[CompY[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[-CompX[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ey1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] ) / (Pinc[]*d) ] ; In rod~{i}   ; Integration Int_1 ; Jacobian JVol ; } } }
-        EndFor
-        { Name Q_subs        ; Value { Integral { [  0.5 * ep0*om0*Fabs[epsr2_im[]      ]        * ( SquNorm[CompY[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[-CompX[{Grad u2d}]*I[]/(om0*ep0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In sub         ; Integration Int_1 ; Jacobian JVol ; } } }
-        { Name Q_rod_out     ; Value { Integral { [  0.5 * ep0*om0*Fabs[epsr_rod_out_im[]]       * ( SquNorm[CompY[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[-CompX[{Grad u2d}]*I[]/(om0*ep0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In rod_out   ; Integration Int_1 ; Jacobian JVol ; } } }
-        { Name Q_layer_dep   ; Value { Integral { [  0.5 * ep0*om0*Fabs[epsr_layer_dep_im[]]     * ( SquNorm[CompY[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[-CompX[{Grad u2d}]*I[]/(om0*ep0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In layer_dep   ; Integration Int_1 ; Jacobian JVol ; } } }
-        { Name Q_layer_cov   ; Value { Integral { [  0.5 * ep0*om0*Fabs[epsr_layer_cov_im[]]     * ( SquNorm[CompY[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[-CompX[{Grad u2d}]*I[]/(om0*ep0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In layer_cov   ; Integration Int_1 ; Jacobian JVol ; } } }
-        { Name Q_tot         ; Value { Integral { [  0.5 * ep0*om0*Fabs[Im[CompZZ[epsilonr[]]]]  * ( SquNorm[CompY[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[-CompX[{Grad u2d}]*I[]/(om0*ep0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In Plot_domain ; Integration Int_1 ; Jacobian JVol ; } } }
-        { Name lambda_step   ; Value { Local { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
-      EndIf
-      If (flag_Hparallel==0)
-        { Name epsr      ; Value { Local { [ CompZZ[epsilonr[]] ]; In Omega; Jacobian JVol; } } }
-        { Name Ez_diff   ; Value { Local { [ {u2d}+u1d[]        ]; In Omega; Jacobian JVol; } } }
-        { Name Ez_tot    ; Value { Local { [ {u2d}+u1[]          ]; In Omega; Jacobian JVol; } } }
-        { Name Ez_totp1  ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 1*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Ez_totp2  ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 2*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Ez_totp3  ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 3*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Ez_totp4  ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 4*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Ez_totm1  ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-1*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Ez_totm2  ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-2*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Ez_totm3  ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-3*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name Ez_totm4  ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-4*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
-        { Name boundary  ; Value { Local { [ bndCol[] ] ; In Plot_bnd ; Jacobian JVol ; } } }
-        
-        For i In {0:2*nb_orders}
-          { Name s_r~{i} ; Value { Integral{ [ expmialpha~{i}[] * ({u2d}+u1d[])/d ] ; In SurfCutSuper1 ; Jacobian JSur ; Integration Int_1 ; } } }
-          { Name s_t~{i} ; Value { Integral{ [ expmialpha~{i}[] * ({u2d}+u1d[])/d ] ; In SurfCutSubs1  ; Jacobian JSur ; Integration Int_1 ; } } }
-          { Name order_t_angle~{i} ; Value { Local{ [-Atan2[Re[alpha~{i}[]],Re[beta_subs~{i}[]]]/deg2rad ] ; In Omega; Jacobian JVol; } } }
-          { Name order_r_angle~{i} ; Value { Local{ [ Atan2[Re[alpha~{i}[]],Re[beta_super~{i}[]]]/deg2rad ] ; In Omega; Jacobian JVol; } } }
-        EndFor
-        For i In {0:2*nb_orders}
-          { Name eff_r~{i} ; Value { Term{ Type Global; [ SquNorm[#i]*beta_super~{i}[]/beta1[] ] ; In SurfCutSuper1 ; } } }
-          { Name eff_t~{i} ; Value { Term{ Type Global; [ SquNorm[#(2*nb_orders+1+i)]*(beta_subs~{i}[]/beta1[])] ; In SurfCutSubs1 ; } } }
-        EndFor
-        For i In {0:N_rods-1:1}
-          { Name Q_rod~{i}  ; Value { Integral { [0.5 * ep0*om0*Fabs[epsr_rods_im[]]          * (SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In rod~{i}    ; Integration Int_1 ; Jacobian JVol ; } } }
-        EndFor
-        { Name Q_subs       ; Value { Integral { [  0.5 * ep0*om0*Fabs[epsr2_im[]]            *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In sub         ; Integration Int_1 ; Jacobian JVol ; } } }
-        { Name Q_rod_out    ; Value { Integral { [  0.5 * ep0*om0*Fabs[epsr_rod_out_im[]]     *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In rod_out   ; Integration Int_1 ; Jacobian JVol ; } } }
-        { Name Q_layer_dep  ; Value { Integral { [  0.5 * ep0*om0*Fabs[epsr_layer_dep_im[]]   *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In layer_dep   ; Integration Int_1 ; Jacobian JVol ; } } }
-        { Name Q_layer_cov  ; Value { Integral { [  0.5 * ep0*om0*Fabs[epsr_layer_cov_im[]]   *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In layer_cov   ; Integration Int_1 ; Jacobian JVol ; } } }
-        { Name Q_tot        ; Value { Integral { [  0.5 * ep0*om0*Fabs[Im[CompXX[epsilonr[]]]]*(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In Plot_domain ; Integration Int_1 ; Jacobian JVol ; } } }
-        { Name lambda_step  ; Value { Local    { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
-      EndIf
-    }
-  }
-}
-
-PostOperation {
-  { Name postop_energy; NameOfPostProcessing postpro_energy ;
-    Operation {
-      If (flag_Hparallel==1)
-        Print[ lambda_step, OnPoint{0,0,0}, Format Table, File > StrCat[myDir, "temp_lambda_step.txt"], SendToServer "GetDP/Lambda_step" ] ;
-        For i In {0:2*nb_orders}
-          Print[ s_r~{i}[SurfCutSuper1], OnGlobal, Store i                , Format Table , File > StrCat[myDir, Sprintf("temp_s_r_%g.txt", i-nb_orders)]];
-          Print[ s_t~{i}[SurfCutSubs1] , OnGlobal, Store (2*nb_orders+1+i), Format Table , File > StrCat[myDir, Sprintf("temp_s_t_%g.txt", i-nb_orders)]];
-        EndFor
-        For i In {1:2*nb_orders-1}
-          Print[ eff_r~{i}[SurfCutSuper1], OnRegion SurfCutSuper1,Store (4*nb_orders+1+i), Format FrequencyTable, File > StrCat[myDir, Sprintf("efficiency_r_%g.txt", i-nb_orders)]];
-          Print[ eff_t~{i}[SurfCutSubs1] , OnRegion SurfCutSubs1 ,Store (6*nb_orders+1+i), Format FrequencyTable, File > StrCat[myDir, Sprintf("efficiency_t_%g.txt", i-nb_orders)]];
-          Print[ order_r_angle~{i}     , OnPoint{0,0,0}, Format Table , File > StrCat[myDir, Sprintf("order_r_angle_%g.txt", i-nb_orders)]];
-          Print[ order_t_angle~{i}     , OnPoint{0,0,0}, Format Table , File > StrCat[myDir, Sprintf("order_t_angle_%g.txt", i-nb_orders)]];          
-        EndFor
-        Print[ eff_r~{nb_orders}[SurfCutSuper1], OnRegion SurfCutSuper1, Format Table, SendToServer "GetDP/R0", File StrCat[myDir, "temp_R0.txt"]];
-        Print[ eff_t~{nb_orders}[SurfCutSubs1] , OnRegion SurfCutSubs1 , Format Table, SendToServer "GetDP/T0", File StrCat[myDir, "temp_T0.txt"]];
-        Print[ Q_tot[Plot_domain]             , OnGlobal  , Format FrequencyTable ,SendToServer "GetDP/total absorption", File > StrCat[myDir, "absorption-Q_tot.txt"]];
-        For i In {0:N_rods-1:1}
-          Print[ Q_rod~{i}[rod~{i}] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, Sprintf("absorption-Q_rod_%g.txt", i+1) ]];           
-        EndFor  
-        Print[ Q_subs[sub]            , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_subs.txt"]];
-        Print[ Q_rod_out[rod_out]     , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_rod_out.txt"]];
-        Print[ Q_layer_dep[layer_dep] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_dep.txt"]];
-        Print[ Q_layer_cov[layer_cov] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_cov.txt"]];
-        Print[ Hz_tot    , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm)];
-        // // debug
-        // Print[ u1    , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("u1%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("u1%.2fnm.pos", lambda0/nm)];
-        // Print[ E1    , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("E1%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("E1%.2fnm.pos", lambda0/nm)];
-        Echo[ Str["For i In {0:2}",
-                  "  View[i].LineWidth = 2;View[i].ColormapNumber = 15;View[0].AxesFormatX = '%.1f';",
-                  "EndFor"], File StrCat[myDir,"tmp0.geo"]] ;
-        // View[i].RangeType = 2;View[i].CustomMin=0.;View[i].CustomMax=1.;
-        If(multiplot)
-          Echo[ Str["For i In {PostProcessing.NbViews-1:0:-1}",
-                    "  If(!StrCmp(View[i].Name, 'boundary') || !StrCmp(View[i].Name, 'boundary_Combine'))",
-                    "    Delete View[i];",
-                    "  EndIf",
-                    "EndFor"], File StrCat[myDir,"tmp1.geo"]] ;
-          Print [ Hz_totp1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_2.pos", lambda0/nm)], ChangeOfCoordinates {$X+1*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-          Print [ Hz_totm1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_3.pos", lambda0/nm)], ChangeOfCoordinates {$X-1*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-          Print [ Hz_totp2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_4.pos", lambda0/nm)], ChangeOfCoordinates {$X+2*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-          Print [ Hz_totm2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_5.pos", lambda0/nm)], ChangeOfCoordinates {$X-2*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-          Print [ Hz_totp3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_6.pos", lambda0/nm)], ChangeOfCoordinates {$X+3*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-          Print [ Hz_totm3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_7.pos", lambda0/nm)], ChangeOfCoordinates {$X-3*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-          Print [ Hz_totp4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_8.pos", lambda0/nm)], ChangeOfCoordinates {$X+4*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-          Print [ Hz_totm4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_9.pos", lambda0/nm)], ChangeOfCoordinates {$X-4*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
-          Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo"]] ;
-          Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary1.pos"], ChangeOfCoordinates {$X+0*d,$Y,$Z} ] ;
-          Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary2.pos"], ChangeOfCoordinates {$X+1*d,$Y,$Z} ] ;
-          Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary3.pos"], ChangeOfCoordinates {$X-1*d,$Y,$Z} ] ;
-          Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary4.pos"], ChangeOfCoordinates {$X+2*d,$Y,$Z} ] ;
-          Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary5.pos"], ChangeOfCoordinates {$X-2*d,$Y,$Z} ] ;
-          Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary6.pos"], ChangeOfCoordinates {$X+3*d,$Y,$Z} ] ;
-          Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary7.pos"], ChangeOfCoordinates {$X-3*d,$Y,$Z} ] ;
-          Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary8.pos"], ChangeOfCoordinates {$X+4*d,$Y,$Z} ] ;
-          Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary9.pos"], ChangeOfCoordinates {$X-4*d,$Y,$Z} ];
-          Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo" ]] ;
-          Echo[ Str["Hide {",
-                    "Point{1,2,7,8,9,10,20,22};",
-                    "Line{1,7,8,9,10,30,32,34,2,3,4,5,6,12,16,20,24,28};",
-                    "Surface{36,48};}",
-                    "Geometry.Color.Lines = {0,0,0};",
-                    "l=PostProcessing.NbViews-1; View[l].ColorTable={Black}; ",
-                    "View[l-1].Visible=1; View[l-1].ShowScale=0;",
-                    "View[l].ShowScale=0; View[l].LineWidth=1.5; View[l].LineType=1;Geometry.LineWidth=0;"],
-                    File StrCat[myDir,"tmp3.geo" ]] ;
-        EndIf
-      EndIf
-      If (flag_Hparallel==0)
-        Print[ lambda_step, OnPoint{0,0,0}, Format Table, File StrCat[myDir, "temp_lambda_step.txt"], SendToServer "GetDP/Lambda_step" ] ;
-        For i In {0:2*nb_orders}
-          Print[ s_r~{i}[SurfCutSuper1], OnGlobal, Store i                , Format Table , File > StrCat[myDir, Sprintf("temp_s_r_%g.txt", i-nb_orders)]];
-          Print[ s_t~{i}[SurfCutSubs1] , OnGlobal, Store (2*nb_orders+1+i), Format Table , File > StrCat[myDir, Sprintf("temp_s_t_%g.txt", i-nb_orders)]];
-        EndFor
-        For i In {1:2*nb_orders-1}
-          Print[ eff_r~{i}[SurfCutSuper1], OnRegion SurfCutSuper1,Store (4*nb_orders+1+i), Format FrequencyTable, File > StrCat[myDir, Sprintf("efficiency_r_%g.txt", i-nb_orders)]];
-          Print[ eff_t~{i}[SurfCutSubs1] , OnRegion SurfCutSubs1 ,Store (6*nb_orders+1+i), Format FrequencyTable, File > StrCat[myDir, Sprintf("efficiency_t_%g.txt", i-nb_orders)]];
-          Print[ order_r_angle~{i}     , OnPoint{0,0,0}, Format Table , File > StrCat[myDir, Sprintf("order_r_angle_%g.txt", i-nb_orders)]];
-          Print[ order_t_angle~{i}     , OnPoint{0,0,0}, Format Table , File > StrCat[myDir, Sprintf("order_t_angle_%g.txt", i-nb_orders)]];          
-        EndFor
-        Print[ eff_r~{nb_orders}[SurfCutSuper1], OnRegion SurfCutSuper1, Format Table, SendToServer "GetDP/R0", File StrCat[myDir, "temp_R0.txt"]];
-        Print[ eff_t~{nb_orders}[SurfCutSubs1] , OnRegion SurfCutSubs1 , Format Table, SendToServer "GetDP/T0", File StrCat[myDir, "temp_T0.txt"]];
-        Print[ Q_tot[Plot_domain]             , OnGlobal  , Format FrequencyTable ,SendToServer "GetDP/total absorption", File > StrCat[myDir, "absorption-Q_tot.txt"]];
-        For i In {0:N_rods-1:1}
-          Print[ Q_rod~{i}[rod~{i}] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, Sprintf("absorption-Q_rod_%g.txt", i+1) ]];           
-        EndFor
-        Print[ Q_subs[sub]             , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_subs.txt"        ] ];
-        Print[ Q_rod_out[rod_out] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_rod_out.txt"] ];
-        Print[ Q_layer_dep[layer_dep] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_dep.txt"] ];
-        Print[ Q_layer_cov[layer_cov] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_cov.txt"] ];
-        Print [ Ez_tot  , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm)];
-        Echo[ Str["For i In {0:2}",
-                  "  View[i].LineWidth = 4;View[i].ColormapNumber = 15;View[0].AxesFormatX = '%.1f';",
-                  "EndFor"], File StrCat[myDir,"tmp0.geo"]] ;
-        If(multiplot)
-        Echo[ Str["For i In {PostProcessing.NbViews-1:0:-1}",
-                  "  If(!StrCmp(View[i].Name, 'boundary') || !StrCmp(View[i].Name, 'boundary_Combine'))",
-                  "    Delete View[i];",
-                  "  EndIf",
-                  "EndFor"], File StrCat[myDir,"tmp1.geo"]];
-        Print [ Ez_totp1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_2.pos", lambda0/nm)], ChangeOfCoordinates {$X+1*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-        Print [ Ez_totm1, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_3.pos", lambda0/nm)], ChangeOfCoordinates {$X-1*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-        Print [ Ez_totp2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_4.pos", lambda0/nm)], ChangeOfCoordinates {$X+2*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-        Print [ Ez_totm2, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_5.pos", lambda0/nm)], ChangeOfCoordinates {$X-2*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-        Print [ Ez_totp3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_6.pos", lambda0/nm)], ChangeOfCoordinates {$X+3*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-        Print [ Ez_totm3, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_7.pos", lambda0/nm)], ChangeOfCoordinates {$X-3*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-        Print [ Ez_totp4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_8.pos", lambda0/nm)], ChangeOfCoordinates {$X+4*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-        Print [ Ez_totm4, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_9.pos", lambda0/nm)], ChangeOfCoordinates {$X-4*d,$Y,$Z}, Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm) ] ;
-        Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo"]] ;
-        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary1.pos"], ChangeOfCoordinates {$X+0*d,$Y,$Z} ] ;
-        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary2.pos"], ChangeOfCoordinates {$X+1*d,$Y,$Z} ] ;
-        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary3.pos"], ChangeOfCoordinates {$X-1*d,$Y,$Z} ] ;
-        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary4.pos"], ChangeOfCoordinates {$X+2*d,$Y,$Z} ] ;
-        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary5.pos"], ChangeOfCoordinates {$X-2*d,$Y,$Z} ] ;
-        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary6.pos"], ChangeOfCoordinates {$X+3*d,$Y,$Z} ] ;
-        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary7.pos"], ChangeOfCoordinates {$X-3*d,$Y,$Z} ] ;
-        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary8.pos"], ChangeOfCoordinates {$X+4*d,$Y,$Z} ] ;
-        Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,"boundary9.pos"], ChangeOfCoordinates {$X-4*d,$Y,$Z} ];
-        Echo[ "Combine ElementsByViewName;", File StrCat[myDir,"tmp2.geo" ]] ;
-        Echo[ Str["Hide {",
-                  "Point{1,2,7,8,9,10,20,22};",
-                  "Line{1,7,8,9,10,30,32,34,2,3,4,5,6,12,16,20,24,28};",
-                  "Surface{36,48};}",
-                  "Geometry.Color.Lines = {0,0,0};",
-                  "l=PostProcessing.NbViews-1; View[l].ColorTable={Black}; ",
-                  "View[l-1].Visible=1; View[l-1].ShowScale=0;",
-                  "View[l].ShowScale=0; View[l].LineWidth=1.5; View[l].LineType=1;Geometry.LineWidth=0;"],
-                  File StrCat[myDir,"tmp3.geo" ]];
-        EndIf
-      EndIf
-    }
-  }
-}
-
-DefineConstant[
-  R_ = {"helmoltz_scalar", Name "GetDP/1ResolutionChoices", Visible 1},
-  C_ = {"-solve -pos -petsc_prealloc 100 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
-  P_ = {"postop_energy", Name "GetDP/2PostOperationChoices", Visible 1}
-];
-
-If(plotRTgraphs)
-  DefineConstant[
-    refl_  = {0, Name "GetDP/R0", ReadOnly 1, Graph "02000000", Visible 1},
-    abs_   = {0, Name "GetDP/total absorption", ReadOnly 1, Graph "00000002", Visible 1},
-    trans_ = {0, Name "GetDP/T0", ReadOnly 1, Graph "000000000002", Visible 1}
-  ];
+If (flag_polar<2)
+    Include "grating2D_scalar.pro";
+Else
+    Include "grating2D_conical.pro";
 EndIf
-- 
GitLab