This example demonstrates how to compute the far-field profile at the focal length of a metasurface lens. The lens design, which is also part of the tutorial, is based on a supercell of binary-grating unit cells. For a review of the binary-grating geometry as well as a demonstration of computing its phasemap, see Tutorial/Mode Decomposition/Phase Map of a Subwavelength Binary Grating. The far-field calculation of the lens contains two separate components: (1) compute the phasemap of the unit cell as a function of a single geometric parameter, the duty cycle, while keeping its height and periodicity fixed (1.8 and 0.3 μm), and (2) form the supercell lens by tuning the local phase of each of a variable number of unit cells according to the quadratic formula for planar wavefront focusing. The design wavelength is 0.5 μm and the focal length is 0.2 mm. The input source is an $E_z$-polarized planewave at normal incidence.
The key to the script is the function grating
with three geometric input arguments (periodicity, height, and list of duty cycles) which performs the two main tasks: (1) for a unit cell, it computes the phase (as well as the transmittance) and then translates this value from the range of [-π,π] of Mode Decomposition to [-2π,0] in order to be consistent with the analytic formula for the local phase and (2) for a supercell, it computes the far-field intensity profile around the focal length of the lens.
import meep as mp
import numpy as np
import matplotlib.pyplot as plt
resolution = 50 # pixels/μm
dpml = 1.0 # PML thickness
dsub = 2.0 # substrate thickness
dpad = 2.0 # padding between grating and PML
lcen = 0.5 # center wavelength
fcen = 1/lcen # center frequency
df = 0.2*fcen # frequency width
focal_length = 200 # focal length of metalens
spot_length = 100 # far field line length
ff_res = 10 # far field resolution (points/μm)
k_point = mp.Vector3(0,0,0)
glass = mp.Medium(index=1.5)
pml_layers = [mp.PML(thickness=dpml,direction=mp.X)]
symmetries=[mp.Mirror(mp.Y)]
def grating(gp,gh,gdc_list):
sx = dpml+dsub+gh+dpad+dpml
src_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub)
mon_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad)
geometry = [mp.Block(material=glass,
size=mp.Vector3(dpml+dsub,mp.inf,mp.inf),
center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub)))]
num_cells = len(gdc_list)
if num_cells == 1:
sy = gp
cell_size = mp.Vector3(sx,sy)
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
component=mp.Ez,
center=src_pt,
size=mp.Vector3(y=sy))]
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
k_point=k_point,
default_material=glass,
sources=sources,
symmetries=symmetries)
flux_obj = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=mon_pt, size=mp.Vector3(y=sy)))
sim.run(until_after_sources=50)
input_flux = mp.get_fluxes(flux_obj)
sim.reset_meep()
geometry.append(mp.Block(material=glass, size=mp.Vector3(gh,gdc_list[0]*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh)))
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
geometry=geometry,
k_point=k_point,
sources=sources,
symmetries=symmetries)
flux_obj = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=mon_pt, size=mp.Vector3(y=sy)))
sim.run(until_after_sources=200)
freqs = mp.get_eigenmode_freqs(flux_obj)
res = sim.get_eigenmode_coefficients(flux_obj, [1], eig_parity=mp.ODD_Z+mp.EVEN_Y)
coeffs = res.alpha
mode_tran = abs(coeffs[0,0,0])**2/input_flux[0]
mode_phase = np.angle(coeffs[0,0,0])
if mode_phase > 0:
mode_phase -= 2*np.pi
return mode_tran, mode_phase
else:
sy = num_cells*gp
cell_size = mp.Vector3(sx,sy)
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
component=mp.Ez,
center=src_pt,
size=mp.Vector3(y=sy))]
for j in range(num_cells):
geometry.append(mp.Block(material=glass,
size=mp.Vector3(gh,gdc_list[j]*gp,mp.inf),
center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh,-0.5*sy+(j+0.5)*gp)))
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
geometry=geometry,
k_point=k_point,
sources=sources,
symmetries=symmetries)
n2f_obj = sim.add_near2far(fcen, 0, 1, mp.Near2FarRegion(center=mon_pt, size=mp.Vector3(y=sy)))
sim.run(until_after_sources=100)
return abs(sim.get_farfields(n2f_obj, ff_res, center=mp.Vector3(-0.5*sx+dpml+dsub+gh+focal_length), size=mp.Vector3(spot_length))['Ez'])**2
In the first of two parts of the calculation, a phasemap of the binary-grating unit cell is generated based on varying the duty cycle from 0.1 to 0.9.
gp = 0.3 # grating periodicity
gh = 1.8 # grating height
gdc = np.linspace(0.1,0.9,30) # grating duty cycle
mode_tran = np.empty((gdc.size))
mode_phase = np.empty((gdc.size))
for n in range(gdc.size):
mode_tran[n], mode_phase[n] = grating(gp,gh,[gdc[n]])
----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00247622 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00647688 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00136495 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.03,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0103281 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000896931 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.0044961 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00125909 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.0382759,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0100789 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00088501 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00477791 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.0013721 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.0465517,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.00952506 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00219703 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00522304 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00136685 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.0548276,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0103021 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00116801 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00504303 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00156713 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.0631034,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.010608 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00124407 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00682497 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000928879 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.0713793,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.012069 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 7 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000815868 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00467396 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000994921 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.0796552,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.00999308 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000801086 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00456905 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00162697 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.087931,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0103731 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.0009408 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00609112 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00112414 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.0962069,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0113389 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00133109 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.004987 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000950098 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.104483,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.00973797 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00119305 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00534606 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00095892 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.112759,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0165179 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 7 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000835896 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00442696 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.001477 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.121034,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0124168 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00136709 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00486684 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000990152 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.12931,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0103049 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.001513 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.005023 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000887156 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.137586,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0100081 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000860929 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.004879 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000782967 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.145862,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0097518 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 7 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000844002 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00591779 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00101209 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.154138,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0102739 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000947952 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00474691 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000870943 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.162414,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.00949287 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000854969 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00593901 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00108814 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.17069,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0101011 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000795126 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00502586 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00118899 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.178966,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0093739 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 9 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00110197 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00501704 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00141692 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.187241,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.013217 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00104403 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00458407 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000869989 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.195517,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.00961399 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00112414 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.0046351 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000927925 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.203793,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.009727 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00181007 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00469708 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00114608 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.212069,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0103791 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 9 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00141215 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00565505 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00124907 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.220345,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.011498 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00082016 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00524902 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000906944 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.228621,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0162721 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 7 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00110412 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00537992 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00147796 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.236897,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.010345 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 7 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00085187 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00499296 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000822783 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.245172,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0116842 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 9 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000720978 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00519085 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000974894 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.253448,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.012033 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 9 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000780106 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00697684 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00137496 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.261724,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0110409 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 7 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000945091 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 time for set_epsilon = 0.00621819 s ----------- run 0 finished at t = 75.0 (7500 timesteps) ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000971079 s Working in 2D dimensions. Computational cell is 7.8 x 0.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,0,0) size (1.8,0.27,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.0130391 s ----------- run 0 finished at t = 225.0 (22500 timesteps) MPB solved for omega_1(2,0,0) = 2 after 8 iters Dominant planewave for band 1: (2.000000,-0.000000,0.000000)
plt.figure(dpi=200)
plt.subplot(1,2,1)
plt.plot(gdc, mode_tran, 'bo-')
plt.xlim(gdc[0],gdc[-1])
plt.xticks([t for t in np.linspace(0.1,0.9,5)])
plt.xlabel("grating duty cycle")
plt.ylim(0.96,1.00)
plt.yticks([t for t in np.linspace(0.96,1.00,5)])
plt.title("transmittance")
plt.subplot(1,2,2)
plt.plot(gdc, mode_phase, 'rs-')
plt.grid(True)
plt.xlim(gdc[0],gdc[-1])
plt.xticks([t for t in np.linspace(0.1,0.9,5)])
plt.xlabel("grating duty cycle")
plt.ylim(-2*np.pi,0)
plt.yticks([t for t in np.linspace(-6,0,7)])
plt.title("phase (radians)")
plt.tight_layout(pad=0.5)
plt.show()
The phasemap is shown above. The left figure shows the transmittance which is nearly unity for all values of the duty cycle; the Fresnel transmittance is 0.96 for the glass-air interface. This is expected since the periodicity is subwavelength. The right figure shows the phase. There is a subregion in the middle of the plot spanning the duty-cycle range of roughly 0.16 to 0.65 in which the phase varies continuously over the full range of -2π to 0. This structural regime is used to design the supercell lens.
In the second part of the calculation, the far-field energy-density profile of three supercell lens designs, comprised of 201, 401, and 801 unit cells, are computed using the quadratic formula for the local phase. Initially, this involves fitting the unit-cell phase data to a finer duty-cycle grid in order to enhance the local-phase interpolation of the supercell. This is important since as the number of unit cells in the lens increases, the local phase via the duty cycle varies more gradually from unit cell to unit cell. However, if the duty cycle becomes too gradual (i.e., less than a tenth of the pixel dimensions), the resolution
may also need to be increased in order to improve the accuracy of subpixel smoothing.
gdc_new = np.linspace(0.16,0.65,500)
mode_phase_interp = np.interp(gdc_new, gdc, mode_phase)
print("phase-range:, {:.6f}".format(mode_phase_interp.max()-mode_phase_interp.min()))
phase_tol = 1e-2
num_cells = [100,200,400]
ff_nc = np.empty((spot_length*ff_res,len(num_cells)))
for k in range(len(num_cells)):
gdc_list = []
for j in range(-num_cells[k],num_cells[k]+1):
phase_local = 2*np.pi/lcen * (focal_length-((j*gp)**2 + focal_length**2)**0.5) # local phase at the center of the j'th unit cell
phase_mod = phase_local % (-2*np.pi) # restrict phase to [-2*pi,0]
if phase_mod > mode_phase_interp.max():
phase_mod = mode_phase_interp.max()
if phase_mod < mode_phase_interp.min():
phase_mod = mode_phase_interp.min()
idx = np.transpose(np.nonzero(np.logical_and(mode_phase_interp > phase_mod-phase_tol, mode_phase_interp < phase_mod+phase_tol))).ravel()
gdc_list.append(gdc_new[idx[0]])
ff_nc[:,k] = grating(gp,gh,gdc_list)
phase-range:, 6.086174 ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00169206 s Working in 2D dimensions. Computational cell is 7.8 x 60.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-30,0) size (1.8,0.109864,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-29.7,0) size (1.8,0.123415,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-29.4,0) size (1.8,0.136671,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-29.1,0) size (1.8,0.152579,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-28.8,0) size (1.8,0.169076,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-28.5,0) size (1.8,0.186752,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-28.2,0) size (1.8,0.0500621,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-27.9,0) size (1.8,0.059489,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-27.6,0) size (1.8,0.0692104,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) ...(+ 192 objects not shown)... time for set_epsilon = 1.79243 s ----------- Meep progress: 12.530000000000001/125.0 = 10.0% done in 4.0s, 35.9s to go on time step 1253 (time=12.53), 0.00319322 s/step Meep progress: 24.92/125.0 = 19.9% done in 8.0s, 32.1s to go on time step 2492 (time=24.92), 0.00322858 s/step Meep progress: 37.88/125.0 = 30.3% done in 12.0s, 27.6s to go on time step 3788 (time=37.88), 0.00308826 s/step Meep progress: 50.44/125.0 = 40.4% done in 16.0s, 23.7s to go on time step 5045 (time=50.45), 0.0031845 s/step Meep progress: 62.02/125.0 = 49.6% done in 20.0s, 20.3s to go on time step 6203 (time=62.03), 0.00345439 s/step Meep progress: 73.87/125.0 = 59.1% done in 24.0s, 16.6s to go on time step 7388 (time=73.88), 0.00337615 s/step Meep progress: 84.91/125.0 = 67.9% done in 28.0s, 13.2s to go on time step 8492 (time=84.92), 0.0036253 s/step Meep progress: 94.47/125.0 = 75.6% done in 32.0s, 10.3s to go on time step 9449 (time=94.49), 0.00418097 s/step Meep progress: 105.32000000000001/125.0 = 84.3% done in 36.0s, 6.7s to go on time step 10534 (time=105.34), 0.00368885 s/step Meep progress: 117.42/125.0 = 93.9% done in 40.0s, 2.6s to go on time step 11744 (time=117.44), 0.00330608 s/step run 0 finished at t = 125.0 (12500 timesteps) get_farfields_array working on point 996 of 1000 (99% done), 0.00401804 s/point ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.00165105 s Working in 2D dimensions. Computational cell is 7.8 x 120.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-60,0) size (1.8,0.09101,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-59.7,0) size (1.8,0.115166,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-59.4,0) size (1.8,0.141974,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-59.1,0) size (1.8,0.174379,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-58.8,0) size (1.8,0.0533026,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-58.5,0) size (1.8,0.0730401,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-58.2,0) size (1.8,0.0933667,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-57.9,0) size (1.8,0.117523,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-57.6,0) size (1.8,0.143741,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) ...(+ 392 objects not shown)... time for set_epsilon = 3.96254 s ----------- Meep progress: 5.39/125.0 = 4.3% done in 4.0s, 88.8s to go on time step 539 (time=5.39), 0.00742406 s/step Meep progress: 11.540000000000001/125.0 = 9.2% done in 8.0s, 78.7s to go on time step 1154 (time=11.54), 0.00651044 s/step Meep progress: 17.45/125.0 = 14.0% done in 12.0s, 74.0s to go on time step 1745 (time=17.45), 0.00677692 s/step Meep progress: 23.16/125.0 = 18.5% done in 16.0s, 70.4s to go on time step 2316 (time=23.16), 0.00700984 s/step Meep progress: 29.25/125.0 = 23.4% done in 20.0s, 65.5s to go on time step 2925 (time=29.25), 0.00657387 s/step Meep progress: 35.27/125.0 = 28.2% done in 24.0s, 61.1s to go on time step 3527 (time=35.27), 0.00664574 s/step Meep progress: 40.980000000000004/125.0 = 32.8% done in 28.0s, 57.5s to go on time step 4098 (time=40.98), 0.00700592 s/step Meep progress: 46.58/125.0 = 37.3% done in 32.0s, 53.9s to go on time step 4658 (time=46.58), 0.00714546 s/step Meep progress: 52.46/125.0 = 42.0% done in 36.0s, 49.8s to go on time step 5246 (time=52.46), 0.00680792 s/step Meep progress: 58.76/125.0 = 47.0% done in 40.0s, 45.1s to go on time step 5876 (time=58.76), 0.00635169 s/step Meep progress: 64.7/125.0 = 51.8% done in 44.0s, 41.0s to go on time step 6470 (time=64.7), 0.00673928 s/step Meep progress: 69.12/125.0 = 55.3% done in 48.0s, 38.8s to go on time step 6912 (time=69.12), 0.00905871 s/step Meep progress: 73.92/125.0 = 59.1% done in 52.0s, 36.0s to go on time step 7392 (time=73.92), 0.00833534 s/step Meep progress: 79.69/125.0 = 63.8% done in 56.0s, 31.9s to go on time step 7969 (time=79.69), 0.00694212 s/step Meep progress: 85.4/125.0 = 68.3% done in 60.0s, 27.8s to go on time step 8540 (time=85.4), 0.00700793 s/step Meep progress: 91.31/125.0 = 73.0% done in 64.0s, 23.6s to go on time step 9131 (time=91.31), 0.0067687 s/step Meep progress: 96.35000000000001/125.0 = 77.1% done in 68.1s, 20.2s to go on time step 9635 (time=96.35), 0.00797192 s/step Meep progress: 99.73/125.0 = 79.8% done in 72.1s, 18.3s to go on time step 9973 (time=99.73), 0.0118546 s/step Meep progress: 104.79/125.0 = 83.8% done in 76.1s, 14.7s to go on time step 10479 (time=104.79), 0.00790892 s/step Meep progress: 109.58/125.0 = 87.7% done in 80.1s, 11.3s to go on time step 10958 (time=109.58), 0.00836883 s/step Meep progress: 114.53/125.0 = 91.6% done in 84.1s, 7.7s to go on time step 11453 (time=114.53), 0.00808508 s/step Meep progress: 120.01/125.0 = 96.0% done in 88.1s, 3.7s to go on time step 12001 (time=120.01), 0.0073038 s/step run 0 finished at t = 125.0 (12500 timesteps) get_farfields_array working on point 446 of 1000 (44% done), 0.00898336 s/point get_farfields_array working on point 897 of 1000 (89% done), 0.00887155 s/point ----------- Initializing structure... Padding y to even number of grid points. Halving computational cell along direction y time for choose_chunkdivision = 0.000861883 s Working in 2D dimensions. Computational cell is 7.8 x 240.3 x 0 with resolution 50 block, center = (-2.4,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-120,0) size (1.8,0.109569,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-119.7,0) size (1.8,0.161122,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-119.4,0) size (1.8,0.0609619,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-119.1,0) size (1.8,0.0983747,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-118.8,0) size (1.8,0.145804,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-118.5,0) size (1.8,0.0518297,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-118.2,0) size (1.8,0.0883587,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-117.9,0) size (1.8,0.132253,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (1.11022e-16,-117.6,0) size (1.8,0.190581,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) ...(+ 792 objects not shown)... time for set_epsilon = 10.2552 s ----------- Meep progress: 2.2800000000000002/125.0 = 1.8% done in 4.0s, 215.4s to go on time step 228 (time=2.28), 0.0175505 s/step Meep progress: 5.0600000000000005/125.0 = 4.0% done in 8.0s, 189.8s to go on time step 506 (time=5.06), 0.0144061 s/step Meep progress: 7.79/125.0 = 6.2% done in 12.0s, 180.8s to go on time step 779 (time=7.79), 0.0146795 s/step Meep progress: 10.61/125.0 = 8.5% done in 16.0s, 172.7s to go on time step 1061 (time=10.61), 0.0142034 s/step Meep progress: 13.370000000000001/125.0 = 10.7% done in 20.0s, 167.2s to go on time step 1337 (time=13.37), 0.0145284 s/step Meep progress: 16.12/125.0 = 12.9% done in 24.0s, 162.4s to go on time step 1612 (time=16.12), 0.014564 s/step Meep progress: 17.88/125.0 = 14.3% done in 28.1s, 168.1s to go on time step 1788 (time=17.88), 0.0228108 s/step Meep progress: 20.5/125.0 = 16.4% done in 32.1s, 163.4s to go on time step 2050 (time=20.5), 0.015286 s/step Meep progress: 23.17/125.0 = 18.5% done in 36.1s, 158.5s to go on time step 2317 (time=23.17), 0.015032 s/step Meep progress: 25.310000000000002/125.0 = 20.2% done in 40.1s, 157.8s to go on time step 2531 (time=25.31), 0.0186941 s/step Meep progress: 28.01/125.0 = 22.4% done in 44.1s, 152.7s to go on time step 2801 (time=28.01), 0.0148529 s/step Meep progress: 30.72/125.0 = 24.6% done in 48.1s, 147.6s to go on time step 3072 (time=30.72), 0.0147973 s/step Meep progress: 33.410000000000004/125.0 = 26.7% done in 52.1s, 142.8s to go on time step 3341 (time=33.41), 0.0148998 s/step Meep progress: 36.11/125.0 = 28.9% done in 56.1s, 138.1s to go on time step 3611 (time=36.11), 0.0148395 s/step Meep progress: 38.31/125.0 = 30.6% done in 60.1s, 136.1s to go on time step 3831 (time=38.31), 0.0182697 s/step Meep progress: 40.82/125.0 = 32.7% done in 64.1s, 132.3s to go on time step 4082 (time=40.82), 0.0159877 s/step Meep progress: 43.5/125.0 = 34.8% done in 68.2s, 127.7s to go on time step 4350 (time=43.5), 0.0149743 s/step Meep progress: 46.17/125.0 = 36.9% done in 72.2s, 123.2s to go on time step 4617 (time=46.17), 0.0150016 s/step Meep progress: 48.86/125.0 = 39.1% done in 76.2s, 118.7s to go on time step 4886 (time=48.86), 0.0149015 s/step Meep progress: 51.57/125.0 = 41.3% done in 80.2s, 114.2s to go on time step 5157 (time=51.57), 0.0148085 s/step Meep progress: 54.26/125.0 = 43.4% done in 84.2s, 109.8s to go on time step 5426 (time=54.26), 0.0148808 s/step Meep progress: 56.97/125.0 = 45.6% done in 88.2s, 105.3s to go on time step 5697 (time=56.97), 0.014799 s/step Meep progress: 59.660000000000004/125.0 = 47.7% done in 92.2s, 101.0s to go on time step 5966 (time=59.66), 0.0149044 s/step Meep progress: 62.36/125.0 = 49.9% done in 96.2s, 96.7s to go on time step 6236 (time=62.36), 0.0148605 s/step Meep progress: 65.07000000000001/125.0 = 52.1% done in 100.2s, 92.3s to go on time step 6507 (time=65.07), 0.014803 s/step Meep progress: 67.77/125.0 = 54.2% done in 104.3s, 88.0s to go on time step 6777 (time=67.77), 0.0148711 s/step Meep progress: 69.79/125.0 = 55.8% done in 108.3s, 85.6s to go on time step 6979 (time=69.79), 0.0198233 s/step Meep progress: 72.57000000000001/125.0 = 58.1% done in 112.3s, 81.1s to go on time step 7257 (time=72.57), 0.0144427 s/step Meep progress: 75.16/125.0 = 60.1% done in 116.3s, 77.1s to go on time step 7516 (time=75.16), 0.0154471 s/step Meep progress: 77.78/125.0 = 62.2% done in 120.3s, 73.0s to go on time step 7778 (time=77.78), 0.0152936 s/step Meep progress: 80.57000000000001/125.0 = 64.5% done in 124.3s, 68.5s to go on time step 8057 (time=80.57), 0.0143637 s/step Meep progress: 83.37/125.0 = 66.7% done in 128.3s, 64.1s to go on time step 8337 (time=83.37), 0.0143218 s/step Meep progress: 86.07000000000001/125.0 = 68.9% done in 132.3s, 59.8s to go on time step 8607 (time=86.07), 0.0148534 s/step Meep progress: 88.82000000000001/125.0 = 71.1% done in 136.3s, 55.5s to go on time step 8882 (time=88.82), 0.0145593 s/step Meep progress: 91.74/125.0 = 73.4% done in 140.3s, 50.9s to go on time step 9174 (time=91.74), 0.0137382 s/step Meep progress: 94.64/125.0 = 75.7% done in 144.3s, 46.3s to go on time step 9464 (time=94.64), 0.0138319 s/step Meep progress: 97.55/125.0 = 78.0% done in 148.3s, 41.7s to go on time step 9755 (time=97.55), 0.0137471 s/step Meep progress: 100.47/125.0 = 80.4% done in 152.4s, 37.2s to go on time step 10047 (time=100.47), 0.0137389 s/step Meep progress: 103.35000000000001/125.0 = 82.7% done in 156.4s, 32.8s to go on time step 10335 (time=103.35), 0.01392 s/step Meep progress: 106.25/125.0 = 85.0% done in 160.4s, 28.3s to go on time step 10625 (time=106.25), 0.0138042 s/step Meep progress: 109.14/125.0 = 87.3% done in 164.4s, 23.9s to go on time step 10914 (time=109.14), 0.0138697 s/step Meep progress: 112.01/125.0 = 89.6% done in 168.4s, 19.5s to go on time step 11201 (time=112.01), 0.0139446 s/step Meep progress: 114.66/125.0 = 91.7% done in 172.4s, 15.5s to go on time step 11466 (time=114.66), 0.01511 s/step Meep progress: 117.28/125.0 = 93.8% done in 176.4s, 11.6s to go on time step 11728 (time=117.28), 0.0153079 s/step Meep progress: 119.9/125.0 = 95.9% done in 180.4s, 7.7s to go on time step 11990 (time=119.9), 0.0152988 s/step Meep progress: 122.79/125.0 = 98.2% done in 184.4s, 3.3s to go on time step 12279 (time=122.79), 0.0138704 s/step run 0 finished at t = 125.0 (12500 timesteps) get_farfields_array working on point 236 of 1000 (23% done), 0.0170414 s/point get_farfields_array working on point 472 of 1000 (47% done), 0.0170053 s/point get_farfields_array working on point 704 of 1000 (70% done), 0.0172601 s/point get_farfields_array working on point 937 of 1000 (93% done), 0.0171973 s/point
Shown below is the supercell lens design involving 201 unit cells. Note that even though periodic boundaries are used in the supercell calculation (via the k_point
), the choice of cell boundaries in the y (or longitudinal) direction is irrelevant given the finite length of the lens. For example, PMLs could also have been used (at the expense of a larger cell). Although add_near2far
does support periodic boundaries (via the nperiods
parameter), it is not necessary for this particular example.
The far-field energy-density profile is shown below for the three lens designs. As the number of unit cells increases, the focal spot becomes sharper and sharper. This is expected since the longer the focal length, the bigger the lens required to demonstrate focusing (which means more unit cells). In this example, the largest lens design contains 801 unit cells which corresponds to 0.24 mm or 1.2X the focal length.
x = np.linspace(focal_length-0.5*spot_length,focal_length+0.5*spot_length,ff_res*spot_length)
plt.figure(dpi=200)
plt.semilogy(x,abs(ff_nc[:,0])**2,'bo-',label='num_cells = {}'.format(2*num_cells[0]+1))
plt.semilogy(x,abs(ff_nc[:,1])**2,'ro-',label='num_cells = {}'.format(2*num_cells[1]+1))
plt.semilogy(x,abs(ff_nc[:,2])**2,'go-',label='num_cells = {}'.format(2*num_cells[2]+1))
plt.xlabel('x coordinate (μm)')
plt.ylabel(r'energy density of far-field electric fields, |E$_z$|$^2$')
plt.title('focusing properties of a binary-grating metasurface lens')
plt.legend(loc='upper right')
plt.tight_layout()
plt.show()