fit d-params cleanup

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

@@ -1,16 +1,14 @@
 from ast import Constant
 from numba import njit, prange, complex128, float64
 import numpy as np
+from scattnlay import mesomie
 
-
-# omega_init = True
 params_count = 4
 
 
 @njit(complex128(float64, float64[:]),
       fastmath=True, cache=True)
 def lorentzian(omega, xvec):
-    # global omega_init
     pc = params_count
     res = 0
     poles = len(xvec)//params_count
@@ -19,31 +17,17 @@ def lorentzian(omega, xvec):
         gamma = xvec[pc*i+0] / (factor**i)
         omega_n = xvec[pc*i+1] / (factor**i)
         f = (xvec[pc*i+2] + 1j*xvec[pc*i+3]) / (factor**i)
-        # if (omega_init):
-        #     print('init', gamma, omega_n, f)
         if (np.abs(f) == 0):
             return res
         res = res + f / (omega * (omega + 1j*gamma) - omega_n**2)
-    # omega_init = False
     return res
 
 
-# @njit(complex128[:](float64[:], float64[:]),
-#       parallel=True, fastmath=True, cache=True)
-def spectra(omega, xvec):
-    # global omega_init
-    # omega_init = True
-    xvec_in = xvec
+@njit(complex128[:](float64[:], float64[:]),
+      fastmath=True, cache=True)
+def multi_lorentzian(omega, xvec):
     poles = len(xvec)//params_count
-    # for i in range(poles):
-    #     for j in range(2):
-    #         if (xvec_in[i+j] < 0):
-    #             xvec_in[i+j] = 0
-    #     # if (xvec_in[i+2] > 0):
-    #         # xvec_in[i+2] = 0
-
     val = np.zeros(omega.size, dtype=np.cdouble)
-    # print('in eval', xvec_in)
     for i in range(omega.size):
-        val[i] = lorentzian(omega[i], xvec_in)
+        val[i] = lorentzian(omega[i], xvec)
     return val

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

@@ -3,16 +3,24 @@
 from matplotlib import pyplot as plt
 import numpy as np
 import cma
-from eval_spectra import spectra, params_count as pc
+from eval_spectra import multi_lorentzian, params_count as pc
+import sys
 
 from_disk = np.loadtxt('rs4-d_perp_interpolated.txt')
+# from_disk = from_disk[:,
+#                       from_disk[0, :] > 0.4
+#                       ]
+# from_disk = from_disk[:,
+#                       from_disk[0, :] < 0.8
+#                       ]
+
 step = 5
 omega_ratio = np.copy(from_disk[0, ::step])
 d_rs4 = from_disk[1, ::step] + 1j*from_disk[2, ::step]
 
 
 def rms(x0):
-    d_fit = spectra(omega_ratio, x0)
+    d_fit = multi_lorentzian(omega_ratio, x0)
     diff_re = np.real(d_rms - d_fit)
     diff_im = np.imag(d_rms - d_fit)
     # rms = np.sqrt(np.sum(np.abs(diff_re)**2))
@@ -31,7 +39,7 @@ x1 = np.copy(x)
 
 
 x0 = np.random.random(pc)
-d_rms = d_rs4 - spectra(omega_ratio, x1)
+d_rms = d_rs4 - multi_lorentzian(omega_ratio, x1)
 x, es = cma.fmin2(rms, x0, sigma0=2)
 # x = np.array([0.0019369902204593222, 0.9978752165162739,
 #              0.05801769075873917, -0.032110084386726336])
@@ -47,7 +55,7 @@ x12 = x
 
 
 x0 = np.random.random(pc)
-d_rms = d_rs4 - spectra(omega_ratio, x12)
+d_rms = d_rs4 - multi_lorentzian(omega_ratio, x12)
 x, es = cma.fmin2(rms, x0, sigma0=0.2)
 # x = np.array([0.1434266344187499, 0.5188802822956783, -
 #              0.00950433613285183, 0.013585684987833985])
@@ -62,13 +70,20 @@ x, es = cma.fmin2(rms, x0, sigma0=0.02)
 #               0.34432793, 0.64182428, -0.0803532, 0.04641341])
 x123 = x
 
+#               [ 0.09349663 -0.8188804  -0.36874431 -0.08079194
+#                 0.25712096  0.56012451 -0.02089828  0.03133694]
+
+#               [ 0.09391149  0.81234806 -0.30526326 -0.01044856
+#                 0.3121473   0.55333457 -0.01684191  0.04727765
+#                 0.11249052  0.88368699 -0.04680872 -0.10982548]
+print(x)
 
-d_fit = spectra(omega_ratio, x)
+d_fit = multi_lorentzian(omega_ratio, x)
 plt.figure('rs4')
 plt.plot(omega_ratio, np.real(d_rms), label='re d')
 plt.plot(omega_ratio, np.imag(d_rms), label='im d')
 plt.plot(omega_ratio, np.real(d_fit), label='re d fit', alpha=0.2, lw=3)
 plt.plot(omega_ratio, np.imag(d_fit), label='im d fit', alpha=0.2, lw=3)
-plt.axhline(y=0.0, color='black', linestyle='-', lw=1)
+plt.axhline(y=0.0, color='black', linestyle='-', lw=1, alpha=0.2)
 plt.legend()
 plt.show()

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

@@ -27,7 +27,7 @@ def eps_m(omega):
     return 1 - omega_p * omega_p / (omega*omega + 1j*omega*gamma)
 
 
-Rs = [2.5, 5, 10, 25]
+Rs = [2.5, 5, 10, 25]  # nm
 y_min = [1e-2, 1e-2, 1e-1, 1e-1]
 y_max = [1e1, 1e1, 5e1, 5e1]
 
@@ -44,7 +44,7 @@ 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,      # R in angstrem
+        mesomie.calc_ab(R*10,      # calc_ab needs R in angstrem
                         x,      # xd
                         x * m,  # xm
                         1,      # eps_d