Tutorial for computing the radiation pattern and total flux of a disc…

… in cylindrical coordinates (#2629)

* tutorial for computing the radiation pattern and total flux of an LED in cylindrical coordinates

* replace infinitely extended LED with finite-sized disc and update figures and results

* various fixes and tweaks to text and script

* add note describing error introduced by finite truncation of near-field surface and possible workaround

doc/docs/Python_Tutorials/Near_to_Far_Field_Spectra.md

@@ -2,7 +2,7 @@
 # Near to Far Field Spectra
 ---
 
-The [near-to-far field transformation](../Python_User_Interface.md#near-to-far-field-spectra) feature is demonstrated using five different examples. There are three steps involved in this type of calculation. First, the "near" surface(s) is defined as a set of surfaces capturing *all* outgoing radiation in *free space* in the desired direction(s). Second, the simulation is run using a pulsed source (or alternatively, a CW source via the [frequency-domain solver](../Python_User_Interface.md#frequency-domain-solver)) to allow Meep to accumulate the DFT fields on the near surface(s). Third, Meep computes the "far" fields at any desired points with the option to save the far fields to an HDF5 file.
+The [near-to-far field transformation](../Python_User_Interface.md#near-to-far-field-spectra) feature in Cartesian (2D/3D) and [cylindrical](../Cylindrical_Coordinates.md) coordinates is demonstrated using six different examples. Generally, there are three steps involved in this type of calculation. First, the "near" surface(s) is defined as a set of surfaces capturing *all* outgoing radiation in *free space* in the desired direction(s). Second, the simulation is run using a pulsed source (or alternatively, a CW source via the [frequency-domain solver](../Python_User_Interface.md#frequency-domain-solver)) to allow Meep to accumulate the DFT fields on the near surface(s). Third, Meep computes the "far" fields at any desired points with the option to save the far fields to an HDF5 file.
 
 [TOC]
 
@@ -357,6 +357,280 @@ if __name__ == '__main__':
     print("norm:, {:.6f}".format(np.linalg.norm(Pr_pec_norm-Pr_theory_norm)))
 ```
 
+Radiation Pattern of a Disc in Cylindrical Coordinates
+------------------------------------------------------
+
+The near-to-far field transformation feature can also be used in [cylindrical coordinates](Cylindrical_Coordinates.md). As a demonstration, we compute the radiation pattern of a dielectric disc and verify Poynting's theorem: the total radiated flux computed from the far fields is equivalent to using the near fields via `add_flux`. (The same result is demonstrated in [Tutorial/Radiation Pattern of an Antenna](#radiation-pattern-of-an-antenna) for 2D Cartesian coordinates.)
+
+The simulation consists of an $E_r$ point-dipole source ($\lambda$ = 1.0 μm) at $r$ = 0.6 μm embedded within a disc (radius of 1.2 μm) of index $n$ = 2.4 above a perfect-metallic ground plane. This is similar to the configuration in [Tutorial/Extraction Efficiency of a Light-Emitting Diode (LED)](Local_Density_of_States.md#extraction-efficiency-of-a-light-emitting-diode-led). Unlike the infinitely extended slab of the LED, a *finite* structure such as the disc ensures that all the power from the dipole emitter is radiated. The LED contains waveguide modes which are more challenging to disentagle from the radiated power.
+
+A schematic of the simulation layout is shown below. The flux and near-field monitors (shown in blue) are overlapping.
+
+![](../images/disc_radiation_layout.png#center)
+
+Obtaining the radiation pattern $P(\theta)$ of the disc involves computing the radial (or "outgoing") flux from the far fields along the circumference of a quarter circle (i.e. angular range of $[0, \pi/2]$). The radius $r$ of the circle needs to be sufficiently large ($\gg \lambda$) to ensure accurate results but is otherwise arbitrary. The total flux is then computed by integrating $P(\theta)$ over the surface of a hemisphere with radius $r$ in [spherical coordinates](https://en.wikipedia.org/wiki/Spherical_coordinate_system):
+
+$$P_{total} = \int_0^{2\pi} \int_0^{\frac{\pi}{2}} P(\theta) r^2 \sin(\theta) d\theta d\phi = 2 \pi r^2 \sum_{n=0}^{N-1} P(\theta_n) \sin(\theta_n) \Delta \theta$$
+
+An angular grid of $N$ equally spaced points in $[0, \pi/2]$ has $\Delta \theta = \frac{\pi}{2(N - 1)}$. Note that the same $r^2 \sin(\theta)$ weighting is necessary for the power in any cone, not just over all angles.
+
+A plot of the radiation pattern in polar coordinates and 3D is shown below. Note regarding the coordinate axes in the polar plot: 0° is in the $+z$ direction which is normal to the ground plane and 90° is in the $+r$ direction which is parallel to the ground plane. This is consistent with the convention for the polar angle $\theta$ used in spherical coordinates.
+
+![](../images/disc_radiation_pattern_polar_vs_3d.png#center)
+
+The total flux computed using the near and far fields is shown to be in close agreement with a relative error of ~7%.
+
+```
+total_flux:, 643.65058 (near), 597.72713 (far), 0.07135 (error)
+```
+
+The error decreases with increasing (1) grid resolution, (2) runtime, and (3) number of angular grid points. However, this only applies to a *closed* near-field surface which is not the case in this example. This is because the ground plane, which extends to infinity, contains $H_r$ and $H_\phi$ fields on its surface which are not zero (unlike the $E_r$ and $E_phi$ fields). These magnetic fields produce equivalent currents which radiate into the far field. The PML in the $r$ direction does not mitigate this effect.
+
+Because the near-field surface actually extends to infinity in the $r$ direction, one approach to reducing the error introduced by its finite truncation would be to simply make the cell size in the $r$ direction larger (the parameter `L` in the script below). Another option which would remove this error entirely would be to simulate the same structure using a closed surface by removing the ground plane and duplicating the structure and source below the $z = 0$ plane. This is known as the method of images. See [Tutorial/Antenna above a Perfect Electric Conductor Ground Plane ](#antenna-above-a-perfect-electric-conductor-ground-plane) for a demonstration of this approach.
+
+The simulation script is in [examples/disc_radiation_pattern.py](https://github.com/NanoComp/meep/blob/master/python/examples/disc_radiation_pattern.py).
+
+```py
+import math
+from typing import Tuple
+
+import matplotlib
+import meep as mp
+import numpy as np
+
+matplotlib.use("agg")
+import matplotlib.pyplot as plt
+
+
+resolution = 100  # pixels/μm
+dpml = 0.5  # thickness of PML
+dair = 1.0  # thickness of air padding
+L = 6.0  # length of non-PML region
+n = 2.4  # refractive index of surrounding medium
+wvl = 1.0  # wavelength (in vacuum)
+
+fcen = 1 / wvl  # center frequency of source/monitor
+
+# field decay threshold for runtime termination criteria
+tol = 1e-8
+
+# number of angular grid points in [0, π/2]
+npts = 100
+
+# grid of polar angles for computing radiated flux in far field
+thetas = np.linspace(0, 0.5 * math.pi, npts)
+
+# radius of quarter circle for computing flux in far field
+r = 1000 * wvl
+
+
+def plot_radiation_pattern_polar(Ptheta: np.ndarray):
+    """Plots the radiation pattern in polar coordinates.
+
+    The angles increase clockwise with zero at the top (+z direction).
+
+    Args:
+        Ptheta: radial flux of the far fields in polar coordinates.
+    """
+    fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(6,6))
+    ax.plot(
+        thetas,
+        Ptheta,
+        "b-",
+    )
+    ax.set_theta_direction(-1)
+    ax.set_theta_offset(0.5 * math.pi)
+    ax.set_thetalim(0, 0.5 * math.pi)
+    ax.grid(True)
+    ax.set_rlabel_position(22)
+    ax.set_title("radiation pattern in polar coordinates")
+
+    if mp.am_master():
+        fig.savefig(
+            "led_radpattern_polar.png",
+            dpi=150,
+            bbox_inches="tight",
+        )
+
+
+def plot_radiation_pattern_3d(Ptheta: np.ndarray):
+    """Plots the radiation pattern in 3d Cartesian coordinates.
+
+    Args:
+        Ptheta: radial flux of the far fields in polar coordinates.
+    """
+    phis = np.linspace(0, 2 * np.pi, npts)
+
+    xs = np.zeros((len(thetas), len(phis)))
+    ys = np.zeros((len(thetas), len(phis)))
+    zs = np.zeros((len(thetas), len(phis)))
+
+    for i, theta in enumerate(thetas):
+        for j, phi in enumerate(phis):
+            xs[i, j] = Ptheta[i] * np.sin(theta) * np.cos(phi)
+            ys[i, j] = Ptheta[i] * np.sin(theta) * np.sin(phi)
+            zs[i, j] = Ptheta[i] * np.cos(theta)
+
+    fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, figsize=(6,6))
+    ax.plot_surface(xs, ys, zs, cmap="inferno")
+    ax.set_title("radiation pattern in 3d")
+
+    if mp.am_master():
+        fig.savefig(
+            "led_radpattern_3d.png",
+            dpi=150,
+            bbox_inches="tight",
+        )
+
+
+def radiation_pattern(sim: mp.Simulation, n2f_mon: mp.DftNear2Far) -> np.ndarray:
+    """Computes the radiation pattern from the far fields.
+
+    Args:
+        sim: a `Simulation` object.
+        n2f_mon: a `DftNear2Far` object returned by `Simulation.add_near2far`.
+
+    Returns:
+        Array of radial Poynting flux, one for each point on the circumference of
+        a quarter circle with angular range of [0, π/2] rad. 0 rad is the +z
+        direction and π/2 is +r.
+    """
+    E = np.zeros((npts, 3), dtype=np.complex128)
+    H = np.zeros((npts, 3), dtype=np.complex128)
+    for n in range(npts):
+        ff = sim.get_farfield(
+            n2f_mon, mp.Vector3(r * math.sin(thetas[n]), 0, r * math.cos(thetas[n]))
+        )
+        E[n, :] = [np.conj(ff[j]) for j in range(3)]
+        H[n, :] = [ff[j + 3] for j in range(3)]
+
+    Pr = np.real(E[:, 1] * H[:, 2] - E[:, 2] * H[:, 1])
+    Pz = np.real(E[:, 0] * H[:, 1] - E[:, 1] * H[:, 0])
+    Prz = np.sqrt(np.square(Pr) + np.square(Pz))
+
+    return Prz
+
+
+def disc_total_flux(dmat: float, h: float) -> Tuple[float, float]:
+    """Computes the total radiated flux from a point dipole embedded
+    within a dielectric disc above a lossless ground plane using
+    its near and far fields as separate calculations.
+
+    Args:
+        dmat: thickness of dielectric disc.
+        h: height of dipole above ground plane as fraction of dmat.
+
+    Returns:
+        A 2-tuple of the total flux computed using the near and far fields,
+        respectively.
+    """
+    sr = L + dpml
+    sz = dmat + dair + dpml
+    cell_size = mp.Vector3(sr, 0, sz)
+
+    boundary_layers = [
+        mp.PML(dpml, direction=mp.R),
+        mp.PML(dpml, direction=mp.Z, side=mp.High),
+    ]
+
+    src_cmpt = mp.Er
+    src_pt = mp.Vector3(0.1 * L, 0, -0.5 * sz + h * dmat)
+    sources = [
+        mp.Source(
+            src=mp.GaussianSource(fcen, fwidth=0.1 * fcen),
+            component=src_cmpt,
+            center=src_pt,
+        )
+    ]
+
+    geometry = [
+        mp.Block(
+            material=mp.Medium(index=n),
+            center=mp.Vector3(0.1 * L, 0, -0.5 * sz + 0.5 * dmat),
+            size=mp.Vector3(0.2 * L, mp.inf, dmat),
+        )
+    ]
+
+    sim = mp.Simulation(
+        resolution=resolution,
+        cell_size=cell_size,
+        dimensions=mp.CYLINDRICAL,
+        m=-1,
+        boundary_layers=boundary_layers,
+        sources=sources,
+        geometry=geometry,
+    )
+
+    # flux monitor
+    flux_mon = sim.add_flux(
+        fcen,
+        0,
+        1,
+        mp.FluxRegion(
+            center=mp.Vector3(0.5 * L, 0, 0.5 * sz - dpml),
+            size=mp.Vector3(L, 0, 0),
+        ),
+        mp.FluxRegion(
+            center=mp.Vector3(L, 0, 0.5 * sz - dpml - 0.5 * (dair + dmat)),
+            size=mp.Vector3(0, 0, dair + dmat),
+        ),
+    )
+
+    # near-field monitor
+    n2f_mon = sim.add_near2far(
+        fcen,
+        0,
+        1,
+        mp.FluxRegion(
+            center=mp.Vector3(0.5 * L, 0, 0.5 * sz - dpml),
+            size=mp.Vector3(L, 0, 0),
+        ),
+        mp.FluxRegion(
+            center=mp.Vector3(L, 0, 0.5 * sz - dpml - 0.5 * (dair + dmat)),
+            size=mp.Vector3(0, 0, dair + dmat),
+        ),
+    )
+
+    fig, ax = plt.subplots()
+    sim.plot2D(ax=ax)
+    if mp.am_master():
+        fig.savefig("disc_simulation_layout.png", dpi=150, bbox_inches="tight")
+
+    sim.run(
+        until_after_sources=mp.stop_when_fields_decayed(
+            50,
+            src_cmpt,
+            src_pt,
+            tol,
+        ),
+    )
+
+    flux_near = mp.get_fluxes(flux_mon)[0]
+
+    Ptheta = radiation_pattern(sim, n2f_mon)
+    plot_radiation_pattern_polar(r * r * Ptheta)
+    plot_radiation_pattern_3d(r * r * Ptheta)
+
+    dtheta = 0.5 * math.pi / (npts - 1)
+    dphi = 2 * math.pi
+    flux_far = np.sum(Ptheta * np.sin(thetas)) * r * r * dtheta * dphi
+
+    return flux_near, flux_far
+
+
+if __name__ == "__main__":
+    disc_thickness = 0.7 * wvl / n
+    dipole_height = 0.5
+
+    near_flux, far_flux = disc_total_flux(disc_thickness, dipole_height)
+
+    err = abs(near_flux - far_flux) / near_flux
+    print(
+        f"total_flux:, {near_flux:.5f} (near), {far_flux:.5f} (far), "
+        f"{err:.5f} (error)"
+    )
+```
+
 Focusing Properties of a Metasurface Lens
 -----------------------------------------