initial commit

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 problems on an isolated arbitrary object:
+* T-matrix computation [1]
+* Quasi-normal modes
+* Plane wave response
+* Green's function (e.g. to compute the Local Density Of States)
+
+It contains various usefull features in electromag:
+* Vector partial waves (AKA vector spherical harmonics)
+* Total field formuation with a point source illumination ("an oriented delta")
+* Linear eigenvalue problems
+* Scattered field formulation
+* Spherical PMLs
+
+
+## Installation
+
+This model requires the following programs:
+* [gmsh](http://www.gmsh.info/)
+* [getdp](http://www.getdp.info/) compiled with python support (see below)
+* python (>3.5.x) with numpy, scipy and matplotlib
+
+## Running the model
+
+Open `scattering.pro` with Gmsh.
+The default parameters are set to compute the T-matrix of a sphere. It retrieves the results from [1].
+
+## Authors
+
+Guillaume Demésy and Brian Stout
+
+## References
+[1] See "[Scattering matrix of arbitrarily shaped objects: Combining Finite Elements and Vector Partial Waves](https://arxiv.org/abs/1802.00596)" for details about T-matrices
+
+## Installation notes
+To enable python support (Python[] function) in GetDP, follow [these instructions (with complex arithmetic)](https://gitlab.onelab.info/getdp/getdp/wikis/GetDP-compilation) and add to the final cmake line:
+
+`-DENABLE_PYTHON=ON -DPYTHON_LIBRARY=/path/to/pythonlib -DPYTHON_INCLUDE_DIR=/path/to/pythoninclude`
+*  On Debian/Ubuntu systems, for python3.6 installed with apt-get,
+    * `/path/to/pythonlib` is `/usr/lib/x86_64-linux-gnu/libpython3.6m.so`
+    * `/path/to/python/include` is `/usr/include/python3.6m`
+* For python versions installed through anaconda in some environment (e.g. env py36 below), a common location is:
+    * `/somepath/anaconda3/envs/py36/lib/libpython3.6m.so`
+    * `/somepath/anaconda3/envs/py36/include/python3.6m`
+

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;

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 Stout's recurrence relations (B_ arrays)
+    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 Stout's recurrence relations (B_ arrays)
+    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
+        print('anm not normalized' , np.trapz(np.trapz((np.sin(theta_sph)*EdotconjMnm).transpose(),theta_sph[:,0]),phi_sph[0,:]))
+        print('normalization     ' , normalize_aM_nm2)
+
+
+
+        # print(np.loadtxt('run_results/a_pe%gpo%g.dat'%(1,int(po))))
+    # print(aM_nm)
+    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')

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)