use data from paper

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

@@ -5,33 +5,28 @@ from matplotlib import markers, pyplot as plt
 import numpy as np
 from scipy import interpolate
 
-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()]
+import scipy.io
+mat = scipy.io.loadmat('d-parameters/rs=4.mat')
+x_mat = mat['omegav'][0]
+d_perp_mat = mat['dperp'][0]*10
 
-xmin_im = np.min(im_d[:, 0])
-xmin_re = np.min(re_d[:, 0])
-xmax_im = np.max(im_d[:, 0])
-xmax_re = np.max(re_d[:, 0])
-x = np.linspace(np.max([xmin_im, xmin_re]), np.min([xmax_im, xmax_re]), 1000)
-im_d_y = interpolate.interp1d(im_d[:, 0],  im_d[:, 1])
-re_d_y = interpolate.interp1d(re_d[:, 0],  re_d[:, 1])
+# 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.txt', data)
+np.savetxt('rs4-d_perp_interpolated.txt', data)
 
-# print(data)
-# plt.plot(im_d[:, 0], im_d[:, 1], marker='o', ls='')
-# plt.plot(x, im_d_y(x))
-# plt.plot(re_d[:, 0], re_d[:, 1], marker='o', ls='')
-# plt.plot(x, re_d_y(x))
-# plt.xlim((0, 1))
-# plt.ylim((-4, 5))
-# plt.show()
 
-from_disk = np.loadtxt('rs4-d_perp.txt')
+from_disk = np.loadtxt('rs4-d_perp_interpolated.txt')
 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()

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

@@ -4,10 +4,15 @@ from matplotlib import pyplot as plt
 import cmath
 from scattnlay import mesomie, mie
 import numpy as np
-
-from_disk = np.loadtxt('rs4-d_perp.txt')
+import scipy.io
 min_lim = 0.4
-max_lim = 0.8
+max_lim = 0.75
+
+# mat = scipy.io.loadmat('d-parameters/rs=4.mat')
+# omega_ratio = mat['omegav'][0]
+# d_perp = mat['dperp'][0]*10
+
+from_disk = np.loadtxt('rs4-d_perp_interpolated.txt')
 omega_ratio = from_disk[0, :]
 d_perp = from_disk[1, :] + 1j*from_disk[2, :]
 
@@ -23,39 +28,45 @@ def eps_m(omega):
 
 
 Rs = [2.5, 5, 10, 25]
-R = 5
+y_min = [1e-2, 1e-2, 1e-1, 1e-1]
+y_max = [1e1, 1e1, 5e1, 5e1]
 
-Qext = []
-Qext_mie = []
-om_rat_plot = []
 # for om_rat in omega_ratio:
-for i in range(len(omega_ratio)):
-    om_rat = omega_ratio[i]
-    if om_rat < min_lim or om_rat > max_lim:
-        continue
-    omega = om_rat*omega_p
-    m = cmath.sqrt(eps_m(omega))
-    x = (omega/c) * R * 1e-9/h_reduced
-    mesomie.calc_ab(R*10,      # 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())
-
-    om_rat_plot.append(om_rat)
-# print(Qext)
-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((0.4, 0.8))
+for fig in range(len(Rs)):
+    R = Rs[fig]
+    Qext = []
+    Qext_mie = []
+    om_rat_plot = []
+    for i in range(len(omega_ratio)):
+        om_rat = omega_ratio[i]
+        if om_rat < min_lim or om_rat > max_lim:
+            continue
+        omega = om_rat*omega_p
+        m = cmath.sqrt(eps_m(omega))
+        x = (omega/c) * R * 1e-9/h_reduced
+        mesomie.calc_ab(R*10,      # 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())
+
+        om_rat_plot.append(om_rat)
+    # 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.legend()
+    plt.yscale('log')
+    plt.xlim((0.4, 0.75))
+    plt.ylim((y_min[fig], y_max[fig]))
+    plt.title("R="+str(R))
 plt.show()

src/mesomie.hpp

@@ -100,8 +100,12 @@ void MesoMie<FloatType>::calc_ab(FloatType R,
       Psi_xd(nmax + 1), Zeta_xd(nmax + 1),  //
       Psi_xm(nmax + 1), Zeta_xm(nmax + 1);
 
-  evalPsiZetaD1D3(std::complex<FloatType>(xd), Psi_xd, Zeta_xd, D1_xd, D3_xd);
-  evalPsiZetaD1D3(std::complex<FloatType>(xm), Psi_xm, Zeta_xm, D1_xm, D3_xm);
+  evalPsiZetaD1D3(std::complex<FloatType>(xd),  //
+                  Psi_xd, Zeta_xd,              //
+                  D1_xd, D3_xd);                //
+  evalPsiZetaD1D3(std::complex<FloatType>(xm),  //
+                  Psi_xm, Zeta_xm,              //
+                  D1_xm, D3_xm);
 
   for (int n = 0; n <= nmax; n++) {
     an_[n] = Psi_xd[n] *