karpov cs works (#68)

examples/NatureComm-2020-Plasmon–emitter_interactions_at_the_nanoscale/evalMie.py

@@ -0,0 +1,49 @@
+from scattnlay import mesomie, mie
+import numpy as np
+
+
+
+def eval_mie_single_spectrum_point(x,m,
+                                   R=None, d_perp=None, d_parl=None):
+    mie.SetLayersSize(x)
+    mie.SetLayersIndex(m)
+    mie.RunMieCalculation()
+    return mie.GetQext(), mie.GetQsca()
+
+
+def eval_mesomie_single_spectrum_point(x, m, n_m, R, d_perp, d_parl):
+    mesomie.calc_ab(R*10,      # calc_ab needs R in angstrem
+                x,      # xd
+                x * m,  # xm
+                n_m * n_m,      # eps_d 
+                m * m,  # eps_m
+                d_parl,      # d_parallel
+                d_perp)      # d_perp
+    mesomie.calc_Q()
+    return mesomie.GetQext(), mesomie.GetQsca()
+
+
+def eval_mie_spectrum(x, m):
+    Qext = []
+    Qsca = []
+    for i in range(len(x)):
+        ext, sca =  eval_mie_single_spectrum_point(x[i], m[i])
+        Qext.append(ext)
+        Qsca.append(sca) 
+    return Qext, Qsca
+
+def eval_mesomie_spectrum(x, m, R, d_perp, d_parl=None, n_m=1):
+    if np.array(d_perp == None).any():
+        d_perp = np.zeros(len(x))
+    if np.array(d_parl == None).any():
+        d_parl = np.zeros(len(x))
+    Qext = []
+    Qsca = []
+    for i in range(len(x)):
+        ext, sca =  eval_mesomie_single_spectrum_point(x[i], m[i], n_m, R,
+                                                d_perp[i], d_parl[i])
+        Qext.append(ext)
+        Qsca.append(sca) 
+    return Qext, Qsca
+
+

examples/NatureComm-2020-Plasmon–emitter_interactions_at_the_nanoscale/fit_d_param.py

@@ -0,0 +1,96 @@
+#!/usr/bin/env python3
+import cmath
+import numpy as np
+from matplotlib import pyplot as plt
+
+from optical_constants import read_refractive_index_from_yaml as read_nk
+
+from evalMie import *
+from karpov_materials import *
+
+WL_min = 325 #nm
+WL_max = 390 #nm
+
+# WL_min = 335 #nm
+# WL_max = 500 #nm
+shell_h = 0.5
+n_m = 1
+
+d_param_filename = ['silver-d_perp_interpolated.txt',
+                    'silver-d_parl_interpolated.txt']
+
+
+def eval_xm(WLs,R, shell_h, n_m):
+    x_cs_in, m_cs_in = [], []
+    x_bulk_in, m_bulk_in = [], []
+
+    w,wp, Gamma = evalDrudeParams(R, WLs)
+    
+    for i in range(len(WLs)):
+        x_const = 2 * n_m *np.pi/ WLs[i]
+        x_cs = [(R-shell_h) * x_const, R * x_const],
+        x_meso = R*x_const 
+
+        n_drude = cmath.sqrt(
+            eps_drude(w[i], wp, Gamma))
+        n_shell = cmath.sqrt(
+            eps_drude(w[i], wp * correction_wp[j], Gamma))
+        m_cs = [n_drude / n_m, n_shell / n_m]
+        m_meso = n_drude / n_m
+
+        x_cs_in.append(x_cs)
+        m_cs_in.append(m_cs)
+        x_bulk_in.append(x_meso)
+        m_bulk_in.append(m_meso)
+    return x_cs_in, m_cs_in, x_bulk_in, m_bulk_in
+
+def plot_Q(WLs, Qext_drude, legend, ls='-'):
+    plt.plot(WLs, Qext_drude, ls=ls)
+    plt.legend(legend, fontsize=14)
+    plt.xlabel("\u03BB, nm", fontsize=14)
+    plt.ylabel("Normalized extintion cross-section, $nm^{2}$", fontsize=14)
+
+
+# correction_wp = [0.97, 0.95, 0.92, 0.89, 0.87, 0.855, 0.84, 0.83]
+
+correction_wp = [0.97, 0.96, 0.95, 0.94, 0.93, 0.92, 0.91, 0.9]
+Rs = [1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5]
+
+correction_wp = correction_wp[:3]
+Rs = Rs[:3]
+
+
+legend_cs = [str(r)+'-'+str(w)+' cs' for r,w in zip(Rs, correction_wp)]
+legend_mie = [str(r)+'-'+str(w) for r,w in zip(Rs, correction_wp)]
+legend_meso = [str(r)+'-'+str(w)+' d' for r,w in zip(Rs, correction_wp)]
+
+WLs_nm, d_perp_in, d_parl_in = get_d_params(WL_min, WL_max, d_param_filename)
+
+index_Ag = read_nk('Ag-Johnson-1972.yml', WLs_nm, units='nm', kind=1)
+m_Johnson = index_Ag[:, 1]
+
+# WLs_nm = np.linspace(WL_min, WL_max, 401)
+Qext = []
+for j in range(len(Rs)):
+
+    x_cs_in, m_cs_in, x_bulk_in, m_bulk_in = eval_xm(WLs_nm,Rs[j], shell_h, n_m)
+
+    Qext_cs_drude, _ = eval_mie_spectrum(x_cs_in, m_cs_in)
+    plt.figure('mie')
+    plot_Q(WLs_nm, Qext_cs_drude, legend_cs, ls='dotted')
+
+    # m_bulk_in = m_Johnson
+
+    Qext_bulk_drude, _ = eval_mie_spectrum(x_bulk_in, m_bulk_in)
+    plt.figure('mie')
+    plot_Q(WLs_nm, Qext_bulk_drude, legend_mie)
+
+    Qext_mesomie, _ = eval_mesomie_spectrum(x_bulk_in, m_bulk_in, Rs[j], d_perp_in, d_parl_in)
+    plt.figure('mie')
+    plot_Q(WLs_nm, Qext_mesomie, legend_meso, ls='--')
+
+    # print(WLs_nm, Qext_drude)
+    # print(max(Qext_drude))
+
+plt.show()
+

examples/NatureComm-2020-Plasmon–emitter_interactions_at_the_nanoscale/karpov_materials.py

@@ -0,0 +1,59 @@
+import numpy as np
+
+
+def evalDrudeParams(R, WLs):
+    GammaB = 3.5 * 10**13
+    Gamma = GammaB + ((1.4 * 10**6) / (R * 10**-9))
+    w = (2 * np.pi * 3 * 1E8) / (WLs * 1E-9)
+    wp = 1.4
+    return w, wp, Gamma
+
+
+def eps_r(w, wp, Gamma):
+    return 5 - (w**2 * (wp * 1E16)**2) / (w**4 - w**2 * Gamma**2)
+
+
+def eps_i(w, wp, Gamma):
+    return (w * (wp * 1E16)**2 * Gamma) / (w**4 - w**2 * Gamma**2)
+
+
+def eps_drude(w, wp, Gamma):
+    return eps_r(w, wp, Gamma) + 1j * eps_i(w, wp, Gamma)
+
+
+def get_d_params(WL_min, WL_max, d_param_filename):
+    from_disk = np.loadtxt(d_param_filename[0])
+    omega_star_ratio = from_disk[0, :]
+    d_perp = from_disk[1, :] + 1j*from_disk[2, :]
+    from_disk = np.loadtxt(d_param_filename[1])
+    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 
+
+    # min_lim_omega_star_ratio = 0.87
+    # max_lim_omega_star_ratio = 0.99
+
+    d_perp_in, d_parl_in = [],[]
+    WLs_nm = []
+    om_rat_plot = []
+    WL_d = 2*np.pi/((omega_star_ratio * omega_p_star/c)/h_reduced)*1e9  # nm
+    for i in range(len(omega_star_ratio)):
+        if WL_d[i] > WL_max: 
+            continue
+        if WL_d[i] < WL_min:
+            continue
+        # 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
+        # om_rat_plot.append(om_star_rat)
+        WLs_nm.append(WL_d[i])
+        d_perp_in.append(d_perp[i])
+        d_parl_in.append(d_parl[i])
+
+    WLs_nm = np.array(WLs_nm)
+    d_perp_in = np.array(d_perp_in)
+    d_parl_in = np.array(d_parl_in)
+    return WLs_nm, d_perp_in, d_parl_in

examples/NatureComm-2020-Plasmon–emitter_interactions_at_the_nanoscale/plot_d_param.py

@@ -21,6 +21,7 @@ d_perp_mat = mat['dperp'][0]*10
 d_parl_mat = mat['dparl'][0]*10
 x = np.linspace(x_mat[0], x_mat[-1], 1001)
 kind = 'linear'
+# kind = 'cubic'
 im_d_y = interpolate.interp1d(x_mat,  np.imag(d_perp_mat), kind=kind)
 re_d_y = interpolate.interp1d(x_mat,  np.real(d_perp_mat))
 data = np.array([x.T, re_d_y(x).T, im_d_y(x).T])

examples/NatureComm-2020-Plasmon–emitter_interactions_at_the_nanoscale/suppFig5.py

@@ -2,9 +2,8 @@
 # -*- coding: UTF-8 -*-
 from matplotlib import pyplot as plt
 import cmath
-from scattnlay import mesomie, mie
+from evalMie import *
 import numpy as np
-import scipy.io
 min_lim = 0.4
 max_lim = 0.75
 
@@ -27,9 +26,11 @@ def eps_m(omega):
     return 1 - omega_p * omega_p / (omega*omega + 1j*omega*gamma)
 
 
+
+
 Rs = [2.5, 5, 10, 25]  # nm
 y_min = [1e-2, 1e-2, 1e-1, 1e-1]
-y_max = [1e1, 1e1, 5e1, 5e1]
+y_max = [1e1, 5e1, 5e1, 5e1]
 
 # for om_rat in omega_ratio:
 for fig in range(len(Rs)):
@@ -37,6 +38,9 @@ for fig in range(len(Rs)):
     Qext = []
     Qext_mie = []
     om_rat_plot = []
+    x_in = []
+    m_in = []
+    d_perp_in = []
     for i in range(len(omega_ratio)):
         om_rat = omega_ratio[i]
         if om_rat < min_lim or om_rat > max_lim:
@@ -44,26 +48,19 @@ for fig in range(len(Rs)):
         omega = om_rat*omega_p
         m = cmath.sqrt(eps_m(omega))
         x = (omega/c) * R * 1e-9/h_reduced
-        mesomie.calc_ab(R*10,      # calc_ab needs R in angstrem
-                        x,      # xd
-                        x * m,  # xm
-                        1,      # eps_d
-                        m * m,  # eps_m
-                        0,      # d_parallel
-                        d_perp[i])      # d_perp
-        mesomie.calc_Q()
-        mie.SetLayersSize(x)
-        mie.SetLayersIndex(m)
-        mie.RunMieCalculation()
-        Qext.append(mesomie.GetQext())
-        Qext_mie.append(mie.GetQext())
-        # print(x, m, Qext[-1] - mie.GetQext())
-
+        x_in.append(x)
+        m_in.append(m)
         om_rat_plot.append(om_rat)
-    # print(Qext)
+        d_perp_in.append(d_perp[i])
+
+    Qext_mie, _ = eval_mie_spectrum(x_in, m_in)
+    Qext, _ = eval_mesomie_spectrum(x_in, m_in, R, d_perp_in)
+
+    print(Qext)
     plt.figure(fig)
     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.plot(om_rat_plot, Qext, label="R="+str(R), lw=1)
     plt.legend()
     plt.yscale('log')
     plt.xlim((0.4, 0.75))

examples/NatureComm-2020-Plasmon–emitter_interactions_at_the_nanoscale/suppFig6.py

@@ -8,6 +8,7 @@ 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, :]
@@ -17,10 +18,10 @@ d_parl = from_disk[1, :] + 1j*from_disk[2, :]
 
 c = 299792458  # m/s
 h_reduced = 6.5821e-16  # eV s
-omega_p = 9.02  # eV
-omega_p_star = 3.81  # eV
-gamma = 0.022  # eV
-eps_d = 1
+omega_p_star = 3.81  # eV 
+# omega_p = 9.02  # eV
+# gamma = 0.022  # eV
+# eps_d = 1
 
 R = 2.5
 y_min = 0
@@ -41,48 +42,49 @@ 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)
+# eps_Ag = index_Ag**2
+# print(index_Ag)
 
 
-def eps_m(omega):
-    return 1 - omega_p * omega_p / (omega*omega + 1j*omega*gamma)
+# 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
+    # 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 = cmath.sqrt(eps_m(omega))
     m = index_Ag[i, 1]
-    print(x, m)
-    mesomie.calc_ab(R*10,      # R in angstrem
-                    x,      # xd
-                    x * m,  # xm
-                    1,      # eps_d
-                    m * m,  # eps_m
-                    d_parl[i],      # d_parallel
-                    d_perp[i])      # d_perp
-    mesomie.calc_Q()
-    mie.SetLayersSize(x)
-    mie.SetLayersIndex(m)
-    mie.RunMieCalculation()
-    Qext.append(mesomie.GetQext())
-    Qext_mie.append(mie.GetQext())
-    # print(x, m, Qext[-1] - mie.GetQext())
+    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,