In this example, we compute the diffraction spectrum of a binary phase grating with finite length. To compute the diffraction spectrum of the infinite periodic structure requires mode decomposition; for a demonstration, see Tutorials/Mode Decomposition/Diffraction Spectrum of a Binary Grating which also describes the grating geometry used in this example (i.e., periodicity of 10 μm, height of 0.5 μm, duty cycle of 0.5, and index 1.5 in air). Note that an infinite periodic structure actually has no spatial separation of the diffracted orders; they are all present at every far-field point. The focus of this tutorial is to demonstrate add_near2far
's support for periodic boundaries.
The simulation involves computing the scattered near fields of a finite-length grating for an Ez-polarized, pulsed planewave source spanning wavelengths of 0.4-0.6 μm at normal incidence. The far fields are then computed for 500 points along a line parallel to the grating axis positioned 100 m away (i.e., ≫ 2D2/λ, the Fraunhofer distance; D=NΛ where N is the number of unit cells and Λ is the grating periodicity, λ is the source wavelength) in the upper half plane of the symmetric finite structure with length corresponding to a 20° cone. The diffraction spectra is computed as the ratio of the energy density of the far fields from two separate runs: (1) an empty cell to obtain the fields from just the incident planewave and (2) a binary-grating unit cell to obtain the scattered fields.
Modeling a finite grating requires specifying the nperiods
parameter of add_near2far
which sums 2*nperiods+1
Bloch-periodic copies of the near fields. However, because of the way in which the edges of the structure are handled, this approach is only an approximation for a finite periodic surface. We will verify that the error from this approximation is O(1/nperiods
) by comparing its result with that of a true finite periodic structure involving multiple periods in a supercell arrangement terminated with a flat surface extending into PML. (There are infinitely many ways to terminate a finite periodic structure, of course, and different choices will have slightly different errors compared to the periodic approximation.)
import meep as mp
import math
import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt
resolution = 25 # pixels/μm
dpml = 1.0 # PML thickness
dsub = 3.0 # substrate thickness
dpad = 3.0 # padding between grating and PML
gp = 10.0 # grating period
gh = 0.5 # grating height
gdc = 0.5 # grating duty cycle
nperiods = 10 # number of unit cells in finite periodic grating
ff_distance = 1e8 # far-field distance from near-field monitor
ff_angle = 20 # far-field cone angle
ff_npts = 500 # number of far-field points
ff_length = ff_distance * math.tan(math.radians(ff_angle))
ff_res = ff_npts / ff_length
sx = dpml + dsub + gh + dpad + dpml
cell_size = mp.Vector3(sx)
pml_layers = [mp.PML(thickness=dpml, direction=mp.X)]
symmetries = [mp.Mirror(mp.Y)]
wvl_min = 0.4 # min wavelength
wvl_max = 0.6 # max wavelength
fmin = 1 / wvl_max # min frequency
fmax = 1 / wvl_min # max frequency
fcen = 0.5 * (fmin + fmax) # center frequency
df = fmax - fmin # frequency width
src_pt = mp.Vector3(-0.5 * sx + dpml + 0.5 * dsub)
sources = [
mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, center=src_pt)
]
k_point = mp.Vector3()
glass = mp.Medium(index=1.5)
sim = mp.Simulation(
resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
k_point=k_point,
default_material=glass,
sources=sources,
)
nfreq = 21
n2f_pt = mp.Vector3(0.5 * sx - dpml - 0.5 * dpad)
n2f_obj = sim.add_near2far(fcen, df, nfreq, mp.Near2FarRegion(center=n2f_pt))
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, n2f_pt, 1e-9))
ff_source = sim.get_farfields(
n2f_obj,
ff_res,
center=mp.Vector3(ff_distance, 0.5 * ff_length),
size=mp.Vector3(y=ff_length),
)
sim.reset_meep()
### unit cell with periodic boundaries
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),
)
]
geometry = [
mp.Block(
material=glass,
size=mp.Vector3(dpml + dsub, mp.inf, mp.inf),
center=mp.Vector3(-0.5 * sx + 0.5 * (dpml + dsub)),
),
mp.Block(
material=glass,
size=mp.Vector3(gh, gdc * gp, mp.inf),
center=mp.Vector3(-0.5 * sx + dpml + dsub + 0.5 * gh),
),
]
sim = mp.Simulation(
resolution=resolution,
split_chunks_evenly=True,
cell_size=cell_size,
boundary_layers=pml_layers,
geometry=geometry,
k_point=k_point,
sources=sources,
symmetries=symmetries,
)
n2f_obj = sim.add_near2far(
fcen,
df,
nfreq,
mp.Near2FarRegion(center=n2f_pt, size=mp.Vector3(y=sy)),
nperiods=nperiods,
)
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, n2f_pt, 1e-9))
ff_unitcell = sim.get_farfields(
n2f_obj,
ff_res,
center=mp.Vector3(ff_distance, 0.5 * ff_length),
size=mp.Vector3(y=ff_length),
)
sim.reset_meep()
### finite periodic grating with flat surface termination extending into PML
num_cells = 2 * nperiods + 1
sy = dpml + num_cells * gp + dpml
cell_size = mp.Vector3(sx, sy)
pml_layers = [mp.PML(thickness=dpml)]
sources = [
mp.Source(
mp.GaussianSource(fcen, fwidth=df),
component=mp.Ez,
center=src_pt,
size=mp.Vector3(y=sy - 2 * dpml),
)
]
geometry = [
mp.Block(
material=glass,
size=mp.Vector3(dpml + dsub, mp.inf, mp.inf),
center=mp.Vector3(-0.5 * sx + 0.5 * (dpml + dsub)),
)
]
for j in range(num_cells):
geometry.append(
mp.Block(
material=glass,
size=mp.Vector3(gh, gdc * gp, mp.inf),
center=mp.Vector3(
-0.5 * sx + dpml + dsub + 0.5 * gh, -0.5 * sy + dpml + (j + 0.5) * gp
),
)
)
sim = mp.Simulation(
resolution=resolution,
split_chunks_evenly=True,
cell_size=cell_size,
boundary_layers=pml_layers,
geometry=geometry,
k_point=k_point,
sources=sources,
symmetries=symmetries,
)
n2f_obj = sim.add_near2far(
fcen, df, nfreq, mp.Near2FarRegion(center=n2f_pt, size=mp.Vector3(y=sy - 2 * dpml))
)
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, n2f_pt, 1e-9))
ff_supercell = sim.get_farfields(
n2f_obj,
ff_res,
center=mp.Vector3(ff_distance, 0.5 * ff_length),
size=mp.Vector3(y=ff_length),
)
norm_err = LA.norm(ff_unitcell["Ez"] - ff_supercell["Ez"]) / nperiods
print("error:, {}, {}".format(nperiods, norm_err))
----------- Initializing structure... Working in 2D dimensions. Computational cell is 8.52 x 0.04 x 0 with resolution 25 time for set_epsilon = 0.00124812 s ----------- field decay(t = 50.02): 0.10083813086471559 / 0.10083813086471559 = 1.0 field decay(t = 100.04): 3.8807013552778427e-16 / 0.10083813086471559 = 3.848446338701171e-15 run 0 finished at t = 100.04 (5002 timesteps) Field time usage: connecting chunks: 0.00231314 s time stepping: 0.169821 s communicating: 0.0808802 s Fourier transforming: 0.0359166 s everything else: 0.156131 s ----------- Initializing structure... Halving computational cell along direction y Working in 2D dimensions. Computational cell is 8.52 x 10 x 0 with resolution 25 block, center = (-2.25,0,0) size (4,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 = (0,0,0) size (0.5,5,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.105644 s ----------- field decay(t = 50.02): 0.09814324428153094 / 0.09814324428153094 = 1.0 field decay(t = 100.04): 7.939880373487791e-06 / 0.09814324428153094 = 8.090093649962981e-05 field decay(t = 150.06): 6.627624476069797e-06 / 0.09814324428153094 = 6.753011401434805e-05 field decay(t = 200.08): 2.146203858444386e-06 / 0.09814324428153094 = 2.1868075323532677e-05 on time step 11243 (time=224.86), 0.00035578 s/step field decay(t = 250.1): 8.032368775670381e-07 / 0.09814324428153094 = 8.184331824846705e-06 field decay(t = 300.12): 2.965177550201506e-07 / 0.09814324428153094 = 3.0212752511988307e-06 field decay(t = 350.14): 1.3086054960500633e-07 / 0.09814324428153094 = 1.3333627858237849e-06 field decay(t = 400.16): 5.285162741753926e-08 / 0.09814324428153094 = 5.385151856803365e-07 field decay(t = 450.18): 1.8336367759809585e-08 / 0.09814324428153094 = 1.8683270452330267e-07 on time step 22548 (time=450.96), 0.000353856 s/step field decay(t = 500.2): 7.33127235003134e-09 / 0.09814324428153094 = 7.469971472515275e-08 field decay(t = 550.22): 2.8980535977675605e-09 / 0.09814324428153094 = 2.952881391870729e-08 field decay(t = 600.24): 9.12195167093687e-10 / 0.09814324428153094 = 9.294528357723631e-09 field decay(t = 650.26): 3.6982320415205274e-10 / 0.09814324428153094 = 3.768198278540582e-09 on time step 33740 (time=674.8), 0.000357411 s/step field decay(t = 700.28): 1.7964764774382362e-10 / 0.09814324428153094 = 1.830463717181505e-09 field decay(t = 750.3000000000001): 1.1056722890352689e-10 / 0.09814324428153094 = 1.1265903192109368e-09 field decay(t = 800.32): 3.1235531072986875e-11 / 0.09814324428153094 = 3.182647089124699e-10 run 0 finished at t = 800.32 (40016 timesteps) Field time usage: connecting chunks: 0.00979877 s time stepping: 9.10302 s communicating: 1.25385 s Fourier transforming: 2.74567 s everything else: 1.15844 s ----------- Initializing structure... Halving computational cell along direction y Working in 2D dimensions. Computational cell is 8.52 x 212 x 0 with resolution 25 block, center = (-2.25,0,0) size (4,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 = (0,-100,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,-90,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,-80,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,-70,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,-60,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,-50,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,-40,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,-30,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,-20,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,-10,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,0,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,10,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,20,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,30,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,40,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,50,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,60,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,70,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,80,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,90,0) size (0.5,5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (0,100,0) size (0.5,5,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 = 3.14788 s ----------- on time step 597 (time=11.94), 0.00670513 s/step on time step 1213 (time=24.26), 0.00650265 s/step on time step 1848 (time=36.96), 0.00630352 s/step on time step 2483 (time=49.66), 0.00630129 s/step field decay(t = 50.02): 0.09814324428153119 / 0.09814324428153119 = 1.0 on time step 3124 (time=62.48), 0.00624995 s/step on time step 3763 (time=75.26), 0.00626324 s/step on time step 4390 (time=87.8), 0.00638796 s/step field decay(t = 100.04): 7.939880373488582e-06 / 0.09814324428153119 = 8.090093649963766e-05 on time step 5016 (time=100.32), 0.00639027 s/step on time step 5655 (time=113.1), 0.00626925 s/step on time step 6286 (time=125.72), 0.00634788 s/step on time step 6913 (time=138.26), 0.00638669 s/step field decay(t = 150.06): 6.627624476069556e-06 / 0.09814324428153119 = 6.753011401434542e-05 on time step 7542 (time=150.84), 0.00636212 s/step on time step 8176 (time=163.52), 0.00631072 s/step on time step 8807 (time=176.14), 0.00633942 s/step on time step 9430 (time=188.6), 0.00642652 s/step field decay(t = 200.08): 3.77258555623287e-08 / 0.09814324428153119 = 3.8439584750336235e-07 on time step 10053 (time=201.06), 0.00642324 s/step on time step 10693 (time=213.86), 0.006254 s/step on time step 11327 (time=226.54), 0.00630954 s/step on time step 11966 (time=239.32), 0.00626743 s/step field decay(t = 250.1): 3.5244910064367445e-10 / 0.09814324428153119 = 3.591170265705177e-09 on time step 12599 (time=251.98), 0.00632365 s/step on time step 13244 (time=264.88), 0.00620512 s/step on time step 13874 (time=277.48), 0.00634975 s/step on time step 14516 (time=290.32), 0.00623725 s/step field decay(t = 300.12): 1.0847045885564748e-11 / 0.09814324428153119 = 1.1052259342934691e-10 run 0 finished at t = 300.12 (15006 timesteps) error:, 10, 5.981324591508142e-05
A plot of (a) the diffraction/far-field spectra and (b) its cross section at a fixed wavelength of 0.5 μm, is generated using the commands below and shown in the accompanying figure for two cases: (1) nperiods = 1
(no tiling; default) and (2) nperiods = 10
(21 copies). Note that because the evenly-spaced points on the line used to compute the far fields are mapped to angles in the plot, the angular data is not evenly spaced. A similar non-uniformity occurs when transforming the far-field data from the frequency to wavelength domain.
freqs = mp.get_near2far_freqs(n2f_obj)
wvl = np.divide(1, freqs)
ff_lengths = np.linspace(0, ff_length, ff_npts)
angles = [math.degrees(math.atan(f)) for f in ff_lengths / ff_distance]
wvl_slice = 0.5
idx_slice = np.where(np.asarray(freqs) == 1 / wvl_slice)[0][0]
rel_enh = np.absolute(ff_unitcell["Ez"]) ** 2 / np.absolute(ff_source["Ez"]) ** 2
plt.figure(dpi=150)
plt.subplot(1, 2, 1)
plt.pcolormesh(wvl, angles, rel_enh, cmap="Blues", shading="flat")
plt.axis([wvl_min, wvl_max, 0, ff_angle])
plt.xlabel("wavelength (μm)")
plt.ylabel("angle (degrees)")
plt.grid(linewidth=0.5, linestyle="--")
plt.xticks([t for t in np.arange(wvl_min, wvl_max + 0.1, 0.1)])
plt.yticks([t for t in range(0, ff_angle + 1, 10)])
plt.title("far-field spectra")
plt.subplot(1, 2, 2)
plt.plot(angles, rel_enh[:, idx_slice], "bo-")
plt.xlim(0, ff_angle)
plt.ylim(0)
plt.xticks([t for t in range(0, ff_angle + 1, 10)])
plt.xlabel("angle (degrees)")
plt.ylabel("relative enhancement")
plt.grid(axis="x", linewidth=0.5, linestyle="--")
plt.title("f.-f. spectra @ λ = {:.1} μm".format(wvl_slice))
plt.tight_layout(pad=0.5)
plt.show()
For the case of nperiods = 1
, three diffraction orders are present in the far-field spectra as broad peaks with finite angular width (a fourth peak/order is also visible). When nperiods = 10
, the diffraction orders become sharp, narrow peaks. The three diffraction orders are labeled in the right inset of the bottom figure as m=1, 3, and 5 corresponding to angles 2.9°, 8.6°, and 14.5° which, along with the diffraction efficiency, can be computed analytically using scalar theory as described in Tutorials/Mode Decomposition/Diffraction Spectrum of a Binary Grating. As an additional validation of the simulation results, the ratio of any two diffraction peaks pa/pb (a,b = 1,3,5,...) is consistent with that of its diffraction efficiencies: b2/a2.
Finally, we verify that the error in add_near2far
— defined as the L2-norm of the difference of the two far-field datasets from the unit- and super-cell calculations normalized by nperiods
— is O(1/nperiods
) by comparing results for three values of nperiods
: 5, 10, and 20. The error values, which are displayed in the output in the line prefixed by error:
, are: 0.0001195599054639075
, 5.981324591508146e-05
, and 2.989829913961854e-05
. The pairwise ratios of these errors is nearly 2 as expected (i.e., doubling nperiods
results in halving the error).
For a single process, the far-field calculation in both runs takes roughly the same amount of time. The wall-clock time is indicated by the getting farfields
category of the Field time usage
statistics displayed as part of the output after time stepping is complete. Time-stepping a supercell, however, which for nperiods=20
is more than 41 times larger than the unit cell (because of the PML termination) results in a total wall-clock time that is more than 40% larger. The slowdown is also due to the requirement of computing 41 times as many Fourier-transformed near fields. Thus, in the case of the unit-cell simulation, the reduced accuracy is a tradeoff for shorter runtime and less storage. In this example which involves multiple output wavelengths, the time for the far-field calculation can be reduced further on a single, shared-memory, multi-core machine via multithreading by compiling Meep with OpenMP and specifying the environment variable OMP_NUM_THREADS
to be an integer greater than one prior to execution.