update

plasmon-modal-efficiency.py

@@ -3,14 +3,15 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import os
-import scipy.special.hankel2 as H2n
+from scipy.special import hankel2 as H2n
 c = 299792458.0
 eps_0 = 8.854187817e-12 # F/m
 pi = np.pi
 verbose = 6
 # r of monitor
 r = 146.513e-9
-debug = True
+#debug = True
+debug = False
 def read_data(dirname):
     data = {}
     WLs = []
@@ -102,8 +103,8 @@ def analyze(data,wl):
     Hx_m = data[5,idx_m][0]
     Hy_m = data[6,idx_m][0]
     Hz_m = data[7,idx_m][0]
-    E_m = np.array([Ex_m,Ey_m,Ez_m])
-    H_m = np.array([Hx_m,Hy_m,Hz_m])
+    E_m = np.transpose(np.array([Ex_m,Ey_m,Ez_m]))
+    H_m = np.transpose(np.array([Hx_m,Hy_m,Hz_m]))
 
     Ex_d = data[2,idx_d][0]
     Ey_d = data[3,idx_d][0]
@@ -111,9 +112,9 @@ def analyze(data,wl):
     Hx_d = data[5,idx_d][0]
     Hy_d = data[6,idx_d][0]
     Hz_d = data[7,idx_d][0]
-    E_d = np.array([Ex_d,Ey_d,Ez_d])
-    H_d = np.array([Hx_d,Hy_d,Hz_d])
-
+    E_d = np.transpose(np.array([Ex_d,Ey_d,Ez_d]))
+    H_d = np.transpose(np.array([Hx_d,Hy_d,Hz_d]))
+    
     k_0 = omega/c #air
     k_sp = k_0*np.sqrt(eps_d*eps_m/(eps_d+eps_m))  # eq5, supmat
 
@@ -128,11 +129,40 @@ def analyze(data,wl):
     e_sp_x_m = -chi_m/(omega*eps_0*eps_m)*np.exp(1j*-chi_m*z_m) # eq6b, supmat
     e_sp_z_m = k_sp/(omega*eps_0*eps_m)*np.exp(1j*-chi_m*z_m) # eq6c, supmat
 
-    if verbose > 5:
-        print("r =",r)
+    zero_m = np.zeros(len(h_sp_m))
+    zero_d = np.zeros(len(h_sp_d))
+
+    E_minus_sp_0_m = np.transpose([1j*H2n(1,k_sp*r)*e_sp_x_m,
+                      zero_m,
+                      H2n(0,k_sp*r)*e_sp_z_m]) # eq11, supmat, replace E_plus to E_minus and H1n to H2n
+    H_minus_sp_0_m = np.transpose([zero_m,
+                      1j*H2n(1,k_sp*r)*h_sp_m,
+                      zero_m]) # eq11, supmat, replace E_plus to E_minus and H1n to H2n
 
-    print(h_sp_m)
+    E_minus_sp_0_d = np.transpose([1j*H2n(1,k_sp*r)*e_sp_x_d,
+                      zero_d,
+                      H2n(0,k_sp*r)*e_sp_z_d]) # eq11, supmat, replace E_plus to E_minus and H1n to H2n
+    H_minus_sp_0_d = np.transpose([zero_d,
+                      1j*H2n(1,k_sp*r)*h_sp_d,
+                      zero_d]) # eq11, supmat, replace E_plus to E_minus and H1n to H2n
 
+    # E_m H_m E_d H_d
+
+    N_sp_0 = (((-1)**0) * (4.0j/(omega*eps_0*k_sp))
+                * (eps_d**2 - eps_m**2) / ((eps_m*eps_d)**(3/2)) )
+    
+    tmp_m = np.cross(E_minus_sp_0_m,H_m) - np.cross(E_m, H_minus_sp_0_m)
+    radail_pojeciton_m = np.transpose(tmp_m)[0]
+    integrand_m = (2*pi/N_sp_0)*radail_pojeciton_m*r 
+
+    tmp_d = np.cross(E_minus_sp_0_d,H_d) - np.cross(E_d, H_minus_sp_0_d)
+    radail_pojeciton_d = np.transpose(tmp_d)[0]
+    integrand_d = (2*pi/N_sp_0)*radail_pojeciton_d*r 
+    
+    A_sp_0_m = np.trapz(integrand_m, z_m)
+    A_sp_0_d = np.trapz(integrand_d, z_d)
+    A_sp_0 = A_sp_0_m + A_sp_0_d
+    return np.absolute(A_sp_0)**2
     # print("S from full field",np.real(np.cross(E,np.conj(H))))
 
     #     print("H0 air  (%5.4g %+5.4gj)"%(np.real(H1_0[wl_idx]), np.imag(H1_0[wl_idx])),
@@ -172,12 +202,44 @@ def analyze(data,wl):
     
 file_ext="pdf"
 
-dirname="bigourdan-Au-sub-Cyl-dipole-W.fsp.1D.monitor_1.results"
-
 def main ():
+    if verbose > 5:
+        print("r =",r)
+    dirname="bigourdan-Au-sub-dipole-W.fsp.1D.monitor_1.results"
     data = read_data(dirname)
+    WLs = []
+    A2 = []
     for wl in data:
-        analyze(data[wl],wl)
+        WLs.append(wl)
+        A2.append(analyze(data[wl],wl))
 
+        #print(WLs)
+    WLs1 = np.array(WLs)
+    A21 = np.array(A2)
+    dirname="bigourdan-Au-sub-Cyl-dipole-W.fsp.1D.monitor_1.results"
+    data = read_data(dirname)
+    WLs = []
+    A2 = []
+    for wl in data:
+        WLs.append(wl)
+        A2.append(analyze(data[wl],wl))
+
+        #print(WLs)
+    WLs2 = np.array(WLs)
+    A22 = np.array(A2)
+    
+    # data = np.vstack((WLs,A2))
+    # print(np.sort(data))
+    
+    plt.plot(WLs1*1e9, A21*275, linestyle='None', marker='o', color="black",label="x 275, no ant.")
+    plt.plot(WLs2*1e9, A22, linestyle='None', marker='*', color="red", label="with antena")
+    plt.legend()
+    plt.xlabel(r'$\lambda$, nm')
+    plt.ylim(0,0.2)
+    plt.ylabel(r'$|A_{sp}|^2$',labelpad=-1)
+    #plt.title(dirname)
+    plt.savefig(dirname+"_A2."+file_ext)
+    plt.clf()
+    plt.close()
 
 main()