In [1]:

    

#------Library loading------

# numpy for matrix computations
import numpy as np; import numpy.ma as ma

# system libraries
import sys

# plotting libraries
%matplotlib inline
import matplotlib.pylab as plt
from matplotlib.patches import Circle, Ellipse

# Generalized Multiparticle Mie import
sys.path.append('../')
import py_gmm

# plotly
from plotly import __version__
import plotly.graph_objs as go
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
init_notebook_mode(connected=True)

print(__version__) # requires version >= 1.9.0


    

    












    




2.0.9

In [2]:

    

# building the optical constant database
eps_db_out=py_gmm.mat.generate_eps_db('../epsilon/',ext='*.edb')
eps_files,eps_names,eps_db=eps_db_out['eps_files'],eps_db_out['eps_names'],eps_db_out['eps_db']

# sphere radius (in nm)
v_r = np.array([ 40.,  40.])

# sphere position (in nm)
m_xyz = np.array([[  -42.5,    0. ,   0. ],
                  [ 42.5,   0. ,   0. ]])

# how many spheres in the target? We guess it from the length of the radius vector
ns = len(v_r)

# sphere composition, calling the names contained in "eps_names", just populated above
target_comp= np.array(['eAuJCSTDf','eAuJCSTDf']) # vector containing the optical constants names

# refractive index of the environment
n_matrix = 1.33  # water


    

In [3]:

    

# Euler angles: (alpha,beta,gamma)=(0,0,0) means a z-directed, x-polarized plane wave
alpha = 0.0   # azimuth
beta = 0.0  # polar
gamma = 0.0  # polarization

# Wavelengths for the specrum computation
wl_min = 300
wl_max = 800
n_wl = 250
v_wl = np.linspace(wl_min,wl_max,n_wl)

# Wavelength for the far field computation
v_wl_lf = [530.0,650] # resonance wavelengths


    

In [4]:

    

n_stop=10 # maximum multipolar expansion order
f_int=0.0; # interaction cutoff (normally set to zero to have full interactions)
lf_ratio=300; # plot sphere local field contribution up to a distance equal to d=lf_ratio*r_sphere
qs_flag='no' # retarded simulation


    

In [5]:

    

# target plot
fig = plt.figure(num=1,figsize=(10,10)) # setting the figure size
ax = fig.add_subplot(1, 1, 1, aspect='equal')  # creating the plotting axis

# plot bounds and eliminating x and y ticks
plt.xlim(-1.1*(-m_xyz[0,0]+v_r[0]),1.1*(m_xyz[1,0]+v_r[ns-1]))
plt.ylim(-1.1*(v_r[0]),1.1*(v_r[0]))
plt.xticks([])
plt.yticks([])

# plotting the target
v_color = ['0.6','0.6']
for c,r,col in zip(m_xyz,v_r,v_color):
    c0=c[0];c1=c[1];
    ax.add_patch(Circle((c0,c1),r,color=col))
    ax.eventson


    

In [6]:

    

# computing the expansion coefficients and cross sections with a loop
m_abcd_ext_sca_abs = []  # list to be filled with the output
for wl in v_wl:
    
    # retrieving optical constants at wl from the database
    e_list=py_gmm.mat.db_to_eps(wl,eps_db,target_comp);
    m_eps=np.column_stack((np.real(e_list),np.imag(e_list)));
    
    # solving the gmm problem (calculating the cross sections and the expansion coefficients)
    out=py_gmm.gmm_py.gmm_f2py_module.expansion_coefficients(m_xyz, # target sphere position in nm
                                                      v_r,   # target sphere radii in nm
                                                      m_eps, # e1 and e2 for each sphere
                                                      f_int,  # interaction coefficient
                                                      n_matrix, # environment refractive index
                                                      wl, # computation wavelength
                                                      alpha,beta,gamma, # euler angles for the incident pw
                                                      n_stop, # maximum number for expansion coefficients
                                                      qs_flag) # quasi static approximation
    m_abcd_ext_sca_abs.append(out)
    
