Skip to content

Commit

Permalink
add example of a nonaxisymmetric dipole with Ey polarization
Browse files Browse the repository at this point in the history
  • Loading branch information
@oskooi
oskooi committed Dec 29, 2024
1 parent 2f8b532 commit 8ba2157
Show file tree
Hide file tree
Showing 5 changed files with 75 additions and 41 deletions.
8 changes: 4 additions & 4 deletions doc/docs/Python_Tutorials/Near_to_Far_Field_Spectra.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
----------------------------------------------------------

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).
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$.

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 (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$.

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$.
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$.

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 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).

Expand Down
Binary file modified doc/docs/images/radiation_pattern_axisymmetric_dipole.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified doc/docs/images/radiation_pattern_nonaxisymmetric_dipole.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
42 changes: 23 additions & 19 deletions python/examples/dipole_in_vacuum_cyl_axisymmetric.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
polar_rad = np.linspace(0, 0.5 * math.pi, NUM_FARFIELD_PTS)


def plot_radiation_pattern(dipole_pol: int, radial_flux: np.ndarray):
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")
Expand All @@ -38,11 +38,11 @@ def plot_radiation_pattern(dipole_pol: int, radial_flux: np.ndarray):
"""
normalized_radial_flux = radial_flux / np.max(radial_flux)

if dipole_pol == mp.X:
if dipole_pol == "x":
dipole_radial_flux = np.square(np.cos(polar_rad))
dipole_radial_flux_label = r"$\cos^2θ$"
dipole_name = "$E_x$"
elif dipole_pol == mp.Z:
else:
dipole_radial_flux = np.square(np.sin(polar_rad))
dipole_radial_flux_label = r"$\sin^2θ$"
dipole_name = "$E_z$"
Expand All @@ -59,7 +59,9 @@ 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 (φ = 0) of an axisymmetric " f"{dipole_name} dipole"
)

if mp.am_master():
fig.savefig(
Expand Down Expand Up @@ -87,8 +89,12 @@ 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] * 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
Expand Down Expand Up @@ -120,13 +126,13 @@ def get_farfields(
FARFIELD_RADIUS_UM * math.cos(polar_rad[n]),
),
)
e_field[n, :] = [np.conj(far_field[j]) for j in range(3)]
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: int, m: int) -> Tuple[np.ndarray, np.ndarray]:
def dipole_in_vacuum(dipole_pol: str, m: int) -> Tuple[np.ndarray, np.ndarray]:
"""Computes the far fields of an axisymmetric point source.
Args:
Expand All @@ -142,7 +148,7 @@ def dipole_in_vacuum(dipole_pol: int, m: int) -> Tuple[np.ndarray, np.ndarray]:

boundary_layers = [mp.PML(thickness=PML_UM)]

