Skip to content

Commit

Permalink
fix bug related to position of near-field monitor in zone-plate examp…
Browse files Browse the repository at this point in the history 
…le (#2783)
  • Loading branch information
@oskooi
oskooi authored Feb 10, 2024
1 parent 995c717 commit 9d93e2b
Show file tree
Hide file tree
Showing 2 changed files with 314 additions and 177 deletions.
277 changes: 191 additions & 86 deletions doc/docs/Python_Tutorials/Cylindrical_Coordinates.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -468,7 +468,7 @@ Focusing Properties of a Binary-Phase Zone Plate

It is also possible to compute a [near-to-far field transformation](../Python_User_Interface.md#near-to-far-field-spectra) in cylindrical coordinates. This is demonstrated in this example for a binary-phase [zone plate](https://en.wikipedia.org/wiki/Zone_plate) which is a rotationally-symmetric diffractive lens used to focus a normally-incident planewave to a single spot.

Using [scalar theory](http://zoneplate.lbl.gov/theory), the radius of the $n$<sup>th</sup> zone can be computed as:
Using [scalar theory](https://en.wikipedia.org/wiki/Zone_plate#Design_and_manufacture), the radius of the $n$<sup>th</sup> zone can be computed as:

$$ r_n^2 = n\lambda (f+\frac{n\lambda}{4})$$

Expand All @@ -481,110 +481,215 @@ where $n$ is the zone index (1,2,3,...,$N$), $f$ is the focal length, and $\lamb
The simulation script is in [examples/zone_plate.py](https://github.com/NanoComp/meep/blob/master/python/examples/zone_plate.py). The notebook is [examples/zone_plate.ipynb](https://nbviewer.jupyter.org/github/NanoComp/meep/blob/master/python/examples/zone_plate.ipynb).

```py
import meep as mp
import numpy as np
import math

import matplotlib.pyplot as plt
import meep as mp
import numpy as np

resolution = 25 # pixels/μm

dpml = 1.0 # PML thickness
dsub = 2.0 # substrate thickness
dpad = 2.0 # padding betweeen zone plate and PML
zh = 0.5 # zone-plate height
zN = 25 # number of zones (odd zones: π phase shift, even zones: none)
focal_length = 200 # focal length of zone plate
spot_length = 100 # far-field line length
ff_res = 10 # far-field resolution
resolution_um = 25

pml_layers = [mp.PML(thickness=dpml)]
pml_um = 1.0
substrate_um = 2.0
padding_um = 2.0
height_um = 0.5
focal_length_um = 200
scan_length_z_um = 100
farfield_resolution_um = 10

wvl_cen = 0.5
frq_cen = 1/wvl_cen
dfrq = 0.2*frq_cen
pml_layers = [mp.PML(thickness=pml_um)]

## radii of zones
## ref: eq. 7 of http://zoneplate.lbl.gov/theory
r = [math.sqrt(n*wvl_cen*(focal_length+n*wvl_cen/4)) for n in range(1,zN+1)]
wavelength_um = 0.5
frequency = 1 / wavelength_um
frequench_width = 0.2 * frequency

sr = r[-1]+dpad+dpml
sz = dpml+dsub+zh+dpad+dpml
cell_size = mp.Vector3(sr,0,sz)
# The number of zones in the zone plate.
# Odd-numbered zones impart a π phase shift and
# even-numbered zones impart no phase shift.
num_zones = 25

sources = [mp.Source(mp.GaussianSource(frq_cen,fwidth=dfrq,is_integrated=True),
component=mp.Er,
center=mp.Vector3(0.5*sr,0,-0.5*sz+dpml),
size=mp.Vector3(sr)),
mp.Source(mp.GaussianSource(frq_cen,fwidth=dfrq,is_integrated=True),
component=mp.Ep,
center=mp.Vector3(0.5*sr,0,-0.5*sz+dpml),
size=mp.Vector3(sr),
amplitude=-1j)]
# Specify the radius of each zone using the equation
# from https://en.wikipedia.org/wiki/Zone_plate.
zone_radius_um = np.zeros(num_zones)
for n in range(1, num_zones + 1):
zone_radius_um[n-1] = math.sqrt(
n * wavelength_um *
(focal_length_um + n * wavelength_um / 4)
)

glass = mp.Medium(index=1.5)
size_r_um = zone_radius_um[-1] + padding_um + pml_um
size_z_um = pml_um + substrate_um + height_um + padding_um + pml_um
cell_size = mp.Vector3(size_r_um, 0, size_z_um)

# Specify a (linearly polarized) planewave at normal incidence.
sources = [
mp.Source(
mp.GaussianSource(
frequency,
fwidth=frequench_width,
is_integrated=True
),
component=mp.Er,
center=mp.Vector3(0.5 * size_r_um, 0, -0.5 * size_z_um + pml_um),
size=mp.Vector3(size_r_um),
),
mp.Source(
mp.GaussianSource(
frequency,
fwidth=frequench_width,
is_integrated=True
),
component=mp.Ep,
center=mp.Vector3(0.5 * size_r_um, 0, -0.5 * size_z_um + pml_um),
size=mp.Vector3(size_r_um),
amplitude=-1j,
),
]

geometry = [mp.Block(material=glass,
size=mp.Vector3(sr,0,dpml+dsub),
center=mp.Vector3(0.5*sr,0,-0.5*sz+0.5*(dpml+dsub)))]
glass = mp.Medium(index=1.5)

for n in range(zN-1,-1,-1):
geometry.append(mp.Block(material=glass if n % 2 == 0 else mp.vacuum,
size=mp.Vector3(r[n],0,zh),
center=mp.Vector3(0.5*r[n],0,-0.5*sz+dpml+dsub+0.5*zh)))
# Add the substrate.
geometry = [
mp.Block(
material=glass,
size=mp.Vector3(size_r_um, 0, pml_um + substrate_um),
center=mp.Vector3(
0.5 * size_r_um,
0,
-0.5 * size_z_um + 0.5 * (pml_um + substrate_um)
),
)
]

sim = mp.Simulation(cell_size=cell_size,
boundary_layers=pml_layers,
resolution=resolution,
sources=sources,
geometry=geometry,
dimensions=mp.CYLINDRICAL,
m=-1)
# Add the zone plates starting with the ones with largest radius.
for n in range(num_zones - 1, -1, -1):
geometry.append(
mp.Block(
material=glass if n % 2 == 0 else mp.vacuum,
size=mp.Vector3(zone_radius_um[n], 0, height_um),
center=mp.Vector3(
0.5 * zone_radius_um[n],
0,
-0.5 * size_z_um + pml_um + substrate_um + 0.5 * height_um
),
)
)

## near-field monitor
n2f_obj = sim.add_near2far(frq_cen,
0,
1,
mp.Near2FarRegion(center=mp.Vector3(0.5*(sr-dpml),0,0.5*sz-dpml),
size=mp.Vector3(sr-dpml)),
mp.Near2FarRegion(center=mp.Vector3(sr-dpml,0,0.5*sz-dpml-0.5*(dsub+zh+dpad)),
size=mp.Vector3(z=dsub+zh+dpad)))
sim = mp.Simulation(
cell_size=cell_size,
boundary_layers=pml_layers,
resolution=resolution_um,
sources=sources,
geometry=geometry,
dimensions=mp.CYLINDRICAL,
m=-1,
)

# Add the near-field monitor (must be entirely in air).
n2f_monitor = sim.add_near2far(
frequency,
0,
1,
mp.Near2FarRegion(
center=mp.Vector3(
0.5 * (size_r_um - pml_um),
0,
0.5 * size_z_um - pml_um
),
size=mp.Vector3(size_r_um - pml_um, 0, 0),
),
mp.Near2FarRegion(
center=mp.Vector3(
size_r_um - pml_um,
0,
0.5 * size_z_um - pml_um - 0.5 * (height_um + padding_um)
),
size=mp.Vector3(0, 0, height_um + padding_um),
),
)

sim.plot2D()
fig, ax = plt.subplots()
sim.plot2D(ax=ax)
if mp.am_master():
plt.savefig("zone_plate_epsilon.png",bbox_inches='tight',dpi=150)

sim.run(until_after_sources=100)

ff_r = sim.get_farfields(n2f_obj,
ff_res,
center=mp.Vector3(0.5*(sr-dpml),0,-0.5*sz+dpml+dsub+zh+focal_length),
size=mp.Vector3(sr-dpml))

ff_z = sim.get_farfields(n2f_obj,
ff_res,
center=mp.Vector3(z=-0.5*sz+dpml+dsub+zh+focal_length),
size=mp.Vector3(z=spot_length))
fig.savefig("zone_plate_layout.png", bbox_inches="tight", dpi=150)

# Timestep the fields until they have sufficiently decayed away.
sim.run(
until_after_sources=mp.stop_when_fields_decayed(
50.0,
mp.Er,
mp.Vector3(0.5 * size_r_um, 0, 0),
1e-6
)
)

E2_r = np.absolute(ff_r['Ex'])**2+np.absolute(ff_r['Ey'])**2+np.absolute(ff_r['Ez'])**2
E2_z = np.absolute(ff_z['Ex'])**2+np.absolute(ff_z['Ey'])**2+np.absolute(ff_z['Ez'])**2
farfields_r = sim.get_farfields(
n2f_monitor,
farfield_resolution_um,
center=mp.Vector3(
0.5 * (size_r_um - pml_um),
0,
-0.5 * size_z_um + pml_um + substrate_um + height_um + focal_length_um
),
size=mp.Vector3(size_r_um - pml_um, 0, 0),
)

farfields_z = sim.get_farfields(
n2f_monitor,
farfield_resolution_um,
center=mp.Vector3(
0,
0,
-0.5 * size_z_um + pml_um + substrate_um + height_um + focal_length_um
),
size=mp.Vector3(0, 0, scan_length_z_um),
)

intensity_r = (
np.absolute(farfields_r["Ex"]) ** 2
+ np.absolute(farfields_r["Ey"]) ** 2
+ np.absolute(farfields_r["Ez"]) ** 2
)
intensity_z = (
np.absolute(farfields_z["Ex"]) ** 2
+ np.absolute(farfields_z["Ey"]) ** 2
+ np.absolute(farfields_z["Ez"]) ** 2
)

# Plot the intensity data and save the result to disk.
fig, ax = plt.subplots(ncols=2)

ax[0].semilogy(
np.linspace(0, size_r_um - pml_um, intensity_r.size),
intensity_r,
"bo-"
)
ax[0].set_xlim(-2, 20)
ax[0].set_xticks(np.arange(0, 25, 5))
ax[0].grid(True, axis="y", which="both", ls="-")
ax[0].set_xlabel(r"$r$ coordinate (μm)")
ax[0].set_ylabel(r"energy density of far fields, |E|$^2$")

ax[1].semilogy(
np.linspace(
focal_length_um - 0.5 * scan_length_z_um,
focal_length_um + 0.5 * scan_length_z_um,
intensity_z.size,
),
intensity_z,
"bo-",
)
ax[1].grid(True, axis="y", which="both", ls="-")
ax[1].set_xlabel(r"$z$ coordinate (μm)")
ax[1].set_ylabel(r"energy density of far fields, |E|$^2$")

fig.suptitle(
f"binary-phase zone plate with focal length $z$ = {focal_length_um} μm"
)

if mp.am_master():
plt.figure(dpi=200)
plt.subplot(1,2,1)
plt.semilogy(np.linspace(0,sr-dpml,len(E2_r)),E2_r,'bo-')
plt.xlim(-2,20)
plt.xticks([t for t in np.arange(0,25,5)])
plt.grid(True,axis="y",which="both",ls="-")
plt.xlabel(r'$r$ coordinate (μm)')
plt.ylabel(r'energy density of far fields, |E|$^2$')
plt.subplot(1,2,2)
plt.semilogy(np.linspace(focal_length-0.5*spot_length,focal_length+0.5*spot_length,len(E2_z)),E2_z,'bo-')
plt.grid(True,axis="y",which="both",ls="-")
plt.xlabel(r'$z$ coordinate (μm)')
plt.ylabel(r'energy density of far fields, |E|$^2$')
plt.suptitle(r"Binary-Phase Zone Plate with Focal Length $z$ = {} μm".format(focal_length))
plt.tight_layout()
plt.savefig("zone_plate_farfields.png")
fig.savefig("zone_plate_farfields.png", dpi=200, bbox_inches="tight")
```

Note that the volume specified in `get_farfields` via `center` and `size` is in cylindrical coordinates. These points must therefore lie in the $\phi = 0$ ($rz = xz$) plane. The fields $E$ and $H$ returned by `get_farfields` can be thought of as either cylindrical ($r$,$\phi$,$z$) or Cartesian ($x$,$y$,$z$) coordinates since these are the same in the $\phi = 0$ plane (i.e., $E_r=E_x$ and $E_\phi=E_y$). Also, `get_farfields` tends to gradually *slow down* as the far-field point gets closer to the near-field monitor. This performance degradation is unavoidable and is due to the larger number of $\phi$ integration points required for accurate convergence of the integral involving the Green's function which diverges as the evaluation point approaches the source point.
Expand Down
Loading

0 comments on commit 9d93e2b

Please sign in to comment.