Konstantin Ladutenko 7 年之前
父節點
當前提交
332d2e26d3
共有 1 個文件被更改,包括 180 次插入0 次删除
  1. 180 0
      plasmon-modal-efficiency.py

+ 180 - 0
plasmon-modal-efficiency.py

@@ -0,0 +1,180 @@
+#!/usr/bin/env python3
+# -*- coding: UTF-8 -*-
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+c = 299792458.0
+eps_0 = 8.854187817e-12 # F/m
+pi = np.pi
+verbose = 6
+# r of monitor
+r = 146.513e-9
+debug = True
+def read_data(dirname):
+    data = {}
+    WLs = []
+    for r,d,f in os.walk(dirname):
+        for fname in f:
+            WLs.append(fname)
+    for fname in WLs:            
+        fdata = np.transpose(
+                    np.genfromtxt(dirname+"/"+fname, delimiter=", ",skip_header=1
+                                      ,dtype=None, encoding = None
+                                      , converters={0: lambda s: complex(s),
+                                                    1: lambda s: complex(s),
+                                                    2: lambda s: complex(s.replace('i', 'j')),
+                                                    3: lambda s: complex(s.replace('i', 'j')),
+                                                    4: lambda s: complex(s.replace('i', 'j')),
+                                                    5: lambda s: complex(s.replace('i', 'j')),
+                                                    6: lambda s: complex(s.replace('i', 'j')),
+                                                    7: lambda s: complex(s.replace('i', 'j')),
+                                                    8: lambda s: complex(s.replace('i', 'j'))
+                                                    }
+                                      )
+                        )
+        data[float(fname[2:-4])]=fdata
+        if debug: break
+    return data
+
+
+def find_nearest(array,value):
+    idx = (np.abs(array-value)).argmin()
+    return array[idx],idx
+
+
+def get_WLs_idx(WLs, data):
+    dist = 1 #mkm
+    mmedia = 1 # vacuum
+    shift = 1 # one mesh step
+    WLs_idx = []
+    for wl in WLs:
+        val, idx = find_nearest(data[dist][mmedia][shift][0,:],wl*1e-9)
+        WLs_idx.append(idx)
+    return WLs_idx
+
+
+# def check_field_match(data_in_air, data_in_gold,wl_idx,z_vec,kappa1,kappa2,eps2): 
+#         z = z_vec[i]*1e-9
+#         if verbose > 8: print("z =",z)
+#         H1_0 = H1[i]/np.exp(-kappa1[wl_idx]*z)
+#         H2_0 = H2[i]/np.exp(-kappa2[wl_idx]*z)
+#         E1_0 = E1[i]/np.exp(-kappa1[wl_idx]*z)
+#         E2_0 = E2[i]/np.exp(-kappa2[wl_idx]*z)
+#         E2_0e = E2[i]/np.exp(-kappa2[wl_idx]*z)*eps2[wl_idx]
+#         if verbose > 8:
+#             print("H0 air  (%5.4g %+5.4gj)"%(np.real(H1_0), np.imag(H1_0)),
+#                   " from H1 (%5.4g %+5.4gj)"%(np.real(H1[i]), np.imag(H1[i])))
+
+
+def analyze(data,wl):
+    # print(data[0,:]) # all z values
+    #data = "z, dip.power, Ex, Ey, Ez, Hx, Hy, Hz, n_Au"
+    #        0, 1        , 2 , 3 , 4 , 5 , 6 , 7 , 8   "
+    lambd = wl
+    omega = 2*pi*c/lambd
+
+    eps_d = complex(1) # air, z>0
+    eps_m = data[8,0]**2 # metal, z<0
+
+    dip_power = data[1,0]
+
+    z = data[0,:]
+
+    idx_d = np.nonzero(z>1e-10)
+    idx_0 = np.nonzero(np.logical_and(z<=1e-10, z>=-1e-10))
+    idx_m = np.nonzero(z<-1e-10)
+
+    z_d = z[idx_d]
+    z_0 = z[idx_0]
+    z_m = z[idx_m]
+
+    if (not np.array_equal(np.hstack((z_m, z_0, z_d)), z)):
+        print("ERROR! loosing z values!")
+        raise
+    
+
+    Ex = data[2,:]
+
+    Ex_m = data[2,idx_m][0]
+    Ey_m = data[3,idx_m][0]
+    Ez_m = data[4,idx_m][0]
+    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])
+
+    Ex_d = data[2,idx_d][0]
+    Ey_d = data[3,idx_d][0]
+    Ez_d = data[4,idx_d][0]
+    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])
+
+    k_0 = omega/c #air
+    k_sp = k_0*np.sqrt(eps_d*eps_m/(eps_d+eps_m))  # eq5, supmat
+
+    chi_d = np.sqrt( eps_d*k_0**2 - k_sp**2 ) # desc. after eq6c, supmat
+    chi_m = np.sqrt( eps_m*k_0**2 - k_sp**2 ) # desc. after eq6c, supmat
+
+    h_sp_d = np.exp(1j*chi_d*z_d) # eq6a, supmat
+    e_sp_x_d = chi_d/(omega*eps_0*eps_d)*np.exp(1j*chi_d*z_d) # eq6b, supmat
+    e_sp_z_d = k_sp/(omega*eps_0*eps_d)*np.exp(1j*chi_d*z_d) # eq6c, supmat
+
+    h_sp_m = np.exp(1j*-chi_m*z_m) # eq6a, supmat
+    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)
+
+    # 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])),
+    #           " from H1 (%5.4g %+5.4gj)"%(np.real(H1[0][wl_idx]), np.imag(H1[0][wl_idx])))
+    # #plasmon_power = 1.0/2.0 * np.real( E1[0] * np.conj(H1[0]))  # TODO check minus sign!!
+    # plasmon_power = -1.0/2.0 * 2.0*np.pi*R * ( # TODO check minus sign!!
+    #     np.real( E1_0 * np.conj(H1_0) )
+    #         / (2.0 * np.real(kappa1))
+    #     +
+    #     np.real( E2_0 * np.conj(H1_0) )
+    #         / (2.0 * np.real(kappa2))        
+    #     )* np.exp( 2.0*np.imag(k_spp)*R )  # TODO check minus sign!!
+    # #print(np.abs(plasmon_power/ dip_power))
+    # eta0 = plasmon_power[0]/ dip_power[0] *100
+    # ppw = plasmon_power[0]
+    # print("\n")
+    # print(dirname)
+    # print("Power: plasmon %4.3g W of dipoles %4.3g W, efficiency %5.3g%%  from:"%(ppw, float(np.abs(dip_power[0])),float(np.abs( eta0))), ppw, eta0)
+    # plt.plot(lambd*1e9, plasmon_power/ dip_power)
+    # plt.ylim(0,0.04)
+    # plt.xlim(550,800)
+
+    # #plt.plot(lambd*1e9, np.real(eps2))
+    # # plt.plot(lambd*1e9, np.real(k_spp))
+    # # plt.plot(lambd*1e9, k_0)
+    # #plt.semilogy(lambd*1e9, np.absolute(plasmon_power/ dip_power))
+    # # # legend = []
+    # # # legend.append(zshift[shift]+"@"+str(WLs[i])+" nm")
+    # # # plt.legend(legend)
+    # # # #plt.xlabel(r'THz')
+    # plt.xlabel(r'$\lambda$, nm')
+    # plt.ylabel(r'$P_{spp}/P_{dipole}$',labelpad=-5)
+    # #plt.title(' R = '+str(core_r)+' nm')
+    # plt.savefig(dirname+"_power_ratio."+file_ext)
+    # plt.clf()
+    # plt.close()
+    
+file_ext="pdf"
+
+dirname="bigourdan-Au-sub-Cyl-dipole-W.fsp.1D.monitor_1.results"
+
+def main ():
+    data = read_data(dirname)
+    for wl in data:
+        analyze(data[wl],wl)
+
+
+main()