QNM expansion

ElectromagneticScattering/scattererTmatrix.geo

@@ -74,3 +74,9 @@ Mesh.VolumeEdges = 0;
 
 Mesh.ElementOrder = 2;
 Mesh.HighOrderOptimize = 1;
+//+
+Hide "*";
+//+
+Show {
+  Point{1}; Point{2}; Curve{1}; Curve{2}; Curve{3}; Surface{1}; Surface{2}; Volume{1}; 
+}

QuasiNormalModeExpansion/QNM_expansion_DtN.geo

+///////////////////////////////
+// Author : Guillaume Demesy //
+///////////////////////////////
+
+Include "QNM_expansion_DtN_data.geo";
+SetFactory("OpenCASCADE");
+lc_bg = lambda0/paramaille;
+lc_sc = lambda0/(paramaille*3);
+
+Rectangle(1) = {-dx_scat/2, -dy_scat/2, 0, dx_scat, dy_scat, 0};
+Disk(2) = {0, 0, 0, r_domain, r_domain};
+If (flag_PML==1)
+  Disk(3) = {0, 0, 0, r_pml, r_pml};
+EndIf
+Coherence;
+
+Physical Point("PrintPoint",100) = {1};
+
+Physical Surface("SCAT",1)       = {1};
+If (flag_PML==1)
+  Physical Surface("BG",2)  = {2};
+  Characteristic Length {6} = lc_bg;
+  Physical Surface("PML",3) = {3};
+  Physical Surface("SRC",4) = {};
+  Characteristic Length {2,3,4,5} = lc_sc;
+  Characteristic Length {1}       = lc_bg;
+  Characteristic Length {7}       = lc_bg;
+  Physical Line("SURF_SM",10)      = {};
+  Else
+  Physical Surface("BG",2)  = {2};
+  Physical Surface("SRC",4) = {};
+  Characteristic Length {2,3,4,5} = lc_sc;
+  Characteristic Length {1}       = lc_bg;
+  Physical Line("SURF_SM",10)      = {1};
+EndIf
+
+Rectangle(10) = {-dx_scat/2+domain_tot*yshift, -dy_scat/2, 0, dx_scat, dy_scat, 0};
+Disk(20) = {domain_tot*yshift, 0, 0, r_domain, r_domain};
+Delete{Surface{10,20};}
+If (flag_PML==1)
+  Disk(30) = {domain_tot*yshift, 0, 0, r_pml, r_pml};
+  Delete{Surface{30};}
+EndIf
+
+If (flag_PML==1)
+  Point(newp) = {-dx_scat/2+domain_tot*yshift, 2*r_pml, 0};
+Else
+  Point(newp) = {-dx_scat/2+domain_tot*yshift, 2*r_domain, 0};
+EndIf
+Geometry.Points = 0;
+

QuasiNormalModeExpansion/QNM_expansion_DtN_data.geo

