make ϕ a parameter with default value of zero

doc/docs/Python_Tutorials/Near_to_Far_Field_Spectra.md

@@ -1706,15 +1706,15 @@ When these two conditions are not met as in the example below involving a small
 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 will be 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$.
+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.
 
-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. 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$.
+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 (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$.
+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$.
 
-The radiation pattern for each case is shown below. The simulation results agree with the analytic formulas.
+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 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).
+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)
 

python/examples/dipole_in_vacuum_cyl_axisymmetric.py

@@ -6,6 +6,7 @@
 """
 
 import argparse
+import cmath
 import math
 from typing import Tuple
 
@@ -19,8 +20,8 @@
 PML_UM = 1.0
 FARFIELD_RADIUS_UM = 1e6 * WAVELENGTH_UM
 NUM_FARFIELD_PTS = 50
+AZIMUTHAL_RAD = 0
 POWER_DECAY_THRESHOLD = 1e-4
-DEBUG_OUTPUT = False
 
 frequency = 1 / WAVELENGTH_UM
 polar_rad = np.linspace(0, 0.5 * math.pi, NUM_FARFIELD_PTS)
@@ -103,17 +104,18 @@ def radiation_pattern(e_field: np.ndarray, h_field: np.ndarray) -> np.ndarray:
 def get_farfields(
     sim: mp.Simulation, n2f_mon: mp.DftNear2Far
 ) -> Tuple[np.ndarray, np.ndarray]:
-    """Computes the far fields from the near fields.
+    """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
-        with six columns for the fields 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").
+        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)
@@ -209,12 +211,6 @@ def dipole_in_vacuum(dipole_pol: str, m: int) -> Tuple[np.ndarray, np.ndarray]:
         )
     )
 
-    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")
-
     e_field, h_field = get_farfields(sim, nearfields_monitor)
 
     return e_field, h_field
@@ -258,15 +254,15 @@ def flux_from_farfields(e_field: np.ndarray, h_field: np.ndarray) -> float:
 
     if args.dipole_pol == "x":
         # An x-polarized dipole can be formed from the superposition
-        # of two left- and right-circularly polarized dipoles.
+        # of left- and right-circularly polarized dipoles.
 
         e_field, h_field = dipole_in_vacuum("x", +1)
-        e_field_total += 0.5 * e_field
-        h_field_total += 0.5 * h_field
+        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
-        h_field_total += 0.5 * h_field
+        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)
@@ -280,6 +276,7 @@ def flux_from_farfields(e_field: np.ndarray, h_field: np.ndarray) -> float:
     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,

python/examples/dipole_in_vacuum_cyl_nonaxisymmetric.py

@@ -6,6 +6,7 @@
 """
 
 import argparse
+import cmath
 import math
 from typing import Tuple
 
@@ -19,8 +20,8 @@
 PML_UM = 1.0
 FARFIELD_RADIUS_UM = 1e6 * WAVELENGTH_UM
 NUM_FARFIELD_PTS = 50
+AZIMUTHAL_RAD = 0
 POWER_DECAY_THRESHOLD = 1e-4
-DEBUG_OUTPUT = False
 
 frequency = 1 / WAVELENGTH_UM
 polar_rad = np.linspace(0, 0.5 * math.pi, NUM_FARFIELD_PTS)
@@ -102,17 +103,18 @@ def radiation_pattern(e_field: np.ndarray, h_field: np.ndarray) -> np.ndarray:
 def get_farfields(
     sim: mp.Simulation, n2f_mon: mp.DftNear2Far
 ) -> Tuple[np.ndarray, np.ndarray]:
-    """Computes the far fields from the near fields.
+    """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
-        with six columns for the fields 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").
+        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)
@@ -191,12 +193,6 @@ def dipole_in_vacuum(
         )
     )
 
-    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")
-
     e_field, h_field = get_farfields(sim, nearfields_monitor)
 
     return e_field, h_field
@@ -250,15 +246,13 @@ def flux_from_farfields(e_field: np.ndarray, h_field: np.ndarray) -> float:
     m = 0
     while True:
         e_field, h_field = dipole_in_vacuum(args.dipole_pol, args.dipole_pos_r, m)
+        e_field_total += e_field * cmath.exp(1j * m * AZIMUTHAL_RAD)
+        h_field_total += h_field * cmath.exp(1j * m * AZIMUTHAL_RAD)
 
-        if args.dipole_pol == "x":
-            e_field_total[:, 0] += e_field[:, 0] * (1 if m == 0 else 2)
-            h_field_total[:, 1] += h_field[:, 1] * (1 if m == 0 else 2)
-            e_field_total[:, 2] += e_field[:, 2] * (1 if m == 0 else 2)
-        else:
-            h_field_total[:, 0] += h_field[:, 0] * (1 if m == 0 else 2)
-            e_field_total[:, 1] += e_field[:, 1] * (1 if m == 0 else 2)
-            h_field_total[:, 2] += h_field[:, 2] * (1 if m == 0 else 2)
+        if m > 0:
+            e_field, h_field = dipole_in_vacuum(args.dipole_pol, args.dipole_pos_r, -m)
+            e_field_total += e_field * cmath.exp(-1j * m * AZIMUTHAL_RAD)
+            h_field_total += h_field * cmath.exp(-1j * m * AZIMUTHAL_RAD)
 
         flux = flux_from_farfields(e_field, h_field)
         if flux > flux_max:
@@ -278,6 +272,7 @@ def flux_from_farfields(e_field: np.ndarray, h_field: np.ndarray) -> float:
     if mp.am_master():
         np.savez(
             "dipole_farfields_nonaxisymmetric.npz",
+            AZIMUTHAL_RAD=AZIMUTHAL_RAD,
             FARFIELD_RADIUS_UM=FARFIELD_RADIUS_UM,
             PML_UM=PML_UM,
             POWER_DECAY_THRESHOLD=POWER_DECAY_THRESHOLD,