dipole-on-surf/plasmon-modal-efficiency.py

#!/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
verbose = 6
# r of monitor
#r = 146.513e-9
r = 800e-9
#debug = True
debug = False

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))
    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)
        ))
        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)
        ))
        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 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.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

    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

    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)
    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])),
    #           " 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 ():
    #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'
    data2 = read_data_mat(filename)
    r,data = data2[0]
    if verbose > 5:
        print("r =",r)
    WLs = []
    A2 = []
    for wl in data:
        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"
    # dirname="bigourdan-Au-sub-Cyl-dipole-W-2mon.fsp.1D.monitor_2.results"
    # data = read_data(dirname)
    filename = 'bigourdan-Au-sub-Cyl-dipole-W-2mon.fsp.1D.mat'
    data2 = read_data_mat(filename)
    r,data = data2[0]
    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, color="black",label="x 275, no ant.")
    plt.plot(WLs2*1e9, A22, color="red", label="with antena")
    plt.legend()
    plt.xlabel(r'$\lambda$, nm')
    plt.ylim(0,0.5)
    plt.ylabel(r'$|A_{sp}|^2$',labelpad=-1)
    #plt.title(dirname)
    plt.savefig(dirname+"_A2."+file_ext)
    plt.clf()
    plt.close()

main()