+///////////////////////////////
+// Author : Guillaume Demesy //
+///////////////////////////////
+
+pp = "Model Options/";
+DefineConstant[
+    flag_DtN         = {2   , Name StrCat[pp, "0DtN"], Choices {0="PML",1="ABC",2="BT2"} , ServerAction "Reset"},
+    flag_vector_case = {0   , Name StrCat[pp, "1polarization case"], Choices {0="(0,0,Ez)",1="(Ex,Ey,0)"} },
+    flag_o2          = {1   , Name StrCat[pp, "2FE order"], Choices {0="order 1",1="order 2"} },
+    neig             = {120  , Name StrCat[pp, "3number of eigenvalues"] }
+];
+  
+nm         = 0.1;
+paramaille = 10;
+
+neig_PML = neig;
+neig_BT  = neig;
+
+a_lat      = 482.5*nm;
+lambda_m   = a_lat;
+
+lambda0   = 500*nm;
+dx_scat   = 200*nm;
+dy_scat   = 400*nm;
+epsr_scat_re = 9;
+epsr_scat_im = 0.5;
+epsr_bg      = 1;
+
+r_domain  = 0.5*Sqrt[dx_scat^2+dy_scat^2]+100*nm;
+r_pml     = r_domain+lambda0;
+
+eps0      = 8.854187817e-3*nm;
+mu0       = 400.*Pi*nm;
+cel       = 1.0/Sqrt[eps0 * mu0];
+
+norm      = dy_scat/(Pi*cel);
+
+omega_target_re  = 0.7/norm;
+omega_target_im  = 1/norm;
+
+target_re_PML = omega_target_re^2 - omega_target_im^2;
+target_im_PML = 2*omega_target_re*omega_target_im;
+target_re_BT  = omega_target_im;
+target_im_BT  = omega_target_re;
+
+eig_min_re = -5/norm;
+eig_max_re =  5/norm;
+eig_min_im = -8/norm;
+eig_max_im =  8/norm;
+
+If (flag_DtN==0)
+  flag_PML = 1;
+  flag_ABC_BAYLISS = 0;
+  slepc_options_rg = Sprintf[" -rg_interval_endpoints %.8lf,%.8lf,%.8lf,%.8lf ",eig_min_re,eig_max_re,eig_min_im,eig_max_im];
+  str_res = "res_modal_PML";
+  str_slepc_solver = " ";
+  domain_tot = 2*r_pml;
+Else
+  flag_PML = 0;
+  flag_ABC_BAYLISS = 1;
+  slepc_options_rg = " ";
+  If (flag_DtN==1)
+    str_res = "res_modal_AR_mur";
+  Else
+    str_res = "res_modal_AR_BT2";
+  EndIf
+  str_slepc_solver = " -nep_type nleigs -nep_rational -nep_target_magnitude -nep_view_pre -nep_monitor_all :temp_eigenvalues.txt ";
+  domain_tot = 2*r_domain;
+EndIf
+
+flag_PRINTVECTOR = 1;
+flag_PRINTFIELDS = 0;
+
+nb_Omegas     = 60;
+RealOmega_min = 0.5 / norm;
+RealOmega_max =   3 / norm;
+theta0        =  30 * Pi/180;
+
+tabrealomegasexpansion = LinSpace[RealOmega_min,RealOmega_max,nb_Omegas];
+// Freq = 0;
+
+// FIXME
+flag_manualexp = 1;
+
+yshift = 1.1;

QuasiNormalModeExpansion/QNM_expansion_dispersive_media.geo

+///////////////////////////////
+// Author : Guillaume Demesy //
+///////////////////////////////
+
+Include "QNM_expansion_dispersive_media_data.geo";
+lc_cell	  = lambda_m/paramaille;
+lc_sq     = lc_cell/2.;
+lc_src    = lc_cell;
+
+Point(1)  = {-a_lat    , -a_lat/2, 0. , lc_cell};
+Point(2)  = {-a_lat    ,  a_lat/2, 0. , lc_cell};
+Point(3)  = {     0    , -a_lat/2, 0. , lc_sq};
+Point(4)  = {     0    ,  a_lat/2, 0. , lc_sq};
+Point(5)  = { a_lat    , -a_lat/2, 0. , lc_cell};
+Point(6)  = { a_lat    ,  a_lat/2, 0. , lc_cell};
+Point(12)  = {-a_lat+0.7 *2*a_lat    , -a_lat/2, 0. , lc_sq};
+Point(13)  = {-a_lat+0.48*2*a_lat    ,  a_lat/2, 0. , lc_sq};
+
+Point(7)  = { cx_src  ,  cy_src , 0.       , lc_src};
+Point(8)  = { cx_src+r_src  ,  cy_src , 0. , lc_src};
+Point(9)  = { cx_src  ,  cy_src+r_src , 0. , lc_src};
+Point(10) = { cx_src-r_src  ,  cy_src , 0. , lc_src};
+Point(11) = { cx_src  ,  cy_src-r_src , 0. , lc_src};
+
+Circle(8)  = {10, 7, 9};
+Circle(9)  = {9, 7, 8};
+Circle(10) = {8, 7, 11};
+Circle(11) = {11, 7, 10};
+
+Line(12) = {7, 8};
+Line(13) = {7, 9};
+Line(14) = {7, 10};
+Line(15) = {7, 11};
+
+Line Loop(1) = {8, -13, 14};
+Plane Surface(1) = {-1};
+Line Loop(2) = {13, 9, -12};
+Plane Surface(2) = {-2};
+Line Loop(3) = {10, -15, 12};
+Plane Surface(3) = {-3};
+Line Loop(4) = {14, -11, -15};
+Plane Surface(4) = {4};
+
+Line(16) = {1, 3};
+Line(17) = {3, 12};
+Line(18) = {12, 5};
+Line(19) = {5, 6};
+Line(20) = {6, 4};
+Line(21) = {4, 13};
+Line(22) = {13, 2};
+Line(23) = {2, 1};
+Line(24) = {3, 13};
+Line(25) = {12, 4};
+Curve Loop(5) = {23, 16, 24, 22};
+Plane Surface(5) = {5};
+Curve Loop(6) = {24, -21, -25, -17};
+Plane Surface(6) = {-6};
+Curve Loop(7) = {25, -20, -19, -18};
+Curve Loop(8) = {8, 9, 10, 11};
+Plane Surface(7) = {-7, 8};
+
+Physical Surface(1) = {6};       // Dispersive medium
+Physical Surface(2) = {5,7};     // Background
+Physical Surface(3) = {1,2,3,4}; // disk green
+Physical Line(10)   = {16,17,18,19,20,21,22,23}; // Outer box
+Physical Line(20)   = {24,25};   // bound lag
+Physical Line(30)   = {12};      // Edge Source Green X
+Physical Line(31)   = {13};      // Edge Source Green Y
+Physical Point(100) = {7};
+
+// dummy lines to see expansion and direct solutions side by side
+pp=newp;
+Point(pp+1)  = {-a_lat+ xshift*2*a_lat    , -a_lat/2 , 0. , lc_cell};
+Point(pp+2)  = {-a_lat+ xshift*2*a_lat    ,  a_lat/2 , 0. , lc_cell};
+Point(pp+3)  = {     0+ xshift*2*a_lat    , -a_lat/2 , 0. , lc_sq};
+Point(pp+4)  = {     0+ xshift*2*a_lat    ,  a_lat/2 , 0. , lc_sq};
+Point(pp+5)  = { a_lat+ xshift*2*a_lat    , -a_lat/2 , 0. , lc_cell};
+Point(pp+6)  = { a_lat+ xshift*2*a_lat    ,  a_lat/2 , 0. , lc_cell};
+Point(pp+12)  = {-a_lat+0.7 *2*a_lat+ xshift*2*a_lat    , -a_lat/2, 0. , lc_sq};
+Point(pp+13)  = {-a_lat+0.48*2*a_lat+ xshift*2*a_lat    ,  a_lat/2, 0. , lc_sq};
+
+c=newl;
+Line(c+16) = {pp+1, pp+3};
+Line(c+17) = {pp+3, pp+12};
+Line(c+18) = {pp+12, pp+5};
+Line(c+19) = {pp+5, pp+6};
+Line(c+20) = {pp+6, pp+4};
+Line(c+21) = {pp+4, pp+13};
+Line(c+22) = {pp+13, pp+2};
+Line(c+23) = {pp+2, pp+1};
+Line(c+24) = {pp+3, pp+13};
+Line(c+25) = {pp+12, pp+4};
+
+// dummy point to enforce window center 
+pp=newp;
+Point(pp+1)  = {-a_lat+ xshift*2*a_lat , 1.5*a_lat , 0. , lc_cell};
+
+Geometry.Points = 0;

