| 12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697 |
- #!/usr/bin/env python3
- # -*- coding: UTF-8 -*-
- import cmath
- import numpy as np
- import scipy.io
- from matplotlib import pyplot as plt
- from optical_constants import read_refractive_index_from_yaml as read_nk
- from scattnlay import mesomie, mie
- from evalMie import *
- from_disk = np.loadtxt('silver-d_perp_interpolated.txt')
- omega_star_ratio = from_disk[0, :]
- d_perp = from_disk[1, :] + 1j*from_disk[2, :]
- from_disk = np.loadtxt('silver-d_parl_interpolated.txt')
- d_parl = from_disk[1, :] + 1j*from_disk[2, :]
- c = 299792458 # m/s
- h_reduced = 6.5821e-16 # eV s
- omega_p_star = 3.81 # eV
- # omega_p = 9.02 # eV
- # gamma = 0.022 # eV
- # eps_d = 1
- R = 2.5
- y_min = 0
- y_max = 2
- min_lim_omega_star_ratio = 0.87
- max_lim_omega_star_ratio = 0.99
- # min_lim_omega_ratio = min_lim_omega_star_ratio * omega_p_star/omega_p
- # max_lim_omega_ratio = max_lim_omega_star_ratio * omega_p_star/omega_p
- # 2 pi / lambda = (omega/c) /h_reduced
- WL = 2*np.pi/((omega_star_ratio * omega_p_star/c)/h_reduced)*1e6 # mkm
- min_WL_available = 0.1879
- max_WL_available = 1.9370
- WL[WL < min_WL_available] = min_WL_available
- WL[WL > max_WL_available] = max_WL_available
- index_Ag = read_nk('Ag-Johnson-1972.yml', WL, kind=1)
- # eps_Ag = index_Ag**2
- # print(index_Ag)
- # def eps_m(omega):
- # return 1 - omega_p * omega_p / (omega*omega + 1j*omega*gamma)
- Qext = []
- Qext_mie = []
- om_rat_plot = []
- x_in = []
- m_in = []
- d_parl_in, d_perp_in = [],[]
- for i in range(len(omega_star_ratio)):
- om_star_rat = omega_star_ratio[i]
- if (om_star_rat < min_lim_omega_star_ratio
- or om_star_rat > max_lim_omega_star_ratio):
- continue
- omega = om_star_rat*omega_p_star
- # WL_mkm = 2*np.pi/((omega/c)/h_reduced)*1e6
- # if WL_mkm < min_WL_available or WL_mkm > max_WL_available:
- # continue
- x = (omega/c) * R * 1e-9/h_reduced
- m = index_Ag[i, 1]
- x_in.append(x)
- m_in.append(m)
- d_parl_in.append(d_parl[i])
- d_perp_in.append(d_perp[i])
- n_m = 1
- om_rat_plot.append(om_star_rat)
- # print('--', om_rat_plot[0], om_rat_plot[-1], len(om_rat_plot))
- # print(x_in[0], x_in[-1], len(x_in))
- # print(m_in[0], m_in[-1], len(m_in))
- # print(R)
- # print(d_perp_in[0],d_perp_in[-1], len(d_perp_in))
- # print(d_parl_in[0],d_parl_in[-1], len(d_parl_in))
- Qext_mie, _ = eval_mie_spectrum(x_in, m_in)
- Qext, _ = eval_mesomie_spectrum(x_in, m_in, R, d_perp_in, d_parl_in)
- plt.plot(om_rat_plot, Qext_mie,
- label='classic', color='gray', lw=4)
- plt.plot(om_rat_plot, Qext,
- label='non-classic', color='red', lw=4)
- plt.legend()
- # plt.yscale('log')
- plt.xlim((min_lim_omega_star_ratio, max_lim_omega_star_ratio))
- plt.ylim((y_min, y_max))
- plt.title("R="+str(R))
- plt.show()
|
