|  | @@ -4,7 +4,7 @@ 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
 | 
	
	
		
			
				|  | @@ -53,13 +53,51 @@ def get_WLs_idx(WLs, data):
 | 
	
		
			
				|  |  |      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,:]
 | 
	
	
		
			
				|  | @@ -70,45 +108,69 @@ def analyze(data, dist, z_vec, wl_idx):
 | 
	
		
			
				|  |  |      kappa1= np.sqrt(k_spp**2 - eps1*k_0**2)
 | 
	
		
			
				|  |  |      kappa2= np.sqrt(k_spp**2 - eps2*k_0**2)
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  | -    # TODO def check_field_match(): 
 | 
	
		
			
				|  |  | -    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
 | 
	
		
			
				|  |  | -        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)*eps2[wl_idx]
 | 
	
		
			
				|  |  | -        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_0), np.imag(E2_0)),
 | 
	
		
			
				|  |  | -              " from E2 (%5.4g %+5.4gj)"%(np.real(E2[i]), np.imag(E2[i])))
 | 
	
		
			
				|  |  | -    # H1_0 = H1/np.exp(-kappa1* 
 | 
	
		
			
				|  |  | -    # print(H1[0], H2[0],H1[0]- H2[0])
 | 
	
		
			
				|  |  | -    # pl_data = (np.absolute(data_gold[:,2,wl_idx]*np.sqrt(dist)))
 | 
	
		
			
				|  |  | -    # plt.semilogy(z_vec, pl_data,marker="o")
 | 
	
		
			
				|  |  | +    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'Z shift, nm')
 | 
	
		
			
				|  |  | -    # plt.ylabel(r'$Abs(E_x) \sqrt{R}$',labelpad=-5)
 | 
	
		
			
				|  |  | -    # # plt.title(' r = '+str(core_r))
 | 
	
		
			
				|  |  | -    # plt.savefig(dirname+"_z."+file_ext)
 | 
	
		
			
				|  |  | -    # plt.clf()
 | 
	
		
			
				|  |  | -    # plt.close()
 | 
	
		
			
				|  |  | +    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="png"
 | 
	
		
			
				|  |  | -dirname="template-dipole-on-sphere-on-surf-z.fsp.results"
 | 
	
		
			
				|  |  | +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]
 | 
	
		
			
				|  |  | +    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]
 | 
	
		
			
				|  |  |  
 | 
	
	
		
			
				|  | @@ -116,7 +178,7 @@ def main ():
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |      #WLs=[300,350,400,450,600,700,800]
 | 
	
		
			
				|  |  |      #WLs=[600,700, 800, 450]
 | 
	
		
			
				|  |  | -    WLs=[600]#, 450]
 | 
	
		
			
				|  |  | +    WLs=[800]#, 450]
 | 
	
		
			
				|  |  |      WLs_idx = get_WLs_idx(WLs, data)
 | 
	
		
			
				|  |  |  
 | 
	
		
			
				|  |  |  
 |