remove c++11 std compile flag (now all major compilers has it as default)

CMakeLists.txt

@@ -9,7 +9,6 @@ message("Hostname:  ${HOSTNAME}")
 message("CMake version:  ${CMAKE_VERSION}")
 
 # Select flags.
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
 set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-O3 -g ")
 set(CMAKE_CXX_FLAGS_RELEASE "-O3 -DNDEBUG")
 set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g")

examples/NatureComm-2020-Plasmon–emitter_interactions_at_the_nanoscale/gold/Au_eps0_imag.csv

@@ -0,0 +1,27 @@
+# from Appl. Phys. B (2017) 123:182 DOI 10.1007/s00340-017-6755-2
+450, 5.24
+451.26398815946527, 5.220008833883355
+457.9631806841046, 5.052174226528707
+466.54284830337957, 4.7262628378284015
+473.4771002148484, 4.477926165899576
+481.82170844729393, 4.064031712684867
+487.6634382566586, 3.8056511056511066
+497.10056859120823, 3.360167196569636
+502.2718751014562, 3.1708653719952062
+514.0248444429287, 2.783804940305224
+519.1961509531767, 2.6374375501703557
+530.9491202946492, 2.254280248883969
+536.1204268048972, 2.1092139022169665
+548.8136336936875, 1.7800499048469973
+556.2180043788153, 1.6460424276568517
+567.1482658663847, 1.4372249510644401
+574.5526365515125, 1.300778017372048
+585.2478386522524, 1.1100125188962693
+592.7697390307949, 1.0020259110634342
+603.8175302117791, 0.8628142600018265
+611.4569602837363, 0.7973717735681927
+622.8573405449647, 0.6988827830507756
+630.379240923507, 0.621181575942142
+641.4270321044913, 0.49461895859589156
+648.9489324830338, 0.40968973222133975
+650, 0.398

examples/NatureComm-2020-Plasmon–emitter_interactions_at_the_nanoscale/gold/Au_eps0_real.csv

@@ -0,0 +1,38 @@
+# from Appl. Phys. B (2017) 123:182 DOI 10.1007/s00340-017-6755-2
+450.0, 9.535
+450.84088126317226, 9.561720990497285
+455.9651758960543, 9.724871841367618
+461.25401209971693, 9.896609579125862
+466.0727295297207, 10.052604690922934
+471.00897665313914, 10.18999488112953
+475.357575309484, 10.311642445374954
+480.52888181973196, 10.404667053327335
+485.8177180233946, 10.483380183133198
+490.8714948402278, 10.528103552341074
+495.4551528834021, 10.56209331293906
+500.27387031340584, 10.594294138768731
+506.3854143709716, 10.626494964598402
+513.7897850560993, 10.6467
+519.5487400334208, 10.658695790428073
+525.1901653173277, 10.6666
+530.5965312144051, 10.6663
+536.0028971114824, 10.658
+541.291733315145, 10.647962181818183
+545.9929210517341, 10.626494964598402
+551.0466978685673, 10.594294138768731
+556.1004746854005, 10.56924905201232
+561.036721808819, 10.533470356646019
+566.6781470927258, 10.476224444059937
+572.4371020700473, 10.397511314254075
+578.5486461276131, 10.333109662594733
+585.2478386522524, 10.297330967228433
+597.1183376871397, 10.283019489081912
+603.4649411315349, 10.247240793715612
+609.9290742693448, 10.21146209834931
+615.5704995532517, 10.18999488112953
+621.9171029976469, 10.16852766390975
+628.3812361354568, 10.15
+634.4927801930226, 10.147060446689968
+640.4867945571735, 10.17568340298301
+645.423041680592, 10.18999488112953
+650.0066997237662, 10.20072848973942

examples/NatureComm-2020-Plasmon–emitter_interactions_at_the_nanoscale/gold/read_eps0.py

@@ -0,0 +1,28 @@
+# -*- coding: utf-8 -*-
+# GPL license from Smuthi project
+
+from scipy.interpolate import interp1d
+import numpy as np
+from matplotlib import pyplot as plt
+
+
+def read_eps0(filename, vacuum_wavelength, kind=1):
+    data_np = np.loadtxt(filename, comments='#', delimiter=',')
+    data_wl = data_np[:, 0]
+    eps0 = data_np[:, 1]
+    f = interp1d(data_wl, eps0, kind=kind)
+    data_out = np.transpose(np.vstack(
+        (vacuum_wavelength, f(vacuum_wavelength))))
+    if len(data_out) == 1:
+        return data_out[0]
+    return data_out
+
+
+WLs = np.linspace(450, 650, 1001)
+eps0_real = read_eps0('Au_eps0_real.csv', WLs, kind=2)
+eps0_real1 = read_eps0('Au_eps0_real.csv', WLs, kind=1)
+eps0_imag = read_eps0('Au_eps0_imag.csv', WLs)
+plt.plot(eps0_real[:, 0], eps0_real[:, 1])
+plt.plot(eps0_real1[:, 0], eps0_real1[:, 1])
+plt.plot(eps0_imag[:, 0], eps0_imag[:, 1])
+plt.show()

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

