# -*- 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')