# extracting coefficients and cross section
v_coeff=[];v_cext=[];v_csca=[];v_cabs=[];
for out in m_abcd_ext_sca_abs:
    v_coeff.append(out[0]);
    v_cext.append(out[1]);
    v_csca.append(out[2]);
    v_cabs.append(out[3]);

# converting the lists to numpy arrays
v_coeff=np.array(v_coeff)
v_cext=np.array(v_cext)
v_csca=np.array(v_csca)
v_cabs=np.array(v_cabs)


    

In [10]:

    

# local field for the first resonance
v_far = []
for wl_lf in v_wl_lf:

    # optical constants
    e_list=py_gmm.mat.db_to_eps(wl_lf,eps_db,target_comp);
    m_eps=np.column_stack((np.real(e_list),np.imag(e_list)));

    # gmm coefficients computation
    out=py_gmm.gmm_py.gmm_f2py_module.expansion_coefficients(m_xyz,  # target sphere position in nm
                                                             v_r,  # target sphere radii in nm
                                                             m_eps,  # e1 and e2 for each sphere
                                                             f_int,  # interaction coefficient
                                                             n_matrix,  # environment refractive index
                                                             wl_lf, # computation wavelength
                                                             alpha,beta,gamma,  # euler angles for the incident pw
                                                             n_stop,  # maximum number for expansion coefficients
                                                             qs_flag)  # quasi static approximation
    v_amnbmn=out[0][:,0] # getting field expansion coefficients
    v_dmncmn=out[0][:,1]

    # local field
    v_emn=py_gmm.gmm_py.gmm_f2py_module.emn(n_stop)[0] # normalization coeffs

    # retrieving the scattered cloud
    k_far = 2.0*np.pi*n_matrix/wl_lf
    m_sc, scatot, error = py_gmm.gmm_py.gmm_f2py_module.efar(k_far,
                                                             n_stop,
                                                             0.0,3.6e2,360,
                                                             180,
                                                             0.0,
                                                             v_amnbmn,
                                                             m_xyz)
    v_far.append(m_sc)


    

In [11]:

    

# cross section plot
f_size=25;
f_size_ticks=20;
plt.figure(1,figsize=(15,10));
plt.plot(v_wl,np.sum(v_cext,axis=1),'k',linewidth=3.0);

# plt title
plt.title('Au dimer',fontsize=f_size)

# axes labels
plt.xlabel(r'wavelength (nm)', fontsize=f_size)
plt.ylabel(r'Cext', fontsize=f_size)

# ticks
plt.xticks(fontsize=f_size_ticks)
plt.yticks(fontsize=f_size_ticks)

# legend
plt.legend((r'Au Integral C$_{ext}$',''),frameon=False,fontsize=f_size-5)

# layout
plt.tight_layout()


    

In [14]:

    

# Plotly scattering plot

# extracting data from output
m_sc_out = v_far[1]
v_theta = m_sc_out[:,0,0]
v_phi = m_sc_out[0,:,1]

theta = m_sc_out[:,:,1].real
phi = m_sc_out[:,:,0].real
r = m_sc_out[:,:,6].real

x = r * np.sin(phi) * np.cos(theta)
y = r * np.cos(phi)
z = r * np.sin(phi) * np.sin(theta)


surface = go.Surface(x=x, y=y, z=z, colorscale='Viridis',surfacecolor=r)
data = [surface]
layout = go.Layout(
    title='Dimer scattering',
    scene=dict(
        xaxis=dict(
            gridcolor='rgb(255, 255, 255)',
            zerolinecolor='rgb(255, 255, 255)',
            showbackground=True,
            backgroundcolor='rgb(230, 230,230)'
        ),
        yaxis=dict(
            gridcolor='rgb(255, 255, 255)',
            zerolinecolor='rgb(255, 255, 255)',
            showbackground=True,
            backgroundcolor='rgb(230, 230,230)'
        ),
        zaxis=dict(
            gridcolor='rgb(255, 255, 255)',
            zerolinecolor='rgb(255, 255, 255)',
            showbackground=True,
            backgroundcolor='rgb(230, 230,230)'
        )
    )
)
fig = go.Figure(data=data, layout=layout)
iplot(fig, filename='parametric-plot-viridis')


    

In [ ]: