fig6

examples/NatureComm-2020-Plasmon–emitter_interactions_at_the_nanoscale/Ag-Johnson-1972.yml

@@ -0,0 +1,58 @@
+# this file is part of refractiveindex.info database
+# refractiveindex.info database is in the public domain
+# copyright and related rights waived via CC0 1.0
+
+REFERENCES: "P. B. Johnson and R. W. Christy. Optical constants of the noble metals, <a href=\"https://doi.org/10.1103/PhysRevB.6.4370\"><i>Phys. Rev. B</i> <b>6</b>, 4370-4379 (1972)</a>"
+COMMENTS: "Room temperature"
+DATA:
+  - type: tabulated nk
+    data: |
+        0.1879 1.07 1.212
+        0.1916 1.10 1.232
+        0.1953 1.12 1.255
+        0.1993 1.14 1.277
+        0.2033 1.15 1.296
+        0.2073 1.18 1.312
+        0.2119 1.20 1.325
+        0.2164 1.22 1.336
+        0.2214 1.25 1.342
+        0.2262 1.26 1.344
+        0.2313 1.28 1.357
+        0.2371 1.28 1.367
+        0.2426 1.30 1.378
+        0.2490 1.31 1.389
+        0.2551 1.33 1.393
+        0.2616 1.35 1.387
+        0.2689 1.38 1.372
+        0.2761 1.41 1.331
+        0.2844 1.41 1.264
+        0.2924 1.39 1.161
+        0.3009 1.34 0.964
+        0.3107 1.13 0.616
+        0.3204 0.81 0.392
+        0.3315 0.17 0.829
+        0.3425 0.14 1.142
+        0.3542 0.10 1.419
+        0.3679 0.07 1.657
+        0.3815 0.05 1.864
+        0.3974 0.05 2.070
+        0.4133 0.05 2.275
+        0.4305 0.04 2.462
+        0.4509 0.04 2.657
+        0.4714 0.05 2.869
+        0.4959 0.05 3.093
+        0.5209 0.05 3.324
+        0.5486 0.06 3.586
+        0.5821 0.05 3.858
+        0.6168 0.06 4.152
+        0.6595 0.05 4.483
+        0.7045 0.04 4.838
+        0.7560 0.03 5.242
+        0.8211 0.04 5.727
+        0.8920 0.04 6.312
+        0.9840 0.04 6.992
+        1.0880 0.04 7.795
+        1.2160 0.09 8.828
+        1.3930 0.13 10.10
+        1.6100 0.15 11.85
+        1.9370 0.24 14.08

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

@@ -0,0 +1,62 @@
+# -*- coding: utf-8 -*-
+"""Provide functionality to read optical constants in format provided by `refractiveindex.info <https://refractiveindex.info/>`_ website"""
+# GPL license from Smuthi project
+
+from scipy.interpolate import interp1d
+import numpy as np
+import yaml
+
+
+def read_refractive_index_from_yaml(filename, vacuum_wavelength, units="mkm", kind=1):
+    """Read optical constants in format provided by refractiveindex.info website.
+
+    Args:
+            filename (str): path and file name for yaml data
+                            downloaded from refractiveindex.info
+            vacuum_wavelength (float or np.array): wavelengths where refractive
+                                                   index data is needed
+            units (str): units for wavelength. currently, microns ('mkm' or 'um')
+                         and nanometers ('nm') can be selected
+            kind (int): order of interpolation
+
+    Returns:
+            A pair (or np.array of pairs) of wavelength and
+            corresponding refractive index (complex)
+    """
+    if units == "nm":
+        factor = 1000
+    elif units in ("mkm", "um"):
+        factor = 1
+    else:
+        raise NotImplementedError("Converting wavelength into '"+units
+                                  + "' units for refractive index data"
+                                  + " was not implemented.")
+
+    with open(filename) as f:
+        the_file = yaml.load(f, yaml.SafeLoader)['DATA'][0]
+    data_type = the_file['type']
+    if data_type != 'tabulated nk':
+        raise NotImplementedError("Input data type '"+data_type
+                                  + "' available in file "+filename
+                                  + " was not implemented.")
+    data = the_file['data'].splitlines()
+    data_split = []
+    for wl in data:
+        data_split.append(wl.split())
+    data_num = []
+    for wl in data_split:
+        record = []
+        for val in wl:
+            record.append(float(val))
+        data_num.append(record)
+    data_np = np.array(data_num)
+    data_wl = data_np[:, 0]*factor
+    index_re = data_np[:, 1]
+    index_im = data_np[:, 2]
+    f_re = interp1d(data_wl, index_re, kind=kind)
+    f_im = interp1d(data_wl, index_im, kind=kind)
+    data_out = np.transpose(np.vstack(
+        (vacuum_wavelength, f_re(vacuum_wavelength)+f_im(vacuum_wavelength)*1j)))
+    if len(data_out) == 1:
+        return data_out[0]
+    return data_out

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

@@ -9,24 +9,45 @@ import scipy.io
 mat = scipy.io.loadmat('d-parameters/rs=4.mat')
 x_mat = mat['omegav'][0]
 d_perp_mat = mat['dperp'][0]*10
-
-# arr2D = np.loadtxt('rs4-im-d_perp.csv', delimiter=',')
-# im_d = arr2D[arr2D[:, 0].argsort()]
-# arr2D = np.loadtxt('rs4-re-d_perp.csv', delimiter=',')
-# re_d = arr2D[arr2D[:, 0].argsort()]
-
 x = np.linspace(x_mat[0], x_mat[-1], 1001)
 im_d_y = interpolate.interp1d(x_mat,  np.imag(d_perp_mat), kind='cubic')
 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])
 np.savetxt('rs4-d_perp_interpolated.txt', data)
 
+mat = scipy.io.loadmat('d-parameters/silver_xd_1.18A.mat')
+x_mat = mat['omegav'][0]
+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'
+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])
+np.savetxt('silver-d_perp_interpolated.txt', data)
+im_d_y = interpolate.interp1d(x_mat,  np.imag(d_parl_mat), kind=kind)
+re_d_y = interpolate.interp1d(x_mat,  np.real(d_parl_mat))
+data = np.array([x.T, re_d_y(x).T, im_d_y(x).T])
+np.savetxt('silver-d_parl_interpolated.txt', data)
+
 
 from_disk = np.loadtxt('rs4-d_perp_interpolated.txt')
+plt.figure('rs4')
 plt.plot(from_disk[0, :], from_disk[1, :], label='re d')
 plt.plot(from_disk[0, :], from_disk[2, :], label='im d')
 plt.plot(x_mat, np.real(d_perp_mat), label='re d mat')
 plt.plot(x_mat, np.imag(d_perp_mat), label='re d mat')
 plt.legend()
 plt.show()
+
+
+from_disk = np.loadtxt('silver-d_perp_interpolated.txt')
+plt.figure('silver')
+plt.plot(from_disk[0, :], from_disk[1, :], label='re d perp')
+plt.plot(from_disk[0, :], from_disk[2, :], label='im d perp')
+from_disk = np.loadtxt('silver-d_parl_interpolated.txt')
+plt.plot(from_disk[0, :], from_disk[1, :], label='re d parl')
+plt.plot(from_disk[0, :], from_disk[2, :], label='im d parl')
+
+plt.legend()
+plt.show()

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

@@ -0,0 +1,93 @@
+#!/usr/bin/env python3
+# -*- coding: UTF-8 -*-
+from matplotlib import pyplot as plt
+import cmath
+from scattnlay import mesomie, mie
+from optical_constants import read_refractive_index_from_yaml as read_nk
+import numpy as np
+import scipy.io
+
+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 = 9.02  # eV
+omega_p_star = 3.81  # 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 = []
+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 = 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]/2,      # 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())
+
+    om_rat_plot.append(om_star_rat)
+
+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()