if dipole_pol == mp.X:
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
Expand All @@ -159,7 +165,7 @@ def dipole_in_vacuum(dipole_pol: int, m: int) -> Tuple[np.ndarray, np.ndarray]:
amplitude=1j if m == 1 else -1j,
),
]
elif dipole_pol == mp.Z:
else:
dipole_pos_r = 0
sources = [
mp.Source(
Expand Down Expand Up @@ -197,7 +203,7 @@ def dipole_in_vacuum(dipole_pol: int, m: int) -> Tuple[np.ndarray, np.ndarray]:
sim.run(
until_after_sources=mp.stop_when_fields_decayed(
20.0,
mp.Er if dipole_pol == mp.X else mp.Ez,
mp.Er if dipole_pol == "x" else mp.Ez,
mp.Vector3(dipole_pos_r, 0, 0),
1e-6,
)
Expand Down Expand Up @@ -251,27 +257,25 @@ def flux_from_farfields(e_field: np.ndarray, h_field: np.ndarray) -> float:
h_field_total = np.zeros((NUM_FARFIELD_PTS, 3), dtype=np.complex128)

if args.dipole_pol == "x":
# A linearly polarized dipole can be formed from the superposition
# An x-polarized dipole can be formed from the superposition
# of two left- and right-circularly polarized dipoles.

e_field, h_field = dipole_in_vacuum(mp.X, +1)
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, h_field = dipole_in_vacuum(mp.X, -1)
e_field, h_field = dipole_in_vacuum("x", -1)
e_field_total += 0.5 * e_field
h_field_total += 0.5 * h_field

elif args.dipole_pol == "z":
e_field, h_field = dipole_in_vacuum(mp.Z, 0)
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(
mp.X if args.dipole_pol == "x" else mp.Z, dipole_radiation_pattern_scaled
)
plot_radiation_pattern(args.dipole_pol, dipole_radiation_pattern_scaled)

if mp.am_master():
np.savez(
Expand Down
66 changes: 48 additions & 18 deletions python/examples/dipole_in_vacuum_cyl_nonaxisymmetric.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,21 +26,29 @@
polar_rad = np.linspace(0, 0.5 * math.pi, NUM_FARFIELD_PTS)


def plot_radiation_pattern(radial_flux: np.ndarray):
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 dipole polarization.
radial_flux: the radial flux in polar coordinates.
"""
normalized_radial_flux = radial_flux / np.max(radial_flux)
dipole_radial_flux = np.square(np.cos(polar_rad))
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.ones(NUM_FARFIELD_PTS)
dipole_radial_flux_label = "constant (1.0)"
dipole_name = "$E_y$"

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=r"$\cos^2θ$")
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)
Expand All @@ -50,7 +58,9 @@ def plot_radiation_pattern(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 a nonaxisymmetric $E_x$ dipole")
ax.set_title(
"radiation pattern (φ = 0) of a nonaxisymmetric " f"{dipole_name} dipole"
)

if mp.am_master():
fig.savefig(
Expand Down Expand Up @@ -78,8 +88,12 @@ 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] * 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
Expand Down Expand Up @@ -111,16 +125,19 @@ def get_farfields(
FARFIELD_RADIUS_UM * math.cos(polar_rad[n]),
),
)
e_field[n, :] = [np.conj(far_field[j]) for j in range(3)]
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_pos_r: mp.Vector3, m: int) -> Tuple[np.ndarray, np.ndarray]:
def dipole_in_vacuum(
dipole_pol: str, dipole_pos_r: mp.Vector3, m: int
) -> Tuple[np.ndarray, np.ndarray]:
"""Computes the far fields of a nonaxisymmetric point source.
Args:
dipole_pol: the dipole polarization.
dipole_pos_r: the radial position of the dipole.
m: angular φ dependence of the fields exp(imφ).
Expand All @@ -133,10 +150,12 @@ def dipole_in_vacuum(dipole_pos_r: mp.Vector3, m: int) -> Tuple[np.ndarray, np.n

boundary_layers = [mp.PML(thickness=PML_UM)]

src_cmpt = mp.Er if dipole_pol == "x" else mp.Ep

sources = [
mp.Source(
src=mp.GaussianSource(frequency, fwidth=0.1 * frequency),
component=mp.Er,
component=src_cmpt,
center=mp.Vector3(dipole_pos_r, 0, 0),
)
]
Expand Down Expand Up @@ -168,7 +187,7 @@ def dipole_in_vacuum(dipole_pos_r: mp.Vector3, m: int) -> Tuple[np.ndarray, np.n

sim.run(
until_after_sources=mp.stop_when_fields_decayed(
20.0, mp.Er, mp.Vector3(dipole_pos_r, 0, 0), 1e-6
20.0, src_cmpt, mp.Vector3(dipole_pos_r, 0, 0), 1e-6
)
)

Expand Down Expand Up @@ -208,11 +227,20 @@ def flux_from_farfields(e_field: np.ndarray, h_field: np.ndarray) -> float:

if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument(
"dipole_pol",
type=str,
choices=["x", "y"],
help="polarization of the electric dipole (x or y)",
)
parser.add_argument(
"dipole_pos_r", type=float, help="radial position of the dipole"
)
args = parser.parse_args()

if args.dipole_pos_r == 0:
raise ValueError(f"dipole_pos_r must be nonzero.")

# Fourier series expansion of the fields from a ring current source
# used to generate a point dipole localized in the azimuthal direction.

Expand All @@ -221,15 +249,16 @@ def flux_from_farfields(e_field: np.ndarray, h_field: np.ndarray) -> float:
flux_max = 0
m = 0
while True:
e_field, h_field = dipole_in_vacuum(args.dipole_pos_r, m)
e_field, h_field = dipole_in_vacuum(args.dipole_pol, args.dipole_pos_r, m)

if m == 0:
e_field_total += e_field
h_field_total += h_field
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:
e_field_total[:, 0] += 2 * e_field[:, 0]
h_field_total[:, 1] += 2 * h_field[:, 1]
e_field_total[:, 2] += 2 * e_field[:, 2]
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)

flux = flux_from_farfields(e_field, h_field)
if flux > flux_max:
Expand All @@ -244,7 +273,7 @@ def flux_from_farfields(e_field: np.ndarray, h_field: np.ndarray) -> float:

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(dipole_radiation_pattern_scaled)
plot_radiation_pattern(args.dipole_pol, dipole_radiation_pattern_scaled)

if mp.am_master():
np.savez(
Expand All @@ -254,6 +283,7 @@ def flux_from_farfields(e_field: np.ndarray, h_field: np.ndarray) -> float:
POWER_DECAY_THRESHOLD=POWER_DECAY_THRESHOLD,
RESOLUTION_UM=RESOLUTION_UM,
WAVELENGTH_UM=WAVELENGTH_UM,
dipole_pol=args.dipole_pol,
dipole_pos_r=args.dipole_pos_r,
dipole_radiation_pattern=dipole_radiation_pattern,
e_field_total=e_field_total,
Expand Down

0 comments on commit 8ba2157

Please sign in to comment.