diff --git a/ElectromagneticScattering/scattererTmatrix.geo b/ElectromagneticScattering/scattererTmatrix.geo
index 3b5cbce7b95bad4fa18bb909aeea638c80222efb..5172c55374ac6b29870f1ba290b37b6ce3bec923 100644
--- a/ElectromagneticScattering/scattererTmatrix.geo
+++ b/ElectromagneticScattering/scattererTmatrix.geo
@@ -6,12 +6,13 @@
 Include "scattererTmatrix_data.geo";
 SetFactory("OpenCASCADE");
 
-In_n =  Sqrt[Fabs[epsr_In_re]];
+In_n  =  Sqrt[Fabs[epsr_In_re]];
+Out_n =  Sqrt[Fabs[epsr_Out_re]];
   
 paramaille_pml = paramaille/1.1;
 
-Out_lc        = lambda_bg/paramaille;
-PML_lc        = lambda_bg/paramaille_pml;
+Out_lc        = lambda0/(paramaille*Out_n);
+PML_lc        = lambda0/(paramaille_pml*Out_n);
 In_lc         = lambda0/(paramaille*In_n*refine_scat);
 CenterScat_lc = lambda0/(paramaille*In_n);
 If(flag_shape==ELL)
diff --git a/ElectromagneticScattering/scattererTmatrix.pro b/ElectromagneticScattering/scattererTmatrix.pro
index 97da573a56ac14c905645cceb6e9eb3358daeca0..3e62be26e38ebfd3abe401a19bf6fc824de60373 100644
--- a/ElectromagneticScattering/scattererTmatrix.pro
+++ b/ElectromagneticScattering/scattererTmatrix.pro
@@ -54,11 +54,11 @@ Function{
 
   epsilonr[Scat_In]   = epsilonr_In[]  * TensorDiag[1.,1.,1.];
   epsilonr[Scat_Out]  = epsilonr_Out[] * TensorDiag[1.,1.,1.];
-  epsilonr[PMLs]      = pml_tens[];
+  epsilonr[PMLs]      = epsilonr_Out[] * pml_tens[];
                           
   epsilonr1[Scat_In]  = epsilonr_Out[] * TensorDiag[1.,1.,1.];
   epsilonr1[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];                
-  epsilonr1[PMLs]     = pml_tens[];
+  epsilonr1[PMLs]     = epsilonr_Out[] * pml_tens[];
                           
   mur[Scat_In]        = TensorDiag[1.,1.,1.];
   mur[Scat_Out]       = TensorDiag[1.,1.,1.];
@@ -204,8 +204,8 @@ PostProcessing {
         { Name Mnm_source~{pe}  ; Value { Local { [   Mnm_source~{pe}[] ]; In All_domains; Jacobian JVol; } } }
 
         For po In {1:p_max}
-          { Name feM~{pe}~{po} ; Value { Integral { [ (1/r_pml_in)^2 *normalize_fenm_Z~{po}[]*  {u}*Conj[Znm~{po}[]]  ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
-          { Name fhM~{pe}~{po} ; Value { Integral { [ (1/r_pml_in)^2 *normalize_fhnm_X~{po}[]*  {u}*Conj[Xnm~{po}[]]  ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
+          { Name feM~{pe}~{po} ; Value { Integral { [ 1/(r_pml_in^2*epsilonr_Out[]) * normalize_fenm_Z~{po}[]*  {u}*Conj[Znm~{po}[]]  ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
+          { Name fhM~{pe}~{po} ; Value { Integral { [ 1/(r_pml_in^2*epsilonr_Out[]) * normalize_fhnm_X~{po}[]*  {u}*Conj[Xnm~{po}[]]  ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
         EndFor
       }
     }
@@ -215,8 +215,8 @@ PostProcessing {
         { Name H_scat            ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
         { Name Nnm_source~{pe}   ; Value {    Local { [   Nnm_source~{pe}[] ]; In All_domains; Jacobian JVol; } } }
         For po In {1:p_max}
-          { Name feN~{pe}~{po}  ; Value { Integral { [ (1/r_pml_in)^2 *normalize_fenm_Z~{po}[]*  {u}*Conj[Znm~{po}[]]  ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
-          { Name fhN~{pe}~{po}  ; Value { Integral { [ (1/r_pml_in)^2 *normalize_fhnm_X~{po}[]*  {u}*Conj[Xnm~{po}[]]  ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
+          { Name feN~{pe}~{po}  ; Value { Integral { [ 1/(r_pml_in^2*epsilonr_Out[]) *normalize_fenm_Z~{po}[]*  {u}*Conj[Znm~{po}[]]  ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
+          { Name fhN~{pe}~{po}  ; Value { Integral { [ 1/(r_pml_in^2*epsilonr_Out[]) *normalize_fhnm_X~{po}[]*  {u}*Conj[Xnm~{po}[]]  ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
         EndFor
       }
     }
diff --git a/ElectromagneticScattering/scattering.pro b/ElectromagneticScattering/scattering.pro
index 789f6bcd253278f7cd422e20c736b11b33b9201a..2aa4078013b3257a76e9266e0d08fad89b889369 100644
--- a/ElectromagneticScattering/scattering.pro
+++ b/ElectromagneticScattering/scattering.pro
@@ -1,709 +1,709 @@
-///////////////////////////////
-// Author : Guillaume Demesy //
-// scattering.pro            //
-///////////////////////////////
-
-Include "scattering_data.geo";
-myDir = "run_results/";
-
-Group {
-  If(flag_cartpml==1)
-    PMLxyz         = Region[1000];
-    PMLxz          = Region[1001];
-    PMLyz          = Region[1002];
-    PMLxy          = Region[1003];
-    PMLz           = Region[1004];
-    PMLy           = Region[1005];
-    PMLx           = Region[1006];
-    Scat_In        = Region[1008];
-    Scat_Out       = Region[1007];
-    // Domains
-    Domain         = Region[{Scat_In,Scat_Out}];
-    PMLs           = Region[{PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}];
-    All_domains    = Region[{Scat_In,Scat_Out,PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}];
-    SurfInt        = Region[{}];
-  EndIf
-  If(flag_cartpml==0)
-    Scat_In        = Region[1];
-    Scat_Out       = Region[2];
-    Domain         = Region[{Scat_In,Scat_Out}];
-    PMLs           = Region[3];
-    All_domains    = Region[{Scat_In,Scat_Out,PMLs}];
-    SurfInt        = Region[{10}];
-  EndIf
-  PrintPoint = Region[2000];
-}
-
-Function{
-  // test={1.4,2.6,3.7};
-  // For i In {0:#test()-1}
-  //   Printf("%g",test(i));
-  // EndFor
-
-  I[]         = Complex[0,1];
-  avoid_sing  = 1.e-14;
-  mu0         = 4*Pi*100.0*nm;
-  epsilon0    = 8.854187817e-3*nm;
-  cel         = 1.0/Sqrt[epsilon0 * mu0];
-  Freq        = cel/lambda0;
-  omega0      = 2.0*Pi*Freq;
-  k_Out       = 2.0*Pi*Sqrt[epsr_Out_re]/lambda0;
-  r3D_sph[]   = Sqrt[X[]^2+Y[]^2+Z[]^2];
-  r2D_sph[]   = Sqrt[X[]^2+Y[]^2];
-  cos_theta[] = Z[]/(r3D_sph[]+avoid_sing);
-  theta[]     = Acos[cos_theta[]];
-  phi[]       = Atan2[-Y[],-X[]]+Pi;
-  
-  For kr In {0:nb_cuts-1}
-    r_sph~{kr}=r_sph_min+kr*(r_sph_max-r_sph_min)/(nb_cuts-1);
-  EndFor
-  
-  If (flag_study==0)
-    Ae     = 1.0;
-    Ah     = Ae*Sqrt[epsilon0*epsr_In_re/mu0];
-    alpha0 = k_Out*Sin[theta0]*Cos[phi0];
-    beta0  = k_Out*Sin[theta0]*Sin[phi0];
-    gamma0 = k_Out*Cos[theta0];
-    Ex0    =  Ae * Cos[psi0]*Cos[theta0]*Cos[phi0] - Ae* Sin[psi0]*Sin[phi0];
-    Ey0    =  Ae * Cos[psi0]*Cos[theta0]*Sin[phi0] + Ae* Sin[psi0]*Cos[phi0];
-    Ez0    = -Ae * Cos[psi0]*Sin[theta0];
-    Hx0    = -1/(omega0*mu0)*(beta0  * Ez0 - gamma0 * Ey0);
-    Hy0    = -1/(omega0*mu0)*(gamma0 * Ex0 - alpha0 * Ez0);
-    Hz0    = -1/(omega0*mu0)*(alpha0 * Ey0 - beta0  * Ex0);
-    Prop[] =  Ae * Complex[ Cos[alpha0*X[]+beta0*Y[]+gamma0*Z[]] , siwt*Sin[alpha0*X[]+beta0*Y[]+gamma0*Z[]] ];
-    Einc[] =  Vector[Ex0*Prop[],Ey0*Prop[],Ez0*Prop[]];
-    Hinc[] =  Vector[Hx0*Prop[],Hy0*Prop[],Hz0*Prop[]];
-    Pinc   =  0.5*Ae*Ae*Sqrt[epsilon0*epsr_In_re/mu0] * Cos[theta0];
-  EndIf
-  If (flag_study==2)
-    // Green Dyadic
-    normrr_p[]  = Sqrt[(X[]-x_p)^2+(Y[]-y_p)^2+(Z[]-z_p)^2];  
-    Gsca[]  = Complex[Cos[-1.*siwt*k_Out*normrr_p[]], 
-                      Sin[-1.*siwt*k_Out*normrr_p[]]]
-                 /(4.*Pi*normrr_p[]);
-    GDfact_diag[] = 1.+(I[]*k_Out*normrr_p[]-1.)/(k_Out*normrr_p[])^2;
-    GDfact_glob[] = (3.-3.*I[]*k_Out*normrr_p[]-(k_Out*normrr_p[])^2)/
-                      (k_Out*normrr_p[]^2)^2;
-    Dyad_Green2[]=
-      Gsca[]*(
-          GDfact_diag[]*TensorDiag[1,1,1]
-        + GDfact_glob[]*
-          Tensor[(X[]-x_p)*(X[]-x_p),(X[]-x_p)*(Y[]-y_p),(X[]-x_p)*(Z[]-z_p),
-                 (Y[]-y_p)*(X[]-x_p),(Y[]-y_p)*(Y[]-y_p),(Y[]-y_p)*(Z[]-z_p),
-                 (Z[]-z_p)*(X[]-x_p),(Z[]-z_p)*(Y[]-y_p),(Z[]-z_p)*(Z[]-z_p)]);
-    // Getdp Func
-    Dyad_Green[]     = DyadGreenHom[x_p,y_p,z_p,XYZ[],k_Out,siwt];
-    Einc_Green_0[]   = Dyad_Green2[]*Vector[1,0,0];
-    Einc_Green_1[]   = Dyad_Green2[]*Vector[0,1,0];
-    Einc_Green_2[]   = Dyad_Green2[]*Vector[0,0,1];
-    CurlDyad_Green[] = CurlDyadGreenHom[x_p,y_p,z_p,XYZ[],k_Out,siwt];
-    Hinc_Green_0[]   = siwt*I[]/(mu0*omega0)*CurlDyad_Green[]*Vector[1,0,0];
-    Hinc_Green_1[]   = siwt*I[]/(mu0*omega0)*CurlDyad_Green[]*Vector[0,1,0];
-    Hinc_Green_2[]   = siwt*I[]/(mu0*omega0)*CurlDyad_Green[]*Vector[0,0,1];
-    // test[] = DyadGreenHom[x_p,y_p,z_p,x_p+1e-1*nm,y_p+1e-1*nm,z_p+1e-1*nm,lambda0,Sqrt[epsr_Out_re]]*Vector[1,0,0];
-  EndIf
-
-  a_pml =  1.;
-  b_pml = -siwt;
-
-  eigtarget = (2.*Pi*cel/lambda0)^2;
-  eigfilter = 1e-4*Sqrt[eigtarget];
-  EigFilter[] = (Norm[$EigenvalueReal] > eigfilter);
-
-  If(flag_cartpml==1)
-    sx[Scat_In]     = 1.;
-    sy[Scat_In]     = 1.;
-    sz[Scat_In]     = 1.;
-    sx[Scat_Out]    = 1.;
-    sy[Scat_Out]    = 1.;
-    sz[Scat_Out]    = 1.;
-    sx[PMLxyz]      = Complex[a_pml,b_pml];
-    sy[PMLxyz]      = Complex[a_pml,b_pml];
-    sz[PMLxyz]      = Complex[a_pml,b_pml];
-    sx[PMLxz]       = Complex[a_pml,b_pml];
-    sy[PMLxz]       = 1.0;
-    sz[PMLxz]       = Complex[a_pml,b_pml];
-    sx[PMLyz]       = 1.0;
-    sy[PMLyz]       = Complex[a_pml,b_pml];
-    sz[PMLyz]       = Complex[a_pml,b_pml];
-    sx[PMLxy]       = Complex[a_pml,b_pml];
-    sy[PMLxy]       = Complex[a_pml,b_pml];
-    sz[PMLxy]       = 1.0;
-    sx[PMLx]        = Complex[a_pml,b_pml];
-    sy[PMLx]        = 1.0;
-    sz[PMLx]        = 1.0;
-    sx[PMLy]        = 1.0;
-    sy[PMLy]        = Complex[a_pml,b_pml];
-    sz[PMLy]        = 1.0;
-    sx[PMLz]        = 1.0;
-    sy[PMLz]        = 1.0;
-    sz[PMLz]        = Complex[a_pml,b_pml];
-    Lxx[]           = sy[]*sz[]/sx[];
-    Lyy[]           = sz[]*sx[]/sy[]; 
-    Lzz[]           = sx[]*sy[]/sz[];
-
-    epsilonr_In[]       = Complex[epsr_In_re  , epsr_In_im];
-    epsilonr_Out[]      = Complex[epsr_Out_re , epsr_Out_im];
-
-    epsilonr[Scat_In]   = epsilonr_In[]  * TensorDiag[1.,1.,1.];
-    epsilonr[Scat_Out]  = epsilonr_Out[] * TensorDiag[1.,1.,1.];
-    epsilonr[PMLs]      = epsr_Out_re * TensorDiag[Lxx[],Lyy[],Lzz[]];
-
-    epsilonr1[Scat_In]  = epsilonr_Out[] * TensorDiag[1.,1.,1.];
-    epsilonr1[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];                
-    epsilonr1[PMLs]     = epsr_Out_re * TensorDiag[Lxx[],Lyy[],Lzz[]];
-
-    mur[Scat_In]        = TensorDiag[1.,1.,1.];
-    mur[Scat_Out]       = TensorDiag[1.,1.,1.];
-    mur[PMLs]           = TensorDiag[Lxx[],Lyy[],Lzz[]];
-  EndIf
-
-  If(flag_cartpml==0)
-    sr[]         = Complex[a_pml,b_pml];
-    sphi[]       = Complex[1,0];
-    stheta[]     = Complex[1,0];
-    r_tild[]     = r_pml_in + (r3D_sph[] - r_pml_in) * sr[];
-    theta_tild[] = theta[];
-    pml_tens_temp1[]  = TensorDiag[(r_tild[]/r3D_sph[])^2 * sphi[]*stheta[]/sr[],
-                                   (r_tild[]/r3D_sph[])   * sr[]*stheta[]/sphi[],
-                                   (r_tild[]/r3D_sph[])   * sphi[]*sr[]/stheta[]];
-    pml_tens_temp2[]  = Rotate[pml_tens_temp1[],0,-theta[]-Pi/2,0];
-    pml_tens[]        = Rotate[pml_tens_temp2[],0,0,-phi[]];
-
-    epsilonr_In[]        = Complex[epsr_In_re  , epsr_In_im];
-    epsilonr_Out[]       = Complex[epsr_Out_re , epsr_Out_im];
-
-    epsilonr[Scat_In]   = epsilonr_In[]  * TensorDiag[1.,1.,1.];
-    epsilonr[Scat_Out]  = epsilonr_Out[] * TensorDiag[1.,1.,1.];
-    epsilonr[PMLs]      = pml_tens[];
-                           
-    epsilonr1[Scat_In]  = epsilonr_Out[] * TensorDiag[1.,1.,1.];
-    epsilonr1[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];                
-    epsilonr1[PMLs]     = pml_tens[];
-                           
-    mur[Scat_In]        = TensorDiag[1.,1.,1.];
-    mur[Scat_Out]       = TensorDiag[1.,1.,1.];
-    mur[PMLs]           = pml_tens[];
-  EndIf
-
-  If(flag_study==RES_PW)
-    source[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc[];
-  EndIf
-  If(flag_study==RES_TMAT)
-    For pe In {1:p_max}
-      ne = Floor[Sqrt[pe]];
-      me = ne*(ne+1) - Floor[pe];
-      Mnm_source~{pe}[] = Mnm[1,ne,me,XYZ[],k_Out];
-      Nnm_source~{pe}[] = Nnm[1,ne,me,XYZ[],k_Out];
-      Mnm_out~{pe}[]    = Mnm[4,ne,me,XYZ[],k_Out];
-      Nnm_out~{pe}[]    = Nnm[4,ne,me,XYZ[],k_Out];
-      source_M~{pe}[]   = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Mnm_source~{pe}[];
-      source_N~{pe}[]   = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Nnm_source~{pe}[];
-    EndFor
-  EndIf
-  
-  If (flag_study==RES_GREEN)
-    sourceG_0[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc_Green_0[];
-    sourceG_1[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc_Green_1[];
-    sourceG_2[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc_Green_2[];
-  EndIf  
-}
-
-Constraint {
-  // {Name Dirichlet; Type Assign;
-  //  Case {
-  //    { Region SurfDirichlet; Value 0.; }
-  //  }
-  // }
-}
-
-Jacobian {
-  { Name JVol ;
-    Case {
-      { Region All ; Jacobian Vol ; }
-    }
-  }
-  { Name JSur ;
-    Case {
-      { Region All ; Jacobian Sur ; }
-    }
-  }
-  { Name JLin ;
-    Case {
-      { Region All ; Jacobian Lin ; }
-    }
-  }
-}
-
-Integration {
-  { Name Int_1 ;
-    Case { 
-      { Type Gauss ;
-        Case { 
-          { GeoElement Point        ; NumberOfPoints   1 ; }
-          { GeoElement Line2        ; NumberOfPoints   4 ; }
-          { GeoElement Triangle2    ; NumberOfPoints  12 ; }
-          { GeoElement Tetrahedron2 ; NumberOfPoints  17 ; }
-        }
-      }
-    }
-  }
-}
-
-FunctionSpace {
-  { Name Hcurl; Type Form1;
-    BasisFunction {
-      { Name sn; NameOfCoef un; Function BF_Edge;
-        Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
-      { Name sn2; NameOfCoef un2; Function BF_Edge_2E;
-        Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
-        If (is_FEM_o2==1)
-          { Name sn3; NameOfCoef un3; Function BF_Edge_3F_b;
-            Support Region[{All_domains,SurfInt}]; Entity FacetsOf[All]; }
-          { Name sn4; NameOfCoef un4; Function BF_Edge_3F_c;
-            Support Region[{All_domains,SurfInt}]; Entity FacetsOf[All]; }
-          { Name sn5; NameOfCoef un5; Function BF_Edge_4E;
-            Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
-         EndIf
-    }
-    Constraint {
-      // { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint Dirichlet; }
-      // { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
-    }
-  }
-}
-
-Formulation {
-  If (flag_study==RES_QNM)
-    {Name QNM_helmholtz_vector; Type FemEquation;
-      Quantity {
-        { Name u; Type Local; NameOfSpace Hcurl;}
-      }
-      Equation {
-        Galerkin {        [      -1/mur[] * Dof{Curl u}       , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1;  }
-        Galerkin { DtDtDof[-1.  * 1/cel^2 * epsilonr[]*Dof{u} , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
-      }
-    }
-  EndIf
-
-  If (flag_study==RES_PW)
-    {Name PW_helmholtz_vector; Type FemEquation;
-      Quantity {
-        { Name u; Type Local; NameOfSpace Hcurl;}
-      }
-      Equation {
-        Galerkin { [-1/mur[]*Dof{Curl u}             , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1;  }
-        Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
-        Galerkin { [source[]                         , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
-      }
-    }
-  EndIf
-  If (flag_study==RES_TMAT)
-    For pe In {1:p_max}
-      {Name VPWM_helmholtz_vector~{pe}; Type FemEquation;
-        Quantity {
-          { Name u; Type Local; NameOfSpace Hcurl;}
-        }
-        Equation {
-          Galerkin { [-1/mur[]*Dof{Curl u}             , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1;  }
-          Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
-          Galerkin { [source_M~{pe}[]                  , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
-        }
-      }
-      {Name VPWN_helmholtz_vector~{pe}; Type FemEquation;
-        Quantity {
-          { Name u; Type Local; NameOfSpace Hcurl;}
-        }
-        Equation {
-          Galerkin { [-1/mur[]*Dof{Curl u}             , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1;  }
-          Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
-          Galerkin { [source_N~{pe}[]                  , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
-        }
-      }
-    EndFor
-    {Name VPWMN_helmholtz_vector; Type FemEquation;
-      Quantity {
-        { Name u; Type Local; NameOfSpace Hcurl;}
-      }
-      Equation {
-        Galerkin { [-1/mur[]*Dof{Curl u}             , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1;  }
-        Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
-        Galerkin { [ (omega0/cel)^2*(epsilonr[]-epsilonr1[])*
-                        ($isN ? Nnm[1,$NE,$ME,XYZ[],k_Out] : Mnm[1,$NE,$ME,XYZ[],k_Out]) 
-                                                     , {u} ]; In Scat_In; Jacobian JVol; Integration Int_1;}
-      }
-    }
-    {Name VPWN_helmholtz_vector_test; Type FemEquation;
-      Quantity {
-        { Name u; Type Local; NameOfSpace Hcurl;}
-      }
-      Equation {
-        Galerkin { [-1/mur[]*Dof{Curl u}             , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1;  }
-        Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
-        Galerkin { [ (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Nnm[1,$NE,$ME,XYZ[],k_Out], {u} ]; In Scat_In; Jacobian JVol; Integration Int_1;}
-      }
-    }
-
-  EndIf
-  If (flag_study==RES_GREEN)
-    For ncomp In {0:2}
-      {Name GreenAll_helmholtz_vector~{ncomp}; Type FemEquation;
-        Quantity {
-          { Name u; Type Local; NameOfSpace Hcurl;}
-        }
-        Equation {
-          Galerkin { [-1/mur[]*Dof{Curl u}             , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1;  }
-          Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
-          Galerkin { [sourceG~{ncomp}[]                , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
-        }
-      }
-    EndFor
-  EndIf
-}
-
-Resolution {
-  If (flag_study==RES_PW)
-    { Name res_PW_helmholtz_vector;
-      System {
-        { Name P; NameOfFormulation PW_helmholtz_vector; Type ComplexValue; Frequency Freq; }
-      }
-      Operation {
-        CreateDir[Str[myDir]];
-        // Evaluate[Python[]{"scattering_init.py"}];
-        Generate[P]; 
-        Solve[P];
-        PostOperation[PW_postop];
-        // Evaluate[Python[]{"scattering_post.py"}];
-      }
-    }
-  EndIf
-  If (flag_study==RES_TMAT)
-      { Name res_VPWall_helmholtz_vector;
-      System {
-        { Name T; NameOfFormulation VPWMN_helmholtz_vector; Type ComplexValue; }
-      }
-      Operation {
-        CreateDir[Str[myDir]];
-        // Evaluate[Python[]{"scattering_init.py"}];
-        
-        For pe In {1:p_max}
-          Evaluate[ $isN = 0 ];
-          Evaluate[ $PE  = pe ];
-          Evaluate[ $NE  = Floor[Sqrt[$PE]] ];
-          Evaluate[ $ME  = $NE*($NE+1) - Floor[$PE] ];
-          If (pe==1)
-            Generate[T];
-            Solve[T];
-          EndIf
-          GenerateRHS[T];
-          SolveAgain[T];
-          PostOperation[VPWM_postop~{pe}];
-          Evaluate[$isN=1];
-          GenerateRHS[T];
-          SolveAgain[T];
-          PostOperation[VPWN_postop~{pe}];
-        EndFor
-        // Evaluate[Python[]{"scattering_post.py"}];
-      }
-    }
-
-  EndIf
-  If (flag_study==RES_GREEN)
-    { Name res_GreenAll_helmholtz_vector;
-      System {
-        For ncomp In {0:2}
-          { Name M~{ncomp}; NameOfFormulation GreenAll_helmholtz_vector~{ncomp}; Type ComplexValue; Frequency Freq; }
-        EndFor
-      }
-      Operation {
-        CreateDir[Str[myDir]];
-        // Evaluate[Python[]{"scattering_init.py"}];
-        For ncomp In {0:2}
-          Generate[M~{ncomp}];
-          Solve[M~{ncomp}];
-          PostOperation[GreenAll_postop~{ncomp}];
-        EndFor
-        // Evaluate[Python[]{"scattering_post.py"}];
-      }
-    }    
-  EndIf
-  If (flag_study==RES_QNM)
-    { Name res_QNM_helmholtz_vector;
-      System{
-        { Name E; NameOfFormulation QNM_helmholtz_vector; Type ComplexValue; }
-      }
-      Operation{
-        CreateDir[Str[myDir]];
-        // Evaluate[Python[]{"scattering_init.py"}];
-        GenerateSeparate[E];
-        EigenSolve[E,neig,eigtarget,0,EigFilter[]];
-        PostOperation[QNM_postop];
-        // Evaluate[Python[]{"scattering_post.py"}];
-      }
-    }
-  EndIf
-
-}
-
-PostProcessing {
-  If (flag_study==RES_QNM)
-    { Name QNM_postpro; NameOfFormulation QNM_helmholtz_vector; NameOfSystem E;
-      Quantity {
-        { Name u ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
-        { Name EigenValuesReal;  Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JVol; } } }
-        { Name EigenValuesImag;  Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JVol; } } }
-      }
-    }
-  EndIf
-
-  If (flag_study==RES_PW)
-    { Name PW_postpro  ; NameOfFormulation PW_helmholtz_vector; NameOfSystem P;
-      Quantity {
-        { Name E_tot             ; Value { Local { [{u}+Einc[]]; In All_domains; Jacobian JVol; } } }
-        { Name normE_tot         ; Value { Local { [Norm[{u}+Einc[]]]; In All_domains; Jacobian JVol; } } }
-        { Name E_scat            ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
-        { Name E_scat_sph        ; Value { Local { [Vector[
-                                    (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
-                                    (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
-                                    (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
-                                   ]];
-                                In All_domains; Jacobian JVol; } } }                                
-        { Name H_scat        ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]         ; In All_domains; Jacobian JVol; } } }
-        { Name H_tot         ; Value { Local { [ siwt*I[]/(mur[]*mu0*omega0)*{Curl u} +Hinc[]]; In All_domains; Jacobian JVol; } } }
-        { Name normalized_losses ; Value { Integral { [ 0.5*omega0*epsilon0*Fabs[Im[epsilonr_In[]]]*(SquNorm[{u}+Einc[]])  ] ; In Scat_In ; Integration Int_1 ; Jacobian JVol ; } } }
-        { Name Poy_dif       ; Value { Local { [  0.5*Re[Cross[{u}        , Conj[       siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]]]  ]; In All_domains; Jacobian JVol; } } }
-        { Name Poy_tot       ; Value { Local { [  0.5*Re[Cross[{u}+Einc[] , Conj[Hinc[]+siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]]]  ]; In All_domains; Jacobian JVol; } } }
-      }
-    }
-  EndIf
-
-  If (flag_study==RES_TMAT)
-    For pe In {1:p_max}
-      { Name VPWM_postpro~{pe}; NameOfFormulation VPWMN_helmholtz_vector; NameOfSystem T;
-        Quantity {
-          { Name E_scat            ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
-          { Name E_scat_sph        ; Value { Local { [Vector[
-                                      (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
-                                      (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
-                                      (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
-                                     ]];
-                                  In All_domains; Jacobian JVol; } } }
-          { Name H_scat            ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
-          { Name Mnm_source~{pe}   ; Value {    Local { [   Mnm_source~{pe}[] ]; In All_domains; Jacobian JVol; } } }
-          For po In {1:p_max}
-            { Name a~{pe}~{po}     ; Value { Integral { [            {u}*Conj[Mnm_out~{po}[]] ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
-            { Name norma~{pe}~{po} ; Value { Integral { [ Mnm_out~{pe}[]*Conj[Mnm_out~{po}[]] ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
-          EndFor
-        }
-      }
-      { Name VPWN_postpro~{pe}; NameOfFormulation VPWMN_helmholtz_vector; NameOfSystem T;
-        Quantity {
-          { Name E_scat            ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
-          { Name E_scat_sph        ; Value { Local { [Vector[
-                                      (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
-                                      (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
-                                      (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
-                                     ]];
-                                  In All_domains; Jacobian JVol; } } }
-          { Name H_scat            ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
-          { Name Nnm_source~{pe}   ; Value {    Local { [   Nnm_source~{pe}[] ]; In All_domains; Jacobian JVol; } } }
-        }
-      }
-    EndFor
-  EndIf
-
-  If (flag_study==RES_GREEN)
-    For ncomp In {0:2}
-      { Name GreenAll_postpro~{ncomp}; NameOfFormulation GreenAll_helmholtz_vector~{ncomp};NameOfSystem M~{ncomp};
-        Quantity {
-          { Name E_tot~{ncomp}     ; Value { Local {[{u}+Einc_Green~{ncomp}[]]; In All_domains; Jacobian JVol; } } }
-          { Name H_tot~{ncomp}     ; Value { Local {[siwt*I[]/(mur[]*mu0*omega0)*{Curl u}+Hinc_Green~{ncomp}[]]; In All_domains; Jacobian JVol; } } }
-          { Name Einc_Green~{ncomp}; Value { Local {[ Einc_Green~{ncomp}[]   ]; In All_domains; Jacobian JVol; } } }
-          { Name Hinc_Green~{ncomp}; Value { Local {[ Hinc_Green~{ncomp}[]   ]; In All_domains; Jacobian JVol; } } }
-          { Name E_scat~{ncomp}    ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
-          { Name E_tot_im~{ncomp}     ; Value { Local {[Im[{u}+Einc_Green~{ncomp}[]]]; In All_domains; Jacobian JVol; } } }
-          { Name Einc_Green_im~{ncomp}; Value { Local {[Im[ Einc_Green~{ncomp}[]   ]]; In All_domains; Jacobian JVol; } } }
-          // { Name Einc_Green_im~{ncomp}; Value { Local {[test[]]; In All_domains; Jacobian JVol; } } }
-          { Name E_scat_im~{ncomp}    ; Value { Local {[Im[{u}                     ]]; In All_domains; Jacobian JVol; } } }
-          { Name E_tot_sph~{ncomp}; Value { Local { [Vector[
-                                      (X[]*CompX[{u}+Einc_Green~{ncomp}[]] 
-                                        + Y[]*CompY[{u}+Einc_Green~{ncomp}[]] 
-                                          + Z[]*CompZ[{u}+Einc_Green~{ncomp}[]])/r3D_sph[],
-                                      (X[]*Z[]*CompX[{u}+Einc_Green~{ncomp}[]] 
-                                        + Y[]*Z[]*CompY[{u}+Einc_Green~{ncomp}[]] 
-                                          - (X[]^2+Y[]^2)*CompZ[{u}+Einc_Green~{ncomp}[]])
-                                            /(r3D_sph[]*r2D_sph[]),
-                                      (-Y[]*CompX[{u}+Einc_Green~{ncomp}[]] 
-                                        + X[]*CompY[{u}+Einc_Green~{ncomp}[]])
-                                          /r2D_sph[] ]];
-                                  In All_domains; Jacobian JVol; } } }
-          { Name E_scat_sph~{ncomp}; Value { Local { [Vector[
-                                      (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
-                                      (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] 
-                                        - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
-                                      (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
-                                     ]];
-                                  In All_domains; Jacobian JVol; } } }
-        { Name H_scat~{ncomp} ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
-      }
-      }
-    EndFor
-  EndIf
-}
-
-
-
-PostOperation {
-  If (flag_study==RES_QNM)
-    { Name QNM_postop; NameOfPostProcessing QNM_postpro ;
-      Operation {
-        Print [EigenValuesReal, OnElementsOf PrintPoint , Format TimeTable, File StrCat[myDir,"EigenValuesReal.txt"]];
-        Print [EigenValuesImag, OnElementsOf PrintPoint , Format TimeTable, File StrCat[myDir,"EigenValuesImag.txt"]];
-        Print [ u             , OnElementsOf All_domains, File StrCat[myDir,"eigenVectors.pos"], EigenvalueLegend];
-      }
-    }
-  EndIf
-
-  If (flag_study==RES_PW)
-    {Name PW_postop; NameOfPostProcessing PW_postpro ;
-      Operation {
-        Print [ E_tot    , OnElementsOf Domain     , File StrCat[myDir,"E_tot.pos"]];
-        Print [ E_tot    , OnElementsOf Domain     , File StrCat[myDir,"E_tot.pos"]];
-        Print [ normE_tot, OnElementsOf Scat_In    , File StrCat[myDir,"normE_tot.pos"]];
-        Print [ H_scat   , OnElementsOf Domain     , File StrCat[myDir,"H_scat.pos"]];
-        Print [ H_tot    , OnElementsOf All_domains, File StrCat[myDir,"H_tot.pos"]];
-        Print [ Poy_dif  , OnElementsOf All_domains, File StrCat[myDir,"Poy_dif.pos"]];
-        Print [ Poy_tot  , OnElementsOf Domain     , File StrCat[myDir,"Poy_tot.pos"]];
-        For kr In {0:nb_cuts-1}
-          Print [ E_scat_sph , OnGrid 
-                {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]} 
-                {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)}, 
-                {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} }, 
-                File StrCat[myDir,StrCat["E_scat_onsphere_sph_PW",Sprintf["_r%g.dat",kr]]], Format Table];
-        EndFor
-      }
-    }
-  EndIf
-
-  If (flag_study==RES_TMAT)
-    For pe In {1:p_max}
-      {Name VPWM_postop~{pe}; NameOfPostProcessing VPWM_postpro~{pe} ;
-        Operation {
-          For kr In {0:nb_cuts-1}
-            Print [ E_scat_sph , OnGrid 
-                  {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]} 
-                  {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)}, 
-                  {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} }, 
-                  File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_sph_M",Sprintf["%g",pe]],Sprintf["_r%g.dat",kr]]], Format Table];
-          EndFor
-          Print [ E_scat , OnGrid
-                {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]} 
-                {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)}, 
-                {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} }, 
-                File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_cart_M",Sprintf["%g",pe]],Sprintf["_r%g.pos",kr]]],
-                Name StrCat["E_scat_onsphere_cart_M",Sprintf["%g",pe]]];
-          For po In {1:p_max}
-            Print[a~{pe}~{po}[SurfInt], OnGlobal  , Format Table, File StrCat[myDir,StrCat[StrCat["a_",Sprintf["pe%g",pe]],Sprintf["po%g.dat",po]]]];
-            Print[norma~{pe}~{po}[SurfInt], OnGlobal  , Format Table, File StrCat[myDir,StrCat[StrCat["norma_",Sprintf["pe%g",pe]],Sprintf["po%g.dat",po]]]];
-          EndFor
-        }
-      }
-      {Name VPWN_postop~{pe}; NameOfPostProcessing VPWN_postpro~{pe} ;
-        Operation {
-          For kr In {0:nb_cuts-1}
-            Print [ E_scat_sph , OnGrid 
-                  {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]} 
-                  {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)}, 
-                  {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} }, 
-                  File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_sph_N",Sprintf["%g",pe]],Sprintf["_r%g.dat",kr]]], Format Table];
-          EndFor
-          Print [ E_scat , OnGrid
-                {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]} 
-                {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)}, 
-                {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} }, 
-                File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_cart_N",Sprintf["%g",pe]],Sprintf["_r%g.pos",kr]]],
-                Name StrCat["E_scat_onsphere_cart_N",Sprintf["%g",pe]]];
-        }
-      }
-    EndFor
-  EndIf
-
-  If (flag_study==RES_GREEN)
-    For ncomp In {0:2}
-      {Name GreenAll_postop~{ncomp}; NameOfPostProcessing GreenAll_postpro~{ncomp} ;
-        Operation {
-          Print [ E_tot_sph~{ncomp} , OnGrid
-                {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]} 
-                {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)}, 
-                {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} }, 
-                File StrCat[myDir,StrCat["E_tot_onsphere_sph_G",Sprintf["%g.dat",ncomp]]] , Format Table];
-          Print [ E_tot~{ncomp} , OnGrid
-                {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]} 
-                {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)}, 
-                {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} }, 
-                File StrCat[myDir,StrCat["E_tot_onsphere_cart_G",Sprintf["%g.pos",ncomp]]]];          
-          Print [ E_scat_im~{ncomp} , OnPoint {x_p,y_p,z_p},
-                File StrCat[myDir,StrCat["E_scat_im_onpt_G",Sprintf["%g.dat",ncomp]]] , Format Table];
-          // Print [ Einc_Green_im~{ncomp} , OnPoint {x_p+0.1*nm,y_p+0.1*nm,z_p+0.1*nm},
-          //       File StrCat[myDir,StrCat["E_inc_im_onpt_G",Sprintf["%g.dat",ncomp]]] , Format Table];
-          // If(flag_FF==1)
-          //   Print [ E_tot~{ncomp}, OnElementsOf Scat_Out , File StrCat[myDir,StrCat["E_tot_",Sprintf["%g.pos",ncomp]]]];
-          //   Print [ H_tot~{ncomp}, OnElementsOf Scat_Out , File StrCat[myDir,StrCat["H_tot_",Sprintf["%g.pos",ncomp]]]];
-          //   Print [ E_tot~{ncomp} , OnGrid
-          //                  {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
-          //                  {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
-          //                  {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
-          //                  File StrCat[myDir,StrCat["E_tot_onsphere_G",Sprintf["%g.pos",ncomp]]]];
-          //   Print [ H_tot~{ncomp} , OnGrid
-          //                   {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
-          //                   {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
-          //                   {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
-          //                   File StrCat[myDir,StrCat["H_tot_onsphere_G",Sprintf["%g.pos",ncomp]]]];
-          //
-          //   Echo[
-          //    Str["Plugin(NearToFarField).Wavenumber        = 2*Pi/lambda0;",
-          //         "Plugin(NearToFarField).PhiStart         = 0;",
-          //         "Plugin(NearToFarField).PhiEnd           = 2*Pi;",
-          //         "Plugin(NearToFarField).NumPointsPhi     = npts_phi;",
-          //         "Plugin(NearToFarField).ThetaStart       = 0;",
-          //         "Plugin(NearToFarField).ThetaEnd         = Pi;",
-          //         "Plugin(NearToFarField).NumPointsTheta   = npts_theta;",
-          //         "Plugin(NearToFarField).EView            = PostProcessing.NbViews-2;",
-          //         "Plugin(NearToFarField).HView            = PostProcessing.NbViews-1;",
-          //         "Plugin(NearToFarField).Normalize        = 0;",
-          //         "Plugin(NearToFarField).dB               = 0;",
-          //         "Plugin(NearToFarField).NegativeTime     = 0;",
-          //         "Plugin(NearToFarField).RFar             = r_pml_in;",
-          //         "Plugin(NearToFarField).Run;"], File StrCat[myDir,"tmp1.geo"]] ;
-          // EndIf
-        }
-      }
-    EndFor
-  EndIf
-}
-
-If (flag_study==RES_QNM)
-  DefineConstant[
-    R_ = {"res_QNM_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
-    C_ = {"-solve -pos -petsc_prealloc 200 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
-    P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
-EndIf
-
-If (flag_study==RES_TMAT)
-  DefineConstant[
-    R_ = {"res_VPWall_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
-    C_ = {"-solve -pos -petsc_prealloc 200 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
-    P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
-EndIf
-If (flag_study==RES_PW)
-  DefineConstant[
-    R_ = {"res_PW_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
-    C_ = {"-solve -pos -petsc_prealloc 200 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
-    P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
-EndIf  
-If (flag_study==RES_GREEN)
-  DefineConstant[
-    R_ = {"res_GreenAll_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
-    C_ = {"-solve -pos -petsc_prealloc 200 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
-    P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
-EndIf
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattering.pro            //
+///////////////////////////////
+
+Include "scattering_data.geo";
+myDir = "run_results/";
+
+Group {
+  If(flag_cartpml==1)
+    PMLxyz         = Region[1000];
+    PMLxz          = Region[1001];
+    PMLyz          = Region[1002];
+    PMLxy          = Region[1003];
+    PMLz           = Region[1004];
+    PMLy           = Region[1005];
+    PMLx           = Region[1006];
+    Scat_In        = Region[1008];
+    Scat_Out       = Region[1007];
+    // Domains
+    Domain         = Region[{Scat_In,Scat_Out}];
+    PMLs           = Region[{PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}];
+    All_domains    = Region[{Scat_In,Scat_Out,PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}];
+    SurfInt        = Region[{}];
+  EndIf
+  If(flag_cartpml==0)
+    Scat_In        = Region[1];
+    Scat_Out       = Region[2];
+    Domain         = Region[{Scat_In,Scat_Out}];
+    PMLs           = Region[3];
+    All_domains    = Region[{Scat_In,Scat_Out,PMLs}];
+    SurfInt        = Region[{10}];
+  EndIf
+  PrintPoint = Region[2000];
+}
+
+Function{
+  // test={1.4,2.6,3.7};
+  // For i In {0:#test()-1}
+  //   Printf("%g",test(i));
+  // EndFor
+
+  I[]         = Complex[0,1];
+  avoid_sing  = 1.e-14;
+  mu0         = 4*Pi*100.0*nm;
+  epsilon0    = 8.854187817e-3*nm;
+  cel         = 1.0/Sqrt[epsilon0 * mu0];
+  Freq        = cel/lambda0;
+  omega0      = 2.0*Pi*Freq;
+  k_Out       = 2.0*Pi*Sqrt[epsr_Out_re]/lambda0;
+  r3D_sph[]   = Sqrt[X[]^2+Y[]^2+Z[]^2];
+  r2D_sph[]   = Sqrt[X[]^2+Y[]^2];
+  cos_theta[] = Z[]/(r3D_sph[]+avoid_sing);
+  theta[]     = Acos[cos_theta[]];
+  phi[]       = Atan2[-Y[],-X[]]+Pi;
+  
+  For kr In {0:nb_cuts-1}
+    r_sph~{kr}=r_sph_min+kr*(r_sph_max-r_sph_min)/(nb_cuts-1);
+  EndFor
+  
+  If (flag_study==0)
+    Ae     = 1.0;
+    Ah     = Ae*Sqrt[epsilon0*epsr_In_re/mu0];
+    alpha0 = k_Out*Sin[theta0]*Cos[phi0];
+    beta0  = k_Out*Sin[theta0]*Sin[phi0];
+    gamma0 = k_Out*Cos[theta0];
+    Ex0    =  Ae * Cos[psi0]*Cos[theta0]*Cos[phi0] - Ae* Sin[psi0]*Sin[phi0];
+    Ey0    =  Ae * Cos[psi0]*Cos[theta0]*Sin[phi0] + Ae* Sin[psi0]*Cos[phi0];
+    Ez0    = -Ae * Cos[psi0]*Sin[theta0];
+    Hx0    = -1/(omega0*mu0)*(beta0  * Ez0 - gamma0 * Ey0);
+    Hy0    = -1/(omega0*mu0)*(gamma0 * Ex0 - alpha0 * Ez0);
+    Hz0    = -1/(omega0*mu0)*(alpha0 * Ey0 - beta0  * Ex0);
+    Prop[] =  Ae * Complex[ Cos[alpha0*X[]+beta0*Y[]+gamma0*Z[]] , siwt*Sin[alpha0*X[]+beta0*Y[]+gamma0*Z[]] ];
+    Einc[] =  Vector[Ex0*Prop[],Ey0*Prop[],Ez0*Prop[]];
+    Hinc[] =  Vector[Hx0*Prop[],Hy0*Prop[],Hz0*Prop[]];
+    Pinc   =  0.5*Ae*Ae*Sqrt[epsilon0*epsr_In_re/mu0] * Cos[theta0];
+  EndIf
+  If (flag_study==2)
+    // Green Dyadic
+    normrr_p[]  = Sqrt[(X[]-x_p)^2+(Y[]-y_p)^2+(Z[]-z_p)^2];  
+    Gsca[]  = Complex[Cos[-1.*siwt*k_Out*normrr_p[]], 
+                      Sin[-1.*siwt*k_Out*normrr_p[]]]
+                 /(4.*Pi*normrr_p[]);
+    GDfact_diag[] = 1.+(I[]*k_Out*normrr_p[]-1.)/(k_Out*normrr_p[])^2;
+    GDfact_glob[] = (3.-3.*I[]*k_Out*normrr_p[]-(k_Out*normrr_p[])^2)/
+                      (k_Out*normrr_p[]^2)^2;
+    Dyad_Green2[]=
+      Gsca[]*(
+          GDfact_diag[]*TensorDiag[1,1,1]
+        + GDfact_glob[]*
+          Tensor[(X[]-x_p)*(X[]-x_p),(X[]-x_p)*(Y[]-y_p),(X[]-x_p)*(Z[]-z_p),
+                 (Y[]-y_p)*(X[]-x_p),(Y[]-y_p)*(Y[]-y_p),(Y[]-y_p)*(Z[]-z_p),
+                 (Z[]-z_p)*(X[]-x_p),(Z[]-z_p)*(Y[]-y_p),(Z[]-z_p)*(Z[]-z_p)]);
+    // Getdp Func
+    Dyad_Green[]     = DyadGreenHom[x_p,y_p,z_p,XYZ[],k_Out,siwt];
+    Einc_Green_0[]   = Dyad_Green2[]*Vector[1,0,0];
+    Einc_Green_1[]   = Dyad_Green2[]*Vector[0,1,0];
+    Einc_Green_2[]   = Dyad_Green2[]*Vector[0,0,1];
+    CurlDyad_Green[] = CurlDyadGreenHom[x_p,y_p,z_p,XYZ[],k_Out,siwt];
+    Hinc_Green_0[]   = siwt*I[]/(mu0*omega0)*CurlDyad_Green[]*Vector[1,0,0];
+    Hinc_Green_1[]   = siwt*I[]/(mu0*omega0)*CurlDyad_Green[]*Vector[0,1,0];
+    Hinc_Green_2[]   = siwt*I[]/(mu0*omega0)*CurlDyad_Green[]*Vector[0,0,1];
+    // test[] = DyadGreenHom[x_p,y_p,z_p,x_p+1e-1*nm,y_p+1e-1*nm,z_p+1e-1*nm,lambda0,Sqrt[epsr_Out_re]]*Vector[1,0,0];
+  EndIf
+
+  a_pml =  1.;
+  b_pml = -siwt;
+
+  eigtarget = (2.*Pi*cel/lambda0)^2;
+  eigfilter = 1e-4*Sqrt[eigtarget];
+  EigFilter[] = (Norm[$EigenvalueReal] > eigfilter);
+
+  If(flag_cartpml==1)
+    sx[Scat_In]     = 1.;
+    sy[Scat_In]     = 1.;
+    sz[Scat_In]     = 1.;
+    sx[Scat_Out]    = 1.;
+    sy[Scat_Out]    = 1.;
+    sz[Scat_Out]    = 1.;
+    sx[PMLxyz]      = Complex[a_pml,b_pml];
+    sy[PMLxyz]      = Complex[a_pml,b_pml];
+    sz[PMLxyz]      = Complex[a_pml,b_pml];
+    sx[PMLxz]       = Complex[a_pml,b_pml];
+    sy[PMLxz]       = 1.0;
+    sz[PMLxz]       = Complex[a_pml,b_pml];
+    sx[PMLyz]       = 1.0;
+    sy[PMLyz]       = Complex[a_pml,b_pml];
+    sz[PMLyz]       = Complex[a_pml,b_pml];
+    sx[PMLxy]       = Complex[a_pml,b_pml];
+    sy[PMLxy]       = Complex[a_pml,b_pml];
+    sz[PMLxy]       = 1.0;
+    sx[PMLx]        = Complex[a_pml,b_pml];
+    sy[PMLx]        = 1.0;
+    sz[PMLx]        = 1.0;
+    sx[PMLy]        = 1.0;
+    sy[PMLy]        = Complex[a_pml,b_pml];
+    sz[PMLy]        = 1.0;
+    sx[PMLz]        = 1.0;
+    sy[PMLz]        = 1.0;
+    sz[PMLz]        = Complex[a_pml,b_pml];
+    Lxx[]           = sy[]*sz[]/sx[];
+    Lyy[]           = sz[]*sx[]/sy[]; 
+    Lzz[]           = sx[]*sy[]/sz[];
+
+    epsilonr_In[]       = Complex[epsr_In_re  , epsr_In_im];
+    epsilonr_Out[]      = Complex[epsr_Out_re , epsr_Out_im];
+
+    epsilonr[Scat_In]   = epsilonr_In[]  * TensorDiag[1.,1.,1.];
+    epsilonr[Scat_Out]  = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+    epsilonr[PMLs]      = epsr_Out_re * TensorDiag[Lxx[],Lyy[],Lzz[]];
+
+    epsilonr1[Scat_In]  = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+    epsilonr1[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];                
+    epsilonr1[PMLs]     = epsr_Out_re * TensorDiag[Lxx[],Lyy[],Lzz[]];
+
+    mur[Scat_In]        = TensorDiag[1.,1.,1.];
+    mur[Scat_Out]       = TensorDiag[1.,1.,1.];
+    mur[PMLs]           = TensorDiag[Lxx[],Lyy[],Lzz[]];
+  EndIf
+
+  If(flag_cartpml==0)
+    sr[]         = Complex[a_pml,b_pml];
+    sphi[]       = Complex[1,0];
+    stheta[]     = Complex[1,0];
+    r_tild[]     = r_pml_in + (r3D_sph[] - r_pml_in) * sr[];
+    theta_tild[] = theta[];
+    pml_tens_temp1[]  = TensorDiag[(r_tild[]/r3D_sph[])^2 * sphi[]*stheta[]/sr[],
+                                   (r_tild[]/r3D_sph[])   * sr[]*stheta[]/sphi[],
+                                   (r_tild[]/r3D_sph[])   * sphi[]*sr[]/stheta[]];
+    pml_tens_temp2[]  = Rotate[pml_tens_temp1[],0,-theta[]-Pi/2,0];
+    pml_tens[]        = Rotate[pml_tens_temp2[],0,0,-phi[]];
+
+    epsilonr_In[]        = Complex[epsr_In_re  , epsr_In_im];
+    epsilonr_Out[]       = Complex[epsr_Out_re , epsr_Out_im];
+
+    epsilonr[Scat_In]   = epsilonr_In[]  * TensorDiag[1.,1.,1.];
+    epsilonr[Scat_Out]  = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+    epsilonr[PMLs]      = epsilonr_Out[] * pml_tens[];
+                           
+    epsilonr1[Scat_In]  = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+    epsilonr1[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];                
+    epsilonr1[PMLs]     = epsilonr_Out[] * pml_tens[];
+                           
+    mur[Scat_In]        = TensorDiag[1.,1.,1.];
+    mur[Scat_Out]       = TensorDiag[1.,1.,1.];
+    mur[PMLs]           = pml_tens[];
+  EndIf
+
+  If(flag_study==RES_PW)
+    source[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc[];
+  EndIf
+  If(flag_study==RES_TMAT)
+    For pe In {1:p_max}
+      ne = Floor[Sqrt[pe]];
+      me = ne*(ne+1) - Floor[pe];
+      Mnm_source~{pe}[] = Mnm[1,ne,me,XYZ[],k_Out];
+      Nnm_source~{pe}[] = Nnm[1,ne,me,XYZ[],k_Out];
+      Mnm_out~{pe}[]    = Mnm[4,ne,me,XYZ[],k_Out];
+      Nnm_out~{pe}[]    = Nnm[4,ne,me,XYZ[],k_Out];
+      source_M~{pe}[]   = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Mnm_source~{pe}[];
+      source_N~{pe}[]   = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Nnm_source~{pe}[];
+    EndFor
+  EndIf
+  
+  If (flag_study==RES_GREEN)
+    sourceG_0[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc_Green_0[];
+    sourceG_1[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc_Green_1[];
+    sourceG_2[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc_Green_2[];
+  EndIf  
+}
+
+Constraint {
+  // {Name Dirichlet; Type Assign;
+  //  Case {
+  //    { Region SurfDirichlet; Value 0.; }
+  //  }
+  // }
+}
+
+Jacobian {
+  { Name JVol ;
+    Case {
+      { Region All ; Jacobian Vol ; }
+    }
+  }
+  { Name JSur ;
+    Case {
+      { Region All ; Jacobian Sur ; }
+    }
+  }
+  { Name JLin ;
+    Case {
+      { Region All ; Jacobian Lin ; }
+    }
+  }
+}
+
+Integration {
+  { Name Int_1 ;
+    Case { 
+      { Type Gauss ;
+        Case { 
+          { GeoElement Point        ; NumberOfPoints   1 ; }
+          { GeoElement Line2        ; NumberOfPoints   4 ; }
+          { GeoElement Triangle2    ; NumberOfPoints  12 ; }
+          { GeoElement Tetrahedron2 ; NumberOfPoints  17 ; }
+        }
+      }
+    }
+  }
+}
+
+FunctionSpace {
+  { Name Hcurl; Type Form1;
+    BasisFunction {
+      { Name sn; NameOfCoef un; Function BF_Edge;
+        Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
+      { Name sn2; NameOfCoef un2; Function BF_Edge_2E;
+        Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
+        If (is_FEM_o2==1)
+          { Name sn3; NameOfCoef un3; Function BF_Edge_3F_b;
+            Support Region[{All_domains,SurfInt}]; Entity FacetsOf[All]; }
+          { Name sn4; NameOfCoef un4; Function BF_Edge_3F_c;
+            Support Region[{All_domains,SurfInt}]; Entity FacetsOf[All]; }
+          { Name sn5; NameOfCoef un5; Function BF_Edge_4E;
+            Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
+         EndIf
+    }
+    Constraint {
+      // { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+      // { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+    }
+  }
+}
+
+Formulation {
+  If (flag_study==RES_QNM)
+    {Name QNM_helmholtz_vector; Type FemEquation;
+      Quantity {
+        { Name u; Type Local; NameOfSpace Hcurl;}
+      }
+      Equation {
+        Galerkin {        [      -1/mur[] * Dof{Curl u}       , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1;  }
+        Galerkin { DtDtDof[-1.  * 1/cel^2 * epsilonr[]*Dof{u} , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
+      }
+    }
+  EndIf
+
+  If (flag_study==RES_PW)
+    {Name PW_helmholtz_vector; Type FemEquation;
+      Quantity {
+        { Name u; Type Local; NameOfSpace Hcurl;}
+      }
+      Equation {
+        Galerkin { [-1/mur[]*Dof{Curl u}             , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1;  }
+        Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
+        Galerkin { [source[]                         , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
+      }
+    }
+  EndIf
+  If (flag_study==RES_TMAT)
+    For pe In {1:p_max}
+      {Name VPWM_helmholtz_vector~{pe}; Type FemEquation;
+        Quantity {
+          { Name u; Type Local; NameOfSpace Hcurl;}
+        }
+        Equation {
+          Galerkin { [-1/mur[]*Dof{Curl u}             , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1;  }
+          Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
+          Galerkin { [source_M~{pe}[]                  , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
+        }
+      }
+      {Name VPWN_helmholtz_vector~{pe}; Type FemEquation;
+        Quantity {
+          { Name u; Type Local; NameOfSpace Hcurl;}
+        }
+        Equation {
+          Galerkin { [-1/mur[]*Dof{Curl u}             , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1;  }
+          Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
+          Galerkin { [source_N~{pe}[]                  , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
+        }
+      }
+    EndFor
+    {Name VPWMN_helmholtz_vector; Type FemEquation;
+      Quantity {
+        { Name u; Type Local; NameOfSpace Hcurl;}
+      }
+      Equation {
+        Galerkin { [-1/mur[]*Dof{Curl u}             , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1;  }
+        Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
+        Galerkin { [ (omega0/cel)^2*(epsilonr[]-epsilonr1[])*
+                        ($isN ? Nnm[1,$NE,$ME,XYZ[],k_Out] : Mnm[1,$NE,$ME,XYZ[],k_Out]) 
+                                                     , {u} ]; In Scat_In; Jacobian JVol; Integration Int_1;}
+      }
+    }
+    {Name VPWN_helmholtz_vector_test; Type FemEquation;
+      Quantity {
+        { Name u; Type Local; NameOfSpace Hcurl;}
+      }
+      Equation {
+        Galerkin { [-1/mur[]*Dof{Curl u}             , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1;  }
+        Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
+        Galerkin { [ (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Nnm[1,$NE,$ME,XYZ[],k_Out], {u} ]; In Scat_In; Jacobian JVol; Integration Int_1;}
+      }
+    }
+
+  EndIf
+  If (flag_study==RES_GREEN)
+    For ncomp In {0:2}
+      {Name GreenAll_helmholtz_vector~{ncomp}; Type FemEquation;
+        Quantity {
+          { Name u; Type Local; NameOfSpace Hcurl;}
+        }
+        Equation {
+          Galerkin { [-1/mur[]*Dof{Curl u}             , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1;  }
+          Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
+          Galerkin { [sourceG~{ncomp}[]                , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
+        }
+      }
+    EndFor
+  EndIf
+}
+
+Resolution {
+  If (flag_study==RES_PW)
+    { Name res_PW_helmholtz_vector;
+      System {
+        { Name P; NameOfFormulation PW_helmholtz_vector; Type ComplexValue; Frequency Freq; }
+      }
+      Operation {
+        CreateDir[Str[myDir]];
+        // Evaluate[Python[]{"scattering_init.py"}];
+        Generate[P]; 
+        Solve[P];
+        PostOperation[PW_postop];
+        // Evaluate[Python[]{"scattering_post.py"}];
+      }
+    }
+  EndIf
+  If (flag_study==RES_TMAT)
+      { Name res_VPWall_helmholtz_vector;
+      System {
+        { Name T; NameOfFormulation VPWMN_helmholtz_vector; Type ComplexValue; }
+      }
+      Operation {
+        CreateDir[Str[myDir]];
+        // Evaluate[Python[]{"scattering_init.py"}];
+        
+        For pe In {1:p_max}
+          Evaluate[ $isN = 0 ];
+          Evaluate[ $PE  = pe ];
+          Evaluate[ $NE  = Floor[Sqrt[$PE]] ];
+          Evaluate[ $ME  = $NE*($NE+1) - Floor[$PE] ];
+          If (pe==1)
+            Generate[T];
+            Solve[T];
+          EndIf
+          GenerateRHS[T];
+          SolveAgain[T];
+          PostOperation[VPWM_postop~{pe}];
+          Evaluate[$isN=1];
+          GenerateRHS[T];
+          SolveAgain[T];
+          PostOperation[VPWN_postop~{pe}];
+        EndFor
+        // Evaluate[Python[]{"scattering_post.py"}];
+      }
+    }
+
+  EndIf
+  If (flag_study==RES_GREEN)
+    { Name res_GreenAll_helmholtz_vector;
+      System {
+        For ncomp In {0:2}
+          { Name M~{ncomp}; NameOfFormulation GreenAll_helmholtz_vector~{ncomp}; Type ComplexValue; Frequency Freq; }
+        EndFor
+      }
+      Operation {
+        CreateDir[Str[myDir]];
+        // Evaluate[Python[]{"scattering_init.py"}];
+        For ncomp In {0:2}
+          Generate[M~{ncomp}];
+          Solve[M~{ncomp}];
+          PostOperation[GreenAll_postop~{ncomp}];
+        EndFor
+        // Evaluate[Python[]{"scattering_post.py"}];
+      }
+    }    
+  EndIf
+  If (flag_study==RES_QNM)
+    { Name res_QNM_helmholtz_vector;
+      System{
+        { Name E; NameOfFormulation QNM_helmholtz_vector; Type ComplexValue; }
+      }
+      Operation{
+        CreateDir[Str[myDir]];
+        // Evaluate[Python[]{"scattering_init.py"}];
+        GenerateSeparate[E];
+        EigenSolve[E,neig,eigtarget,0,EigFilter[]];
+        PostOperation[QNM_postop];
+        // Evaluate[Python[]{"scattering_post.py"}];
+      }
+    }
+  EndIf
+
+}
+
+PostProcessing {
+  If (flag_study==RES_QNM)
+    { Name QNM_postpro; NameOfFormulation QNM_helmholtz_vector; NameOfSystem E;
+      Quantity {
+        { Name u ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+        { Name EigenValuesReal;  Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JVol; } } }
+        { Name EigenValuesImag;  Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JVol; } } }
+      }
+    }
+  EndIf
+
+  If (flag_study==RES_PW)
+    { Name PW_postpro  ; NameOfFormulation PW_helmholtz_vector; NameOfSystem P;
+      Quantity {
+        { Name E_tot             ; Value { Local { [{u}+Einc[]]; In All_domains; Jacobian JVol; } } }
+        { Name normE_tot         ; Value { Local { [Norm[{u}+Einc[]]]; In All_domains; Jacobian JVol; } } }
+        { Name E_scat            ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+        { Name E_scat_sph        ; Value { Local { [Vector[
+                                    (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+                                    (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+                                    (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+                                   ]];
+                                In All_domains; Jacobian JVol; } } }                                
+        { Name H_scat        ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]         ; In All_domains; Jacobian JVol; } } }
+        { Name H_tot         ; Value { Local { [ siwt*I[]/(mur[]*mu0*omega0)*{Curl u} +Hinc[]]; In All_domains; Jacobian JVol; } } }
+        { Name normalized_losses ; Value { Integral { [ 0.5*omega0*epsilon0*Fabs[Im[epsilonr_In[]]]*(SquNorm[{u}+Einc[]])  ] ; In Scat_In ; Integration Int_1 ; Jacobian JVol ; } } }
+        { Name Poy_dif       ; Value { Local { [  0.5*Re[Cross[{u}        , Conj[       siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]]]  ]; In All_domains; Jacobian JVol; } } }
+        { Name Poy_tot       ; Value { Local { [  0.5*Re[Cross[{u}+Einc[] , Conj[Hinc[]+siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]]]  ]; In All_domains; Jacobian JVol; } } }
+      }
+    }
+  EndIf
+
+  If (flag_study==RES_TMAT)
+    For pe In {1:p_max}
+      { Name VPWM_postpro~{pe}; NameOfFormulation VPWMN_helmholtz_vector; NameOfSystem T;
+        Quantity {
+          { Name E_scat            ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+          { Name E_scat_sph        ; Value { Local { [Vector[
+                                      (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+                                      (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+                                      (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+                                     ]];
+                                  In All_domains; Jacobian JVol; } } }
+          { Name H_scat            ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
+          { Name Mnm_source~{pe}   ; Value {    Local { [   Mnm_source~{pe}[] ]; In All_domains; Jacobian JVol; } } }
+          For po In {1:p_max}
+            { Name a~{pe}~{po}     ; Value { Integral { [            {u}*Conj[Mnm_out~{po}[]] ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
+            { Name norma~{pe}~{po} ; Value { Integral { [ Mnm_out~{pe}[]*Conj[Mnm_out~{po}[]] ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
+          EndFor
+        }
+      }
+      { Name VPWN_postpro~{pe}; NameOfFormulation VPWMN_helmholtz_vector; NameOfSystem T;
+        Quantity {
+          { Name E_scat            ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+          { Name E_scat_sph        ; Value { Local { [Vector[
+                                      (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+                                      (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+                                      (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+                                     ]];
+                                  In All_domains; Jacobian JVol; } } }
+          { Name H_scat            ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
+          { Name Nnm_source~{pe}   ; Value {    Local { [   Nnm_source~{pe}[] ]; In All_domains; Jacobian JVol; } } }
+        }
+      }
+    EndFor
+  EndIf
+
+  If (flag_study==RES_GREEN)
+    For ncomp In {0:2}
+      { Name GreenAll_postpro~{ncomp}; NameOfFormulation GreenAll_helmholtz_vector~{ncomp};NameOfSystem M~{ncomp};
+        Quantity {
+          { Name E_tot~{ncomp}     ; Value { Local {[{u}+Einc_Green~{ncomp}[]]; In All_domains; Jacobian JVol; } } }
+          { Name H_tot~{ncomp}     ; Value { Local {[siwt*I[]/(mur[]*mu0*omega0)*{Curl u}+Hinc_Green~{ncomp}[]]; In All_domains; Jacobian JVol; } } }
+          { Name Einc_Green~{ncomp}; Value { Local {[ Einc_Green~{ncomp}[]   ]; In All_domains; Jacobian JVol; } } }
+          { Name Hinc_Green~{ncomp}; Value { Local {[ Hinc_Green~{ncomp}[]   ]; In All_domains; Jacobian JVol; } } }
+          { Name E_scat~{ncomp}    ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+          { Name E_tot_im~{ncomp}     ; Value { Local {[Im[{u}+Einc_Green~{ncomp}[]]]; In All_domains; Jacobian JVol; } } }
+          { Name Einc_Green_im~{ncomp}; Value { Local {[Im[ Einc_Green~{ncomp}[]   ]]; In All_domains; Jacobian JVol; } } }
+          // { Name Einc_Green_im~{ncomp}; Value { Local {[test[]]; In All_domains; Jacobian JVol; } } }
+          { Name E_scat_im~{ncomp}    ; Value { Local {[Im[{u}                     ]]; In All_domains; Jacobian JVol; } } }
+          { Name E_tot_sph~{ncomp}; Value { Local { [Vector[
+                                      (X[]*CompX[{u}+Einc_Green~{ncomp}[]] 
+                                        + Y[]*CompY[{u}+Einc_Green~{ncomp}[]] 
+                                          + Z[]*CompZ[{u}+Einc_Green~{ncomp}[]])/r3D_sph[],
+                                      (X[]*Z[]*CompX[{u}+Einc_Green~{ncomp}[]] 
+                                        + Y[]*Z[]*CompY[{u}+Einc_Green~{ncomp}[]] 
+                                          - (X[]^2+Y[]^2)*CompZ[{u}+Einc_Green~{ncomp}[]])
+                                            /(r3D_sph[]*r2D_sph[]),
+                                      (-Y[]*CompX[{u}+Einc_Green~{ncomp}[]] 
+                                        + X[]*CompY[{u}+Einc_Green~{ncomp}[]])
+                                          /r2D_sph[] ]];
+                                  In All_domains; Jacobian JVol; } } }
+          { Name E_scat_sph~{ncomp}; Value { Local { [Vector[
+                                      (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+                                      (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] 
+                                        - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+                                      (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+                                     ]];
+                                  In All_domains; Jacobian JVol; } } }
+        { Name H_scat~{ncomp} ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
+      }
+      }
+    EndFor
+  EndIf
+}
+
+
+
+PostOperation {
+  If (flag_study==RES_QNM)
+    { Name QNM_postop; NameOfPostProcessing QNM_postpro ;
+      Operation {
+        Print [EigenValuesReal, OnElementsOf PrintPoint , Format TimeTable, File StrCat[myDir,"EigenValuesReal.txt"]];
+        Print [EigenValuesImag, OnElementsOf PrintPoint , Format TimeTable, File StrCat[myDir,"EigenValuesImag.txt"]];
+        Print [ u             , OnElementsOf All_domains, File StrCat[myDir,"eigenVectors.pos"], EigenvalueLegend];
+      }
+    }
+  EndIf
+
+  If (flag_study==RES_PW)
+    {Name PW_postop; NameOfPostProcessing PW_postpro ;
+      Operation {
+        Print [ E_tot    , OnElementsOf Domain     , File StrCat[myDir,"E_tot.pos"]];
+        Print [ E_tot    , OnElementsOf Domain     , File StrCat[myDir,"E_tot.pos"]];
+        Print [ normE_tot, OnElementsOf Scat_In    , File StrCat[myDir,"normE_tot.pos"]];
+        Print [ H_scat   , OnElementsOf Domain     , File StrCat[myDir,"H_scat.pos"]];
+        Print [ H_tot    , OnElementsOf All_domains, File StrCat[myDir,"H_tot.pos"]];
+        Print [ Poy_dif  , OnElementsOf All_domains, File StrCat[myDir,"Poy_dif.pos"]];
+        Print [ Poy_tot  , OnElementsOf Domain     , File StrCat[myDir,"Poy_tot.pos"]];
+        For kr In {0:nb_cuts-1}
+          Print [ E_scat_sph , OnGrid 
+                {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]} 
+                {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)}, 
+                {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} }, 
+                File StrCat[myDir,StrCat["E_scat_onsphere_sph_PW",Sprintf["_r%g.dat",kr]]], Format Table];
+        EndFor
+      }
+    }
+  EndIf
+
+  If (flag_study==RES_TMAT)
+    For pe In {1:p_max}
+      {Name VPWM_postop~{pe}; NameOfPostProcessing VPWM_postpro~{pe} ;
+        Operation {
+          For kr In {0:nb_cuts-1}
+            Print [ E_scat_sph , OnGrid 
+                  {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]} 
+                  {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)}, 
+                  {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} }, 
+                  File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_sph_M",Sprintf["%g",pe]],Sprintf["_r%g.dat",kr]]], Format Table];
+          EndFor
+          Print [ E_scat , OnGrid
+                {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]} 
+                {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)}, 
+                {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} }, 
+                File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_cart_M",Sprintf["%g",pe]],Sprintf["_r%g.pos",kr]]],
+                Name StrCat["E_scat_onsphere_cart_M",Sprintf["%g",pe]]];
+          For po In {1:p_max}
+            Print[a~{pe}~{po}[SurfInt], OnGlobal  , Format Table, File StrCat[myDir,StrCat[StrCat["a_",Sprintf["pe%g",pe]],Sprintf["po%g.dat",po]]]];
+            Print[norma~{pe}~{po}[SurfInt], OnGlobal  , Format Table, File StrCat[myDir,StrCat[StrCat["norma_",Sprintf["pe%g",pe]],Sprintf["po%g.dat",po]]]];
+          EndFor
+        }
+      }
+      {Name VPWN_postop~{pe}; NameOfPostProcessing VPWN_postpro~{pe} ;
+        Operation {
+          For kr In {0:nb_cuts-1}
+            Print [ E_scat_sph , OnGrid 
+                  {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]} 
+                  {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)}, 
+                  {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} }, 
+                  File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_sph_N",Sprintf["%g",pe]],Sprintf["_r%g.dat",kr]]], Format Table];
+          EndFor
+          Print [ E_scat , OnGrid
+                {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]} 
+                {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)}, 
+                {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} }, 
+                File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_cart_N",Sprintf["%g",pe]],Sprintf["_r%g.pos",kr]]],
+                Name StrCat["E_scat_onsphere_cart_N",Sprintf["%g",pe]]];
+        }
+      }
+    EndFor
+  EndIf
+
+  If (flag_study==RES_GREEN)
+    For ncomp In {0:2}
+      {Name GreenAll_postop~{ncomp}; NameOfPostProcessing GreenAll_postpro~{ncomp} ;
+        Operation {
+          Print [ E_tot_sph~{ncomp} , OnGrid
+                {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]} 
+                {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)}, 
+                {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} }, 
+                File StrCat[myDir,StrCat["E_tot_onsphere_sph_G",Sprintf["%g.dat",ncomp]]] , Format Table];
+          Print [ E_tot~{ncomp} , OnGrid
+                {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]} 
+                {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)}, 
+                {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} }, 
+                File StrCat[myDir,StrCat["E_tot_onsphere_cart_G",Sprintf["%g.pos",ncomp]]]];          
+          Print [ E_scat_im~{ncomp} , OnPoint {x_p,y_p,z_p},
+                File StrCat[myDir,StrCat["E_scat_im_onpt_G",Sprintf["%g.dat",ncomp]]] , Format Table];
+          // Print [ Einc_Green_im~{ncomp} , OnPoint {x_p+0.1*nm,y_p+0.1*nm,z_p+0.1*nm},
+          //       File StrCat[myDir,StrCat["E_inc_im_onpt_G",Sprintf["%g.dat",ncomp]]] , Format Table];
+          // If(flag_FF==1)
+          //   Print [ E_tot~{ncomp}, OnElementsOf Scat_Out , File StrCat[myDir,StrCat["E_tot_",Sprintf["%g.pos",ncomp]]]];
+          //   Print [ H_tot~{ncomp}, OnElementsOf Scat_Out , File StrCat[myDir,StrCat["H_tot_",Sprintf["%g.pos",ncomp]]]];
+          //   Print [ E_tot~{ncomp} , OnGrid
+          //                  {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
+          //                  {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+          //                  {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+          //                  File StrCat[myDir,StrCat["E_tot_onsphere_G",Sprintf["%g.pos",ncomp]]]];
+          //   Print [ H_tot~{ncomp} , OnGrid
+          //                   {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
+          //                   {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+          //                   {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+          //                   File StrCat[myDir,StrCat["H_tot_onsphere_G",Sprintf["%g.pos",ncomp]]]];
+          //
+          //   Echo[
+          //    Str["Plugin(NearToFarField).Wavenumber        = 2*Pi/lambda0;",
+          //         "Plugin(NearToFarField).PhiStart         = 0;",
+          //         "Plugin(NearToFarField).PhiEnd           = 2*Pi;",
+          //         "Plugin(NearToFarField).NumPointsPhi     = npts_phi;",
+          //         "Plugin(NearToFarField).ThetaStart       = 0;",
+          //         "Plugin(NearToFarField).ThetaEnd         = Pi;",
+          //         "Plugin(NearToFarField).NumPointsTheta   = npts_theta;",
+          //         "Plugin(NearToFarField).EView            = PostProcessing.NbViews-2;",
+          //         "Plugin(NearToFarField).HView            = PostProcessing.NbViews-1;",
+          //         "Plugin(NearToFarField).Normalize        = 0;",
+          //         "Plugin(NearToFarField).dB               = 0;",
+          //         "Plugin(NearToFarField).NegativeTime     = 0;",
+          //         "Plugin(NearToFarField).RFar             = r_pml_in;",
+          //         "Plugin(NearToFarField).Run;"], File StrCat[myDir,"tmp1.geo"]] ;
+          // EndIf
+        }
+      }
+    EndFor
+  EndIf
+}
+
+If (flag_study==RES_QNM)
+  DefineConstant[
+    R_ = {"res_QNM_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+    C_ = {"-solve -pos -petsc_prealloc 200 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
+    P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf
+
+If (flag_study==RES_TMAT)
+  DefineConstant[
+    R_ = {"res_VPWall_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+    C_ = {"-solve -pos -petsc_prealloc 200 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
+    P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf
+If (flag_study==RES_PW)
+  DefineConstant[
+    R_ = {"res_PW_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+    C_ = {"-solve -pos -petsc_prealloc 200 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
+    P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf  
+If (flag_study==RES_GREEN)
+  DefineConstant[
+    R_ = {"res_GreenAll_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+    C_ = {"-solve -pos -petsc_prealloc 200 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
+    P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf