Tutorial example for computing the radiation pattern of axisymmetric and nonaxisymmetric linearly polarized dipoles in cylindrical coordinates

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 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.
+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]
 
@@ -1702,3 +1702,20 @@ In this example, to demonstrate agreement between the far fields and DFT fields,
 When these two conditions are not met as in the example below involving a small `dpad` and large `d2`, the error from the finite truncation and numerical dispersion can be large and therefore result in a significant mismatch between the far fields computed using the near-to-far field transformation versus the actual DFT fields at the same location.
 
 ![Comparison of the far fields from the near-to-far field transformation dominated by errors and the DFT fields at the same location for a holey-waveguide cavity.](../images/farfields_vs_DFTfields_holeycavity_mismatch.png#center)
+
+Radiation Pattern of an Antenna in Cylindrical Coordinates
+----------------------------------------------------------
+
+Earlier we showed how to compute the [radiation pattern of an antenna with linear polarization in vacuum](#radiation-pattern-of-an-antenna) using a simulation in 2D Cartesian coordinates. The same calculation can also be performed using [cylindrical coordinates](../Exploiting_Symmetry.md#cylindrical-symmetry). We will compute the radiation pattern for two different cases in which the dipole is (1) axisymmetric (i.e., at $r = 0$) or (2) nonaxisymmetric (i.e., at $r > 0$). The radiation pattern is validated using the analytic result from antenna theory.
+
+In this example, the radiation pattern is computed for $\phi = 0$ (i.e., the $rz$ or $xz$ plane). Note that the radiation pattern for an $\hat{x}$ or $\hat{y}$ polarized dipole is *not* rotationally invariant about the $z$ axis. This means that the radiation pattern depends on the choice of $\phi$. (This is different than the computation of the [extraction efficiency for nonaxisymmetric dipoles](Cylindrical_Coordinates.md#nonaxisymmetric-dipole-sources) for which the radiation pattern is rotationally invariant because the dipoles are arranged along the circumference of a circle.)
+
+For (1), there are two dipole configurations: $E_x$ and $E_z$. An $E_z$ dipole is positioned at $r = 0$ with $m = 0$. This involves a single simulation. An $E_x$ dipole at $r = 0$, however, involves the superposition of left- and right-circularly polarized dipoles ($E_r \pm iE_\phi$) as described in [Tutorial/Scattering Cross Section of a Finite Dielectric Cylinder](Cylindrical_Coordinates.md#scattering-cross-section-of-a-finite-dielectric-cylinder). This requires *two* simulations. The computation of the radiation pattern of an $E_x$ dipole at $r = 0$ is different from the [computation of its extraction efficiency](Local_Density_of_States.md#extraction-efficiency-of-a-light-emitting-diode-led) which involves a *single* $E_r$ source with either $m = +1$ or $m = -1$. This is because the latter calculation involves a circularly polarized source which emits exactly half the power as a linearly polarized source even though their radiation patterns are different: $\frac{1}{2}(1 + \cos^2\theta)$ vs. $\cos^2\theta$.
+
+For (2), an $E_x$ (or equivalently an $E_r$) dipole positioned at $r > 0$ requires a [Fourier-series expansion of the fields](Cylindrical_Coordinates.md#nonaxisymmetric-dipole-sources) from an $E_r$ "ring" current source with azimuthal dependence $\exp(im\phi)$. The $m = -1$ fields can be obtained directly from the $m = +1$ fields using the formulas $(E_r, E_\phi, E_z)_m = (E_r, -E_\phi, E_z)_{-m}$ and $(H_r, H_\phi, H_z)_m = (-H_r, H_\phi, -H_z)_{-m}$. These formulas can be used to simplify the expressions for the Fourier-series expansion of the fields at $\phi = 0$: $\vec{E}_{tot}(\theta) = (E_r, 0, E_z)_{m=0} + 2\sum_{m=1}^M (E_r, 0, E_z)_m$ and $\vec{H}_{tot}(\theta) = (0, H_\phi, 0)_{m=0} + 2\sum_{m=1}^M (0, H_\phi, 0)_m$. An $E_y$ (or equivalently an $E_\phi$) dipole involves flipping the sign of the $-m$ fields resulting in: $\vec{E}_{tot}(\theta) = (0, E_\phi, 0)_{m=0} + 2\sum_{m=1}^M (0, E_\phi, 0)_m$ and $\vec{H}_{tot}(\theta) = (H_r, 0, H_z)_{m=0} + 2\sum_{m=1}^M (H_r, 0, H_z)_m$. We will compute the radiation pattern for $E_x$ and $E_y$ dipoles at $r = 0.1$ μm.
+
+The radiation pattern for each case is shown below. The simulation results agree with the analytic formulas. The simulation scripts are in [examples/dipole_in_vacuum_cyl_axisymmetric.py](https://github.com/NanoComp/meep/blob/master/python/examples/dipole_in_vacuum_cyl_axisymmetric.py) and [examples/dipole_in_vacuum_cyl_nonaxisymmetric.py](https://github.com/NanoComp/meep/blob/master/python/examples/dipole_in_vacuum_cyl_nonaxisymmetric.py).
+
+![](../images/radiation_pattern_axisymmetric_dipole.png#center)
+
+![](../images/radiation_pattern_nonaxisymmetric_dipole.png#center)

python/examples/dipole_in_vacuum_cyl_axisymmetric.py

@@ -0,0 +1,290 @@
+"""Radiation pattern of an axisymmetric dipole in cylindrical coordinates.
+
+Tutorial Reference:
+
+https://meep.readthedocs.io/en/latest/Python_Tutorials/Near_to_Far_Field_Spectra/#radiation-pattern-of-an-antenna-in-cylindrical-coordinates
+"""
+
+import argparse
+import cmath
+import math
+from typing import Tuple
+
+import matplotlib.pyplot as plt
+import meep as mp
+import numpy as np
+
+
+RESOLUTION_UM = 50
+WAVELENGTH_UM = 1.0
+PML_UM = 1.0
+FARFIELD_RADIUS_UM = 1e6 * WAVELENGTH_UM
+NUM_FARFIELD_PTS = 50
+AZIMUTHAL_RAD = 0
+POWER_DECAY_THRESHOLD = 1e-4
+
+frequency = 1 / WAVELENGTH_UM
+polar_rad = np.linspace(0, 0.5 * math.pi, NUM_FARFIELD_PTS)
+
+
+def plot_radiation_pattern(dipole_pol: str, radial_flux: np.ndarray):
+    """Plots the radiation pattern in polar coordinates.
+
+    The angles increase clockwise with zero in the +z direction (the "pole")
+    and π/2 in the +r direction (the "equator").
+
+    Args:
+        dipole_pol: the polarization the electric dipole. Either x or z.
+        radial_flux: the radial flux in polar coordinates.
+    """
+    normalized_radial_flux = radial_flux / np.max(radial_flux)
+
+    if dipole_pol == "x":
+        dipole_radial_flux = np.square(np.cos(polar_rad))
+        dipole_radial_flux_label = r"$\cos^2θ$"
+        dipole_name = "$E_x$"
+    else:
+        dipole_radial_flux = np.square(np.sin(polar_rad))
+        dipole_radial_flux_label = r"$\sin^2θ$"
+        dipole_name = "$E_z$"
+
+    fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(6, 6))
+    ax.plot(polar_rad, normalized_radial_flux, "b-", label="Meep")
+    ax.plot(polar_rad, dipole_radial_flux, "r--", label=dipole_radial_flux_label)
+    ax.legend()
+    ax.set_theta_direction(-1)
+    ax.set_theta_offset(0.5 * math.pi)
+    ax.set_thetalim(0, 0.5 * math.pi)
+    ax.set_rmax(1)
+    ax.set_rticks([0, 0.5, 1])
+    ax.grid(True)
+    ax.set_rlabel_position(22)
+    ax.set_ylabel("radial flux (a.u.)")
+    ax.set_title(
+        "radiation pattern (φ = 0) of an axisymmetric " f"{dipole_name} dipole"
+    )
+
+    if mp.am_master():
+        fig.savefig(
+            "dipole_radiation_pattern_axisymmetric.png",
+            dpi=150,
+            bbox_inches="tight",
+        )
+
+    relative_error = np.linalg.norm(
+        normalized_radial_flux - dipole_radial_flux
+    ) / np.linalg.norm(dipole_radial_flux)
+    print(f"relative error in radiation pattern:, {relative_error}")
+
+
+def radiation_pattern(e_field: np.ndarray, h_field: np.ndarray) -> np.ndarray:
+    """Computes the radiation pattern from the far fields.
+
+    Args:
+        e_field, h_field: the electric (Er, Ep, Ez) and magnetic (Hr, Hp, Hz)
+          far fields, respectively.
+
+    Returns:
+        The radial Poynting flux as a 1D array. One element for each point on
+        the circumference of a quarter circle with angular range of
+        [0, π/2] rad. 0 radians is the +z direction (the "pole") and π/2 is
+        the +r direction (the "equator").
+    """
+    flux_x = np.real(
+        e_field[:, 1] * np.conj(h_field[:, 2]) - e_field[:, 2] * np.conj(h_field[:, 1])
+    )
+    flux_z = np.real(
+        e_field[:, 0] * np.conj(h_field[:, 1]) - e_field[:, 1] * np.conj(h_field[:, 0])
+    )
+    flux_r = np.sqrt(np.square(flux_x) + np.square(flux_z))
+
+    return flux_r
+
+
+def get_farfields(
+    sim: mp.Simulation, n2f_mon: mp.DftNear2Far
+) -> Tuple[np.ndarray, np.ndarray]:
+    """Computes the far fields from the near fields for φ = 0 (rz plane).
+
+    Args:
+        sim: a `Simulation` object.
+        n2f_mon: a `DftNear2Far` object returned by `Simulation.add_near2far`.
+
+    Returns:
+        The electric (Er, Ep, Ez) and magnetic (Hr, Hp, Hz) far fields. One row
+        for each point on the circumference of a quarter circle with angular
+        range of [0, π/2] rad. Each row has six columns for the fields.
+        0 radians is the +z direction (the "pole") and π/2 is the +r direction
+        (the "equator").
+    """
+    e_field = np.zeros((NUM_FARFIELD_PTS, 3), dtype=np.complex128)
+    h_field = np.zeros((NUM_FARFIELD_PTS, 3), dtype=np.complex128)
+    for n in range(NUM_FARFIELD_PTS):
+        far_field = sim.get_farfield(
+            n2f_mon,
+            mp.Vector3(
+                FARFIELD_RADIUS_UM * math.sin(polar_rad[n]),
+                0,
+                FARFIELD_RADIUS_UM * math.cos(polar_rad[n]),
+            ),
+        )
+        e_field[n, :] = [far_field[j] for j in range(3)]
+        h_field[n, :] = [far_field[j + 3] for j in range(3)]
+
+    return e_field, h_field
+
+
+def dipole_in_vacuum(dipole_pol: str, m: int) -> Tuple[np.ndarray, np.ndarray]:
+    """Computes the far fields of an axisymmetric point source.
+
+    Args:
+        dipole_pol: the polarization of the electric dipole. Either x or z.
+        m: angular φ dependence of the fields exp(imφ).
+
+    Returns:
+        A 2-tuple containing the electric and magnetic far fields as 1D arrays.
+    """
+    sr = 2.0
+    sz = 4.0
+    cell_size = mp.Vector3(sr + PML_UM, 0, sz + 2 * PML_UM)
+
+    boundary_layers = [mp.PML(thickness=PML_UM)]
+
+    if dipole_pol == "x":
+        # An Er source at r = 0 needs to be slighty offset due to a bug.
+        # https://github.com/NanoComp/meep/issues/2704
+        dipole_pos_r = 1.5 / RESOLUTION_UM
+        sources = [
+            mp.Source(
+                src=mp.GaussianSource(frequency, fwidth=0.1 * frequency),
+                component=mp.Er,
+                center=mp.Vector3(dipole_pos_r, 0, 0),
+            ),
+            mp.Source(
+                src=mp.GaussianSource(frequency, fwidth=0.1 * frequency),
+                component=mp.Ep,
+                center=mp.Vector3(dipole_pos_r, 0, 0),
+                amplitude=1j if m == 1 else -1j,
+            ),
+        ]
+    else:
+        dipole_pos_r = 0
+        sources = [
+            mp.Source(
+                src=mp.GaussianSource(frequency, fwidth=0.1 * frequency),
+                component=mp.Ez,
+                center=mp.Vector3(dipole_pos_r, 0, 0),
+            )
+        ]
+
+    sim = mp.Simulation(
+        resolution=RESOLUTION_UM,
+        cell_size=cell_size,
+        dimensions=mp.CYLINDRICAL,
+        m=m,
+        boundary_layers=boundary_layers,
+        sources=sources,
+        force_complex_fields=True,
+    )
+
+    nearfields_monitor = sim.add_near2far(
+        frequency,
+        0,
+        1,
+        mp.FluxRegion(
+            center=mp.Vector3(0.5 * sr, 0, 0.5 * sz), size=mp.Vector3(sr, 0, 0)
+        ),
+        mp.FluxRegion(center=mp.Vector3(sr, 0, 0), size=mp.Vector3(0, 0, sz)),
+        mp.FluxRegion(
+            center=mp.Vector3(0.5 * sr, 0, -0.5 * sz),
+            size=mp.Vector3(sr, 0, 0),
+            weight=-1.0,
+        ),
+    )
+
+    sim.run(
+        until_after_sources=mp.stop_when_fields_decayed(
+            20.0,
+            mp.Er if dipole_pol == "x" else mp.Ez,
+            mp.Vector3(dipole_pos_r, 0, 0),
+            1e-6,
+        )
+    )
+
+    e_field, h_field = get_farfields(sim, nearfields_monitor)
+
+    return e_field, h_field
+
+
+def flux_from_farfields(e_field: np.ndarray, h_field: np.ndarray) -> float:
+    """Computes the flux from the far fields.
+
+    Args:
+        e_field, h_field: the electric (Er, Ep, Ez) and magnetic (Hr, Hp, Hz)
+          far fields, respectively.
+
+    Returns:
+        The Poynting flux obtained from the far fields.
+    """
+    dphi = 2 * math.pi
+    dtheta = 0.5 * math.pi / (NUM_FARFIELD_PTS - 1)
+    dipole_radiation_pattern = radiation_pattern(e_field, h_field)
+    flux = (
+        np.sum(dipole_radiation_pattern * np.sin(polar_rad))
+        * FARFIELD_RADIUS_UM**2
+        * dtheta
+        * dphi
+    )
+
+    return flux
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "dipole_pol",
+        type=str,
+        choices=["x", "z"],
+        help="polarization of the electric dipole (x or z)",
+    )
+    args = parser.parse_args()
+
+    e_field_total = np.zeros((NUM_FARFIELD_PTS, 3), dtype=np.complex128)
+    h_field_total = np.zeros((NUM_FARFIELD_PTS, 3), dtype=np.complex128)
+
+    if args.dipole_pol == "x":
+        # An x-polarized dipole can be formed from the superposition
+        # of left- and right-circularly polarized dipoles.
+
+        e_field, h_field = dipole_in_vacuum("x", +1)
+        e_field_total += 0.5 * e_field * cmath.exp(1j * AZIMUTHAL_RAD)
+        h_field_total += 0.5 * h_field * cmath.exp(1j * AZIMUTHAL_RAD)
+
+        e_field, h_field = dipole_in_vacuum("x", -1)
+        e_field_total += 0.5 * e_field * cmath.exp(-1j * AZIMUTHAL_RAD)
+        h_field_total += 0.5 * h_field * cmath.exp(-1j * AZIMUTHAL_RAD)
+
+    else:
+        e_field, h_field = dipole_in_vacuum("z", 0)
+        e_field_total += e_field
+        h_field_total += h_field
+
+    dipole_radiation_pattern = radiation_pattern(e_field_total, h_field_total)
+    dipole_radiation_pattern_scaled = dipole_radiation_pattern * FARFIELD_RADIUS_UM**2
+    plot_radiation_pattern(args.dipole_pol, dipole_radiation_pattern_scaled)
+
+    if mp.am_master():
+        np.savez(
+            "dipole_farfields_axisymmetric.npz",
+            AZIMUTHAL_RAD=AZIMUTHAL_RAD,
+            FARFIELD_RADIUS_UM=FARFIELD_RADIUS_UM,
+            PML_UM=PML_UM,
+            POWER_DECAY_THRESHOLD=POWER_DECAY_THRESHOLD,
+            RESOLUTION_UM=RESOLUTION_UM,
+            WAVELENGTH_UM=WAVELENGTH_UM,
+            dipole_pol=args.dipole_pol,
+            dipole_radiation_pattern=dipole_radiation_pattern,
+            e_field_total=e_field_total,
+            h_field_total=h_field_total,
+            polar_rad=polar_rad,
+        )