#!/usr/bin/env python3 # -*- coding: UTF-8 -*- from scipy.special import hankel2 as H2n import matplotlib.pyplot as plt import numpy as np import os import scipy.io c = 299792458.0 eps_0 = 8.854187817e-12 # F/m pi = np.pi # r of monitor #r = 146.513e-9 r = 800e-9 # debug = True # verbose = 7 debug = False verbose = 6 def read_data_mat(fname): #data = "z, dip.power, Ex, Ey, Ez, Hx, Hy, Hz, n_Au" # 0, 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 " mat = scipy.io.loadmat(fname) lambd = np.reshape(mat["lambda"],(-1)) dippower = np.reshape(mat["dippower"],(-1)) vacpower = np.reshape(mat["vacpower"],(-1)) z = np.reshape(mat["z"],(-1)) n_Au = np.reshape(mat["n_fdtd"],(-1)) onez = np.ones((len(z))) data1 = {} r1 = mat["mon1_E"][0,0][4][0,0] E1 = mat["mon1_E"][0,0][0] H1 = mat["mon1_H"][0,0][0] data2 = {} r2 = mat["mon2_E"][0,0][4][0,0] E2 = mat["mon2_E"][0,0][0] H2 = mat["mon2_H"][0,0][0] for i in range(len(lambd)): fdata = np.vstack(( z.astype(np.complex128), dippower[i]*onez.astype(np.complex128) ,E1[:,0,i], E1[:,1,i], E1[:,2,i] ,H1[:,0,i], H1[:,1,i], H1[:,2,i] ,n_Au[i]*onez.astype(np.complex128), dippower[i]*onez.astype(np.complex128) )) data1[lambd[i]]=fdata fdata = np.vstack(( z.astype(np.complex128), dippower[i]*onez.astype(np.complex128) ,E2[:,0,i], E2[:,1,i], E2[:,2,i] ,H2[:,0,i], H2[:,1,i], H2[:,2,i] ,n_Au[i]*onez.astype(np.complex128), dippower[i]*onez.astype(np.complex128) )) data2[lambd[i]]=fdata if debug: break return ((r1,data1),(r2,data2)) 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 cross(a, b): c = [a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]] return c def analyze(data,wl): # print(data[0,:]) # all z values #data = "z, dip.power, Ex, Ey, Ez, Hx, Hy, Hz, n_Au, vac. power" # 0, 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9" 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] vac_power = data[9,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_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.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] 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.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 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 # # TODO!!! alt # chi_d = np.sqrt( k_sp**2 - eps_d*k_0**2 ) # desc. after eq6c, supmat # chi_m = np.sqrt( k_sp**2 - eps_m*k_0**2 ) # desc. after eq6c, supmat mul = -1 #TODO!!! minus??? h_sp_d = np.exp(mul*1j*chi_d*z_d) # eq6a, supmat e_sp_x_d = chi_d/(omega*eps_0*eps_d)*np.exp(mul*1j*chi_d*z_d) # eq6b, supmat e_sp_z_d = -k_sp/(omega*eps_0*eps_d)*np.exp(mul*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 > 8: print("WL =",1e9*wl) # if 1e9*wl > 763 and 1e9*wl < 766: print(1j*1e9/chi_d) print(1j*1e9/chi_m) print("_d") print(1e9*z_d) print(h_sp_d) print("_m") print(1e9*z_m) print(h_sp_m) 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 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) #tmp_m = np.cross(E_minus_sp_0_m,np.conj(H_m)) + np.cross(np.conj(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 # a = E_minus_sp_0_m[0]*1e9 # print("\n\nE_minus:", a ) # b = H_m[0]*1e9 # print("H_m:",b) # print("\ncross:",np.cross(E_minus_sp_0_m*1e9,H_m*1e9)[0]) # print("cross:",np.array(cross(a,b))) tmp_d = np.cross(E_minus_sp_0_d,H_d) - np.cross(E_d, H_minus_sp_0_d) #tmp_d = np.cross(E_minus_sp_0_d,np.conj(H_d)) + np.cross(np.conj(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])), # " 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" def main (monitor_index): #dirname="bigourdan-Au-sub-dipole-W.fsp.1D.monitor_1.results" #dirname="bigourdan-Au-sub-dipole-W-2mon.fsp.1D.monitor_2.results" #data = read_data(dirname) #filename = 'bigourdan-Au-sub-dipole-W-2mon.fsp.1D.mat' filename="bg-Au-sub-dipole-Au.fsp.1D.mat" data2 = read_data_mat(filename) r,data = data2[monitor_index] if verbose > 5: print("r =",r) WLs = [] A2 = [] for wl in data: WLs.append(wl) A2.append(analyze(data[wl],wl)) WLs1 = np.array(WLs) A21 = np.array(A2) #return # #dirname="bigourdan-Au-sub-Cyl-dipole-W.fsp.1D.monitor_1.results" # dirname="bigourdan-Au-sub-Cyl-dipole-W-2mon.fsp.1D.monitor_2.results" # data = read_data(dirname) #filename = 'bg-Au-sub-Au-dipole-Au.fsp.1D.mat' filename = 'bg-Au-sub-Si-dipole-Au.fsp.1D.mat' data2 = read_data_mat(filename) r,data = data2[monitor_index] 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, A22/A21, color="black",label="eff.") # plt.plot(WLs1*1e9, A21, color="black",label="x 275, no ant.") # plt.plot(WLs2*1e9, A22, color="red", label="with antena") plt.legend() plt.xlabel(r'$\lambda$, nm') # if monitor_index == 0: # plt.ylim(0,2) # else: # plt.ylim(0,0.32) #plt.ylim(0,2) plt.ylabel(r'$^{|A^{ant}_{sp}|^2}/_{|A^0_{sp}|^2}$',labelpad=1,size=14) plt.title("r = %4.1f nm"%(r*1e9)) plt.savefig(filename+str(monitor_index)+"_A2."+file_ext) plt.clf() plt.close() main(0) main(1)