diff --git a/doc/docs/Python_Tutorials/Cylindrical_Coordinates.md b/doc/docs/Python_Tutorials/Cylindrical_Coordinates.md
index ad0c02549..b6bb5aff9 100644
--- a/doc/docs/Python_Tutorials/Cylindrical_Coordinates.md
+++ b/doc/docs/Python_Tutorials/Cylindrical_Coordinates.md
@@ -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})$$
 
@@ -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.
diff --git a/python/examples/zone_plate.py b/python/examples/zone_plate.py
index 17db0e0b9..97daa0c34 100644
--- a/python/examples/zone_plate.py
+++ b/python/examples/zone_plate.py
@@ -1,153 +1,185 @@
+"""Computes the diffraction spectra of a zone plate in cylindrical coords."""
+
 import math
 
 import matplotlib.pyplot as plt
+import meep as mp
 import numpy as np
 
-import meep as mp
 
-resolution = 25  # pixels/μm
+resolution_um = 25
 
-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
+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
 
-pml_layers = [mp.PML(thickness=dpml)]
+pml_layers = [mp.PML(thickness=pml_um)]
 
-wvl_cen = 0.5
-frq_cen = 1 / wvl_cen
-dfrq = 0.2 * frq_cen
+wavelength_um = 0.5
+frequency = 1 / wavelength_um
+frequench_width = 0.2 * frequency
 
-## 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)
-]
+# 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
+
+# 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)
+    )
 
-sr = r[-1] + dpad + dpml
-sz = dpml + dsub + zh + dpad + dpml
-cell_size = mp.Vector3(sr, 0, sz)
+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(frq_cen, fwidth=dfrq, is_integrated=True),
+        mp.GaussianSource(frequency, fwidth=frequench_width, is_integrated=True),
         component=mp.Er,
-        center=mp.Vector3(0.5 * sr, 0, -0.5 * sz + dpml),
-        size=mp.Vector3(sr),
+        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(frq_cen, fwidth=dfrq, is_integrated=True),
+        mp.GaussianSource(frequency, fwidth=frequench_width, is_integrated=True),
         component=mp.Ep,
-        center=mp.Vector3(0.5 * sr, 0, -0.5 * sz + dpml),
-        size=mp.Vector3(sr),
+        center=mp.Vector3(0.5 * size_r_um, 0, -0.5 * size_z_um + pml_um),
+        size=mp.Vector3(size_r_um),
         amplitude=-1j,
     ),
 ]
 
 glass = mp.Medium(index=1.5)
 
+# Add the substrate.
 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)),
+        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)
+        ),
     )
 ]
 
-geometry.extend(
-    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 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,
+            ),
+        )
     )
-    for n in range(zN - 1, -1, -1)
-)
 
 sim = mp.Simulation(
     cell_size=cell_size,
     boundary_layers=pml_layers,
-    resolution=resolution,
+    resolution=resolution_um,
     sources=sources,
     geometry=geometry,
     dimensions=mp.CYLINDRICAL,
     m=-1,
 )
 
-## near-field monitor
-n2f_obj = sim.add_near2far(
-    frq_cen,
+# 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 * (sr - dpml), 0, 0.5 * sz - dpml),
-        size=mp.Vector3(sr - dpml),
+        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(sr - dpml, 0, 0.5 * sz - dpml - 0.5 * (dsub + zh + dpad)),
-        size=mp.Vector3(z=dsub + zh + dpad),
+        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)
+    fig.savefig("zone_plate_layout.png", bbox_inches="tight", dpi=150)
 
-sim.run(until_after_sources=100)
+# 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
+    )
+)
 
-ff_r = sim.get_farfields(
-    n2f_obj,
-    ff_res,
+farfields_r = sim.get_farfields(
+    n2f_monitor,
+    farfield_resolution_um,
     center=mp.Vector3(
-        0.5 * (sr - dpml), 0, -0.5 * sz + dpml + dsub + zh + focal_length
+        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(sr - dpml),
+    size=mp.Vector3(size_r_um - pml_um, 0, 0),
 )
 
-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),
+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),
 )
 
-E2_r = (
-    np.absolute(ff_r["Ex"]) ** 2
-    + np.absolute(ff_r["Ey"]) ** 2
-    + np.absolute(ff_r["Ez"]) ** 2
+intensity_r = (
+    np.absolute(farfields_r["Ex"]) ** 2
+    + np.absolute(farfields_r["Ey"]) ** 2
+    + np.absolute(farfields_r["Ez"]) ** 2
 )
-E2_z = (
-    np.absolute(ff_z["Ex"]) ** 2
-    + np.absolute(ff_z["Ey"]) ** 2
-    + np.absolute(ff_z["Ez"]) ** 2
+intensity_z = (
+    np.absolute(farfields_z["Ex"]) ** 2
+    + np.absolute(farfields_z["Ey"]) ** 2
+    + np.absolute(farfields_z["Ez"]) ** 2
 )
 
-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(list(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(f"binary-phase zone plate with focal length $z$ = {focal_length} μm")
+# 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$")
 
-    plt.tight_layout()
-    plt.savefig("zone_plate_farfields.png")
+fig.suptitle(f"binary-phase zone plate with focal length $z$ = {focal_length_um} μm")
+
+if mp.am_master():
+    fig.savefig("zone_plate_farfields.png", dpi=200, bbox_inches="tight")