#!/usr/bin/env python3 # -*- coding: UTF-8 -*- # # Copyright (C) 2021 Konstantin Ladutenko # # This file is part of python-scattnlay # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # The only additional remark is that we expect that all publications # describing work using this software, or all commercial products # using it, cite the following reference: # [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by # a multilayered sphere," Computer Physics Communications, # vol. 180, Nov. 2009, pp. 2348-2354. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . from scattnlay import mie import matplotlib.pyplot as plt import numpy as np from optical_constants import read_refractive_index_from_yaml as get_index def gauss(x, mu, sigma): return 1/(sigma * np.sqrt(2 * np.pi)) * np.exp( - (x - mu)**2 / (2 * sigma**2)) # fill_factor = 0.8 from_WL = 300 to_WL = 1100 WL_points= 100 WLs = np.linspace(from_WL, to_WL, WL_points) r_points = 200 from_r = 40/2. to_r =80/2. all_r = np.linspace(from_r, to_r, r_points) r_mean = 58.3/2. r_std = 6.3/2. # from_r = 1. # to_r =15. # all_r = np.linspace(from_r, to_r, r_points) # r_mean = 6 # r_std = 2 core_r = 30 inner_shell_h = 8 outer_shell_h = 10 host_media = 1.33 r_weights = gauss(all_r, r_mean,r_std)/len(all_r) plt.plot(all_r, r_weights ) plt.xlabel("R, nm") plt.ylabel("amount") index_SiO2 = get_index("refractiveindex_info/SiO2-Gao.yml", WLs, units='nm') # index_Au = get_index("refractiveindex_info/Au-McPeak.yml", WLs, units='nm') index_Au = get_index("refractiveindex_info/Au-Johnson.yml", WLs, units='nm') index_TiO2 = get_index("refractiveindex_info/TiO2-Sarkar.yml", WLs, units='nm') # index_TiO2[:,1] += 0.5j # index_Au[:,1] = index_Au[:,1]* fill_factor + (1-fill_factor) # print("Au index before correction, max = ", np.max(np.imag(index_Au[:,1]))) # # Taking into account gold free electrons damping # # contributed by surface scattering and bulk dumping # # See eq 1 in [1] -> doi: https://doi.org/10.1186/s11671-018-2670-7 # eps_exp = index_Au[:,1]**2 # c=299792458 # speed of light # h= 4.135667516e-15 # eV*s, Planck constant # w = h*c/(WLs*1e-9) # eV, frequency # w_p =8.55 # eV, gold plasmon frequency # g_b = 18.4e-3 # eV, bulk dumping # v_f = 1.4e6 # m/s, Fermi velocity of electrons in gold # A = 1.33 # fit parameter for 16 nm gold shell, see Table 2 in [1] # L_b = (4.*((core_r+inner_shell_h)**3 - core_r**3)/ # (3.*((core_r+inner_shell_h)**2 + core_r**2))) # # g_s = v_f/L_b # eq 2 in [1] # g_s = h*A*v_f/(inner_shell_h*1e-9) # eq 4 in [1] # # g_b *= 0.2 # # g_s *= 0.5 # eps_Au = (eps_exp # + # w_p**2 / ( w * (w + 1j*g_b) ) # - # w_p**2 / ( w * (w + 1j*(g_b+g_s)) ) # ) # index_Au[:,1] = np.sqrt(eps_Au) index_Au[:,1] += 1.6j # print(f"L_b={L_b}") # # print(w) # print(f"g_s={g_s}, g_b={g_b} ") # print("Au index after, max = ", np.max(np.imag(index_Au[:,1]))) x = np.ones((3), dtype = np.float64) m = np.ones((3), dtype = np.complex128) Qext_core_shell = np.zeros(len(WLs)) Qext_3l = np.zeros(len(WLs)) for i in range(len(WLs)): WL = WLs[i] for j in range(len(all_r)): core_r = all_r[j] # inner_shell_h = all_r[j] weight = r_weights[j] # print(core_r) x = host_media*2.0*np.pi/WL*np.array([core_r, core_r+inner_shell_h, core_r+inner_shell_h+outer_shell_h]) m = np.array([index_SiO2[i][1], index_Au[i][1], index_TiO2[i][1]] )/host_media # print(x, m) mie.SetLayersSize(x) mie.SetLayersIndex(m) mie.RunMieCalculation() Qext_3l[i] += mie.GetQext()*weight x = host_media*2.0*np.pi/WL*np.array([core_r, core_r+inner_shell_h]) m = np.array([index_SiO2[i][1], index_Au[i][1]])/host_media mie.SetLayersSize(x) mie.SetLayersIndex(m) mie.RunMieCalculation() Qext_core_shell[i] += mie.GetQext()*weight fig, axs2 = plt.subplots(1,1)#, sharey=True, sharex=True) axs2.plot(WLs, Qext_3l, color="purple") axs2.plot(WLs, Qext_core_shell, color="lime") axs2.set_xlabel("WL, nm") axs2.set_ylabel("Extinction, a.u.") # axs2 = axs.twinx() # axs2.plot(np.array(core_r_vec)*2,an_vec[:,0],"b.",lw=0.8, markersize=1.9,label="$a_0$") # axs2.plot(np.array(core_r_vec)*2,bn_vec[:,0],"b-", markersize=1.9,label="$b_0$") # axs2.plot(np.array(core_r_vec)*2,an_vec[:,1],"g.",lw=0.8, markersize=1.9,label="$a_1$") # axs2.plot(np.array(core_r_vec)*2,bn_vec[:,1],"g-", markersize=1.9,label="$b_1$") # axs2.legend(loc="upper right") # axs2.tick_params('y', colors='black') # axs2.set_ylim(0,1) # axs2.set_ylabel("Mie",color="black") plt.savefig("spectra.pdf",pad_inches=0.02, bbox_inches='tight') plt.show() plt.clf() plt.close()