simplify the nonaxisymmetric computation by obtaining the -m fields d…

…irectly from the +m fields using the analytic expression

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]
 
@@ -1706,13 +1706,13 @@ When these two conditions are not met as in the example below involving a small
 Radiation Pattern of an Antenna in Cylindrical Coordinates
 ----------------------------------------------------------
 
-In an earlier section, we showed how to compute the [radiation pattern of an antennna](#radiation-pattern-of-an-antenna) in vacuum with **linear** polarization using 2D Cartesian coordinates. The same calculation can also be performed using [cylindrical coordinates](Exploiting_Symmetry.md#cylindrical-symmetry). We will demonstrate this for two cases in which the dipole is (1) axisymmetric (i.e., at $r = 0$) or (2) nonaxisymmetric (i.e., at $r > 0$). In this example, the radiation pattern is computed for $\phi = 0$ (i.e., the $rz$ or $xz$ plane).
+Earlier we showed how to compute the [radiation pattern of an antennna in vacuum with linear polarization](#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 demonstrate this for two different cases in which the dipole is (1) axisymmetric (i.e., at $r = 0$) or (2) nonaxisymmetric (i.e., at $r > 0$). In this example, the radiation pattern is computed for $\phi = 0$ (i.e., the $rz$ or $xz$ plane).
 
-For (1), an $E_z$ dipole is positioned exactly 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 as described in [Tutorial/Scattering Cross Section of a Finite Dielectric Cylinder](Cylindrical_Coordinates.md#scattering-cross-section-of-a-finite-dielectric-cylinder). This involves *two* simulations. Note that 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 in the latter a circularly polarized source emits the same power as a linearly polarized source.
+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 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. Note that 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 $\pm m$ fields are computed separately (even though they are complex conjugates) and combined as part of the calculation of the total fields of the expansion.
+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) = \vec{E}_{m=0}(\theta) + 2\sum_{m=1}^M (E_r, 0, E_z)_m$ and $\vec{H}_{tot}(\theta) = \vec{H}_{m=0}(\theta) + 2\sum_{m=1}^M (0, H_\phi, 0)_m$.
 
-The figures show the radiation pattern for each of these cases obtained using the scripts. The simulation and analytic results show good agreement.
+The radiation pattern for each of these cases obtained using the corresponding scripts is shown below. The simulation results show agreement 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).
 

python/examples/dipole_in_vacuum_cyl_axisymmetric.py

@@ -1,4 +1,9 @@
-"""Radiation pattern of an axisymmetric dipole in cylindrical coordinates."""
+"""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 math
@@ -35,11 +40,11 @@ def plot_radiation_pattern(dipole_pol: int, radial_flux: np.ndarray):
 
     if dipole_pol == mp.X:
         dipole_radial_flux = np.square(np.cos(polar_rad))
-        dipole_radial_flux_label = "$\cos^2θ$"
+        dipole_radial_flux_label = r"$\cos^2θ$"
         dipole_name = "$E_x$"
     elif dipole_pol == mp.Z:
         dipole_radial_flux = np.square(np.sin(polar_rad))
-        dipole_radial_flux_label = "$\sin^2θ$"
+        dipole_radial_flux_label = r"$\sin^2θ$"
         dipole_name = "$E_z$"
 
     fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(6, 6))
@@ -54,20 +59,18 @@ def plot_radiation_pattern(dipole_pol: int, radial_flux: np.ndarray):
     ax.grid(True)
     ax.set_rlabel_position(22)
     ax.set_ylabel("radial flux (a.u.)")
-    ax.set_title("radiation pattern of an axisymmetric " + dipole_name +
-                 " dipole")
+    ax.set_title("radiation pattern of an axisymmetric " + dipole_name + " dipole")
 
     if mp.am_master():
         fig.savefig(
             "dipole_radiation_pattern_axisymmetric.png",
-            dpi=100,
+            dpi=150,
             bbox_inches="tight",
         )
 
-    relative_error = (
-        np.linalg.norm(normalized_radial_flux - dipole_radial_flux) /
-        np.linalg.norm(dipole_radial_flux)
-    )
+    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}")
 
 
@@ -84,18 +87,15 @@ def radiation_pattern(e_field: np.ndarray, h_field: np.ndarray) -> np.ndarray:
         [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] * 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_x = 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_r = np.sqrt(np.square(flux_x) + np.square(flux_z))
 
     return flux_r
 
 
 def get_farfields(
-        sim: mp.Simulation,
-        n2f_mon: mp.DftNear2Far
+    sim: mp.Simulation, n2f_mon: mp.DftNear2Far
 ) -> Tuple[np.ndarray, np.ndarray]:
     """Computes the far fields from the near fields.
 
@@ -117,19 +117,16 @@ def get_farfields(
             mp.Vector3(
                 FARFIELD_RADIUS_UM * math.sin(polar_rad[n]),
                 0,
-                FARFIELD_RADIUS_UM * math.cos(polar_rad[n])
-            )
+                FARFIELD_RADIUS_UM * math.cos(polar_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)]
 
     return e_field, h_field
 
 
-def dipole_in_vacuum(
-        dipole_pol: int,
-        m: int
-) -> Tuple[np.ndarray, np.ndarray]:
+def dipole_in_vacuum(dipole_pol: int, m: int) -> Tuple[np.ndarray, np.ndarray]:
     """Computes the far fields of an axisymmetric point source.
 
     Args:
@@ -153,22 +150,22 @@ def dipole_in_vacuum(
             mp.Source(
                 src=mp.GaussianSource(frequency, fwidth=0.1 * frequency),
                 component=mp.Er,
-                center=mp.Vector3(dipole_pos_r, 0, 0)
+                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
-            )
+                amplitude=1j if m == 1 else -1j,
+            ),
         ]
     elif dipole_pol == mp.Z:
         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)
+                center=mp.Vector3(dipole_pos_r, 0, 0),
             )
         ]
 
@@ -179,46 +176,38 @@ def dipole_in_vacuum(
         m=m,
         boundary_layers=boundary_layers,
         sources=sources,
-        force_complex_fields=True
+        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)
+            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
-        )
+            weight=-1.0,
+        ),
     )
 
     sim.run(
         until_after_sources=mp.stop_when_fields_decayed(
             20.0,
             mp.Er if dipole_pol == mp.X else mp.Ez,
             mp.Vector3(dipole_pos_r, 0, 0),
-            1e-6
+            1e-6,
         )
     )
 
     if DEBUG_OUTPUT:
         fig, ax = plt.subplots()
         sim.plot2D(ax=ax, show_monitors=True)
         if mp.am_master():
-            fig.savefig(
-                "dipole_in_vacuum_cyl_layout.png",
-                dpi=150,
-                bbox_inches="tight"
-            )
+            fig.savefig("dipole_in_vacuum_cyl_layout.png", dpi=150, bbox_inches="tight")
 
     e_field, h_field = get_farfields(sim, nearfields_monitor)
 
@@ -239,10 +228,10 @@ def flux_from_farfields(e_field: np.ndarray, h_field: np.ndarray) -> float:
     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
+        np.sum(dipole_radiation_pattern * np.sin(polar_rad))
+        * FARFIELD_RADIUS_UM**2
+        * dtheta
+        * dphi
     )
 
     return flux
@@ -254,7 +243,7 @@ def flux_from_farfields(e_field: np.ndarray, h_field: np.ndarray) -> float:
         "dipole_pol",
         type=str,
         choices=["x", "z"],
-        help="polarization of the electric dipole (x or z)"
+        help="polarization of the electric dipole (x or z)",
     )
     args = parser.parse_args()
 
@@ -278,15 +267,10 @@ def flux_from_farfields(e_field: np.ndarray, h_field: np.ndarray) -> float:
         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
-    )
+    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(
-        mp.X if args.dipole_pol == "x" else mp.Z,
-        dipole_radiation_pattern_scaled
+        mp.X if args.dipole_pol == "x" else mp.Z, dipole_radiation_pattern_scaled
     )
 
     if mp.am_master():