Improving quality of imaging through microcylinder

[Jontza]

I am imaging 5 metal rods through a microcylinder. The rods are placed under the microcylinder, and the rods are on top of a metal plate. Plane wave propagates through the microcylinder from above, and the light reflected back from the rods is gathered at all the boundaries of the domain and the near-to-far-field transformation is performed on it.

Below there is images of the simulation domain, and the image that can be formed from the simulation.
In the image it can be seen that there is some kind of interference in the vertical direction.

Is there something that could be done to reduce this interference and to obtain clearer image?

Here is the code I am using:

  import meep as mp
  import numpy as np
  import os
  import math
  import matplotlib.pyplot as plt
  import argparse
  from PTP import PTP2D
  from scipy.interpolate import RectBivariateSpline
  
  print(mp.__version__)
  
  time_steps = 600
  
  radius = 5.5 #radius of the sphere or cylinder
  
  n_rod = 5   #number of rods
  r_rod = 0.1   #radius of the rods
  pitch = 0.1+2*r_rod   #pitch of the rods
  depth = 2*r_rod   #depth of the grooves
  angle = 0
  h_bot = r_rod*2
  
  gap_z = 0.05   #size of the gap between the rods and the cylinder/sphere
  
  index = 1.47   #refractive index of the cylinder/sphere
  wl = 0.51   #center wavelength
  wl_min = 0.4   #minimum wavelength
  wl_max = 0.7   #maximum wavelength
  nfreq = 5   #number of wavelengths sampled
  
  src_cmpt = 1   #[1->Ex / 2->Ey / 3->Ez]
  beam_x0 = mp.Vector3(0,0,0)    # beam focus (relative to source center)
  rot_angle = 5  # CCW rotation angle about y axis (0: +y axis)
  beam_kdir = mp.Vector3(0,0,-1).rotate(mp.Vector3(0,1,0),math.radians(rot_angle))  # beam propagation direction
  beam_w0 = mp.inf  # beam waist radius
  beam_E0 = mp.Vector3(1,0,0)
  
  nres = 80   #resolution in the computational cell, points/um
  
  z_proj = 20   #location of the far-field projection plane
  sxy_proj = 30    #size of the far-field projection plane, default: 40
  fres = 5   #resolution in the far field, default: 5
  
  z_bp = -5 #location of the backward projection plane, now the zf2 below
  zf2 = np.arange(-12.7, -12.8, -0.1)
  
  res_upsample = 501 #number of how many points the far-field data will be upsampled to, default: 2001
  
  def field_BivariateSpline(field, xw, yw, xnum, ynum):
      # Use RectBivariateSpline to interpolate data on a uniform grid
      x = np.linspace(-xw/2, xw/2, field.shape[0])
      y = np.linspace(-yw/2, yw/2, field.shape[1])
      x_new = np.linspace(-xw/2, xw/2, xnum)
      y_new = np.linspace(-yw/2, yw/2, ynum)
      fieldr_intp = RectBivariateSpline(y, x, field.real)
      fieldi_intp = RectBivariateSpline(y, x, field.imag)
      fieldr_new = fieldr_intp(y_new, x_new)
      fieldi_new = fieldi_intp(y_new, x_new)
      field_new = fieldr_new + 1j * fieldi_new
      return field_new

  
  def main(args):
      
  
      resolution = args.near_res     # pixels/um
  
      lambdas = args.wl   # wavelength
      fcen = 1/lambdas    # center frequency
      f_max = 1/wl_min    # maximum frequency
      f_min = 1/wl_max    # minimum frequency
      df = f_max-f_min           # bandwidth for pulsed source, default: 1
  
      gap = args.gap      # distance between source and lens (s')
      ms_r = args.R       # radius of microsphere
  
      nLens = args.nLens  # refractive index of microsphere
      ms = mp.Medium(index=nLens)
      ms_pec = mp.metal
  
      ## define the simulation region
      dpml = wl_max/2				# thickness of PML, default: 0.5
      pad_ff = wl_max			# pad between the microsphere and the NFTFF boundaries, default: 0.5
      pad_pml = wl_max/4				# pad between the NFTFF boundaries and the PML boundaries
      
      if ms_r != 0:
          sxy = ms_r*2+2*pad_ff		# Width of X and Y
  
      else:
          sxy = xsrc*2+pad_ff*2
      
      if ms_r != 0:
          sz = ms_r*2+gap+2*pad_ff+r_rod*2+h_bot	# Width of Z
      else:
          sz = gap+pad_ff*2 #sxy #
  
      sy = sxy
      sx = sxy
      
      pixel_s = 1/resolution
      
      print(pixel_s)
      
      sx_ff = sxy+2*pad_pml
      sy_ff = sy+2*pad_pml
      sz_ff = sz+2*pad_pml
      
      z_proj = args.z_proj                # The z coordinate of the projection plane
      x_width = args.x_width              # The x width of the projection plane
      y_width = args.y_width              # The y width of the projection plane
      far_resolution = args.far_res       # pixels/um
      
      cell_size = mp.Vector3(sx_ff+2*dpml, sy_ff+2*dpml, sz_ff+2*dpml) #sy+2*pad_pml+2*dpml
      cell_shift = pad_ff/2 #(ms_r*2+gap)/2
      geometry_center = mp.Vector3(z=cell_shift)
      boundary_layers = [mp.PML(dpml)]# [mp.PML(dpml, direction=mp.X),mp.PML(dpml, direction=mp.Z),mp.Absorber(dpml*2, direction=mp.Y)]
  
      src_z = 0.5*sz_ff+cell_shift#-(gap+ms_r)			# z coordinate of the source
      geometry = []
      
      bar_loc = np.arange(-pitch*(np.ceil(n_rod/2)-1),pitch*(np.ceil(n_rod/2)),pitch)
      bar_z = -ms_r-gap-r_rod
      
  
      bottom_z = -ms_r-gap-h_bot/2-depth
  
      rod_ax=np.tan(np.deg2rad(angle))
  
      ## define the microsphere and the rods
      
      ###   Bottom of the blu-ray   ###
      
      geometry.append(
              mp.Cylinder(center=mp.Vector3(x=0,y=0,z=bottom_z),
                    height=h_bot, 
                    radius=0.5*sx,
                    axis=mp.Vector3(0,0,1),
                    material=ms_pec)
          )
      
      ###   Rods   ###
          
          
      for i in bar_loc:
          geometry.append(
              mp.Cylinder(center=mp.Vector3(x=i,y=0,z=bar_z),
                    height=sy, 
                    radius=r_rod,
                    axis=mp.Vector3(rod_ax,1,0), # mp.Vector3(0.267949,1,0),
                    material=ms_pec)
    	  )
      
      ###   Cylinder on the rods   ###
      
      if ms_r != 0:
          geometry.append(
                  mp.Cylinder(center=mp.Vector3(x=0,y=0,z=0),
            	        height=sy, 
                        radius=ms_r,
                        axis=mp.Vector3(0,1,0),
                	    material=ms)
          )
          geom_suffix = 'cylinder_y'
      
      output_dir = 'micro{}-out'.format(geom_suffix)
      
      ## define the orientation of dipole source and symmetries
      symmetries = []
      if args.src_cmpt == 1:
          src_cmpt = mp.Ex
          suffix = 'x'
          output_func = mp.output_efield_x
          symmetries =  [mp.Rotate2(mp.Z, phase=-1)] # [mp.Mirror(mp.X, phase=-1),mp.Mirror(mp.Y, phase=1)] #
      elif args.src_cmpt == 2:
          src_cmpt = mp.Ey
          suffix = 'y'
          output_func = mp.output_efield_y
          symmetries = [mp.Rotate2(mp.Z)]
      elif args.src_cmpt == 3:
          src_cmpt = mp.Ez
          suffix = 'z'
          output_func = mp.output_efield_z
          symmetries = [mp.Rotate2(mp.Z)]
      
      sources = [
	  mp.GaussianBeamSource(src=mp.GaussianSource(fcen, fwidth=df, is_integrated=True),
                  center=mp.Vector3(0, 0, src_z),
                  size=mp.Vector3(x=sxy+2*pad_pml+2*dpml,y=sy+2*pad_pml+2*dpml),#mp.Vector3(x=mp.inf,y=mp.inf),#
                  beam_x0=beam_x0,
                  beam_kdir=beam_kdir,
                  beam_w0=beam_w0,
                  beam_E0=beam_E0
                )
      ]
  
      ## construct the FDTD simulation
      
      sim = mp.Simulation(cell_size=cell_size,
                          geometry_center = geometry_center,
                          resolution=resolution,
                          geometry=geometry,
                          sources=sources,
                          dimensions=3,
                          boundary_layers=boundary_layers,
                          symmetries=symmetries,
                          loop_tile_base_db = 10000,
                          loop_tile_base_eh = 10000,
                          force_complex_fields=True,
                          split_chunks_evenly=True
                         )
      
      mp.verbosity(2)
      
      dpx = 4*pixel_s    
  
      ## Define the field monitors for NTFF
      ## The Near2FarRegion should not lie in the PML regions
      
      nearfield = sim.add_near2far(fcen, df, nfreq,
                                   mp.Near2FarRegion(center=mp.Vector3(z=0.5*sz_ff)+geometry_center,
                                                     size=mp.Vector3(x=sx_ff, y=sy_ff),
                                                     weight=+1),
                                   mp.Near2FarRegion(center=mp.Vector3(z=-0.5*sz_ff)+geometry_center,
                                                     size=mp.Vector3(x=sx_ff, y=sy_ff),
                                                     weight=-1),
                                   mp.Near2FarRegion(center=mp.Vector3(x=0.5*sx_ff)+geometry_center,
                                                     size=mp.Vector3(z=sz_ff, y=sy_ff),
                                                     weight=+1),
                                   mp.Near2FarRegion(center=mp.Vector3(x=-0.5*sx_ff)+geometry_center,
                                                     size=mp.Vector3(z=sz_ff, y=sy_ff),
                                                     weight=-1),
                                   mp.Near2FarRegion(center=mp.Vector3(y=0.5*sy_ff)+geometry_center,
                                                     size=mp.Vector3(z=sz_ff, x=sy_ff),
                                                     weight=+1),
                                   mp.Near2FarRegion(center=mp.Vector3(y=-0.5*sy_ff)+geometry_center,
                                                     size=mp.Vector3(z=sz_ff, x=sx_ff),
                                                     weight=-1)
                                   )
  
      
      ## Define DFT monitor to record the fields in the simulation region
      
      vol = mp.Volume(center=geometry_center, size=mp.Vector3(sx_ff, sy_ff, sz)) #default: (sxy, 0, sz)
      vol_xz = mp.Volume(center=geometry_center, size=mp.Vector3(sx_ff, 0, sz_ff)) 
      vol_yz = mp.Volume(center=geometry_center, size=mp.Vector3(0, sy_ff, sz_ff)) 
      vol_xy = mp.Volume(center=mp.Vector3(z=bar_z-r_rod), size=mp.Vector3(sx_ff, sy_ff, 0)) 
      
      dfts = sim.add_dft_fields([src_cmpt], fcen, df, nfreq, where=vol) #default: ([src_cmpt], fcen, fcen, 1, where=vol)
      dfts_xz = sim.add_dft_fields([src_cmpt], fcen, df, nfreq, where=vol_xz) #default: ([src_cmpt], fcen, fcen, 1, where=vol)
      dfts_yz = sim.add_dft_fields([src_cmpt], fcen, df, nfreq, where=vol_yz) #default: ([src_cmpt], fcen, fcen, 1, where=vol)
      dfts_xy = sim.add_dft_fields([src_cmpt], fcen, df, nfreq, where=vol_xy) #default: ([src_cmpt], fcen, fcen, 1, where=vol)
  
  
      ## Output
      
      sub_dir = 'ms={}um'.format(ms_r)
      directory = '{}/{}'.format(output_dir, sub_dir)
      file_suffix = '{}_r={:.1f}_n_{:.2f}_t={:.1f}_J{}_wl={:.3f}_{}_rods_angle_{:.0f}degree_nr_{:.0f}_sz_{:.2f}_fr_{:.0f}_nfreq_{:.0f}_z_proj_{:.0f}_sxy_proj_{:.0f}'.format(geom_suffix, ms_r, nLens, time_steps, suffix, lambdas, n_rod, angle, resolution, sz, far_resolution, nfreq, z_proj, sxy_proj)
      if not os.path.exists(directory):
          os.makedirs(directory)
      sim.use_output_directory(directory)
  
  
      ## Run simulation
      ## If using until_after_sources, the presence of WGM will result in very slow convergence.
      
      sim.run(mp.in_volume(mp.Volume(center=geometry_center, size=mp.Vector3(sxy, sy, sz))),
                           until=time_steps)
      
      ## Visualize the structure in the simulation region to check if the setup is correct
      
      x, y, z, w = sim.get_array_metadata(vol_xz)
      eps_data = sim.get_array(vol= vol_xz, component=mp.Dielectric)
      e_data = np.absolute(sim.get_dft_array(dfts_xz, src_cmpt, 0))
      extent_xz = [min(x), max(x), min(z), max(z)]
      fig = plt.figure()
      plt.imshow(eps_data.transpose()[::-1], interpolation='spline36', cmap='binary', extent=extent_xz)
      plt.imshow(e_data.transpose()[::-1]**2, interpolation='spline36', cmap='jet', alpha=0.9, extent=extent_xz)
      title = '$|E{}|^2$'.format(suffix)
      plt.xlabel('x (um)')
      plt.ylabel('z (um)')
      plt.title(title, fontsize=16)
      fig.savefig('{}/structure_xz_{}.png'.format(directory,file_suffix))
      plt.close(fig)
      
      x, y, z, w = sim.get_array_metadata(vol_yz)
      eps_data = sim.get_array(vol= vol_yz, component=mp.Dielectric)
      e_data = np.absolute(sim.get_dft_array(dfts_yz, src_cmpt, nfreq-1))
      extent_yz = [min(y), max(y), min(z), max(z)]
      fig = plt.figure()
      plt.imshow(eps_data.transpose()[::-1], interpolation='spline36', cmap='binary', extent=extent_yz)
      plt.imshow(e_data.transpose()[::-1]**2, interpolation='spline36', cmap='jet', alpha=0.8, extent=extent_yz)
      title = '$|E{}|^2$'.format(suffix)
      plt.xlabel('y (um)')
      plt.ylabel('z (um)')
      plt.title(title, fontsize=16)
      fig.savefig('{}/structure_yz_{}.png'.format(directory,file_suffix))
      plt.close(fig)
      
      x, y, z, w = sim.get_array_metadata(vol_xy)
      eps_data = sim.get_array(vol= vol_xy, component=mp.Dielectric)
      e_data = np.absolute(sim.get_dft_array(dfts_xy, src_cmpt, nfreq-1))
      extent_xy = [min(x), max(x), min(y), max(y)]
      fig = plt.figure()
      plt.imshow(eps_data.transpose()[::-1], interpolation='spline36', cmap='binary', extent=extent_xy)
      plt.imshow(e_data.transpose()[::-1]**2, interpolation='spline36', cmap='jet', alpha=0.8, extent=extent_xy)
      title = '$|E{}|^2$'.format(suffix)
      plt.xlabel('x (um)')
      plt.ylabel('y (um)')
      plt.title(title, fontsize=16)
      fig.savefig('{}/structure_xy_{}.png'.format(directory,file_suffix))
      plt.close(fig)
      
      ## Near-to-far-field transformation 1
      print('\n')
      print('#############################')
      print('Extract Farfield Information 1')
      print('#############################')
      print('\n')
      
      dx = 1/far_resolution
      dy = 1/far_resolution
      nx = x_width*far_resolution
      ny = y_width*far_resolution
      xs = np.linspace(-x_width/2, x_width/2, nx, endpoint=True)
      ys = np.linspace(-y_width/2, y_width/2, ny, endpoint=True)
      
      #Projection 1
  
      counter = 0
      E_ff = []
      H_ff = []
      xy_pos = []
      xy_mesh = np.meshgrid(xs,ys, indexing='ij')
      #S_ff = []
      
      x_pos=xy_mesh[0].ravel()
      y_pos=xy_mesh[1].ravel()
      
      xy_pos=[list(a) for a in zip(x_pos, y_pos)]     
      
      ff = sim.get_farfields(nearfield,
                              far_resolution,
                              center = mp.Vector3(0, 0, z_proj),
                              size = mp.Vector3(x_width, y_width, 0)
                             )
      
      
      print('#############################')
      print('[Finished] Extract Farfield Information')
      print('#############################')
      print('\n')
      
      E_ff=[]
      Ex = ff['Ex'].ravel()
      Ey = ff['Ey'].ravel()
      Ez = ff['Ez'].ravel()
      
      Hx = ff['Hx'].ravel()
      Hy = ff['Hy'].ravel()
      Hz = ff['Hz'].ravel()
      
      E_ff=[list(a) for a in zip(Ex, Ey, Ez)]
      H_ff=[list(a) for a in zip(Hx, Hy, Hz)]
      
      print('#############################')
      print('Get Wavelengths')
      print('#############################')
      print('\n')
      
      frequencies = mp.get_near2far_freqs(nearfield)
      wavelengths=[]
      for f in frequencies:
          wavelengths.append(f**(-1))
          print('Wavelength: {}'.format(f**(-1)))
          
      print('#############################')
      print('[Finished] Get Wavelengths')
      print('#############################')
      print('\n')
    
      E_ff_new = []
      H_ff_new = []
      xy_pos_new = []  
      
      print('#############################')
      print('Reshape Electric and Magnetic Fields')
      print('#############################')
      print('\n')
      
      E_ff = np.asarray(E_ff).reshape(nx, ny, 3*nfreq)
      H_ff = np.asarray(H_ff).reshape(nx, ny, 3*nfreq)
      xy_pos = np.asarray(xy_pos).reshape(nx, ny, 2)
      
      print('#############################')
      print('[Finished] Reshape Electric and Magnetic Fields')
      print('#############################')
      print('\n')
      
      ## Save npz file
      
      print('#############################')
      print('Save .npz file')
      print('#############################')
      print('\n')
  
      np.savez(file='{}/farfield_{}'.format(directory, file_suffix), E=E_ff, H=H_ff, xy_pos=xy_pos, z_proj=z_proj, lambdas=wavelengths)
  
      print('#############################')
      print('[Finished] Save .npz file')
      print('#############################')
      print('\n')
      
      ## Up-sampling
      xnum = res_upsample #default: 2001
      ynum = res_upsample #default: 2001
      x_new = np.linspace(-x_width/2, x_width/2, xnum, endpoint=True)
      y_new = np.linspace(-y_width/2, y_width/2, ynum, endpoint=True)
      dx_new = x_width / (xnum-1)
      dy_new = y_width / (ynum-1)
      extent = [-x_width/2, x_width/2, -y_width/2, y_width/2]
  
      print('#############################')
      print('Backward propagation')
      print('#############################')
      print('\n')
  
  
      # BivariateSpline interpolation
  
      e1_img_bp = np.zeros((xnum,ynum))
      e2_img_bp = np.zeros((xnum,ynum))
      
      #zf = args.zf
      for zf in zf2:
          for i in range(nfreq):
          
              #print('[Progress] {}/{}'.format(i, nfreq))
              Ex_new = field_BivariateSpline(E_ff[:, :, 0+i*3], x_width, y_width, xnum, ynum)
              Ey_new = field_BivariateSpline(E_ff[:, :, 1+i*3], x_width, y_width, xnum, ynum)
              Ez_new = field_BivariateSpline(E_ff[:, :, 2+i*3], x_width, y_width, xnum, ynum)
              
              print(wavelengths[i])
              ## Backward propagation
  
              ex_PTP = PTP2D(Ex_new, wavelengths[i], x_width, y_width)
              ey_PTP = PTP2D(Ey_new, wavelengths[i], x_width, y_width)
              ez_PTP = PTP2D(Ez_new, wavelengths[i], x_width, y_width)
  
              ex_PTP.lp_filter()
              ey_PTP.lp_filter()
              ez_PTP.lp_filter()
  
              ex_img_bp = ex_PTP.aspw(zf=zf-z_proj)
              ey_img_bp = ey_PTP.aspw(zf=zf-z_proj)
              ez_img_bp = ez_PTP.aspw(zf=zf-z_proj)
              e1_img_bp = np.abs(ex_img_bp) ** 2 + np.abs(ey_img_bp) ** 2 + np.abs(ez_img_bp) ** 2
              e2_img_bp = e2_img_bp + np.abs(ex_img_bp) ** 2 + np.abs(ey_img_bp) ** 2 + np.abs(ez_img_bp) ** 2
  
              crop_width = sxy
              crop_indx = np.where(np.abs(x_new) <= crop_width / 2)
              crop_indy = np.where(np.abs(y_new) <= crop_width / 2)
              fig = plt.figure()
              plt.imshow(np.transpose(e1_img_bp[np.min(crop_indx):(np.max(crop_indx)+1), np.min(crop_indy):(np.max(crop_indy)+1)]),
                         cmap='hot', extent=[-crop_width/2, crop_width/2, -crop_width/2, crop_width/2], interpolation='spline36')
              plt.title('\u03BB = {:.3f} um / z = {:.1f} um'.format(wavelengths[i], zf), fontsize=16)
              plt.colorbar()
              plt.xticks(fontsize=12)
              plt.yticks(fontsize=12)
              plt.xlabel('x (um)', fontsize=14)
              plt.ylabel('y (um)', fontsize=14)
              plt.savefig('{}/image_bp_wl_{}um_{}_{}.png'.format(directory, wavelengths[i], zf, file_suffix))
              plt.close(fig)
      
  
  
          ## Visualize the virtual image
          crop_width = sxy
          crop_indx = np.where(np.abs(x_new) <= crop_width / 2)
          crop_indy = np.where(np.abs(y_new) <= crop_width / 2)
          fig = plt.figure()
          plt.imshow(np.transpose(e2_img_bp[np.min(crop_indx):(np.max(crop_indx)+1), np.min(crop_indy):(np.max(crop_indy)+1)]),
                     cmap='hot', extent=[-crop_width/2, crop_width/2, -crop_width/2, crop_width/2], interpolation='spline36')
          plt.title('white light, center \u03BB = {:.3f} um / z = {:.1f} um'.format(lambdas, zf), fontsize=16)
          plt.colorbar()
          plt.xticks(fontsize=12)
          plt.yticks(fontsize=12)
          plt.xlabel('x (um)', fontsize=14)
          plt.ylabel('y (um)', fontsize=14)
          plt.savefig('{}/image_bp_tot_zf_{}um_{}.png'.format(directory, zf, file_suffix))
          plt.close(fig)
      
  
  
      print('#############################')
      print('[Finished] Backward propagation')
      print('#############################')
      print('\n')
    	  
  
      # BivariateSpline interpolation
      
  
      Ex_new = []
      Ey_new = []
      Ez_new = []
      
      Hx_new = []
      Hy_new = []
      Hz_new = []
      
      E_ff_img = np.zeros((xnum,ynum))
      
      for i in range(nfreq):
          print('[Progress] {}/{}'.format(i, nfreq))
          Ex_new.append(field_BivariateSpline(E_ff[:, :, 0+i*3], x_width, y_width, xnum, ynum))
          Ey_new.append(field_BivariateSpline(E_ff[:, :, 1+i*3], x_width, y_width, xnum, ynum))
          Ez_new.append(field_BivariateSpline(E_ff[:, :, 2+i*3], x_width, y_width, xnum, ynum))
      
          Hx_new.append(field_BivariateSpline(H_ff[:, :, 0+i*3], x_width, y_width, xnum, ynum))
          Hy_new.append(field_BivariateSpline(H_ff[:, :, 1+i*3], x_width, y_width, xnum, ynum))
          Hz_new.append(field_BivariateSpline(H_ff[:, :, 2+i*3], x_width, y_width, xnum, ynum))
          
          E_ff_img = E_ff_img + np.abs(field_BivariateSpline(E_ff[:, :, 0+i*3], x_width, y_width, xnum, ynum)) ** 2 + np.abs(field_BivariateSpline(E_ff[:, :, 1+i*3], x_width, y_width, xnum, ynum)) ** 2 + np.abs(field_BivariateSpline(E_ff[:, :, 2+i*3], x_width, y_width, xnum, ynum)) ** 2
      
          counter=0
      
          '''
          for x in x_new:
              for y in y_new:
                  counter = counter + 1
                  print('[Progress] {}/{}'.format(counter, xnum * ynum))
                  xy_pos_new.append([x, y])
          '''
          xy_mesh = np.meshgrid(x_new,x_new, indexing='ij')
          #S_ff = []
  
          x_pos=xy_mesh[0].ravel()
          y_pos=xy_mesh[1].ravel()
  
          xy_pos_new=[list(a) for a in zip(x_pos, y_pos)]
      
      for i,j,k in zip(Ex_new,Ey_new,Ez_new):
          E_ff_new.append((i,j,k))
          
      for i,j,k in zip(Hx_new,Hy_new,Hz_new):
          H_ff_new.append((i,j,k))
      
      
      ## Visualize the field at the far field plane
      
      E_ff_img2 = np.abs(field_BivariateSpline(E_ff[:, :, 0], x_width, y_width, xnum, ynum)) ** 2 + np.abs(field_BivariateSpline(E_ff[:, :, 1], x_width, y_width, xnum, ynum)) ** 2 + np.abs(field_BivariateSpline(E_ff[:, :, 2], x_width, y_width, xnum, ynum)) ** 2
      
      crop_width = sxy_proj
      crop_indx = np.where(np.abs(x_new) <= crop_width / 2)
      crop_indy = np.where(np.abs(y_new) <= crop_width / 2)
      fig = plt.figure()
      plt.imshow(np.transpose(E_ff_img[np.min(crop_indx):(np.max(crop_indx)+1), np.min(crop_indy):(np.max(crop_indy)+1)]),
                 cmap='hot', extent=[-crop_width/2, crop_width/2, -crop_width/2, crop_width/2], interpolation='spline36')
      plt.title('\u03BB = {:.3f} um / z = {:.1f} um'.format(lambdas, zf), fontsize=16)
      plt.colorbar()
      plt.xticks(fontsize=12)
      plt.yticks(fontsize=12)
      plt.xlabel('x (um)', fontsize=14)
      plt.ylabel('y (um)', fontsize=14)
      plt.savefig('{}/image_ff_{}.png'.format(directory, file_suffix))
      plt.close(fig)
      
      crop_width = sxy_proj
      crop_indx = np.where(np.abs(x_new) <= crop_width / 2)
      crop_indy = np.where(np.abs(y_new) <= crop_width / 2)
      fig = plt.figure()
      plt.imshow(np.transpose(E_ff_img2[np.min(crop_indx):(np.max(crop_indx)+1), np.min(crop_indy):(np.max(crop_indy)+1)]),
                 cmap='hot', extent=[-crop_width/2, crop_width/2, -crop_width/2, crop_width/2], interpolation='spline36')
      plt.title('\u03BB = {:.3f} um / z = {:.1f} um'.format(lambdas, zf), fontsize=16)
      plt.colorbar()
      plt.xticks(fontsize=12)
      plt.yticks(fontsize=12)
      plt.xlabel('x (um)', fontsize=14)
      plt.ylabel('y (um)', fontsize=14)
      plt.savefig('{}/image_ff2_{}.png'.format(directory, file_suffix))
      plt.close(fig)
  
  if __name__ == '__main__':
      parser = argparse.ArgumentParser(prog='PROG')
      
      parser.add_argument('-R', type=float, default=radius, dest='R',
                          help='radius of microsphere (default: 2.0 um)')
      parser.add_argument('-gap', type=float, default=gap_z, dest='gap',
                          help='gap between source and hemisphere (default: 0.05 um)')
      parser.add_argument('-nLens', type=float, default=index, dest='nLens',
                          help='refractive index of hemisphere (default: 1.46)')
      parser.add_argument('-wl', type=float, default=wl, dest='wl',
                          help='working wavelength (default: 0.405 um)')
      parser.add_argument('--src_cmpt', type=int, default=src_cmpt, dest='src_cmpt',
                          help='dipole orientation: [1->Ex / 2->Ey / 3->Ez] (default: 1)')
      parser.add_argument('-near_res', type=float, default=nres, dest='near_res',
                          help='the resolution of mesh grid (default: 50 pts/um)')
      parser.add_argument('-z_proj', type=float, default=z_proj, dest='z_proj',
                          help='the position of the projection plane (default: 7 um)')
      parser.add_argument('-x_width', type=float, default=sxy_proj, dest='x_width',
                          help='the x width of the projection plane (default: 40 um)')
      parser.add_argument('-y_width', type=float, default=sxy_proj, dest='y_width',
                          help='the y width of the projection plane (default: 40 um)')
      parser.add_argument('-far_res', type=float, default=fres, dest='far_res',
                          help='the sampling rate on the projection plane (default: 5 pts/um)')
      parser.add_argument('-zf', type=float, default=z_bp, dest='zf',
                          help='the final plane of backward propagation (default: -4.5 um)')
      
      args, unknown = parser.parse_known_args([])
      print(args)
      main(args)