@@ -0,0 +1,226 @@
+#!/usr/bin/env python3
+# -*- coding: UTF-8 -*-
+from glob import glob
+from operator import truediv
+from matplotlib import pyplot as plt
+import cmath
+from scattnlay import mesomie, mie
+import numpy as np
+import scipy.io
+from functools import wraps
+import time
+from eval_spectra import multi_lorentzian, params_count as pc
+import cma
+import sys
+from optical_constants import read_refractive_index_from_yaml as read_nk
+
+isPlot = True
+# isPlot = False
+
+
+def timeit(func):
+    @wraps(func)
+    def timeit_wrapper(*args, **kwargs):
+        start_time = time.perf_counter()
+        result = func(*args, **kwargs)
+        end_time = time.perf_counter()
+        total_time = end_time - start_time
+        print(
+            f'Function {func.__name__}{args} {kwargs} Took {total_time:.4f} seconds')
+        return result
+    return timeit_wrapper
+
+
+min_lim = 0.4
+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')
+from_disk = from_disk[:,
+                      from_disk[0, :] > min_lim
+                      ]
+from_disk = from_disk[:,
+                      from_disk[0, :] < max_lim
+                      ]
+
+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
+
+# omega_p_star = 3.81  # eV #orig
+omega_p_star = 5.81  # eV
+omega = omega_ratio*omega_p_star
+# WL_mkm = 2*np.pi/((omega/c)/h_reduced)*1e6
+
+WL = 2*np.pi/((omega_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)
+# plt.figure('index')
+# plt.plot(index_Ag[:, 0], np.real(index_Ag[:, 1]))
+# plt.plot(index_Ag[:, 0], np.imag(index_Ag[:, 1]))
+# plt.show()
+# sys.exit()
+
+
+def eps_m(omega):
+    return 1 - omega_p * omega_p / (omega*omega + 1j*omega*gamma)
+
+
+n = 1 if isPlot else 1
+Rs = [2.5]*n  # nm
+y_min = [1e-2]*n
+y_max = [1e1]*n
+
+# Rs = [2.5, 5, 10, 25]  # nm
+# y_min = [1e-2, 1e-2, 1e-1, 1e-1]
+# y_max = [1e1, 1e1, 5e1, 5e1]
+
+# for om_rat in omega_ratio:
+R = 0
+
+
+@ timeit
+def run():
+    for fig in range(len(Rs)):
+        global R
+        R = Rs[fig]
+        Qext = []
+        Qext_fit = []
+        om_rat_plot = []
+        x_fit_from_d_perp = np.array([0.09391149,  0.81234806, -0.30526326, -0.01044856,
+                                      0.3121473,   0.55333457, -0.01684191,  0.04727765,
+                                      0.11249052,  0.88368699, -0.04680872, -0.10982548])
+        omega = omega_ratio*omega_p
+        x = (omega/c) * R * 1e-9/h_reduced
+        m_orig = np.sqrt(eps_m(omega))
+        # m = np.sqrt(eps_m(omega))
+
+        # m = index_Ag
+        m = index_Ag[:, 1]
+        # plt.figure('index2')
+        # plt.plot(x, np.real(m), label='orig')
+        # plt.plot(x, np.imag(m), label='orig')
+        # plt.plot(x, np.real(m), label='mod')
+        # plt.plot(x, np.imag(m), label='mod')
+        # plt.legend()
+        # plt.show()
+        # sys.exit()
+
+        Qext = getMesoQext(R, x, m_orig, d_perp)
+
+        def getQfit(x_fit):
+            d_fit = multi_lorentzian(omega_ratio, x_fit)
+            return getMesoQext(R, x, m, d_fit)
+
+        def rms_base(Q_rms, Q_fit):
+            diff_re = np.real(Q_rms - Q_fit)
+            diff_im = np.imag(Q_rms - Q_fit)
+            rms = (np.sum(np.abs(diff_re)**2))
+            rms += (np.sum(np.abs(diff_im)**2))
+            return rms
+
+        def rms(x0):
+            Q_fit = getQfit(x0)
+            return rms_base(Q_rms, Q_fit)
+
+        def rms_x12(x0):
+            x0[:pc] = x0[:pc]+x1
+            Q_fit = getQfit(x0)
+            return rms_base(Q_rms, Q_fit)
+
+        def rms_x1(x0):
+            Q_fit = getQfit(np.hstack((x1, x0)))
+            return rms_base(Q_rms, Q_fit)
+
+        x_fit_from_d_perp = np.array([0.09391149,  0.81234806, -0.30526326, -0.01044856,
+                                      0.3121473,   0.55333457, -0.01684191,  0.04727765,
+                                      0.11249052,  0.88368699, -0.04680872, -0.10982548])
+
+        x0 = np.random.random(pc)
+        Q_rms = Qext
+
+        xfit, es = cma.fmin2(rms, x0, sigma0=2)
+        # xfit = np.array([0.32740616895990493, -0.844352181880013, -
+        #                  0.5682552466755, 0.015925488550861087])
+        x1 = np.copy(xfit)
+        Qext_fit = getQfit(xfit)
+        # print(x1)
+
+        # x0 = np.random.random(pc)
+        # # x0 = np.copy(x1)
+        # x0 = np.hstack((x1, x0))
+        # xfit, es = cma.fmin2(rms, x0, sigma0=0.2)
+        # # xfit = np.array([-9.03308419e-01, -1.77584180e+00, -2.36259518e+01, 3.41837008e+00,
+        #                  5.75724471e-06, , 2.64230613e+00, , 4.72267415e+01, -1.61624064e+00])
+        # # xfit = np.array([0.32740616895990493, -0.844352181880013, -
+        # #  0.5682552466755, 0.015925488550861087])
+        # x2 = np.copy(xfit)
+        # Qext_fit = getQfit(np.hstack((x1, x2)))
+
+        print(xfit)
+
+        om_rat_plot = omega_ratio
+
+        if isPlot:
+            Qext_fit = getQfit(xfit)
+            Qext_no_d = getMesoQext(R, x, m_orig, np.zeros(len(x)))
+            Qext_no_d_mod = getMesoQext(R, x, m, np.zeros(len(x)))
+            d_fit = multi_lorentzian(omega_ratio, xfit)
+            plt.figure(fig)
+            plt.plot(om_rat_plot, Qext_fit,
+                     label='fitted orig using mod', color='gray', lw=4)
+            plt.plot(om_rat_plot, Qext_no_d,
+                     label='d=0 orig', color='gray', lw=4)
+            plt.plot(om_rat_plot, Qext_no_d_mod,
+                     label='d=0 mod', color='gray', lw=4)
+            plt.plot(om_rat_plot, Qext, label='direct orig', color='red', lw=1)
+            plt.legend()
+            # plt.yscale('log')
+            # plt.xlim((0.4, 0.404))
+            # plt.ylim((y_min[fig]*1.9, y_min[fig]*2.3))
+            plt.xlim((0.4, 0.75))
+            # plt.xlim((0.5, 0.6))
+            # plt.ylim((y_min[fig], y_max[fig]))
+            plt.title("Qext for R="+str(R))
+
+            plt.figure("d_param "+str(fig))
+            plt.plot(om_rat_plot, np.real(d_fit), label='re fit', lw=4)
+            plt.plot(om_rat_plot, np.real(d_perp), label='re d_perp')
+            plt.plot(om_rat_plot, np.imag(d_fit), label='im fit', lw=4)
+            plt.plot(om_rat_plot, np.imag(d_perp), label='im d_perp')
+            plt.legend()
+
+    if isPlot:
+        plt.show()
+
+
+# @timeit
+def getMesoQext(R, x, m, d_perp):
+    Qext = np.zeros((len(x)))
+    for i in range(len(x)):
+        mesomie.calc_ab(R*10,      # calc_ab needs R in angstrem
+                        x[i],      # xd
+                        x[i] * m[i],  # xm
+                        1,      # eps_d
+                        m[i] * m[i],  # eps_m
+                        0,      # d_parallel
+                        d_perp[i])      # d_perp
+        mesomie.calc_Q()
+        Qext[i] = mesomie.GetQext()
+    return Qext
+
+
+run()

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

@@ -0,0 +1,223 @@
+#!/usr/bin/env python3
+# -*- coding: UTF-8 -*-
+from glob import glob
+from operator import truediv
+from matplotlib import pyplot as plt
+import cmath
+from scattnlay import mesomie, mie
+import numpy as np
+import scipy.io
+from functools import wraps
+import time
+from eval_spectra import multi_lorentzian, params_count as pc
+import cma
+import sys
+from optical_constants import read_refractive_index_from_yaml as read_nk
+
+isPlot = True
+# isPlot = False
+
+
+def timeit(func):
+    @wraps(func)
+    def timeit_wrapper(*args, **kwargs):
+        start_time = time.perf_counter()
+        result = func(*args, **kwargs)
+        end_time = time.perf_counter()
+        total_time = end_time - start_time
+        print(
+            f'Function {func.__name__}{args} {kwargs} Took {total_time:.4f} seconds')
+        return result
+    return timeit_wrapper
+
+
+min_lim = 0.4
+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')
+from_disk = from_disk[:,
+                      from_disk[0, :] > min_lim
+                      ]
+from_disk = from_disk[:,
+                      from_disk[0, :] < max_lim
+                      ]
+
+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
+
+# omega_p_star = 3.81  # eV #orig
+omega_p_star = 5.81  # eV
+omega = omega_ratio*omega_p_star
+# WL_mkm = 2*np.pi/((omega/c)/h_reduced)*1e6
+
+WL = 2*np.pi/((omega_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)
+# plt.figure('index')
+# plt.plot(index_Ag[:, 0], np.real(index_Ag[:, 1]))
+# plt.plot(index_Ag[:, 0], np.imag(index_Ag[:, 1]))
+# plt.show()
+# sys.exit()
+
+
+def eps_m(omega):
+    return 1 - omega_p * omega_p / (omega*omega + 1j*omega*gamma)
+
+
+n = 1 if isPlot else 1
+Rs = [2.5]*n  # nm
+y_min = [1e-2]*n
+y_max = [1e1]*n
+
+# Rs = [2.5, 5, 10, 25]  # nm
+# y_min = [1e-2, 1e-2, 1e-1, 1e-1]
+# y_max = [1e1, 1e1, 5e1, 5e1]
+
+# for om_rat in omega_ratio:
+R = 0
+
+
+@ timeit
+def run():
+    for fig in range(len(Rs)):
+        global R
+        R = Rs[fig]
+        Qext = []
+        Qext_fit = []
+        om_rat_plot = []
+        x_fit_from_d_perp = np.array([0.09391149,  0.81234806, -0.30526326, -0.01044856,
+                                      0.3121473,   0.55333457, -0.01684191,  0.04727765,
+                                      0.11249052,  0.88368699, -0.04680872, -0.10982548])
+        omega = omega_ratio*omega_p
+        x = (omega/c) * R * 1e-9/h_reduced
+        # m_orig = np.sqrt(eps_m(omega))
+        m = np.sqrt(eps_m(omega))
+
+        # m = index_Ag[:, 1]
+
+        # plt.figure('index2')
+        # plt.plot(x, np.real(m), label='orig')
+        # plt.plot(x, np.imag(m), label='orig')
+        # plt.plot(x, np.real(m), label='mod')
+        # plt.plot(x, np.imag(m), label='mod')
+        # plt.legend()
+        # plt.show()
+        # sys.exit()
+
+        Qext = getMesoQext(R, x, m, d_perp)
+
+        def getQfit(x_fit):
+            d_fit = multi_lorentzian(omega_ratio, x_fit)
+            return getMesoQext(R, x, m, d_fit)
+
+        def rms_base(Q_rms, Q_fit):
+            diff_re = np.real(Q_rms - Q_fit)
+            diff_im = np.imag(Q_rms - Q_fit)
+            rms = (np.sum(np.abs(diff_re)**2))
+            rms += (np.sum(np.abs(diff_im)**2))
+            return rms
+
+        def rms(x0):
+            Q_fit = getQfit(x0)
+            return rms_base(Q_rms, Q_fit)
+
+        def rms_x12(x0):
+            x0[:pc] = x0[:pc]+x1
+            Q_fit = getQfit(x0)
+            return rms_base(Q_rms, Q_fit)
+
+        def rms_x1(x0):
+            Q_fit = getQfit(np.hstack((x1, x0)))
+            return rms_base(Q_rms, Q_fit)
+
+        x_fit_from_d_perp = np.array([0.09391149,  0.81234806, -0.30526326, -0.01044856,
+                                      0.3121473,   0.55333457, -0.01684191,  0.04727765,
+                                      0.11249052,  0.88368699, -0.04680872, -0.10982548])
+
+        x0 = np.random.random(pc)
+        Q_rms = Qext
+
+        xfit, es = cma.fmin2(rms, x0, sigma0=2)
+        # xfit = np.array([0.32740616895990493, -0.844352181880013, -
+        #                  0.5682552466755, 0.015925488550861087])
+        x1 = np.copy(xfit)
+        Qext_fit = getQfit(xfit)
+        # print(x1)
+
+        # x0 = np.random.random(pc)
+        # # x0 = np.copy(x1)
+        # x0 = np.hstack((x1, x0))
+        # xfit, es = cma.fmin2(rms, x0, sigma0=0.2)
+        # # xfit = np.array([-9.03308419e-01, -1.77584180e+00, -2.36259518e+01, 3.41837008e+00,
+        #                  5.75724471e-06, , 2.64230613e+00, , 4.72267415e+01, -1.61624064e+00])
+        # # xfit = np.array([0.32740616895990493, -0.844352181880013, -
+        # #  0.5682552466755, 0.015925488550861087])
+        # x2 = np.copy(xfit)
+        # Qext_fit = getQfit(np.hstack((x1, x2)))
+
+        print(xfit)
+
+        om_rat_plot = omega_ratio
+
+        if isPlot:
+            Qext_fit = getQfit(xfit)
+            Qext_no_d = getMesoQext(R, x, m, np.zeros(len(x)))
+            d_fit = multi_lorentzian(omega_ratio, xfit)
+            plt.figure(fig)
+            plt.plot(om_rat_plot, Qext_fit,
+                     label='fitted', color='gray', lw=4)
+            plt.plot(om_rat_plot, Qext_no_d,
+                     label='d=0', color='gray', lw=4)
+            plt.plot(om_rat_plot, Qext, label='direct', color='red', lw=1)
+            plt.legend()
+            # plt.yscale('log')
+            # plt.xlim((0.4, 0.404))
+            # plt.ylim((y_min[fig]*1.9, y_min[fig]*2.3))
+            plt.xlim((0.4, 0.75))
+            # plt.xlim((0.5, 0.6))
+            # plt.ylim((y_min[fig], y_max[fig]))
+            plt.title("Qext for R="+str(R))
+
+            plt.figure("d_param "+str(fig))
+            plt.plot(om_rat_plot, np.real(d_fit), label='re fit', lw=4)
+            plt.plot(om_rat_plot, np.real(d_perp), label='re d_perp')
+            plt.plot(om_rat_plot, np.imag(d_fit), label='im fit', lw=4)
+            plt.plot(om_rat_plot, np.imag(d_perp), label='im d_perp')
+            plt.legend()
+
+    if isPlot:
+        plt.show()
+
+
+# @timeit
+def getMesoQext(R, x, m, d_perp):
+    Qext = np.zeros((len(x)))
+    for i in range(len(x)):
+        mesomie.calc_ab(R*10,      # calc_ab needs R in angstrem
+                        x[i],      # xd
+                        x[i] * m[i],  # xm
+                        1,      # eps_d
+                        m[i] * m[i],  # eps_m
+                        0,      # d_parallel
+                        d_perp[i])      # d_perp
+        mesomie.calc_Q()
+        Qext[i] = mesomie.GetQext()
+    return Qext
+
+
+run()

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

@@ -50,7 +50,11 @@ index_Ag = read_nk('Ag-Johnson-1972.yml', WL, kind=1)
 eps_Ag = index_Ag**2
 # print(index_Ag)
 
-factor = 1
+factor = 1  # TODO see app phys b 2017 fig 3, Karpov
+
+# 1. Нужна правильная eps (см выше) + размерный фактор в gamma
+# 2. Варьируем концентрацию
+# 3. На выходе спектр с красным сдвигом. jpcc
 
 
 def eps_m(omega, eps_inf, omega_p_local):

setup.py

@@ -49,7 +49,8 @@ ext_dp = Extension("scattnlay_dp",
                    ["src/pb11-wrapper.cc"],
                    language="c++",
                    include_dirs=[np.get_include(), pb.get_include()],
-                   extra_compile_args=['-std=c++11'])
+                   #    extra_compile_args=['-std=c++11']
+                   )
 # extra_compile_args=['-std=c++11', '-O3',
 #                     '-mavx2', '-mfma',
 #                     '-finline-limit=1000000', '-ffp-contract=fast']),
@@ -57,7 +58,8 @@ ext_mp = Extension("scattnlay_mp",
                    ["src/pb11-wrapper.cc"],
                    language="c++",
                    include_dirs=[np.get_include(), pb.get_include()],
-                   extra_compile_args=['-std=c++11', '-DMULTI_PRECISION=100'])
+                   extra_compile_args=[  # '-std=c++11',
+                       '-DMULTI_PRECISION=100'])
 # extra_compile_args=['-std=c++11', '-O3',
 #                     '-mavx2', '-mfma',
 #                     '-finline-limit=1000000', '-ffp-contract=fast',