initial match of supp fig 5 and simulation results

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

@@ -14,7 +14,7 @@ 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]), 100)
+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])
 

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

@@ -1,37 +1,61 @@
 #!/usr/bin/env python3
 # -*- coding: UTF-8 -*-
-from cProfile import label
-from matplotlib import markers, pyplot as plt
+from matplotlib import pyplot as plt
+import cmath
+from scattnlay import mesomie, mie
 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()]
-
-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]), 100)
-im_d_y = interpolate.interp1d(im_d[:, 0],  im_d[:, 1])
-re_d_y = interpolate.interp1d(re_d[:, 0],  re_d[:, 1])
-
-data = np.array([x.T, re_d_y(x).T, im_d_y(x).T])
-np.savetxt('rs4-d_perp.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')
-plt.plot(from_disk[0, :], from_disk[1, :], label='re d')
-plt.plot(from_disk[0, :], from_disk[2, :], label='im d')
+min_lim = 0.4
+max_lim = 0.8
+omega_ratio = from_disk[0, :]
+d_perp = from_disk[1, :] + 1j*from_disk[2, :]
+
+c = 299792458  # m/s
+h_reduced = 6.5821e-16  # eV s
+omega_p = 5.9  # eV
+gamma = 0.1  # eV
+eps_d = 1
+
+
+def eps_m(omega):
+    return 1 - omega_p * omega_p / (omega*omega + 1j*omega*gamma)
+
+
+Rs = [2.5, 5, 10, 25]
+R = 5
+
+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))
 plt.show()

scattnlay/__init__.py

@@ -33,3 +33,5 @@
 from scattnlay.main import scattcoeffs, expancoeffs, scattnlay, fieldnlay
 from scattnlay.main import mie, mie_mp
 from scattnlay.main import mesomie  # , mesomie_mp
+# import scattnlay_dp
+# print('using scattnlay from ', scattnlay_dp.__file__)

src/mesomie.hpp

@@ -82,13 +82,14 @@ void MesoMie<FloatType>::calc_Q() {
 //******************************************************************************
 template <typename FloatType>
 void MesoMie<FloatType>::calc_ab(FloatType R,
-                                 std::complex<FloatType> xd,
+                                 FloatType xd,
                                  std::complex<FloatType> xm,
                                  std::complex<FloatType> eps_d,
                                  std::complex<FloatType> eps_m,
                                  std::complex<FloatType> d_parallel,
                                  std::complex<FloatType> d_perp) {
-  x_ = R;
+  x_ = xd;
+  //   std::cout << "xd: " << xd << "  R: " << R << std::endl;
   double xx = static_cast<double>(x_);
   int nmax = std::round(xx + 11 * std::pow(xx, (1.0 / 3.0)) + 1);
   an_.resize(nmax + 1, static_cast<FloatType>(0.0));
@@ -106,36 +107,40 @@ void MesoMie<FloatType>::calc_ab(FloatType R,
     an_[n] = Psi_xd[n] *
              (                                                               //
                  eps_m * xd * D1_xd[n] - eps_d * xm * D1_xm[n] +             //
-                 (eps_m - eps_d) *                                           //
+                 (                                                           //
+                     (eps_m - eps_d) *                                       //
                      (                                                       //
                          static_cast<FloatType>(n * (n + 1)) * d_perp +      //
                          xd * D1_xd[n] * xm * D1_xm[n] * d_parallel          //
                          ) /                                                 //
                      R                                                       //
+                     )                                                       //
                  ) /                                                         //
              (                                                               //
                  Zeta_xd[n] *                                                //
                  (                                                           //
                      eps_m * xd * D3_xd[n] - eps_d * xm * D1_xm[n] +         //
-                     (eps_m - eps_d) *                                       //
+                     (                                                       //
+                         (eps_m - eps_d) *                                   //
                          (                                                   //
                              static_cast<FloatType>(n * (n + 1)) * d_perp +  //
                              xd * D3_xd[n] * xm * D1_xm[n] * d_parallel      //
                              ) /                                             //
                          R                                                   //
+                         )                                                   //
                      )                                                       //
              );
     bn_[n] = Psi_xd[n] *
-             (                                                          //
-                 xd * D1_xd[n] - xm * D1_xm[n] +                        //
-                 (xm * xm - xd * xd) * d_parallel / R                   //
-                 ) /                                                    //
-             (                                                          //
-                 Zeta_xd[n] * (                                         //
-                                  xd * D3_xd[n] - xm * D1_xm[n] +       //
-                                  (xm * xm - xd * xd) * d_parallel / R  //
-                                  )                                     //
-             );                                                         //
+             (                                                            //
+                 xd * D1_xd[n] - xm * D1_xm[n] +                          //
+                 ((xm * xm - xd * xd) * d_parallel / R)                   //
+                 ) /                                                      //
+             (                                                            //
+                 Zeta_xd[n] * (                                           //
+                                  xd * D3_xd[n] - xm * D1_xm[n] +         //
+                                  ((xm * xm - xd * xd) * d_parallel / R)  //
+                                  )                                       //
+             );                                                           //
   }
 }
 

src/nmie.hpp

@@ -574,7 +574,7 @@ class MesoMie {
   }
 
   void calc_ab(FloatType R,
-               std::complex<FloatType> xd,
+               FloatType xd,
                std::complex<FloatType> xm,
                std::complex<FloatType> eps_d,
                std::complex<FloatType> eps_m,

tests/test_bulk_sphere.cc

@@ -87,8 +87,8 @@ TEST(BulkSphere, MesoMieDu) {
   for (const auto& data : parameters_and_results) {
     auto x = std::get<0>(data);
     auto m = std::get<1>(data);
-    mesomie.calc_ab(x,        // R
-                    {x, 0},   // xd
+    mesomie.calc_ab(1e-10,    // R
+                    x,        // xd
                     x * m,    // xm
                     {1, 0},   // eps_d
                     m * m,    // eps_m

tests/test_py.py

@@ -21,6 +21,7 @@ test_cases = [
     [10000, 1.33+1e-5j, 2.004089, 1.723857],  # ,'f'],
     [10000, 1.5+1j, 2.004368, 1.236574],  # ,'j'],
     [10000, 10+10j, 2.005914, 1.795393],  # ,'m'],
+    # [1.8263846985116234, 0.02867488311561525+1.2957040351341687j, 3, 3]
 ]
 
 
@@ -29,8 +30,6 @@ class TestStringMethods(unittest.TestCase):
     def test_bulk_mesomie(self):
         tol = 3e-7
         for solver in [mesomie]:
-            if solver is None:
-                continue
             print('Using solver: ', solver)
             for case in test_cases:
                 x = case[0]
@@ -45,6 +44,7 @@ class TestStringMethods(unittest.TestCase):
                 solver.calc_Q()
                 Qext = solver.GetQext()
                 Qsca = solver.GetQsca()
+                # print(x, m, Qext)
                 self.assertTrue((case[2]-Qext)/Qext < tol)
                 self.assertTrue((case[3]-Qsca)/Qsca < tol)