Merge branch 'master' of https://gitlab.onelab.info/doc/models

ElectromagneticScattering/README.md

+A Onelab model for 3D scattering problems in nanophotonics.
+
+## Synopsis
+
+This project contains a [Onelab](http://onelab.info/wiki/ONELAB) model for solving various 3D electromagnetic scattering problems on an isolated object :
+* T-matrix computation: See https://arxiv.org/abs/1802.00596 for details
+* Quasi-normal modes
+* Plane wave response
+* Green's function (LDOS)
+
+## Installation
+
+This model requires the following (open-source) programs:
+* [onelab](http://onelab.info/wiki/ONELAB), which bundles both [gmsh](http://www.gmsh.info/) and [getdp](http://www.getdp.info/)
+* python (v2.7.x or v3.6.x) with numpy, scipy and matplotlib
+
+## Running the model
+
+Open `scattering.pro` with Gmsh.
+
+## Authors
+
+Guillaume Demésy and Brian Stout
\ No newline at end of file

ElectromagneticScattering/scattering_data.geo

+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattering_data.geo       //
+///////////////////////////////
+
+nm       = 1.;
+epsilon0 = 8.854187817e-3*nm;
+mu0      = 400.*Pi*nm;
+cel      = 1.0/(Sqrt[epsilon0 * mu0]);
+deg2rad  = Pi/180.;
+
+pp0        = "1Geometry/0";
+pp1        = "2Study Type/0";
+pp2        = "3Electromagnetic parameters/0";
+pp3        = "4Mesh size and PMLs parameters/0";
+pp4        = "5Postpro options/0";
+pp5        = "6Post plot options/0";
+close_menu = 0;
+colorro    = "LightGrey";
+colorppOK   = "Ivory";
+colorppWA   = "LightSalmon1";
+colorppNO   = "LightSalmon4";
+
+ELL    = 0;
+PARALL = 1;
+CYL    = 2;
+CONE   = 3;
+TOR    = 4;
+
+RES_PW    = 0;
+RES_TMAT  = 1;
+RES_GREEN = 2;
+RES_QNM   = 3;
+
+DefineConstant[
+  flag_shape = {ELL , Name StrCat[pp0, "0Scatterer shape"], 
+    Choices {ELL="ellispoid",PARALL="parallelepiped",CYL="cylinder",CONE="cone",TOR="split torus"},Closed 0,GmshOption "Reset", Autocheck 0}];
+
+If (flag_shape==ELL)
+  DefineConstant[
+    ell_rx = {73.45 , Name StrCat[pp0  , "1ellipsoid X-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+    ell_ry = {73.45 , Name StrCat[pp0  , "2ellipsoid Y-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+    ell_rz = {73.45 , Name StrCat[pp0  , "3ellipsoid Z-radius [nm]"], Highlight Str[colorppOK] , Closed 0}
+  ];
+  // FIXME gmsh Max?
+  If (ell_rx>ell_ry)
+    rbb_tmp=ell_rx;
+  Else rbb_tmp=ell_ry;
+  EndIf
+  If (ell_rz>rbb_tmp)
+    rbb=ell_rz;
+  Else rbb=rbb_tmp;
+  EndIf
+EndIf
+If (flag_shape==PARALL)
+  DefineConstant[
+    par_ax = {300 , Name StrCat[pp0  , "1cube X-edge size [nm]"], Highlight Str[colorppOK]  , Closed 0},
+    par_ay = {100 , Name StrCat[pp0  , "2cube Y-edge size [nm]"], Highlight Str[colorppOK]  , Closed 0},
+    par_az = {600 , Name StrCat[pp0  , "3cube Z-edge size [nm]"], Highlight Str[colorppOK]  , Closed 0}];
+  rbb = 0.5*Sqrt[par_ax^2+par_ay^2+par_az^2];
+EndIf
+If (flag_shape==CYL)
+  DefineConstant[
+    cyl_rx = {300 , Name StrCat[pp0  , "1cylinder X-radius [nm]"], Highlight Str[colorppOK]  , Closed 0},
+    cyl_ry = {300 , Name StrCat[pp0  , "2cylinder Y-radius [nm]"], Highlight Str[colorppOK]  , Closed 0},
+    cyl_h  = {300 , Name StrCat[pp0  , "3cylinder height [nm]"]  , Highlight Str[colorppOK]  , Closed 0}];
+  // FIXME gmsh Max?
+  If (cyl_rx>cyl_ry)
+    rbb_tmp=cyl_rx;
+  Else rbb_tmp=cyl_ry;
+  EndIf
+  rbb = 0.5*Sqrt[rbb_tmp^2+cyl_h^2];
+EndIf
+If (flag_shape==CONE)
+  DefineConstant[
+    cone_rx = {300 , Name StrCat[pp0  , "1cone basis X-radius [nm]"], Highlight Str[colorppOK]  , Closed 0},
+    cone_ry = {300 , Name StrCat[pp0  , "1cone basis Y-radius [nm]"], Highlight Str[colorppOK]  , Closed 0},
+    cone_h  = {300 , Name StrCat[pp0  , "2cone height [nm]"]        , Highlight Str[colorppOK]  , Closed 0}];
+  // FIXME gmsh Max?
+  If (cone_rx>cone_ry)
+    rbb_tmp=cone_rx;
+  Else rbb_tmp=cone_ry;
+  EndIf    
+  If (2.0*cone_h/3.0>Sqrt[rbb_tmp^2+cone_h^2/9.0])
+    rbb=2.0*cone_h/3.0;
+  Else rbb=Sqrt[rbb_tmp^2+cone_h^2/9.0];
+  EndIf
+EndIf
+If (flag_shape==TOR)
+  DefineConstant[
+    tor_r1    = { 300 , Name StrCat[pp0 , "1torus radius 1  [nm]"], Highlight Str[colorppOK] , Closed 0},
+    tor_r2x   = { 100 , Name StrCat[pp0 , "2torus radius 2x [nm]"], Highlight Str[colorppOK] , Closed 0},
+    tor_r2z   = { 50  , Name StrCat[pp0 , "3torus radius 2z [nm]"], Highlight Str[colorppOK] , Closed 0},
+    tor_angle = { 340 , Name StrCat[pp0 , "4torus angle [deg]"]   , Highlight Str[colorppOK] , Closed 0, Min 5, Max 355}];
+  rbb = tor_r1+tor_r2x;
+EndIf
+DefineConstant[
+  rot_theta = {0 , Name StrCat[pp0  , "5rotate scatterer (polar) [deg]"] , Highlight Str[colorppOK]  , Closed 0, Min 0, Max 180},
+  rot_phi   = {0 , Name StrCat[pp0  , "6rotate scatterer (azimut) [deg]"], Highlight Str[colorppOK]  , Closed 0, Min 0, Max 360}];
+
+DefineConstant[
+  flag_study = { RES_TMAT , Name StrCat[pp1, "0Study type"], 
+        Choices {RES_PW="Plane Wave Response",
+                 RES_TMAT="T-Matrix", 
+                 RES_GREEN="Green's Tensor",
+                 RES_QNM="Quasi-Normal Modes"},Closed 0,GmshOption "Reset", Autocheck 0}];
+
+DefineConstant[
+  epsr_In_re  = { 9., Name StrCat[pp2  , "0scatterer permittivity (real) []"] , Highlight Str[colorppOK] , Closed 0},
+  epsr_In_im  = { 0., Name StrCat[pp2  , "1scatterer permittivity (imag) []"] , Highlight Str[colorppOK] , Closed 0},
+  epsr_Out_re = { 1., Name StrCat[pp2  , "2background permittivity (real) []"], Highlight Str[colorppOK] , Closed 0},
+  epsr_Out_im = { 0., Name StrCat[pp2  , "3background permittivity (imag) []"], Highlight Str[colorppOK] , Closed 0}
+];
+
+If (flag_study!=RES_QNM)
+  DefineConstant[
+    lambda0 = {587.6 , Name StrCat[pp2  , "4wavelength [nm]"] , Highlight Str[colorppOK]  , Closed 0,GmshOption "Reset", Autocheck 0}];
+Else
+  DefineConstant[
+    lambda0 = {1000 , Name StrCat[pp2  , "4wavelength target   [nm]"] , Highlight Str[colorppOK]  , Closed 0,GmshOption "Reset", Autocheck 0},
+    neig    = {30   , Name StrCat[pp2  , "5number of eigenvalues []"] , Highlight Str[colorppOK]  , Closed 0,GmshOption "Reset", Autocheck 0}
+    ];
+EndIf
+lambda_bg = lambda0/Sqrt[epsr_Out_re];
+  
+If (flag_study==RES_PW)
+  DefineConstant[
+    theta0 = {0 , Name StrCat[pp2  , "5plane wave theta [deg]"], Highlight Str[colorppOK]  , Closed 0, Min 0, Max 180},
+    phi0   = {0 , Name StrCat[pp2  , "6plane wave phi [deg]"]  , Highlight Str[colorppOK]  , Closed 0, Min 0, Max 360},
+    psi0   = {0 , Name StrCat[pp2  , "7plane wave psi [deg]"]  , Highlight Str[colorppOK]  , Closed 0, Min 0, Max 180}];
+    theta0 = theta0 *deg2rad;
+    phi0 = phi0 * deg2rad;
+    psi0 = psi0 * deg2rad;
+EndIf
+If (flag_study==RES_GREEN)
+  DefineConstant[
+    x_p = {0    , Name StrCat[pp2  , "5Green X-point [nm]"]  , Highlight Str[colorppOK]  , Closed 0},
+    y_p = {0    , Name StrCat[pp2  , "6Green Y-point [nm]"]  , Highlight Str[colorppOK]  , Closed 0},
+    z_p = {-rbb*1.1 , Name StrCat[pp2  , "7Green Z-point [nm]"]  , Highlight Str[colorppOK]  , Closed 0}];
+    x_p=x_p*nm;
+    y_p=y_p*nm;
+    z_p=z_p*nm;
+EndIf
+
+DefineConstant[
+  n_max = {1 , Name StrCat[pp2  , "8n_max integer"], Highlight Str[colorppOK]  , Closed 0, Min 0, Max 5},
+  siwt  = {0 , Name StrCat[pp2  , "9Time sign e^(+|-iwt)"], Choices {0="e^(-iwt)",2="e^(+iwt)"}  , Closed 0}
+];
+DefineConstant[
+  flag_cartpml = {0 , Name StrCat[pp3  , "0PML type"] , Choices {0="spherical",1="cartesian"}, Closed 0},
+  pml_size     = {lambda_bg/1.5, Name StrCat[pp3  , "1PML thickness [nm]"] , Highlight Str[colorppWA]  , Closed 0, Min (lambda_bg/10), Max (3*lambda_bg)},
+  paramaille   = {4. , Name StrCat[pp3  , "2mesh size"], Highlight Str[colorppWA]  , Closed 0},
+  refine_scat  = {2. , Name StrCat[pp3  , "3scatterer mesh refinement"], Highlight Str[colorppWA]  , Closed 0}
+  is_FEM_o2    = {1 , Name StrCat[pp3   , "4Interpolation order "] , Choices {0="order 1",1="order 2"}, Closed 0}
+];
+DefineConstant[
+  space2pml  = {lambda_bg/10.*4 , Name StrCat[pp4  , "0space around scatterer [nm]"] , Highlight Str[colorppNO]    , Closed 1,Min (lambda_bg/10.), Max (3*lambda_bg)},
+  dist_rcut  = {lambda_bg/10.   , Name StrCat[pp4  , "1min dist from object & PML"] , Highlight Str[colorppNO]    , Closed 1},
+  nb_cuts    = {3               , Name StrCat[pp4  , "2number of cuts [integer]"] , Highlight Str[colorppNO]    , Closed 1},
+  npts_theta = {50              , Name StrCat[pp4  , "3polar sampling [integer]"], Highlight Str[colorppNO]     , Closed 0,Min 50, Max 300},
+  npts_phi   = {100             , Name StrCat[pp4  , "4azimuthal sampling [integer]"], Highlight Str[colorppNO] , Closed 0,Min 100, Max 600}
+];
+DefineConstant[
+  flag_plotcuts = {0, Choices{0,1}, Name StrCat[pp5, "Plot radial cuts?"]},
+  flag_FF       = {1, Choices{0,1}, Name StrCat[pp5, "Plot far field?"]}
+];
+
+
+
+// FIXME conditional for normalization
+If (flag_shape==ELL)
+  ell_rx = ell_rx*nm;
+  ell_ry = ell_ry*nm;
+  ell_rz = ell_rz*nm;
+EndIf
+If (flag_shape==PARALL)
+    par_ax = par_ax*nm;
+    par_ay = par_ay*nm;
+    par_az = par_az*nm;
+EndIf
+If (flag_shape==CYL)
+  cyl_rx = cyl_rx*nm;
+  cyl_ry = cyl_ry*nm;
+  cyl_h  = cyl_h *nm;
+EndIf
+If (flag_shape==CONE)
+  cone_rx = cone_rx*nm;
+  cone_ry = cone_ry*nm;
+  cone_h = cone_h*nm;
+EndIf
+If (flag_shape==TOR)
+  tor_r1 = tor_r1*nm;
+  tor_r2x = tor_r2x*nm;
+  tor_r2z = tor_r2z*nm;
+EndIf
+
+lambda0   = lambda0*nm;
+lambda_bg = lambda_bg*nm;
+pml_size  = pml_size*nm;
+dist_rcut = dist_rcut*nm;
+space2pml = space2pml*nm;
+rbb       = rbb*nm;
+
+r_pml_in  = space2pml+rbb;
+r_pml_out = space2pml+rbb+pml_size;
+r_sph_min = rbb+dist_rcut;
+r_sph_max = space2pml+rbb-dist_rcut;
+
+sph_scan = 0.00001;
+
+npts_plot_theta = 25;
+npts_plot_phi   = 50;
+
+p_max = n_max*n_max+2*n_max;
+siwt=siwt-1;

ElectromagneticScattering/scattering_init.py

+# -*- coding: utf-8 -*-
+
+"""
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattering_init.py        //
+///////////////////////////////
+"""
+
+import sys,os
+import numpy as np
+from scipy.special import jv, yv, hankel1, hankel2
+sys.path.append(os.getcwd())
+from matplotlib import cm
+import pylab as pl
+from scattering_tmp import *
+np.set_printoptions(precision=2)
+pi = np.pi
+
+ps       = np.linspace(1,p_max,p_max)
+r_sphs   = np.linspace(r_sph_min,r_sph_max,nb_cuts)
+k_Out    = 2*pi*np.sqrt(epsr_Out_re)/lambda0
+epsr_In  = epsr_In_re+1j*epsr_In_im
+epsr_Out = epsr_Out_re+1j*epsr_Out_im
+
+phi_range   = np.linspace(sph_scan,2*pi-sph_scan,npts_phi)
+theta_range = np.linspace(sph_scan,  pi-sph_scan,npts_theta)
+[phi_sph,theta_sph]=np.meshgrid(phi_range,theta_range)
+cos_theta = np.cos(theta_range)
+sin_theta = np.sin(phi_range)
+
+def field_VSH_expansion(post_filename):
+    m_max = n_max
+    p_max = n_max**2 +2*n_max
+    FF_Xnm_t        = np.zeros((npts_theta,npts_phi,p_max),dtype=complex)
+    FF_Xnm_p        = np.zeros((npts_theta,npts_phi,p_max),dtype=complex)
+    FF_erCrossXnm_t = np.zeros((npts_theta,npts_phi,p_max),dtype=complex)
+    FF_erCrossXnm_p = np.zeros((npts_theta,npts_phi,p_max),dtype=complex)
+    #####################################
+    ##### sn, pn ,un
+    ##### Brian's recurrence relations
+    B_Pnpms        = np.zeros((npts_theta,m_max+1,n_max+1))
+    B_unpms        = np.zeros((npts_theta,m_max+1,n_max+1))
+    B_snpms        = np.zeros((npts_theta,m_max+1,n_max+1))
+    ### Init
+    B_Pnpms[:,0,0] = np.sqrt(1./(4.*pi))
+    B_unpms[:,0,0] =  0.
+    B_unpms[:,1,1] = -0.25*np.sqrt(3./pi)
+    for k in range(npts_theta):
+        u=cos_theta[k]
+        for n in range(1,n_max+1):
+            B_Pnpms[k,n,n]   = -np.sqrt((2.*float(n)+1.)/(2.*float(n))) * np.sqrt(1.-u**2) * B_Pnpms[k,n-1,n-1]
+            B_Pnpms[k,n-1,n] =  u*np.sqrt(2.*float(n)+1.) * B_Pnpms[k,n-1,n-1]
+        for n in range(2,n_max+1):
+            B_unpms[k,n,n]   = -np.sqrt( (float(n)*(2.*float(n)+1.)) / (2.*(float(n)+1.)*(float(n)-1.)) ) * np.sqrt(1.-u**2) * B_unpms[k,n-1,n-1]
+            B_unpms[k,n-1,n] =  np.sqrt( (2.*float(n)+1.)*(float(n)-1.)/(float(n)+1.) ) * u * B_unpms[k,n-1,n-1]
+        for n in range(2,n_max+1):
+            for m in range(n-2+1):
+                B_Pnpms[k,m,n] = np.sqrt((4.*(float(n))**2-1.)/((float(n))**2-(float(m))**2)) * u * B_Pnpms[k,m,n-1] \
+                                        - np.sqrt( ( (2.*float(n)+1.)*((float(n)-1.)**2-(float(m))**2) ) \
+                                                            / ( (2.*float(n)-3.)*((float(n))**2-(float(m))**2) ) )*B_Pnpms[k,m,n-2]
+
+                B_unpms[k,m,n] = np.sqrt( ((4.*(float(n))**2-1.)*(float(n)-1.))/((float(n)**2-float(m)**2)*(float(n)+1.)) ) * u * B_unpms[k,m,n-1] \
+                                        - np.sqrt( ((2.*float(n)+1.) * (float(n)-1.) * (float(n)-2.) * (float(n-m)-1.) *(float(n+m)-1.)) \
+                                                 / ((2.*float(n)-3.) * (float(n)**2-float(m)**2)*float(n)*(float(n)+1.)))*B_unpms[k,m,n-2]
+        for n in range(0,n_max+1):
+            m=0
+            B_snpms[k,m,n]  = 1./float(m+1)*np.sqrt((float(n+m)+1.)*(float(n-m))) * np.sqrt(1.-u**2) *\
+                              B_unpms[k,m+1,n] + u*B_unpms[k,m,n]
+        for n in range(1,n_max+1):
+            for m in range(1,n+1):
+                B_snpms[k,m,n]  = float(n)/float(m) * u * B_unpms[k,m,n] - float(m+n)/float(m) * \
+                                    np.sqrt( ( (2.*float(n)+1.)*(float(n)-float(m))*(float(n)-1.) ) / \
+                                     ( (2.*float(n)-1.)*(float(n)+float(m))*(float(n)+1.) ) )*B_unpms[k,m,n-1]
+    B_Pnmms = np.zeros_like(B_Pnpms)
+    B_unmms = np.zeros_like(B_unpms)
+    B_snmms = np.zeros_like(B_snpms)
+    for m in range(m_max+1):
+        B_Pnmms[:,m,:]       = (-1.0)**m     * B_Pnpms[:,m,:]
+        B_unmms[:,m,:]       = (-1.0)**(m+1) * B_unpms[:,m,:]
+        B_snmms[:,m,:]       = (-1.0)**(m)   * B_snpms[:,m,:]
+    B_Pnmms=B_Pnmms[:,::-1,:]
+    B_unmms=B_unmms[:,::-1,:]
+    B_snmms=B_snmms[:,::-1,:]
+    B_Pnms       = np.concatenate((B_Pnmms,B_Pnpms[:,1::,:]),axis=1)
+    B_unms       = np.concatenate((B_unmms,B_unpms[:,1::,:]),axis=1)
+    B_snms       = np.concatenate((B_snmms,B_snpms[:,1::,:]),axis=1)
+    
+    #####################################
+    ##### sn, pn ,un
+    ##### Brian's recurrence relations
+    m_max = n_max
+    p_max = n_max*n_max+2*n_max
+    aM_nm   = np.zeros(p_max,dtype=complex)
+    bN_nm   = np.zeros(p_max,dtype=complex)
+    fenm_Y  = np.zeros(p_max,dtype=complex)
+    fenm_Z  = np.zeros(p_max,dtype=complex)
+    fhnm_X  = np.zeros(p_max,dtype=complex)
+
+    B_Pnpms        = np.zeros((npts_theta,m_max+1,n_max+1))
+    B_unpms        = np.zeros((npts_theta,m_max+1,n_max+1))
+    B_snpms        = np.zeros((npts_theta,m_max+1,n_max+1))
+    ### Init
+    B_Pnpms[:,0,0] = np.sqrt(1./(4.*pi))
+    B_unpms[:,0,0] =  0.
+    B_unpms[:,1,1] = -0.25*np.sqrt(3./pi)
+    for k in range(npts_theta):
+        u=cos_theta[k]
+        for n in range(1,n_max+1):
+            B_Pnpms[k,n,n]   = -np.sqrt((2.*float(n)+1.)/(2.*float(n))) * np.sqrt(1.-u**2) * B_Pnpms[k,n-1,n-1]
+            B_Pnpms[k,n-1,n] =  u*np.sqrt(2.*float(n)+1.) * B_Pnpms[k,n-1,n-1]
+        for n in range(2,n_max+1):
+            B_unpms[k,n,n]   = -np.sqrt( (float(n)*(2.*float(n)+1.)) / (2.*(float(n)+1.)*(float(n)-1.)) ) * np.sqrt(1.-u**2) * B_unpms[k,n-1,n-1]
+            B_unpms[k,n-1,n] =  np.sqrt( (2.*float(n)+1.)*(float(n)-1.)/(float(n)+1.) ) * u * B_unpms[k,n-1,n-1]
+        for n in range(2,n_max+1):
+            for m in range(n-2+1):
+                B_Pnpms[k,m,n] = np.sqrt((4.*(float(n))**2-1.)/((float(n))**2-(float(m))**2)) * u * B_Pnpms[k,m,n-1] \
+                                        - np.sqrt( ( (2.*float(n)+1.)*((float(n)-1.)**2-(float(m))**2) ) \
+                                                            / ( (2.*float(n)-3.)*((float(n))**2-(float(m))**2) ) )*B_Pnpms[k,m,n-2]
+
+                B_unpms[k,m,n] = np.sqrt( ((4.*(float(n))**2-1.)*(float(n)-1.))/((float(n)**2-float(m)**2)*(float(n)+1.)) ) * u * B_unpms[k,m,n-1] \
+                                        - np.sqrt( ((2.*float(n)+1.) * (float(n)-1.) * (float(n)-2.) * (float(n-m)-1.) *(float(n+m)-1.)) \
+                                                 / ((2.*float(n)-3.) * (float(n)**2-float(m)**2)*float(n)*(float(n)+1.)))*B_unpms[k,m,n-2]
+        for n in range(0,n_max+1):
+            m=0
+            B_snpms[k,m,n]  = 1./float(m+1)*np.sqrt((float(n+m)+1.)*(float(n-m))) * np.sqrt(1.-u**2) *  B_unpms[k,m+1,n] + u*B_unpms[k,m,n]
+        for n in range(1,n_max+1):
+            for m in range(1,n+1):
+                B_snpms[k,m,n]  = float(n)/float(m) * u * B_unpms[k,m,n] - float(m+n)/float(m) * \
+                                    np.sqrt( ( (2.*float(n)+1.)*(float(n)-float(m))*(float(n)-1.) ) / ( (2.*float(n)-1.)*(float(n)+float(m))*(float(n)+1.) ) )*B_unpms[k,m,n-1]
+    B_Pnmms = np.zeros_like(B_Pnpms)
+    B_unmms = np.zeros_like(B_unpms)
+    B_snmms = np.zeros_like(B_snpms)
+    for m in range(m_max+1):
+        B_Pnmms[:,m,:]       = (-1.0)**m     * B_Pnpms[:,m,:]
+        B_unmms[:,m,:]       = (-1.0)**(m+1) * B_unpms[:,m,:]
+        B_snmms[:,m,:]       = (-1.0)**(m)   * B_snpms[:,m,:]
+    B_Pnmms=B_Pnmms[:,::-1,:]
+    B_unmms=B_unmms[:,::-1,:]
+    B_snmms=B_snmms[:,::-1,:]
+    B_Pnms       = np.concatenate((B_Pnmms,B_Pnpms[:,1::,:]),axis=1)
+    B_unms       = np.concatenate((B_unmms,B_unpms[:,1::,:]),axis=1)
+    B_snms       = np.concatenate((B_snmms,B_snpms[:,1::,:]),axis=1)
+    E_scat_onsphere_sph = np.array( [np.loadtxt(post_filename,usecols=[8])
+                                        + 1j*np.loadtxt(post_filename,usecols=[11]),
+                                        np.loadtxt(post_filename,usecols=[9])
+                                        + 1j*np.loadtxt(post_filename,usecols=[12]),
+                                        np.loadtxt(post_filename,usecols=[10])
+                                        + 1j*np.loadtxt(post_filename,usecols=[13])])
+    E_scat_onsphere_sph_r  = E_scat_onsphere_sph[0,:].reshape(npts_phi,npts_theta,order='F').transpose()
+    E_scat_onsphere_sph_t  = E_scat_onsphere_sph[1,:].reshape(npts_phi,npts_theta,order='F').transpose()
+    E_scat_onsphere_sph_p  = E_scat_onsphere_sph[2,:].reshape(npts_phi,npts_theta,order='F').transpose()
+    for ko in range(p_max):
+        po = ps[ko]
+        n  = int(np.sqrt(po))
+        m  = n*(n+1) - int(po)
+        # print('=========>> po',po,'n',n,'m',m)
+        B_PnmN_costheta   = np.tile( B_Pnms[:,m+m_max,n],(npts_phi,1)).transpose()
+        B_UnmN_costheta   = np.tile( B_unms[:,m+m_max,n],(npts_phi,1)).transpose()
+        B_SnmN_costheta   = np.tile( B_snms[:,m+m_max,n],(npts_phi,1)).transpose()
+        B_Ynm_r           =       B_PnmN_costheta * np.exp(1j*float(m)*phi_sph)
+        B_Ynm_t           =          np.zeros_like(phi_sph)
+        B_Ynm_p           =          np.zeros_like(phi_sph)
+        B_Xnm_r           =          np.zeros_like(phi_sph)
+        B_Xnm_t           = 1j  * B_UnmN_costheta * np.exp(1j*float(m)*phi_sph)
+        B_Xnm_p           = -1. * B_SnmN_costheta * np.exp(1j*float(m)*phi_sph)
+        B_Znm_r           =          np.zeros_like(phi_sph)
+        B_Znm_t           =  1. * B_SnmN_costheta * np.exp(1j*float(m)*phi_sph)
+        B_Znm_p           = 1j  * B_UnmN_costheta * np.exp(1j*float(m)*phi_sph)
+        B_erCrossXnm_r    = np.zeros_like(phi_sph,dtype=complex)
+        B_erCrossXnm_t    = -B_Xnm_p
+        B_erCrossXnm_p    =  B_Xnm_t
+        FF_Xnm_t[:,:,ko]  = B_Xnm_t
+        FF_Xnm_p[:,:,ko]  = B_Xnm_p
+        FF_erCrossXnm_t[:,:,ko] = B_erCrossXnm_t
+        FF_erCrossXnm_p[:,:,ko] = B_erCrossXnm_p
+        
+        
+        sph_bessel_n_ofkr       = np.sqrt(pi/(2.*k_Out*r_sph))*outgoing_sph_hankel(float(n  )+0.5,k_Out*r_sph)
+        sph_bessel_nminus1_ofkr = np.sqrt(pi/(2.*k_Out*r_sph))*outgoing_sph_hankel(float(n-1)+0.5,k_Out*r_sph)
+        dRicatti_dx_ofkr        = (k_Out * r_sph * (sph_bessel_nminus1_ofkr-(float(n+1)/((k_Out*r_sph))) * sph_bessel_n_ofkr) + sph_bessel_n_ofkr)
+        B_Mnm_r = 0.
+        B_Mnm_t = sph_bessel_n_ofkr * B_Xnm_t
+        B_Mnm_p = sph_bessel_n_ofkr * B_Xnm_p
+        B_Nnm_r = 1./(k_Out*r_sph) * np.sqrt(float(n*(n+1))) * sph_bessel_n_ofkr * B_Ynm_r
+        B_Nnm_t = 1./(k_Out*r_sph) * dRicatti_dx_ofkr * B_Znm_t
+        B_Nnm_p = 1./(k_Out*r_sph) * dRicatti_dx_ofkr * B_Znm_p
+
+        B_EdotconjYnm           = E_scat_onsphere_sph_r*B_Ynm_r.conjugate()
+        B_EdotconjZnm           = E_scat_onsphere_sph_t*B_Znm_t.conjugate() + E_scat_onsphere_sph_p*B_Znm_p.conjugate()
+        B_EdotconjXnm           = E_scat_onsphere_sph_t*B_Xnm_t.conjugate() + E_scat_onsphere_sph_p*B_Xnm_p.conjugate()
+        normalize_fhnm_X        = 1./sph_bessel_n_ofkr
+        normalize_fenm_Y        = k_Out*r_sph/(sph_bessel_n_ofkr*np.sqrt(float(n)*(float(n)+1.)) )
+        normalize_fenm_Z        = k_Out*r_sph/dRicatti_dx_ofkr
+        fenm_Y[int(po)-1]     = np.trapz(np.trapz((np.sin(theta_sph)*B_EdotconjYnm).transpose(),theta_sph[:,0]),phi_sph[0,:])*normalize_fenm_Y
+        fenm_Z[int(po)-1]     = np.trapz(np.trapz((np.sin(theta_sph)*B_EdotconjZnm).transpose(),theta_sph[:,0]),phi_sph[0,:])*normalize_fenm_Z
+        fhnm_X[int(po)-1]     = np.trapz(np.trapz((np.sin(theta_sph)*B_EdotconjXnm).transpose(),theta_sph[:,0]),phi_sph[0,:])*normalize_fhnm_X
+
+        EdotconjMnm = E_scat_onsphere_sph_r*B_Mnm_r.conjugate() + E_scat_onsphere_sph_t*B_Mnm_t.conjugate() + E_scat_onsphere_sph_p*B_Mnm_p.conjugate()
+        EdotconjNnm = E_scat_onsphere_sph_r*B_Nnm_r.conjugate() + E_scat_onsphere_sph_t*B_Nnm_t.conjugate() + E_scat_onsphere_sph_p*B_Nnm_p.conjugate()
+        MnmconjMnm  = B_Mnm_r*B_Mnm_r.conjugate() + B_Mnm_t*B_Mnm_t.conjugate() + B_Mnm_p*B_Mnm_p.conjugate()
+        NnmconjNnm  = B_Nnm_r*B_Nnm_r.conjugate() + B_Nnm_t*B_Nnm_t.conjugate() + B_Nnm_p*B_Nnm_p.conjugate()
+        MnmconjNnm  = B_Mnm_r*B_Nnm_r.conjugate() + B_Mnm_t*B_Nnm_t.conjugate() + B_Mnm_p*B_Nnm_p.conjugate()
+        XnmconjXnm  = B_Xnm_r*B_Xnm_r.conjugate() + B_Xnm_t*B_Xnm_t.conjugate() + B_Xnm_p*B_Xnm_p.conjugate()
+        YnmconjYnm  = B_Ynm_r*B_Ynm_r.conjugate() + B_Ynm_t*B_Ynm_t.conjugate() + B_Ynm_p*B_Ynm_p.conjugate()
+        ZnmconjZnm  = B_Znm_r*B_Znm_r.conjugate() + B_Znm_t*B_Znm_t.conjugate() + B_Znm_p*B_Znm_p.conjugate()
+        XnmconjYnm  = B_Xnm_r*B_Ynm_r.conjugate() + B_Xnm_t*B_Ynm_t.conjugate() + B_Xnm_p*B_Ynm_p.conjugate()
+        YnmconjZnm  = B_Ynm_r*B_Znm_r.conjugate() + B_Ynm_t*B_Znm_t.conjugate() + B_Ynm_p*B_Znm_p.conjugate()
+        ZnmconjXnm  = B_Znm_r*B_Xnm_r.conjugate() + B_Znm_t*B_Xnm_t.conjugate() + B_Znm_p*B_Xnm_p.conjugate()
+        normalize_aM_nm2 = np.trapz(np.trapz((np.sin(theta_sph)*MnmconjMnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+        normalize_bN_nm2 = np.trapz(np.trapz((np.sin(theta_sph)*NnmconjNnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+        orth             = np.trapz(np.trapz((np.sin(theta_sph)*MnmconjNnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+        orthX            = np.trapz(np.trapz((np.sin(theta_sph)*XnmconjXnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+        orthY            = np.trapz(np.trapz((np.sin(theta_sph)*YnmconjYnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+        orthZ            = np.trapz(np.trapz((np.sin(theta_sph)*ZnmconjZnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+        orth1            = np.trapz(np.trapz((np.sin(theta_sph)*XnmconjYnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+        orth2            = np.trapz(np.trapz((np.sin(theta_sph)*YnmconjZnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+        orth3            = np.trapz(np.trapz((np.sin(theta_sph)*ZnmconjXnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+        aM_nm[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*EdotconjMnm).transpose(),theta_sph[:,0]),phi_sph[0,:]) / normalize_aM_nm2
+        bN_nm[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*EdotconjNnm).transpose(),theta_sph[:,0]),phi_sph[0,:]) / normalize_bN_nm2
+    return [fhnm_X,fenm_Y,fenm_Z,aM_nm,bN_nm,FF_Xnm_t,FF_Xnm_p,FF_erCrossXnm_t,FF_erCrossXnm_p]
+
+###########################
+def plot_farfield(far_field_sph,filename):
+    vmax_sph = np.max(far_field_sph)
+    vmin_sph = np.min(far_field_sph)
+    x_sph = far_field_sph/vmax_sph*np.sin(theta_sph) * np.cos(phi_sph)
+    y_sph = far_field_sph/vmax_sph*np.sin(theta_sph) * np.sin(phi_sph)
+    z_sph = far_field_sph/vmax_sph*np.cos(theta_sph)
+    fig = pl.figure()
+    ax = fig.add_subplot(111, projection='3d')
+    ax.set_aspect('equal')
+    surf=ax.plot_surface(x_sph,y_sph,z_sph,\
+                    facecolors=cm.viridis(far_field_sph/vmax_sph),\
+                    rstride=2,\
+                    cstride=2,\
+                    linewidth=0,\
+                    vmin = vmin_sph,\
+                    vmax = vmax_sph,\
+                    shade=True,\
+                    alpha=0.5,\
+                    antialiased=False)
+    cset = ax.contourf(x_sph, y_sph, z_sph,1,fc='k', zdir='z', offset=-1)
+    cset = ax.contourf(x_sph, y_sph, z_sph,1,fc='k', zdir='x', offset=-1)
+    cset = ax.contourf(x_sph, y_sph, z_sph,1,fc='k', zdir='y', offset=1 )
+    surf.set_edgecolor('k')
+    max_range = 0.5*np.max(np.array([x_sph.max()-x_sph.min(), y_sph.max()-y_sph.min(), z_sph.max()-z_sph.min()]))
+    mid_x = (x_sph.max()+x_sph.min())*0.5
+    mid_y = (y_sph.max()+y_sph.min())*0.5
+    mid_z = (z_sph.max()+z_sph.min())*0.5
+    ax.set_xlim(mid_x-max_range, mid_x+max_range)
+    ax.set_ylim(mid_y-max_range, mid_y+max_range)
+    ax.set_zlim(mid_z-max_range, mid_z+max_range)
+    ax.xaxis.set_ticklabels([]);ax.xaxis.set_label('x')
+    ax.yaxis.set_ticklabels([]);ax.yaxis.set_label('y')
+    ax.zaxis.set_ticklabels([]);ax.zaxis.set_label('z')
+    pl.savefig(filename,bbox_inches='tight')
+    pl.close('all')

ElectromagneticScattering/scattering_post.py

+"""
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattering_post.py        //
+///////////////////////////////
+"""
+
+Tmatrix_rfull_ab  = np.zeros((2*p_max,2*p_max,nb_cuts),dtype=complex)
+Tmatrix_rfull_fy  = np.zeros((2*p_max,2*p_max,nb_cuts),dtype=complex)
+Tmatrix_rfull_fz  = np.zeros((2*p_max,2*p_max,nb_cuts),dtype=complex)
+tab_E_NTF_t_G = np.zeros((npts_theta,npts_phi,3),dtype=complex)
+tab_E_NTF_p_G = np.zeros((npts_theta,npts_phi,3),dtype=complex)
+tab_E_NTF_t_PW = np.zeros((npts_theta,npts_phi),dtype=complex)
+tab_E_NTF_p_PW = np.zeros((npts_theta,npts_phi),dtype=complex)
+
+my_dir='./run_results/'
+if flag_study==0:
+    nbNM = 1
+    p_max_in=1
+elif flag_study==1:
+    nbNM = 2
+    p_max_in=p_max
+elif flag_study==2:
+    nbNM = 1
+    p_max_in=p_max
+    nb_cuts = 3
+elif flag_study==3:
+    p_max_in=0
+
+if siwt==1:
+    outgoing_sph_hankel = hankel2
+else:
+    outgoing_sph_hankel = hankel1
+
+print('Postprocessing...')
+for ke in range(p_max_in):
+    for isN in range(nbNM):
+        pe = ps[ke]
+        ne =int(np.sqrt(pe))
+        me = ne*(ne+1) - int(pe)
+        aM_nm   = np.zeros((nb_cuts,p_max),dtype=complex)
+        bN_nm   = np.zeros((nb_cuts,p_max),dtype=complex)
+        fenm_Y  = np.zeros((nb_cuts,p_max),dtype=complex)
+        fenm_Z  = np.zeros((nb_cuts,p_max),dtype=complex)
+        fhnm_X  = np.zeros((nb_cuts,p_max),dtype=complex)
+        for nr in range(nb_cuts):
+            r_sph      = r_sphs[nr]
+            # print('======>> postprocessing r_sph',r_sph)
+            par_gmsh_getdp=open('parameters_gmsh_getdp.dat', 'a')
+            par_gmsh_getdp.write('r_sph = %3.15e;\n'  %(r_sph) )
+            par_gmsh_getdp.close()
+            if flag_study==1:            
+                if isN==0:post_filename = my_dir+'E_scat_onsphere_sph_M%g_r%g.dat'%((pe,nr))
+                else:     post_filename = my_dir+'E_scat_onsphere_sph_N%g_r%g.dat'%((pe,nr))
+            elif flag_study==0:
+                post_filename = my_dir+'E_scat_onsphere_sph_PW_r%g.dat'%(nr)
+            elif flag_study==2:
+                post_filename = my_dir+'E_tot_onsphere_sph_G%g.dat'%(nr)
+            res=field_VSH_expansion(post_filename)
+            fhnm_X[nr,:]    = res[0]
+            fenm_Y[nr,:]    = res[1]
+            fenm_Z[nr,:]    = res[2]
+            aM_nm[nr,:]     = res[3]
+            bN_nm[nr,:]     = res[4]
+            FF_Xnm_t        = res[5]
+            FF_Xnm_p        = res[6]
+            FF_erCrossXnm_t = res[7]
+            FF_erCrossXnm_p = res[8]
+        if flag_study==1:
+            if isN==0:
+                Tmatrix_rfull_ab[ke,0:p_max,:]      =  aM_nm.transpose()
+                Tmatrix_rfull_ab[ke,p_max:2*p_max,:] =  bN_nm.transpose()
+                Tmatrix_rfull_fy[ke,0:p_max,:]      = fhnm_X.transpose()
+                Tmatrix_rfull_fy[ke,p_max:2*p_max,:] = fenm_Y.transpose()
+                Tmatrix_rfull_fz[ke,0:p_max,:]      = fhnm_X.transpose()
+                Tmatrix_rfull_fz[ke,p_max:2*p_max,:] = fenm_Z.transpose()
+            if isN==1:
+                Tmatrix_rfull_ab[p_max+ke,0:p_max,:]      =  aM_nm.transpose()
+                Tmatrix_rfull_ab[p_max+ke,p_max:2*p_max,:] =  bN_nm.transpose()
+                Tmatrix_rfull_fy[p_max+ke,0:p_max,:]      = fhnm_X.transpose()
+                Tmatrix_rfull_fy[p_max+ke,p_max:2*p_max,:] = fenm_Y.transpose()
+                Tmatrix_rfull_fz[p_max+ke,0:p_max,:]      = fhnm_X.transpose()
+                Tmatrix_rfull_fz[p_max+ke,p_max:2*p_max,:] = fenm_Z.transpose()
+            Tmatrix_ab     = np.mean(Tmatrix_rfull_ab,axis=2)
+            Tmatrix_fy     = np.mean(Tmatrix_rfull_fy,axis=2)
+            Tmatrix_fz     = np.mean(Tmatrix_rfull_fz,axis=2)
+            std_Tmatrix_ab = np.std(Tmatrix_rfull_ab,axis=2)
+            std_Tmatrix_fy = np.std(Tmatrix_rfull_fy,axis=2)
+            std_Tmatrix_fz = np.std(Tmatrix_rfull_fz,axis=2)
+            np.savez('T_matrix_ellips_sphPML_epsre_%.2lf_eps_im_%.2lf.npz'%(epsr_In.real,epsr_In.imag), \
+                    Tmatrix_rfull_ab = Tmatrix_rfull_ab, \
+                    Tmatrix_rfull_fy = Tmatrix_rfull_fy, \
+                    Tmatrix_rfull_fz = Tmatrix_rfull_fz,\
+                    r_sphs           = r_sphs,\
+                    n_max            = n_max,\
+                    ps               = ps,\
+                    p_max            = p_max,\
+                    lambda0          = lambda0, \
+                    Tmatrix_ab       = Tmatrix_ab, \
+                    Tmatrix_fy       = Tmatrix_fy, \
+                    Tmatrix_fz       = Tmatrix_fz, \
+                    std_Tmatrix_ab   = std_Tmatrix_ab,\
+                    std_Tmatrix_fy   = std_Tmatrix_fy,\
+                    std_Tmatrix_fz   = std_Tmatrix_fz)
+        if flag_study==0:
+            mean_aM_nm  = np.mean(aM_nm ,axis=0)
+            mean_bN_nm  = np.mean(bN_nm ,axis=0)
+            mean_fenm_Y = np.mean(fenm_Y,axis=0)
+            mean_fenm_Z = np.mean(fenm_Z,axis=0)
+            mean_fhnm_X = np.mean(fhnm_X,axis=0)
+            std_aM_nm   = np.std( aM_nm ,axis=0)
+            std_bN_nm   = np.std( bN_nm ,axis=0)
+            std_fenm_Y  = np.std( fenm_Y,axis=0)
+            std_fenm_Z  = np.std( fenm_Z,axis=0)
+            std_fhnm_X  = np.std( fhnm_X,axis=0)
+            for ko in range(p_max):
+                tab_E_NTF_t_PW += (1j)**ko * mean_aM_nm[ko]*1j*FF_Xnm_t[:,:,ko] + (1j)**ko * mean_bN_nm[ko]*FF_erCrossXnm_t[:,:,ko]
+                tab_E_NTF_p_PW += (1j)**ko * mean_aM_nm[ko]*1j*FF_Xnm_p[:,:,ko] + (1j)**ko * mean_bN_nm[ko]*FF_erCrossXnm_p[:,:,ko]
+            I_NTF_PW = (np.abs(tab_E_NTF_t_PW))**2 + (np.abs(tab_E_NTF_p_PW))**2
+            plot_farfield(I_NTF_PW[:,:],my_dir+'farfield_PW.pdf')
+        if flag_study==2:
+            tens_Green_scat=np.zeros((3,3))
+            for nr in range(nb_cuts):
+                tens_Green_scat[:,nr] = np.loadtxt(my_dir+'E_scat_im_onpt_G%g.dat'%(nr))[8:11]
+                for ko in range(p_max):
+                    tab_E_NTF_t_G[:,:,nr] += (1j)**ko * aM_nm[nr,ko]*1j*FF_Xnm_t[:,:,ko] + (1j)**ko * bN_nm[nr,ko]*FF_erCrossXnm_t[:,:,ko]
+                    tab_E_NTF_p_G[:,:,nr] += (1j)**ko * aM_nm[nr,ko]*1j*FF_Xnm_p[:,:,ko] + (1j)**ko * bN_nm[nr,ko]*FF_erCrossXnm_p[:,:,ko]
+            I_NTF_G = (np.abs(tab_E_NTF_t_G))**2 + (np.abs(tab_E_NTF_p_G))**2
+            tens_Green_inc = -siwt*np.eye(3)*k_Out/(6.*pi)
+            tens_Green_tot = tens_Green_scat+tens_Green_inc
+            for k in range(3):plot_farfield(I_NTF_G[:,:,k],my_dir+'farfield_green%g.pdf'%(k))
+
+if flag_study==1:
+    print('Tmatrix_ab:')
+    print(Tmatrix_ab)
+    print('std_Tmatrix_ab:')
+    print(std_Tmatrix_ab)
+    print('Tmatrix_fy:')
+    print(Tmatrix_fy)
+    print('np.abs(2*Tmatrix_fy+1):')
+    print(np.abs(2*Tmatrix_fy+1))
+    print('Tmatrix_fz:')
+    print(Tmatrix_fz)
+    print('np.abs(2*Tmatrix_fz+1):')
+    print(np.abs(2*Tmatrix_fz+1))
+    print('np.abs(2*Tmatrix_ab+1):')
+    print(np.abs(2*Tmatrix_ab+1))
+elif flag_study==0:
+    print('anm:')
+    print(mean_aM_nm)
+    print('bnm')
+    print(mean_bN_nm)
+elif flag_study==2:
+    print('tens_Green_tot')
+    print(tens_Green_tot)
+    print('tens_Green_tot/G0')
+    print(tens_Green_tot/tens_Green_inc[0,0])
+    
+elif flag_study==3:
+    eigs = np.loadtxt(my_dir+'EigenValuesReal.txt',usecols=[5]) + 1j*np.loadtxt(my_dir+'EigenValuesImag.txt',usecols=[5])
+    eigs *= (nm/1e-9)**2
+    pl.savez(my_dir+'eigs.npz',eigs=eigs)
+    pl.figure(figsize=(20,15))
+    pl.plot( eigs.real  ,   eigs.imag,'+',mfc='none')
+    pl.grid()
+    pl.savefig(my_dir+'eigenvalues.pdf')
\ No newline at end of file

NonLinearEVP/NonLinearEVP.geo

+///////////////////////////////////
+//// Author : Guillaume Demesy ////
+////        .geo               ////
+///////////////////////////////////
+
+Include "NonLinearEVP_data.geo";
+lc_cell = a_lat/paramaille;
+lc_sq   = lc_cell/4;
+lc_pml  = lc_cell*1.2;
+lc_sqa  = lc_sq;
+
+lc_sq_inside = lc_sq*1.4;
+epsc         = lc_sq*0.8;
+
+
+If (flag_Tmesh==0)
+  Point(1)  = {-a_lat/2. ,-d_sq/2.-space2pml, 0. , lc_cell};
+  Point(2)  = { a_lat/2. ,-d_sq/2.-space2pml, 0. , lc_cell};
+  Point(3)  = { a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell};
+  Point(4)  = {-a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell};
+
+  Point(5)  = {-d_sq/2. ,-d_sq/2., 0. , lc_sq};
+  Point(6)  = { d_sq/2. ,-d_sq/2., 0. , lc_sq};
+  Point(7)  = { d_sq/2. , d_sq/2., 0. , lc_sq};
+  Point(8)  = {-d_sq/2. , d_sq/2., 0. , lc_sq};
+
+  Point(9)  = {-a_lat/2. ,-d_sq/2.-space2pml-pmlsize, 0. , lc_pml};
+  Point(10) = { a_lat/2. ,-d_sq/2.-space2pml-pmlsize, 0. , lc_pml};
+  Point(11) = { a_lat/2. , d_sq/2.+space2pml+pmlsize, 0. , lc_pml};
+  Point(12) = {-a_lat/2. , d_sq/2.+space2pml+pmlsize, 0. , lc_pml};
+
+  Point(13) = {-a_lat/2. , a_lat/2., 0. , lc_sqa};
+  Point(14) = {-a_lat/2. ,-a_lat/2., 0. , lc_sqa};
+  Point(15) = { a_lat/2. , a_lat/2., 0. , lc_sqa};
+  Point(16) = { a_lat/2. ,-a_lat/2., 0. , lc_sqa};
+
+  Line(1) = {1, 2};
+  Line(3) = {3, 4};
+  Line(5) = {5, 6};
+  Line(6) = {6, 7};
+  Line(7) = {7, 8};
+  Line(8) = {8, 5};
+
+  Line(15) = {1, 14};
+  Line(16) = {14, 13};
+  Line(17) = {13, 4};
+  Line(18) = {2, 16};
+  Line(19) = {16, 15};
+  Line(20) = {15, 3};
+
+  Line(9) = {4, 12};
+  Line(10) = {12, 11};
+  Line(11) = {11, 3};
+  Line(12) = {1, 9};
+  Line(13) = {9, 10};
+  Line(14) = {10, 2};
+
+  Line Loop(1) = {12, 13, 14, -1};
+  Plane Surface(1) = {1};
+  Line Loop(2) = {5, 6, 7, 8};
+  Plane Surface(2) = {2};
+  Line Loop(3) = {17, -3, -20, -19, -18, -1, 15, 16};
+  Plane Surface(3) = {2, -3};
+  Line Loop(4) = {9, 10, 11, 3};
+  Plane Surface(4) = {-4};
+
+  Periodic Line { 14,18,19,20,11 } = {12,15,16,17,9 } Translate {a_lat,0,0} ; 
+
+  Physical Surface(100) = {2};   // 1 dom in
+  Physical Surface(101) = {3};  // 2 dom out
+  Physical Surface(102) = {1};  // PML bot
+  Physical Surface(103) = {4};  // PML top
+  Physical Line(1001)   = {12,15,16,17,9}; // bloch x left
+  Physical Line(1002)   = {14,18,19,20,11}; // bloch x right
+  Physical Line(1003)   = {10}; // top bound
+  Physical Line(1004)   = {13}; // bot bound
+  Physical Line(1005)   = {5,6,7,8}; // bound for lag
+  Physical Point(10000) = {1};   // Printpoint
+EndIf
+If (flag_Tmesh==1)
+
+  Point(1)  = {-a_lat/2. ,-d_sq/2.-space2pml, 0. , lc_cell};
+  Point(2)  = { a_lat/2. ,-d_sq/2.-space2pml, 0. , lc_cell};
+  Point(3)  = { a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell};
+  Point(4)  = {-a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell};
+
+  Point(5)  = {-d_sq/2. ,-d_sq/2., 0. , lc_sq*15};
+  Point(6)  = { d_sq/2. ,-d_sq/2., 0. , lc_sq*15};
+  Point(7)  = { d_sq/2. , d_sq/2., 0. , lc_sq*15};
+  Point(8)  = {-d_sq/2. , d_sq/2., 0. , lc_sq*15};
+
+  Point(9)  = {-a_lat/2. ,-d_sq/2.-space2pml-pmlsize, 0. , lc_pml};
+  Point(10) = { a_lat/2. ,-d_sq/2.-space2pml-pmlsize, 0. , lc_pml};
+  Point(11) = { a_lat/2. , d_sq/2.+space2pml+pmlsize, 0. , lc_pml};
+  Point(12) = {-a_lat/2. , d_sq/2.+space2pml+pmlsize, 0. , lc_pml};
+
+  Point(13) = {-d_sq/2.-epsc , d_sq/2.+epsc, 0. , lc_sqa};
+  Point(14) = {-d_sq/2.      , d_sq/2.+epsc, 0. , lc_sqa};
+  Point(15) = {-d_sq/2.+epsc , d_sq/2.+epsc, 0. , lc_sqa};
+  Point(16) = { d_sq/2.-epsc , d_sq/2.+epsc, 0. , lc_sqa};
+  Point(17) = { d_sq/2.      , d_sq/2.+epsc, 0. , lc_sqa};
+  Point(18) = { d_sq/2.+epsc , d_sq/2.+epsc, 0. , lc_sqa};
+  Point(19) = { d_sq/2.+epsc , d_sq/2.     , 0. , lc_sqa};
+  Point(20) = { d_sq/2.+epsc , d_sq/2.-epsc, 0. , lc_sqa};
+  Point(21) = { d_sq/2.+epsc ,-d_sq/2.+epsc, 0. , lc_sqa};
+  Point(22) = { d_sq/2.+epsc ,-d_sq/2.     , 0. , lc_sqa};
+  Point(23) = { d_sq/2.+epsc ,-d_sq/2.-epsc, 0. , lc_sqa};
+  Point(24) = { d_sq/2.      ,-d_sq/2.-epsc, 0. , lc_sqa};
+  Point(25) = { d_sq/2.-epsc ,-d_sq/2.-epsc, 0. , lc_sqa};
+  Point(26) = {-d_sq/2.+epsc ,-d_sq/2.-epsc, 0. , lc_sqa};
+  Point(27) = {-d_sq/2.      ,-d_sq/2.-epsc, 0. , lc_sqa};
+  Point(28) = {-d_sq/2.-epsc ,-d_sq/2.-epsc, 0. , lc_sqa};
+  Point(29) = {-d_sq/2.-epsc ,-d_sq/2.     , 0. , lc_sqa};
+  Point(30) = {-d_sq/2.-epsc ,-d_sq/2.+epsc, 0. , lc_sqa};
+  Point(31) = {-d_sq/2.-epsc , d_sq/2.-epsc, 0. , lc_sqa};
+  Point(32) = {-d_sq/2.-epsc , d_sq/2.     , 0. , lc_sqa};
+
+  Point(33)  = {-d_sq/2.+epsc ,-d_sq/2.+epsc, 0. , lc_sqa};
+  Point(34)  = { d_sq/2.-epsc ,-d_sq/2.+epsc, 0. , lc_sqa};
+  Point(35)  = { d_sq/2.-epsc , d_sq/2.-epsc, 0. , lc_sqa};
+  Point(36)  = {-d_sq/2.+epsc , d_sq/2.-epsc, 0. , lc_sqa};
+
+  Point(37)  = {-d_sq/2.+epsc ,-d_sq/2., 0. , lc_sqa};
+  Point(38)  = { d_sq/2.-epsc ,-d_sq/2., 0. , lc_sqa};
+  Point(39)  = { d_sq/2.-epsc , d_sq/2., 0. , lc_sqa};
+  Point(45)  = {-d_sq/2.+epsc , d_sq/2., 0. , lc_sqa};
+  Point(46)  = {-d_sq/2. ,-d_sq/2.+epsc, 0. , lc_sqa};
+  Point(47)  = { d_sq/2. ,-d_sq/2.+epsc, 0. , lc_sqa};
+  Point(48)  = { d_sq/2. , d_sq/2.-epsc, 0. , lc_sqa};
+  Point(49)  = {-d_sq/2. , d_sq/2.-epsc, 0. , lc_sqa};
+
+  Point(54) = {-a_lat/2. , a_lat/2., 0. , lc_sqa};
+  Point(55) = {-a_lat/2. ,-a_lat/2., 0. , lc_sqa};
+  Point(56) = { a_lat/2. , a_lat/2., 0. , lc_sqa};
+  Point(58) = { a_lat/2. ,-a_lat/2., 0. , lc_sqa};
+
+  Line(1) = {1, 2};
+  Line(3) = {3, 4};
+  Line(9) = {4, 12};
+  Line(10) = {12, 11};
+  Line(11) = {11, 3};
+  Line(12) = {1, 9};
+  Line(13) = {9, 10};
+  Line(14) = {10, 2};
+
+  Line(35) = {8, 45};
+  Line(36) = {45, 45};
+  Line(37) = {39, 39};
+  Line(38) = {45, 39};
+  Line(39) = {39, 7};
+  Line(40) = {7, 48};
+  Line(41) = {48, 47};
+  Line(42) = {47, 6};
+  Line(43) = {6, 38};
+  Line(44) = {38, 37};
+  Line(45) = {37, 5};
+  Line(46) = {5, 46};
+  Line(47) = {46, 49};
+  Line(48) = {49, 8};
+  Line(49) = {31, 49};
+  Line(50) = {49, 36};
+  Line(51) = {36, 45};
+  Line(52) = {45, 15};
+  Line(53) = {32, 8};
+  Line(54) = {8, 14};
+  Line(55) = {13, 8};
+  Line(56) = {8, 36};
+  Line(57) = {31, 8};
+  Line(58) = {8, 15};
+  Line(59) = {16, 39};
+  Line(60) = {39, 35};
+  Line(61) = {35, 48};
+  Line(62) = {48, 20};
+  Line(63) = {19, 7};
+  Line(64) = {7, 17};
+  Line(65) = {16, 7};
+  Line(66) = {7, 20};
+  Line(67) = {35, 7};
+  Line(68) = {7, 18};
+  Line(69) = {34, 47};
+  Line(70) = {47, 21};
+  Line(71) = {34, 38};
+  Line(72) = {38, 25};
+  Line(73) = {25, 6};
+  Line(74) = {6, 24};
+  Line(75) = {6, 22};
+  Line(76) = {6, 21};
+  Line(77) = {23, 6};
+  Line(78) = {6, 34};
+  Line(79) = {30, 46};
+  Line(80) = {46, 33};
+  Line(81) = {33, 37};
+  Line(82) = {37, 26};
+  Line(83) = {26, 5};
+  Line(84) = {5, 27};
+  Line(85) = {5, 28};
+  Line(86) = {5, 29};
+  Line(87) = {5, 30};
+  Line(88) = {5, 33};
+  Line(89) = {33, 36};
+  Line(90) = {36, 35};
+  Line(91) = {35, 34};
+  Line(92) = {34, 33};
+  Line(93) = {13, 14};
+  Line(94) = {14, 15};
+  Line(95) = {15, 16};
+  Line(96) = {16, 17};
+  Line(97) = {17, 18};
+  Line(98) = {18, 19};
+  Line(99) = {19, 20};
+  Line(100) = {20, 21};
+  Line(101) = {21, 22};
+  Line(102) = {22, 23};
+  Line(103) = {23, 24};
+  Line(104) = {24, 25};
+  Line(105) = {25, 26};
+  Line(106) = {26, 27};
+  Line(107) = {27, 28};
+  Line(108) = {28, 29};
+  Line(109) = {29, 30};
+  Line(110) = {30, 31};
+  Line(111) = {31, 32};
+  Line(112) = {32, 13};
+  Line(113) = {1, 55};
+  Line(114) = {55, 54};
+  Line(115) = {54, 4};
+  Line(116) = {2, 58};
+  Line(117) = {58, 56};
+  Line(118) = {56, 3};
+
+  Line Loop(1) = {53, -55, -112};
+  Plane Surface(1) = {1};
+  Line Loop(2) = {55, 54, -93};
+  Plane Surface(2) = {2};
+  Line Loop(3) = {54, 94, -58};
+  Plane Surface(3) = {-3};
+  Line Loop(4) = {35, 52, -58};
+  Plane Surface(4) = {4};
+  Line Loop(5) = {111, 53, -57};
+  Plane Surface(5) = {-5};
+  Line Loop(6) = {49, 48, -57};
+  Plane Surface(6) = {6};
+  Line Loop(7) = {48, 56, -50};
+  Plane Surface(7) = {-7};
+  Line Loop(8) = {56, 51, -35};
+  Plane Surface(8) = {8};
+  Line Loop(9) = {51, 38, 60, -90};
+  Plane Surface(9) = {-9};
+  Line Loop(10) = {38, -59, -95, -52};
+  Plane Surface(10) = {10};
+  Line Loop(11) = {59, 39, -65};
+  Plane Surface(11) = {11};
+  Line Loop(12) = {65, 64, -96};
+  Plane Surface(12) = {12};
+  Line Loop(13) = {64, 97, -68};
+  Plane Surface(13) = {-13};
+  Line Loop(14) = {68, 98, 63};
+  Plane Surface(14) = {-14};
+  Line Loop(15) = {63, 66, -99};
+  Plane Surface(15) = {15};
+  Line Loop(16) = {66, -62, -40};
+  Plane Surface(16) = {-16};
+  Line Loop(17) = {40, -61, 67};
+  Plane Surface(17) = {-17};
+  Line Loop(18) = {67, -39, 60};
+  Plane Surface(18) = {18};
+  Line Loop(19) = {91, 69, -41, -61};
+  Plane Surface(19) = {19};
+  Line Loop(20) = {70, -100, -62, 41};
+  Plane Surface(20) = {20};
+  Line Loop(21) = {42, 76, -70};
+  Plane Surface(21) = {21};
+  Line Loop(22) = {76, 101, -75};
+  Plane Surface(22) = {-22};
+  Line Loop(23) = {42, 78, 69};
+  Plane Surface(23) = {-23};
+  Line Loop(24) = {71, -43, 78};
+  Plane Surface(24) = {24};
+  Line Loop(25) = {72, 73, 43};
+  Plane Surface(25) = {25};
+  Line Loop(26) = {73, 74, 104};
+  Plane Surface(26) = {-26};
+  Line Loop(27) = {74, -103, 77};
+  Plane Surface(27) = {27};
+  Line Loop(28) = {77, 75, 102};
+  Plane Surface(28) = {-28};
+  Line Loop(29) = {92, 81, -44, -71};
+  Plane Surface(29) = {29};
+  Line Loop(30) = {82, -105, -72, 44};
+  Plane Surface(30) = {30};
+  Line Loop(31) = {83, -45, 82};
+  Plane Surface(31) = {-31};
+  Line Loop(32) = {83, 84, -106};
+  Plane Surface(32) = {32};
+  Line Loop(33) = {84, 107, -85};
+  Plane Surface(33) = {-33};
+  Line Loop(34) = {85, 108, -86};
+  Plane Surface(34) = {-34};
+  Line Loop(35) = {86, 109, -87};
+  Plane Surface(35) = {-35};
+  Line Loop(36) = {87, 79, -46};
+  Plane Surface(36) = {-36};
+  Line Loop(37) = {46, 80, -88};
+  Plane Surface(37) = {-37};
+  Line Loop(38) = {88, 81, 45};
+  Plane Surface(38) = {-38};
+  Line Loop(39) = {110, 49, -47, -79};
+  Plane Surface(39) = {-39};
+  Line Loop(40) = {50, -89, -80, 47};
+  Plane Surface(40) = {-40};
+  Line Loop(43) = {3, 9, 10, 11};
+  Plane Surface(43) = {-43};
+  Line Loop(44) = {1, -14, -13, -12};
+  Plane Surface(44) = {-44};
+  Line Loop(45) = {113, 114, 115, -3, -118, -117, -116, -1};
+  Line Loop(46) = {106, 107, 108, 109, 110, 111, 112, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105};
+  Plane Surface(45) = {-45,-46};
+  Line Loop(47) = {89, 90, 91, 92};
+  Plane Surface(46) = {-47};
+
+
+  Periodic Line { 14,116,117,118,11 } = {12,113, 114, 115,9 } Translate {a_lat,0,0} ;
+
+  Periodic Line {47}  = {110} Translate {epsc,0,0} ;
+  Periodic Line {89}  = {47}  Translate {epsc,0,0} ;
+  Periodic Line {41}  = {91}  Translate {epsc,0,0} ;
+  Periodic Line {100} = {41}  Translate {epsc,0,0} ;
+
+  Periodic Line {44} = {105} Translate {0,epsc,0} ;
+  Periodic Line {92} = {44}  Translate {0,epsc,0} ;
+  Periodic Line {38} = {90}  Translate {0,epsc,0} ;
+  Periodic Line {95} = {38}  Translate {0,epsc,0} ;
+
+  // Transfinite Surface { expression-list } | "*" < = { expression-list } > < Left | Right | Alternate | AlternateRight | AlternateLeft > ;
+  Transfinite Surface { 9 } Left;
+  Transfinite Surface { 10 } Right;
+
+  Transfinite Surface { 19 } Left;
+  Transfinite Surface { 20 } Left;
+
+  Transfinite Surface { 29 } Left;
+  Transfinite Surface { 30 } Left;
+
+  Transfinite Surface { 39 } Left;
+  Transfinite Surface { 40 } Right;
+
+
+  Physical Surface(100) = {7,8,9,17,18,19,23,24,29,37,38,40,46};   // 1 dom in
+  Physical Surface(101) = {45, 1, 2, 3, 4, 10, 5, 6, 39, 36, 35, 34, 33, 32, 31, 30, 25, 26, 27, 28, 22, 21, 20, 16, 15, 14, 13, 12, 11}; // out
+  Physical Surface(102) = {43};  // PML top
+  Physical Surface(103) = {44};  // PML bot
+  Physical Line(1001)   = {12,113,114,115,9}; // bloch x left
+  Physical Line(1002)   = {14,116,117,118,11}; // bloch x right
+  Physical Line(1003)   = {10}; // top bound
+  Physical Line(1004)   = {13}; // bot bound
+  Physical Line(1005)   = {38,39,40,41,42,43,44,45,46,47,48,35}; // bound for lag
+  Physical Point(10000) = {1};   // Printpoint
+EndIf
+

NonLinearEVP/NonLinearEVP_data.geo

+///////////////////////////////////
+//// Author : Guillaume Demesy ////
+////   _data.geo               ////
+///////////////////////////////////
+
+nm       = 1.;
+epsilon0 = 8.854187817e-3*nm;
+mu0      = 400.*Pi*nm;
+cel      = 1.0/(Sqrt[epsilon0 * mu0]);
+deg2rad  = Pi/180;
+
+pp0        = "1Geometry/0";
+pp1        = "2Polarization-Bloch/0";
+pp3        = "3Eigenvalue problem parameters/0";
+pp2        = "4Drude Model parameters/0";
+pp4        = "5Discretization/0";
+pp5        = "6Formulations/0";
+pp6        = "7Output/0";
+close_menu = 0;
+colorpp0   = "MintCream";
+colorpp1   = "Ivory";
+colorpp2   = "PaleGoldenRod";
+colorpp3   = "Ivory";
+colorpp4   = "LightSalmon";
+colorpp5   = "Ivory";
+
+DefineConstant[ a_lat = {50     , Name StrCat[pp0  , "1grating period d [nm]"]          , Highlight Str[colorpp0]  , Closed close_menu} ];
+
+// normalization factor
+norm = a_lat/(2.*Pi*cel);
+
+DefineConstant[
+  d_sq           = {0.806 , Name StrCat[pp0 , "2sq [d]"] , Highlight Str[colorpp0]  , Closed close_menu} , 
+  space2pml      = {1     , Name StrCat[pp0 , "3PML distance to object [d]"]  , Highlight Str[colorpp0]  , Closed close_menu} , 
+  pmlsize        = {5     , Name StrCat[pp0 , "4PML thickness [d]"] , Highlight Str[colorpp0]  , Closed close_menu} , 
+  
+  flag_Hparallel = {1     , Name StrCat[pp1 , "1polarization case"] , Choices {0="E //",1="H //"} },
+  kx             = {0.75  , Name StrCat[pp1 , "2kx [Pi\a]"] , Highlight Str[colorpp1]  , Closed close_menu} , 
+
+  eps_oo_1       = {1     , Name StrCat[pp2 , "0eps_oo_1 [ - ]"] , Highlight Str[colorpp2]  , Closed close_menu} , 
+  om_d_1         = {1.1   , Name StrCat[pp2 , "1om_d_1 [2Pic\a]"] , Highlight Str[colorpp2]  , Closed close_menu} , 
+  gam_1          = {0.05  , Name StrCat[pp2 , "2gam_1 [2Pic\a]"] , Highlight Str[colorpp2]  , Closed close_menu} , 
+  
+  neig          =  {1     , Name StrCat[pp3 , "0Number of eigenvalues [int]"] , Highlight Str[colorpp3]  , Closed close_menu} ,
+  eig_target_re  = {0.0077, Name StrCat[pp3 , "1EV real part target [2Pic\a]"] , Highlight Str[colorpp3]  , Closed close_menu} , 
+  eig_target_im  = {0.2598, Name StrCat[pp3 , "2EV imag part target [2Pic\a]"] , Highlight Str[colorpp3]  , Closed close_menu} , 
+  eig_min_re     = {0.007 , Name StrCat[pp3 , "3EV real min [2Pic\a]"] , Highlight Str[colorpp3]  , Closed close_menu} , 
+  eig_max_re     = {0.009 , Name StrCat[pp3 , "4EV real max [2Pic\a]"] , Highlight Str[colorpp3]  , Closed close_menu} , 
+  eig_min_im     = {0.25  , Name StrCat[pp3 , "5EV imag min [2Pic\a]"] , Highlight Str[colorpp3]  , Closed close_menu} , 
+  eig_max_im     = {0.27  , Name StrCat[pp3 , "6EV imag max [2Pic\a]"] , Highlight Str[colorpp3]  , Closed close_menu} , 
+
+  paramaille     = {4     , Name StrCat[pp4 , "0number of mesh elements per period []"]  , Highlight Str[colorpp4]  , Closed close_menu} , 
+  flag_Tmesh     = {0     , Name StrCat[pp4 , "2locally structured mesh?"] , Choices {0="unstruct",1="struct"} },
+  flag_o2        = {1     , Name StrCat[pp4 , "3FEM order"] , Choices {0="o1",1="o2"} },
+  
+  flag_res       = {2     , Name StrCat[pp5  , "0resolution type"], 
+                      Choices {0="Aux_E" ,1="PEP_E" ,2="NEP_E" ,3="Lag_E" ,4="PEP_h", 5="all"},ServerAction "ResetDatabase"},
+  flag_outEigvec = {1     , Name StrCat[pp4, "output eigenvector?"], Choices{0,1}}
+];
+
+// Normalization
+d_sq           = d_sq          * a_lat;
+space2pml      = space2pml     * a_lat;
+pmlsize        = pmlsize       * a_lat;
+kx             = kx            * Pi/a_lat;
+eig_target_re  = eig_target_re / norm;
+eig_target_im  = eig_target_im / norm;
+eig_min_re     = eig_min_re    / norm;
+eig_max_re     = eig_max_re    / norm;
+eig_min_im     = eig_min_im    / norm;
+eig_max_im     = eig_max_im    / norm;
+om_d_1         = om_d_1        / norm;
+gam_1          = gam_1         / norm;
+
+eps_oo_2       = 1;
+om_d_2         = 0;
+gam_2          = 0;
+
+slepc_options_rg = StrCat(" -rg_interval_endpoints ",
+                     Sprintf("%.8lf,",eig_min_re),
+                     Sprintf("%.8lf,",eig_max_re),
+                     Sprintf("%.8lf,",eig_min_im),
+                     Sprintf("%.8lf",eig_max_im));
\ No newline at end of file

NonLinearEVP/README.TXT

+A demo of non-linear eigenvalue problems in GetDP.
+
+This model computes one selected eigenfrequency of a grating made of a 
+frequency-dispersive material. Five different formulations are given,
+calling the polynomial and (rational) non-linear solvers of the SLEPc library
+thanks to the unified Eig operator.
+
+Quick start
+-----------
+
+Open `NonLinearEVP.pro' with Gmsh.
+
+Additional info
+---------------
+
+Some documentation is available here: https://arxiv.org/abs/1802.02363
\ No newline at end of file