Adding ElectromagneticScattering model

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 general 3D electromagnetic scattering problems:
+* 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)
+
+## Running the model
+
+Open `scattering.pro` with Gmsh.
+
+## Authors
+
+Guillaume Demésy
\ 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