fix two errors in Near_to_Far_Fields_Spectra.md

doc/docs/Python_Tutorials/Near_to_Far_Field_Spectra.md

@@ -173,7 +173,7 @@ plt.show()
 
 As a second example, we compute the radiation pattern of an antenna positioned a given height $h$ above a perfect-electric conductor (PEC) ground plane. Depending on the wavelength and height of the antenna, self-interference effects due to reflections from the ground plane will produce well-defined lobes in the radiation pattern. The challenge in setting up this calculation is that because the ground plane is infinitely extended, it is not possible to enclose the antenna by a near-field surface. A non-closed near-field surface unfortunately gives rise to truncation errors which is described in more detail in the [section below](#truncation-errors-from-a-non-closed-near-field-surface).
 
-A workaround is to transform this problem into radiation in free space by making use of the fact that the effect of the ground plane can be exactly reproduced by two antennas of *opposite* phase separated by a distance of $2h$. This is known as the method of images. Additionally, the odd-mirror symmetry plane can be used to divide the cell in half in order to reduce the computional cost.
+A workaround is to transform this problem into radiation in free space by making use of the fact that the effect of the ground plane can be exactly reproduced by two antennas of *opposite* phase separated by a distance of $2h$. This is known as the method of images. Additionally, the odd-mirror symmetry plane can be used to divide the cell in half in order to reduce the computational cost.
 
 We can validate the radiation pattern computed by Meep using analytic theory. The radiation pattern of a two-element antenna array is equivalent to the radiation pattern of a single antenna multiplied by its "array factor" (AF) as derived in Section 6.2 "Two-Element Array" of [Antenna Theory: Analysis and Design, Fourth Edition (2016)](https://www.amazon.com/Antenna-Theory-Analysis-Constantine-Balanis/dp/1118642066) by C.A. Balanis. In this example, we consider an $E_z$-polarized antenna at a vacuum wavelength ($\lambda$) of 0.65 μm embedded within a medium with $n$ of 1.2 and positioned 1.25 μm above the ground plane. The outgoing (radial) flux is computed along the circumference of a circle with radius 1000$\lambda$ (or 650 μm) centered at the midpoint between the two antennas. The angular range is [0,90] degrees with 0° being the direction normal to the ground plane. A schematic showing the simulation layout and a plot of the radiation pattern computed by Meep and analytic theory are shown in the figure below. There is good agreement between the two results.
 
@@ -384,7 +384,7 @@ The total flux computed using the near and far fields is shown to be in close ag
 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.
+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.
 
@@ -635,6 +635,350 @@ if __name__ == "__main__":
     )
 ```
 
+### Extraction Efficiency of a Collection of Dipoles in a Disc
+
+[Tutorial/Radiation Pattern of a Disc in Cylindrical Coordinates](Near_to_Far_Field_Spectra.md#radiation-pattern-of-a-disc-in-cylindrical-coordinates) demonstrated the procedure for computing the radiation pattern of a *single* dipole (actually a "ring" current source with angular dependence $exp(im\phi)$). [Tutorial/Nonaxisymmetric Dipole Sources](Cylindrical_Coordinates.md#nonaxisymmetric-dipole-sources) described the method for modeling a point dipole at $r > 0$ in cylindrical coordinates using a Fourier-series expansion of the fields in $\phi$. [Tutorial/Extraction Efficiency of a Light-Emitting Diode](Local_Density_of_States.md#extraction-efficiency-of-a-light-emitting-diode-led) described the procedure for computing the extraction efficiency of a dipole at $r = 0$. These three results can be combined to compute the extraction efficiency for a point dipole *anywhere* in the cylindrical cell. Computing the extraction efficiency of an actual light-emitting diode (LED), however, involves a collection of spatially incoherent dipole emitters. [Tutorial/Stochastic Dipole Emission in Light Emitting Diodes](Custom_Source.md#stochastic-dipole-emission-in-light-emitting-diodes) described a method for computing the emission of a collection of dipoles using a series of single-dipole calculations and averaging the emission profiles in post processing. The example used a 2D simulation involving a 1D binary grating (or photonic crystal). This tutorial demonstrates how this approach for modeling spatially incoherent dipoles can be extended to cylindrical coordinates for structures with rotational symmetry.
+
+The example uses the same setup as the [previous tutorial](#radiation-pattern-of-a-disc-in-cylindrical-coordinates) involving a dielectric disc above a lossless-reflector ground plane. The dipoles are arranged on a line extending from $r = 0$ to $r = R$ where $R$ is the disc radius. The height of the dipoles ($z$ coordinate) within the disc is fixed. The radiation pattern $P(r,\theta)$ for a dipole at $r$ is computed using a Fourier-series expansion in $\phi$. The *total* radiation pattern $P(\theta)$ for $N$ dipoles with equal spacing $\Delta r$ is the average of the per-dipole results:
+
+$$P(\theta) = \frac{1}{R} \int_0^R P(r,\theta) s(r) 2\pi rdr = \frac{1}{N} \sum_{n=0}^{N-1} P(r_n,\theta) s(r_n) 2\pi r_n \Delta r$$,
+
+where $s(r)$ is a weighting function necessary for ensuring equal contribution from all dipoles relative to the dipole at $r = 0$. Note: an $E_r$ dipole at $r = 0$ must be placed at $r_0 = 1.5\Delta r$ due to an [interpolation bug](https://github.com/NanoComp/meep/issues/2704). $s(r)$ can be determined empirically by computing the emission in vacuum for a set of dipoles at different radial positions. The emission profile of a dipole in vacuum is a constant independent of its position. This criteria is used to obtain: $s(r) = \frac{1}{2(r/r_0)^2}$. This weighting function is also used to average the flux emitted by each dipole (obtained using using the LDOS feature). This quantity is the denominator in the expression for the extraction efficiency.
+
+This figure shows the radiation pattern from $N=11$ dipoles with $\lambda$ of 1.0 $\mu$m in the middle of a disc of height 0.29 $\mu$m, radius 1.2 $\mu$m, and refractive index 2.4.
+
+![](../images/disc_dipoles_radiation_pattern.png#center)
+
+The simulation script is in [examples/disc_extraction_efficiency.py](https://github.com/NanoComp/meep/blob/master/python/examples/disc_extraction_efficiency.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
+
+
+N_DISC = 2.4  # refractive index of disc
+DISC_RADIUS_UM = 1.2  # radius of disc
+WAVELENGTH_UM = 1.0  # wavelength (in vacuum)
+
+# radius of quarter circle for computing flux in far field
+FARFIELD_RADIUS_UM = 1000 * WAVELENGTH_UM
+
+RESOLUTION = 50  # pixels/μm
+
+
+def plot_radiation_pattern_polar(theta_rad: np.ndarray, radial_flux: np.ndarray):
+    """Plots the radiation pattern in polar coordinates.
+
+    Args:
+      theta_rad: angles of the radiation pattern. The angles increase clockwise
+        with zero at the pole (+z direction) and π/2 at the equator (+r
+        direction).
+      radial_flux: radial flux of the far fields in polar coordinates.
+    """
+    fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(6, 6))
+    ax.plot(
+        theta_rad,
+        radial_flux,
+        "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_ylabel("radial flux (a.u.)")
+    ax.set_title("radiation pattern in polar coordinates")
+
+    if mp.am_master():
+        fig.savefig(
+            "disc_radiation_pattern_polar.png",
+            dpi=150,
+            bbox_inches="tight",
+        )
+
+
+def plot_radiation_pattern_3d(theta_rad: np.ndarray, radial_flux: np.ndarray):
+    """Plots the radiation pattern in 3d Cartesian coordinates.
+
+    Args:
+      theta_rad: angles of the radiation pattern.
+      radial_flux: radial flux of the far fields in polar coordinates.
+    """
+    phis = np.linspace(0, 2 * np.pi, num_angles)
+
+    xs = np.zeros((len(theta_rad), len(phis)))
+    ys = np.zeros((len(theta_rad), len(phis)))
+    zs = np.zeros((len(theta_rad), len(phis)))
+
+    for i, theta in enumerate(theta_rad):
+        for j, phi in enumerate(phis):
+            xs[i, j] = radial_flux[i] * np.sin(theta) * np.cos(phi)
+            ys[i, j] = radial_flux[i] * np.sin(theta) * np.sin(phi)
+            zs[i, j] = radial_flux[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")
+    ax.set_box_aspect((np.amax(xs), np.amax(ys), np.amax(zs)))
+    ax.set_zlabel("radial flux (a.u.)")
+    ax.set(xticklabels=[], yticklabels=[])
+
+    if mp.am_master():
+        fig.savefig(
+            "disc_radiation_pattern_3d.png",
+            dpi=150,
+            bbox_inches="tight",
+        )
+
+
+def radiation_pattern(
+    theta_rad: np.ndarray, sim: mp.Simulation, n2f_mon: mp.DftNear2Far
+) -> np.ndarray:
+    """Computes the radiation pattern from the near fields.
+
+    Args:
+      theta_rad: angles of the radiation pattern.
+      sim: a `Simulation` object.
+      n2f_mon: a `DftNear2Far` object returned by `Simulation.add_near2far`.
+    """
+    e_field = np.zeros((theta_rad.shape[0], 3), dtype=np.complex128)
+    h_field = np.zeros((theta_rad.shape[0], 3), dtype=np.complex128)
+    for n in range(num_angles):
+        far_field = sim.get_farfield(
+            n2f_mon,
+            mp.Vector3(
+                FARFIELD_RADIUS_UM * math.sin(theta_rad[n]),
+                0,
+                FARFIELD_RADIUS_UM * math.cos(theta_rad[n]),
+            ),
+        )
+        e_field[n, :] = [np.conj(far_field[j]) for j in range(3)]
+        h_field[n, :] = [far_field[j + 3] for j in range(3)]
+
+    flux_r = np.real(e_field[:, 1] * h_field[:, 2] - e_field[:, 2] * h_field[:, 1])
+    flux_z = np.real(e_field[:, 0] * h_field[:, 1] - e_field[:, 1] * h_field[:, 0])
+    flux_rz = np.sqrt(np.square(flux_r) + np.square(flux_z))
+
+    return flux_rz
+
+
+def radiation_pattern_flux(theta_rad: np.ndarray, radial_flux: np.ndarray) -> float:
+    """Computes the total flux from the radiation pattern.
+
+    Based on integrating the radiation pattern over solid angles spanned by
+    polar angles in the range of [0, π/2].
+
+    Args:
+      theta_rad: angles of the radiation pattern.
+      radial_flux: radial flux of the far fields in polar coordinates.
+    """
+    dphi = 2 * math.pi
+    dtheta_rad = theta_rad[1] - theta_rad[0]
+
+    total_flux = (
+        np.sum(radial_flux * np.sin(theta_rad))
+        * FARFIELD_RADIUS_UM**2
+        * dtheta_rad
+        * dphi
+    )
+
+    return total_flux
+
+
+def dipole_in_disc(
+    t_disc_um: float, h_disc: float, rpos: float, m: int, theta_rad: np.ndarray
+) -> Tuple[float, np.ndarray]:
+    """Computes the total flux and radiation pattern of a dipole in a disc.
+
+    Args:
+      t_disc_um: thickness of disc.
+      h_disc: height of dipole above ground plane as fraction of t_disc_um.
+      rpos: radial position of dipole.
+      m: angular φ dependence of the fields exp(imφ).
+      theta_rad: angles of the radiation pattern.
+
+    Returns:
+      A 2-tuple of the total flux and the radiation pattern.
+    """
+    t_pml_um = 0.5  # thickness of PML
+    t_air_um = 1.0  # thickness of air padding above disc
+    length_r_um = 5.0  # length of cell in r
+
+    frequency = 1 / WAVELENGTH_UM  # center frequency of source/monitor
+
+    # field decay threshold for runtime termination criteria
+    decay_tol = 1e-6
+
+    size_r = length_r_um + t_pml_um
+    size_z = t_disc_um + t_air_um + t_pml_um
+    cell_size = mp.Vector3(size_r, 0, size_z)
+
+    boundary_layers = [
+        mp.PML(t_pml_um, direction=mp.R),
+        mp.PML(t_pml_um, direction=mp.Z, side=mp.High),
+    ]
+
+    src_cmpt = mp.Er
+    src_pt = mp.Vector3(rpos, 0, -0.5 * size_z + h_disc * t_disc_um)
+    sources = [
+        mp.Source(
+            src=mp.GaussianSource(frequency, fwidth=0.1 * frequency),
+            component=src_cmpt,
+            center=src_pt,
+        )
+    ]
+
+    geometry = [
+        mp.Block(
+            material=mp.Medium(index=N_DISC),
+            center=mp.Vector3(0.5 * DISC_RADIUS_UM, 0, -0.5 * size_z + 0.5 * t_disc_um),
+            size=mp.Vector3(DISC_RADIUS_UM, mp.inf, t_disc_um),
+        )
+    ]
+
+    sim = mp.Simulation(
+        resolution=RESOLUTION,
+        cell_size=cell_size,
+        dimensions=mp.CYLINDRICAL,
+        m=m,
+        boundary_layers=boundary_layers,
+        sources=sources,
+        geometry=geometry,
+    )
+
+    n2f_mon = sim.add_near2far(
+        frequency,
+        0,
+        1,
+        mp.FluxRegion(
+            center=mp.Vector3(0.5 * length_r_um, 0, 0.5 * size_z - t_pml_um),
+            size=mp.Vector3(length_r_um, 0, 0),
+        ),
+        mp.FluxRegion(
+            center=mp.Vector3(
+                length_r_um, 0, 0.5 * size_z - t_pml_um - 0.5 * (t_air_um + t_disc_um)
+            ),
+            size=mp.Vector3(0, 0, t_air_um + t_disc_um),
+        ),
+    )
+
+    sim.run(
+        mp.dft_ldos(frequency, 0, 1),
+        until_after_sources=mp.stop_when_fields_decayed(
+            50,
+            src_cmpt,
+            src_pt,
+            decay_tol,
+        ),
+    )
+
+    delta_vol = 2 * np.pi * rpos / (RESOLUTION**2)
+    dipole_flux = -np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) * delta_vol
+
+    dipole_radiation_pattern = radiation_pattern(theta_rad, sim, n2f_mon)
+
+    return dipole_flux, dipole_radiation_pattern
+
+
+if __name__ == "__main__":
+    disc_thickness = 0.7 * WAVELENGTH_UM / N_DISC
+    dipole_height = 0.5
+    num_dipoles = 11
+    dipole_rpos = np.linspace(0, DISC_RADIUS_UM, num_dipoles)
+
+    # An Er source at r=0 needs to be slighty offset.
+    # https://github.com/NanoComp/meep/issues/2704
+    dipole_rpos[0] = 1.5 / RESOLUTION
+
+    delta_rpos = dipole_rpos[2] - dipole_rpos[1]
+
+    # number of angular grid points in [0, π/2]
+    num_angles = 100
+
+    # grid of polar angles for computing radiated flux in far field
+    theta_rad = np.linspace(0, 0.5 * math.pi, num_angles)
+
+    # r = 0 source requires a single simulation with m = ±1.
+    m = -1
+    dipole_flux, dipole_radiation_pattern = dipole_in_disc(
+        disc_thickness, dipole_height, dipole_rpos[0], m, theta_rad
+    )
+
+    flux_total = dipole_flux * dipole_rpos[0] * delta_rpos
+    radiation_pattern_total = dipole_radiation_pattern * dipole_rpos[0] * delta_rpos
+
+    print(
+        f"dipole:, {dipole_rpos[0]:.4f}, "
+        f"{radiation_pattern_flux(theta_rad, dipole_radiation_pattern):.6f}"
+    )
+
+    # r > 0 source requires Fourier-series expansion of φ.
+    flux_tol = 1e-3  # threshold flux to determine when to truncate expansion
+    for rpos in dipole_rpos[1:]:
+        dipole_flux_total = 0
+        dipole_radiation_pattern_total = np.zeros((num_angles,))
+        dipole_radiation_pattern_flux_max = 0
+        m = 0
+        while True:
+            dipole_flux, dipole_radiation_pattern = dipole_in_disc(
+                disc_thickness, dipole_height, rpos, m, theta_rad
+            )
+            dipole_flux_total += dipole_flux if m == 0 else 2 * dipole_flux
+            dipole_radiation_pattern_total += (
+                dipole_radiation_pattern if m == 0 else 2 * dipole_radiation_pattern
+            )
+
+            dipole_radiation_pattern_flux = radiation_pattern_flux(
+                theta_rad, dipole_radiation_pattern
+            )
+            if dipole_radiation_pattern_flux > dipole_radiation_pattern_flux_max:
+                dipole_radiation_pattern_flux_max = dipole_radiation_pattern_flux
+
+            if m > 0 and (
+                (dipole_radiation_pattern_flux / dipole_radiation_pattern_flux_max)
+                < flux_tol
+            ):
+                break
+
+            print(
+                f"dipole-m:, {rpos:.4f}, {m}, " f"{dipole_radiation_pattern_flux:.6f}"
+            )
+            m += 1
+
+        scale_factor = 1 / (2 * (rpos / dipole_rpos[0]) ** 2)
+        flux_total += dipole_flux_total * scale_factor * rpos * delta_rpos
+        radiation_pattern_total += (
+            dipole_radiation_pattern_total * scale_factor * rpos * delta_rpos
+        )
+
+        dipole_radiation_pattern_total_flux = radiation_pattern_flux(
+            theta_rad, dipole_radiation_pattern_total
+        )
+        print(
+            f"dipole:, {rpos:.4f}, {m}, " f"{dipole_radiation_pattern_total_flux:.6f}"
+        )
+
+    flux_total /= num_dipoles
+    radiation_pattern_total /= num_dipoles
+
+    radiation_pattern_total_flux = radiation_pattern_flux(
+        theta_rad, radiation_pattern_total
+    )
+    extraction_efficiency = radiation_pattern_total_flux / flux_total
+    print(f"exteff:, {extraction_efficiency:.6f}")
+
+    plot_radiation_pattern_polar(radiation_pattern_total * FARFIELD_RADIUS_UM**2)
+    plot_radiation_pattern_3d(radiation_pattern_total * FARFIELD_RADIUS_UM**2)
+```
+
 Focusing Properties of a Metasurface Lens
 -----------------------------------------