QuasiNormalModeExpansion/QNM_expansion_dispersive_media.pro

+///////////////////////////////
+// Author : Guillaume Demesy //
+///////////////////////////////
+
+Include "QNM_expansion_dispersive_media_data.geo";
+
+Group {
+  // Domains
+  Om_1		   = Region[1];
+  Om_src     = Region[3];
+  Om_2		   = Region[{2,3}];
+  Om	  	   = Region[{Om_1,Om_2}];
+  Om_but_src = Region[{Om,-Om_src}];
+
+  Surfdiri     = Region[10];
+  Bound        = Region[20];
+  pt_src       = Region[100];
+  green_x_src  = Region[30];
+  green_y_src  = Region[31];
+  PrintPoint   = Region[100];
+  Om_plot     = Region[{Om,-Om_src}];
+}
+
+Function{
+  I[]               = Complex[0,1];
+  k0[]              = $realomega/cel;
+  mur[]             = TensorDiag[1,1,1];
+  epsr_nod[]        = TensorDiag[1,1,1];
+  om_lorentz_1[]    = Complex[om_lorentz_1_re,om_lorentz_1_im];
+  mur_direct[]      = 1;
+  epsr_direct[Om_1] = 1 + a_lorentz_1/($realomega-om_lorentz_1[]) - a_lorentz_1/($realomega+Conj[om_lorentz_1[]]);
+  epsr_direct[Om_2] = 1;
+  epsr_mau1[]       = 1 + I[]*a_lorentz_1/($Eigenvalue-I[]*om_lorentz_1[]) - I[]*a_lorentz_1/($Eigenvalue+I[]*Conj[om_lorentz_1[]]);
+  
+  If (flag_vector_case==0) 
+      source_RHS[Om_src]     =  Vector[0,0,1];
+      source_RHS[Om_but_src] =  Vector[0,0,0];
+  Else
+      source_RHS[Om_src]     =  Vector[1,0,0];
+      source_RHS[Om_but_src] =  Vector[0,0,0];
+  EndIf
+}
+
+Constraint {
+	{Name Dirichlet; Type Assign;
+		Case {
+			{ Region Surfdiri ; Value 0.; }
+		}
+	}
+  
+  { Name Dirichlet_src;
+    Case {
+      { Region green_x_src;
+        Type Assign; 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 Line        ; NumberOfPoints  4 ; }
+	        { GeoElement Triangle    ; NumberOfPoints  12 ; }
+			  }
+      }
+    }
+  }
+}
+
+FunctionSpace {
+  { Name Hgrad_perp; Type Form1P;
+    BasisFunction {
+      { Name un;  NameOfCoef un;  Function BF_PerpendicularEdge_1N; Support Region[Om]; Entity NodesOf[All]; }
+      { Name un2; NameOfCoef un2; Function BF_PerpendicularEdge_2E; Support Region[Om]; Entity EdgesOf[All]; }
+      If(flag_o2==1)
+        { Name un3; NameOfCoef un3; Function BF_PerpendicularEdge_2F; Support Region[Om]; Entity FacetsOf[All]; }
+        { Name un4; NameOfCoef un4; Function BF_PerpendicularEdge_3E; Support Region[Om]; Entity EdgesOf[All]; }
+        { Name un5; NameOfCoef un5; Function BF_PerpendicularEdge_3F; Support Region[Om]; Entity FacetsOf[All]; }
+      EndIf
+     }
+    Constraint {
+      { NameOfCoef un;  EntityType NodesOf ; NameOfConstraint Dirichlet; }
+      { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+      If(flag_o2==1)
+        { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+        { NameOfCoef un4; EntityType EdgesOf  ; NameOfConstraint Dirichlet; }
+        { NameOfCoef un5; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+      EndIf
+    }
+  }
+
+  { Name Hcurl; Type Form1;
+    BasisFunction {
+      { Name sn;  NameOfCoef un;  Function BF_Edge   ; Support Region[Om]; Entity EdgesOf[All]; }
+      { Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[Om]; Entity EdgesOf[All]; }
+      If(flag_o2==1)
+        { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[Om]; Entity FacetsOf[All]; }
+        { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[Om]; Entity FacetsOf[All]; }
+        { Name sn5; NameOfCoef un5; Function BF_Edge_4E  ; Support Region[Om]; Entity EdgesOf[All]; }
+      EndIf
+    }
+    Constraint {
+      { NameOfCoef un ; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+      { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+      { NameOfCoef un3; EntityType FacetsOf; NameOfConstraint Dirichlet; }
+      { NameOfCoef un4; EntityType FacetsOf; NameOfConstraint Dirichlet; }
+      { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+
+    }
+  }
+  
+  { Name Hcurl_direct_green_lin; Type Form1;
+    BasisFunction {
+      { Name sn;  NameOfCoef un ; Function BF_Edge   ; Support Region[{Om,green_x_src}]; Entity EdgesOf[All]; }
+      { Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[{Om,green_x_src}]; Entity EdgesOf[All]; }
+      If(flag_o2==1)
+        { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[Om]; Entity FacetsOf[All]; }
+        { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[Om]; Entity FacetsOf[All]; }
+        { Name sn5; NameOfCoef un5; Function BF_Edge_4E  ; Support Region[Om]; Entity EdgesOf[All]; }
+      EndIf
+    }
+    Constraint{
+      { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet_src; }
+      { NameOfCoef un3; EntityType FacetsOf; NameOfConstraint Dirichlet_src; }
+      { NameOfCoef un4; EntityType FacetsOf; NameOfConstraint Dirichlet_src; }
+      { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint Dirichlet_src; }
+    }
+  }
+}
+
+Formulation {
+  {Name modal_helmholtz; Type FemEquation;
+    Quantity{
+      If (flag_vector_case==0) 
+        { Name u; Type Local; NameOfSpace Hgrad_perp;}
+      Else
+        { Name u; Type Local; NameOfSpace Hcurl;}
+      EndIf
+    }
+    Equation {
+      Galerkin { Eig[  1/mur[] *    Dof{Curl u}, {Curl u} ]; Rational 1; In Om  ; Jacobian JSur; Integration Int_1;}
+      Galerkin { Eig[ epsr_nod[]/cel^2 * Dof{u}, {u}      ]; Rational 2; In Om_1; Jacobian JSur; Integration Int_1;}
+      Galerkin { Eig[ epsr_nod[]/cel^2 * Dof{u}, {u}      ]; Rational 3; In Om_2; Jacobian JSur; Integration Int_1;}
+    }
+  }
+
+  // if other matrices, fails to print eigs (conflict between M1 matrices)
+  {Name form_RHS; Type FemEquation;
+    Quantity{
+      If (flag_vector_case==0) 
+        { Name u; Type Local; NameOfSpace Hgrad_perp;}
+      Else
+        { Name u; Type Local; NameOfSpace Hcurl;}
+      EndIf
+    }
+    Equation{ 
+      Galerkin { [   1. * Dof{u} , {u} ]; In Om; Jacobian JVol; Integration Int_1;}
+      Galerkin { [ source_RHS[] , {u} ]; In Om; Jacobian JVol; Integration Int_1;  }
+    }
+  }
+  
+  {Name direct_helmholtz; Type FemEquation;
+    Quantity{
+      If (flag_vector_case==0) 
+        { Name u; Type Local; NameOfSpace Hgrad_perp;}
+      Else
+        { Name u; Type Local; NameOfSpace Hcurl;}
+      EndIf
+    }
+    Equation {
+      Galerkin { [ -1/mur_direct[]*Dof{Curl u}  , {Curl u}]; In Om           ; Jacobian JSur; Integration Int_1; }
+      Galerkin { [ k0[]^2*epsr_direct[]*Dof{u}    , {u}]   ; In Om           ; Jacobian JSur; Integration Int_1; }
+      Galerkin { [ -source_RHS[]                , {u}]     ; In Om_src; Jacobian JVol; Integration Int_1;  }
+    }
+  }
+}
+
+Resolution {
+  { Name res_nleig_rat_exp;
+    System{
+      { Name V ; NameOfFormulation form_RHS   ; Type ComplexValue; }
+      { Name M1; NameOfFormulation modal_helmholtz ; Type ComplexValue; }
+      { Name M2; NameOfFormulation direct_helmholtz ; Type ComplexValue; }
+    }
+    Operation{
+      CreateDir["out"];
+      For k In {0:#tabrealomegasexpansion()-1}
+        Evaluate[$realomega=tabrealomegasexpansion(k)];
+        GenerateListOfRHS[V,Om,#tabrealomegasexpansion()];
+      EndFor
+      GenerateSeparate[M1];
+      EigenSolveAndExpand[M1,neig,eig_target_re,eig_target_im,
+          { {1}, {-1,-2*om_lorentz_1_im,2*a_lorentz_1*om_lorentz_1_re-om_lorentz_1_mod^2,0,0}, {1,0,0} } ,
+          { {1}, {-1,-2*om_lorentz_1_im,-om_lorentz_1_mod^2},                                  {1} } ,
+          tabrealomegasexpansion(),
+          V];
+      PostOperation[postop_lorentz_nleig_rat_exp];
+      For k In {0:#tabrealomegasexpansion()-1}
+        Evaluate[$realomega=tabrealomegasexpansion(k)];
+        Generate[M2];
+        Solve[M2];
+        PostOperation[postop_direct_helmholtz];
+      EndFor
+      For k In {0:#tabrealomegasexpansion()-1}
+        PostOperation[postop_gmsh~{k}];
+        Sleep[0.3];
+      EndFor
+    }
+  }
+}
+
+PostProcessing {
+  { Name postpro_lorentz_nleig_rat_exp; NameOfFormulation modal_helmholtz;
+    Quantity {
+      { Name intnormu       ; Value { Integral { [ Norm[{u}]    ] ; In Om    ; Jacobian JSur; Integration Int_1 ; } } }
+      { Name intnormu_rod   ; Value { Integral { [ Norm[{u}]    ] ; In Om_1    ; Jacobian JSur; Integration Int_1 ;} } }
+      { Name u   ; Value { Local { [ {u}    ] ; In Om; Jacobian JSur; } } }
+      { Name u_im; Value { Local { [ Im[{u}]    ] ; In Om; Jacobian JSur; } } }
+      { Name EigenValuesReal;  Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JSur; } } }
+      { Name EigenValuesImag;  Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JSur; } } }
+      { Name vol ; Value { Integral { [1]   ; In PrintPoint       ; Jacobian JVol; Integration Int_1 ; } } }
+      { Name eig_re ; Value { Integral { [$EigenvalueReal/$vol] ; In PrintPoint ; Jacobian JVol; Integration Int_1 ; } } }
+      { Name eig_im ; Value { Integral { [$EigenvalueImag/$vol] ; In PrintPoint ; Jacobian JVol; Integration Int_1 ; } } }
+    }
+  }
+
+  { Name postpro_direct_helmholtz; NameOfFormulation direct_helmholtz;
+    Quantity {
+      { Name eps          ; Value { Local { [ epsr_direct[]] ; In Om   ; Jacobian JSur; } } }
+      { Name source_RHS   ; Value { Local { [ source_RHS[] ] ; In Om; Jacobian JSur; } } }
+      { Name u_direct     ; Value { Local { [ {u}          ] ; In Om; Jacobian JSur; } } }
+      { Name u_direct_im  ; Value { Local { [ Im[{u}]      ] ; In Om; Jacobian JSur; } } }
+      { Name intnormu_rod ; Value { Integral { [ Norm[{u}] ] ; In Om_1 ; Jacobian JSur; Integration Int_1 ;} } }
+	  }
+  }
+}
+
+
+PostOperation { 
+  { Name postop_lorentz_nleig_rat_exp; NameOfPostProcessing postpro_lorentz_nleig_rat_exp ;
+    Operation {
+      Print [vol[PrintPoint] , OnGlobal, TimeStep 0, StoreInVariable $vol,Format Table];
+      For j In {0:(neig-1)}
+        Print [eig_re[PrintPoint] , OnGlobal, TimeStep j , Format Table , SendToServer "GetDP/eig_re"];
+        Print [eig_im[PrintPoint] , OnGlobal, TimeStep j , Format Table , SendToServer "GetDP/eig_im"];
+      EndFor
+      Print [EigenValuesReal, OnElementsOf PrintPoint, Format TimeTable, File "out/EigenValuesReal.txt"];
+      Print [EigenValuesImag, OnElementsOf PrintPoint, Format TimeTable, File "out/EigenValuesImag.txt"];
+      Print [ u  , OnElementsOf Om, File Sprintf["out/u.pos"], EigenvalueLegend];
+      Print [ u  , OnPoint {cx_src,-cy_src,0}, Format TimeTable, File "out/u_expand_on_probe_point.txt"];
+      For k In {0:#tabrealomegasexpansion()-1}
+        // Print [ u_im, OnElementsOf Om, File Sprintf["out/u_exp_%03g.out",k], Name Sprintf["u_expand_imag_%03g (left)",k] , EigenvalueLegend, TimeStep {neig+k}];
+        Print [ u_im, OnGrid {$A,$B, 0} {-a_lat:a_lat:a_lat/30, -a_lat/2:a_lat/2:a_lat/15, 0}, File Sprintf["out/u_exp_%03g.out",k], Name Sprintf["u_expand_imag_%03g (left)",k] , EigenvalueLegend, TimeStep {neig+k}];
+      EndFor
+    }
+  }
+  { Name postop_direct_helmholtz; NameOfPostProcessing postpro_direct_helmholtz ;
+    Operation {
+      Print[ u_direct_im, OnGrid {$A,$B, 0} {-a_lat:a_lat:a_lat/30, -a_lat/2:a_lat/2:a_lat/15, 0}, File "out/u_direct.pos" , ChangeOfCoordinates {$X+xshift*2*a_lat,$Y,$Z} , Name "u_direct_imag (right)"];      
+      Print[ u_direct , OnPoint {cx_src,-cy_src,0}, Format TimeTable, File >> "out/u_direct_on_probe_point.txt"];
+      Print [intnormu_rod[Om_1], OnGlobal, Format TimeTable, File "out/intnormu_direct_on_src_point.out" ];
+		}
+  }
+  For k In {0:#tabrealomegasexpansion()-1}
+    { Name postop_gmsh~{k}; NameOfPostProcessing postpro_direct_helmholtz ;
+      Operation {
+        Echo[Str["Geometry.Points = 0;","Geometry.Color.Curves = {0,0,0};","Mesh.SurfaceEdges = 0;",
+        "k=PostProcessing.NbViews-1;","View[k].Visible = 0;",
+        "View.Visible = 0;","View[0].Visible = 1;", "View[0].ColormapNumber = 7;",
+        Sprintf["View[%g].Visible = 0;",k],
+        Sprintf["View[%g].Visible = 0;",nb_Omegas+k],
+        Sprintf["View[%g].ColormapNumber = 7;",k+1],
+        Sprintf["View[%g].ColormapNumber = 7;",nb_Omegas+k+1],
+        Sprintf["View[%g].Visible = 1;",k+1],
+        Sprintf["View[%g].Visible = 1;",nb_Omegas+k+1]], File "tmp0.geo"];
+        // Echo[Str["Plugin(Annotate).Text= "TOP : DIRECT -- BOTTOM : EXPANSION\";" , "Plugin(Annotate).Run;"]];
+      }
+    }  
+  EndFor
+}
+
+
+slepc_options_rg  = Sprintf[" -rg_interval_endpoints %.8lf,%.8lf,%.8lf,%.8lf ",eig_min_re,eig_max_re,eig_min_im,eig_max_im];
+DefineConstant[
+  R_ = {"res_nleig_rat_exp", Name "GetDP/1ResolutionChoices", Visible 1},
+  C_ = {StrCat["-pre res_nleig_rat_exp -cal -slepc -nep_type nleigs -nep_rational -nep_target_magnitude -nep_view_pre -nep_ncv 300 -nep_monitor_all :temp_eigenvalues.txt",slepc_options_rg] , Name "GetDP/9ComputeCommand", Visible 1}
+];
+
+DefineConstant[
+  eig_re_  = {0, Name "GetDP/eig_re", ReadOnly 1, Graph "10", Visible 0},
+  eig_im_  = {0, Name "GetDP/eig_im", ReadOnly 1, Graph "01", Visible 0}
+];
+
+// DefineConstant[
+//   R_ = {"res_nleig_rat_exp", Name "GetDP/1ResolutionChoices", Visible 1},
+//   C_ = {"-pre res_nleig_rat_exp -cal -slepc -nep_max_it 1000 -nep_type nleigs -nep_rational -nep_target_magnitude -nep_ncv 100 -rg_interval_endpoints 0.00000000,0.05763407,0.00390394,3.12315286 -nep_view_pre -nep_monitor_all :temp_eigenvalues.txt", Name "GetDP/9ComputeCommand", Visible 1}
+// ];

QuasiNormalModeExpansion/QNM_expansion_dispersive_media_data.geo

+///////////////////////////////
+// Author : Guillaume Demesy //
+///////////////////////////////
+
+pp = "Model Options/";
+DefineConstant[
+  flag_vector_case = {1   , Name StrCat[pp, "0polarization case"], Choices {0="(0,0,Ez)",1="(Ex,Ey,0)"} },
+  flag_o2          = {1   , Name StrCat[pp, "1FE order"], Choices {0="order 1",1="order 2"} },
+  neig             = {100 , Name StrCat[pp, "0number of eigenvalues"]}
+];
+  
+nm            = 0.1;
+N_true_poles  = 1;
+paramaille    = 5;
+
+a_lat         = 482.5*nm;
+lambda_m      = a_lat;
+
+eps0      = 8.854187817e-3*nm;
+mu0       = 400.*Pi*nm;
+cel       = 1.0/Sqrt[eps0 * mu0];
+norm      = a_lat/(2.*Pi*cel);
+
+om_p_josab    = .4*norm;
+om_0_josab    = 0.3*norm;
+gama_josab    = 0.05*norm;
+om_lorentz_1_im  = -gama_josab/2;
+om_lorentz_1_re  = Sqrt[om_0_josab^2-om_lorentz_1_im^2];
+a_lorentz_1      = -om_p_josab^2/(2*om_lorentz_1_re);
+om_lorentz_1_mod = Sqrt[om_lorentz_1_re^2+om_lorentz_1_im^2];
+
+eig_target_re =  om_lorentz_1_im/4;
+eig_target_im =  1/norm;
+eig_max_im    =  8/norm;
+eig_min_im    =  0.01/norm;
+eig_min_re    =  0;
+eig_max_re    =  0.8*Fabs[om_lorentz_1_im];
+
+cx_src = 0.75*2*a_lat - a_lat;
+cy_src = 0.75*a_lat - a_lat/2;
+r_src  = a_lat/100;
+
+nb_Omegas     = 50;
+RealOmega_min = 0.1;
+RealOmega_max = 1.1;
+tabrealomegasexpansion = LinSpace[RealOmega_min,RealOmega_max,nb_Omegas];
+Freq = 0;
+
+xshift = 1.1;