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

#!/usr/bin/env python3
# -*- coding: UTF-8 -*-
import numpy as np
import matplotlib.pyplot as plt
c = 299792458
pi = np.pi
verbose = 6
def read_data(dirname, distance, zshift):
    media = [1,2] # 1 - positive zshift, 2 - negative (need to add a minus sign for real shift).
    #min_mesh_step = 2.5 #nm
    data = []
    data.append([])
    for x in distance:
        data.append([])
        data[x].append([])
        for m in media:
            data[x].append([])
            for z in zshift:            
                monitor_name = "mon_x"+str(x)+"mkm_media"+str(m)+"_zshift"+z+"nm"
                data[x][m].append(
                    np.transpose(
                    np.genfromtxt(dirname+"/"+monitor_name+".txt", 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'))
                                                    }
                                      )
                        )
                    )
    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): 
    H1 = data_in_air[:,6,wl_idx]
    H2 = data_in_gold[:,6,wl_idx]
    E1 = data_in_air[:,4,wl_idx]
    E2 = data_in_gold[:,4,wl_idx]
    for i in range(len(z_vec)):
        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])))
            print("H0 gold (%5.4g %+5.4gj)"%(np.real(H2_0), np.imag(H2_0)),
                  " from H2 (%5.4g %+5.4gj)"%(np.real(H2[i]), np.imag(H2[i])))
            print("E0 air  (%5.4g %+5.4gj)"%(np.real(E1_0), np.imag(E1_0)),
                  " from E1 (%5.4g %+5.4gj)"%(np.real(E1[i]), np.imag(E1[i])))
            print("E0*eps2 (%5.4g %+5.4gj)"%(np.real(E2_0e), np.imag(E2_0e)),
                  " from E2 (%5.4g %+5.4gj)"%(np.real(E2[i]), np.imag(E2[i])))
            print("E0 gold (%5.4g %+5.4gj)"%(np.real(E2_0), np.imag(E2_0)))


def analyze(data, dist, z_vec, wl_idx):
    ''' dist in mkm!!!
    '''
    #data = [dist][mmedia][shift] "lambda, dip.power, Ex, Ey, Ez, Hx, Hy, Hz, n_Au"
    #                              0     , 1        , 2 , 3 , 4 , 5 , 6 , 7 , 8   "
    data_in_air = np.array(data[dist][1])
    data_in_gold = np.array(data[dist][2])
    lambd = data_in_air[0][0,:]
    omega = 2*pi*c/lambd
    dip_power = data_in_air[0][1,:]

    Ex = data_in_air[0,2,0]
    Ey = data_in_air[0,3,0]
    Ez = data_in_air[0,4,0]
    Hx = data_in_air[0,5,0]
    Hy = data_in_air[0,6,0]
    Hz = data_in_air[0,7,0]
    E = np.array([Ex,Ey,Ez])
    H = np.array([Hx,Hy,Hz])
    print("S from full field",np.real(np.cross(E,np.conj(H))))
    
    eps1 = complex(1)
    n_Au = data_in_air[0][8,:]
    eps2 = n_Au**2

    k_0 = omega/c #air
    k_spp = k_0*np.sqrt(eps1*eps2/(eps1+eps2))
    kappa1= np.sqrt(k_spp**2 - eps1*k_0**2)
    kappa2= np.sqrt(k_spp**2 - eps2*k_0**2)

    check_field_match(data_in_air, data_in_gold,wl_idx,z_vec,kappa1,kappa2,eps2)

    H1 = data_in_air[:,6]
    E1 = data_in_air[:,4]
    
    z = z_vec[0]*1e-9

    if verbose > 5: print("Using data from air monitor at z =",z)
    H1_0 = H1[0]/np.exp(-kappa1*z)
    E1_0 = E1[0]/np.exp(-kappa1*z)
    E2_0 = E1[0]/eps2
    if verbose > 5: 
        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])))
        print("E0 air  (%5.4g %+5.4gj)"%(np.real(E1_0[wl_idx]), np.imag(E1_0[wl_idx])),
              " from E1 (%5.4g %+5.4gj)"%(np.real(E1[0][wl_idx]), np.imag(E1[0][wl_idx])))
        print("E0 gold (%5.4g %+5.4gj)"%(np.real(E2_0[wl_idx]), np.imag(E2_0[wl_idx])), " from E1")

    R = dist*1e-6
    print("R =",R)
    #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,1.0)

    #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="template-dipole-on-sphere-on-surf-z.fsp.results"
#dirname="Au-JC-R100-Au-JC.fsp.results"
#dirname="Au-McPeak-R100-Si-Green.fsp.results"
#dirname="Au-McPeak-R100-Au-McPeak.fsp.results"
#dirname="Au-McPeak-R0.fsp.results"
#dirname="Au-McPeak-R100-Si-Green-1500.fsp.results"
#dirname="Au-McPeak-R100-Si-Green-1500-l.fsp.results"
dirname="Au-McPeak-R50-Si-Green-1500-l.fsp.results"
def main ():
    distance = [1,2,3,4,5,6,7,8,9,10] #mkm
    zshift = ["5","20","200","400","600"]
    z_vec = [int(val) for val in zshift]

    data = read_data(dirname, distance, zshift)

    #WLs=[300,350,400,450,600,700,800]
    #WLs=[600,700, 800, 450]
    WLs=[800]#, 450]
    WLs_idx = get_WLs_idx(WLs, data)


    dist = 10 #mkm
    wl_idx = WLs_idx[0]
    
    analyze(data, dist, z_vec, wl_idx)


    # legend = []
    # for shift in range(len(zshift)):
    #     for i in range(len(WLs)):
    #         pl_data = []
    #         idx = WLs_idx[i]
    #         legend.append(zshift[shift]+"@"+str(WLs[i])+" nm")
    #         for dist in distance:
    #             pl_data.append(np.absolute(data[dist][mmedia][shift][2,idx]*np.sqrt(dist)))
    #         print(len(pl_data))
    #         plt.semilogy(distance, pl_data,marker="o")
    # plt.legend(legend)
    # # #plt.xlabel(r'THz')
    # plt.xlabel(r'Monitor R, $\mu$m')
    # plt.ylabel(r'$Abs(E_x) \sqrt{R}$',labelpad=-5)
    # # plt.title(' r = '+str(core_r))
    # plt.savefig(dirname+"_WLs."+file_ext)
    # plt.clf()
    # plt.close()
main()