Tutorial example for radiation pattern of a dipole using 1D calculation and Brillouin-zone integration

[oskooi]

Closes #2866.

This is still a work in progress. This implementation is based on the procedure described in #2866 (comment) for computing the radiation pattern of a dipole in vacuum using a 1D calculation and Brillouin-zone integration.

The radiation pattern of a dipole should be similar to a "donut" with rotational symmetry around the axis defined by the dipole orientation. However, as shown in the image below, the radiation pattern generated by this script is rotationally symmetric around the $z$ axis. This seems to be incorrect because the dipole is oriented in the $x$ direction.

For reference, the script used to generate the radiation pattern from the farfield data produced as the output from the script is provided below.

from typing import Tuple

import matplotlib.pyplot as plt
import numpy as np


def spherical_to_cartesian(
	polar_rad,
	azimuthal_rad
) -> Tuple[float, float, float]:
    """Converts a point in spherical to Cartesian coordinates."""

    rx = np.sin(polar_rad) * np.cos(azimuthal_rad)
    ry = np.sin(polar_rad) * np.sin(azimuthal_rad)
    rz = np.cos(polar_rad)

    return rx, ry, rz

if __name__ == "__main__":
    data = np.load("dipole_farfields.npz")

    farfield_radius_um = data['FARFIELD_RADIUS_UM']
    polar_rad = data['polar_rad']
    azimuthal_rad = data['azimuthal_rad']
    num_polar = data['NUM_POLAR']
    num_azimuthal = data['NUM_AZIMUTHAL']

    flux_z = np.zeros((num_polar, num_azimuthal))
    xs = np.zeros((num_polar, num_azimuthal))
    ys = np.zeros((num_polar, num_azimuthal))
    zs = np.zeros((num_polar, num_azimuthal))

    for i in range(num_polar):
	for j in range(num_azimuthal):

            flux_z[i, j] = np.real(
		np.conj(data['Ex'][i, j]) * data['Hy'][i, j] -
		np.conj(data['Ey'][i, j]) * data['Hx'][i, j]
            )

            xs[i, j], ys[i, j], zs[i, j] = spherical_to_cartesian(
		polar_rad[i], azimuthal_rad[j]
            )

            xs[i, j] *= flux_z[i, j]
            ys[i, j] *= flux_z[i, j]
            zs[i, j] *= flux_z[i, j]

    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("z flux (a.u.)", labelpad=30)
    ax.tick_params(axis="z", which="major", pad=15)
    ax.set(xticklabels=[], yticklabels=[])

    fig.savefig(
	"radiation_pattern_3d.png",
	dpi=150,
	bbox_inches="tight",
    )

[codecov-commenter]

⚠️ Please install the to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 73.71%. Comparing base (f29a8c7) to head (eef5dd1).
Report is 36 commits behind head on master.

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2917      +/-   ##
==========================================
- Coverage   73.81%   73.71%   -0.10%     
==========================================
  Files          18       18              
  Lines        5423     5449      +26     
==========================================
+ Hits         4003     4017      +14     
- Misses       1420     1432      +12     

---- 🚨 Try these New Features:

• Flaky Tests Detection - Detect and resolve failed and flaky tests

[stevengj]

You need to take into account $k_x, k_y$ in the near-to-far calculation, since the current sheets are oscillating in $(x,y)$. (Note that this will also give you nonzero $z$ components in the far field.)

[stevengj]

(Note that you might run into an issue with the 1d PML, because when you are right on the light line the fields won't decay.)

[oskooi]

You need to take into account $k_x, k_y$ in the near-to-far calculation, since the current sheets are oscillating in $(x,y)$. (Note that this will also give you nonzero $z$ components in the far field.)

Here are the equations for the S- and P-polarized farfields given a current sheet oscillating in $(x,y)$.

[oskooi]

After updating the script to include the corrected formulas for obtaining the farfields given a current sheet, the farfields given an $x$-directed dipole is not rotationally symmetric but is still not quite "donut" like.

[oskooi]

from typing import Tuple

import matplotlib.pyplot as plt
import numpy as np


def spherical_to_cartesian(
       polar_rad,
       azimuthal_rad
) -> Tuple[float, float, float]:
   """Converts a point in spherical to Cartesian coordinates."""

   rx = np.sin(polar_rad) * np.cos(azimuthal_rad)
   ry = np.sin(polar_rad) * np.sin(azimuthal_rad)
   rz = np.cos(polar_rad)

   return rx, ry, rz


if __name__ == "__main__":
   data = np.load("dipole_farfields.npz")

   farfield_radius_um = data['FARFIELD_RADIUS_UM']
   polar_rad = data['polar_rad']
   azimuthal_rad = data['azimuthal_rad']
   num_polar = data['NUM_POLAR']
   num_azimuthal = data['NUM_AZIMUTHAL']

   flux_x = np.zeros((num_polar, num_azimuthal))
   flux_y = np.zeros((num_polar, num_azimuthal))
   flux_z = np.zeros((num_polar, num_azimuthal))
   xs = np.zeros((num_polar, num_azimuthal))
   ys = np.zeros((num_polar, num_azimuthal))
   zs = np.zeros((num_polar, num_azimuthal))

   for i in range(num_polar):
       for j in range(num_azimuthal):

           flux_x[i, j] = np.real(
               np.conj(data['Ey'][i, j]) * data['Hz'][i, j] -
               np.conj(data['Ez'][i, j]) * data['Hy'][i, j]
           )

           flux_y[i, j] = np.real(
               np.conj(data['Ez'][i, j]) * data['Hx'][i, j] -
               np.conj(data['Ex'][i, j]) * data['Hz'][i, j]
           )

           flux_z[i, j] = np.real(
               np.conj(data['Ex'][i, j]) * data['Hy'][i, j] -
               np.conj(data['Ey'][i, j]) * data['Hx'][i, j]
           )

           flux_r = np.sqrt(flux_x[i, j]**2 + flux_y[i, j]**2 + flux_z[i, j]**2)

           rx, ry, rz = spherical_to_cartesian(polar_rad[i], azimuthal_rad[j])

           xs[i, j] = rx * flux_x[i, j]
           ys[i, j] = ry * flux_y[i, j]
           zs[i, j] = rz * flux_z[i, j]

   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("z flux (a.u.)", labelpad=30)
   ax.tick_params(axis="z", which="major", pad=15)
   ax.set(xticklabels=[], yticklabels=[])

   fig.savefig(
       "radiation_pattern_3d.png",
       dpi=150,
       bbox_inches="tight",
   )

[stevengj]

I think you want use

xs[i, j] = rx
ys[i, j] = ry
zs[i, j] = rz
flux_r[i,j] = rx * flux_x[i, j] + ry * flux_y[i, j] + rz * flux_z[i, j]

with the surface given by (xs,ys,zs) but the color given by flux_r. Here is an example in matplotlib — you have to use the facecolors argument to plot_surface. Something like:

c = matplotlib.cm.magma(flux_r)
plot_surface(xs,ys,zs, facecolors=c)

or just

pcolor(xs,ys, flux_r, cmap="magma")

[oskooi]

I fixed a couple of bugs in the calculation of the far fields. However, even with these fixes the radiation pattern for $\phi = 0$ (i.e., a slice of $\phi$ in the range $[0, 2\pi]$) is quite different than the expected $\cos^2(\theta)$. (A fixed $\phi$ is intended to simplify the analysis from a 3D contour which can be difficult to visualize correctly to a 2D polar plot.) The script used to generate the 2D radiation pattern (polar plot of the radial flux) using the simulation data is provided below. For reference, the radiation pattern for a dipole in 3D resembling a "donut" is shown as a separate figure below.

At this point, it might be worth investigating whether setting up the 3D simulation using a 1D grid itself is bug free (e.g., see #2916).

import math
from typing import Tuple

import matplotlib
import matplotlib.pyplot as plt
import numpy as np


if __name__ == "__main__":
    data = np.load("dipole_farfields.npz")

    polar_rad = data['polar_rad']
    num_polar = data['NUM_POLAR']
    flux_x = np.zeros(num_polar)
    flux_z = np.zeros(num_polar)

    for i in range(num_polar):
        flux_x[i] = np.real(
            np.conj(data['Ey'][i, 0]) * data['Hz'][i, 0] -
            np.conj(data['Ez'][i, 0]) * data['Hy'][i, 0]
	)

	flux_z[i] = np.real(
            np.conj(data['Ex'][i, 0]) * data['Hy'][i, 0] -
            np.conj(data['Ey'][i, 0]) * data['Hx'][i, 0]
        )

    flux_r = np.sqrt(np.square(flux_x) + np.square(flux_z))

    np.savez(
        "dipole_radiation_pattern.npz",
        flux_r=flux_r,
        polar_rad=polar_rad,
    )

    fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(6, 6))
    ax.plot(polar_rad, flux_r / np.max(flux_r), 'b-', label="Meep")
    ax.plot(polar_rad, np.square(np.cos(polar_rad)), 'r--', label="cos$^2$(θ)")
    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.legend()

    fig.savefig(
        "radiation_pattern_jx.png",
        dpi=150,
        bbox_inches="tight",
    )

[stevengj]

Suggested change

	ex0 = frequency * current_amplitude / (2 * kz)
	ex0 = kz * current_amplitude / (2 * frequency)

[stevengj]

shifting the index seems wrong. the point index hasn't changed in this process, you only rotated the fields and r.

[oskooi]

The results for the radiation pattern of a dipole in vacuum are improved after making these two modifications to the script:

1. Correcting the analytic formulas for the far fields from a sheet current. These errors were discovered rederiving the formulas.
2. Specifying a radius of 1 rather than 1e6 for the spherical surface for the far fields. It's not obvious why the radius of the far field surface should affet the results but it does so in a significant way.

The radiation pattern of the dipole computed using the Brillouin-zone integration is now similar to the analytic result. However, the results are not yet identical. The discrepancy should go away when using a finer and finer grid for the Brillouin-zone integration. This is something I am currently investigating.

[stevengj]

1. Probably for large R the far fields E(kx,ky) become very oscillatory functions of (kx,ky), so that your integral is not being computed accurately. Would be nice to plot E(kx,ky) for a few different R to confirm this.
2. For small R, the angle is not being computed correctly — the angle should be from the location of your dipole source to the far-field point, not from the near-field monitor to the far-field point.

[oskooi]

1. Probably for large R the far fields E(kx,ky) become very oscillatory functions of (kx,ky), so that your integral is not being computed accurately. Would be nice to plot E(kx,ky) for a few different R to confirm this.

Increasing $R$ does indeed lead to increasing oscillations in the integrand. The figure below shows the real part of the far field $E_x$ at a fixed position on the hemispherical surface over the ($k_x$, $k_y$) grid for two different $R$: 1.0 vs. 1e6.

This explains why $R$ had to be made small (< $\lambda$ = 1 μm) in order to be able generate a smooth radiation profile.

[stevengj]

There are numerical methods (Levin, 1994) to integrate fast-oscillating integrals like this accurately and efficiently, though I can't find any implementation in Python. However, in this case we really only care about the $r \to \infty$ limit, which suggests that we really should be able to get an exact asymptotic expression (e.g. by doing a saddle-point integration and taking the limit).

I suspect that the result in the limit will be proportional to just integrating the powers from each $k$ separately, though I don't know what the proportionality constant will be.

[oskooi]

I suspect that the result in the limit will be proportional to just integrating the powers from each $k$ separately, though I don't know what the proportionality constant will be.

After reorganizing the calculation to compute the radial Poynting flux $P_r = \sum_k P_k$ via $P_k = E_k^* \times H_k$ rather than $P_r = E_{total}^* \times H_{total}$ where $E_{total} = \sum_k E_k$ and $H_{total} = \sum_k H_k$ fields at each grid point $(\theta, \phi)$ in the far field via Brillouin-zone integration, the resulting radiation profile at all $\phi$ grid points matches $\cos(\theta)$ rather than the expected $\cos^2(\theta)$. The radiation profile at $\phi = 0$ is shown below. Note that $E_k = E^0_k e^{ikr}$ where $E^0_k$ is the amplitude obtained analytically and $r$ is the far-field radius.

The key feature of computing $P_r = \sum_k P_k$ as we discussed is that it does not involve choosing an (arbitrary) radius $r$ of the hemispherical surface for the far-field grid because it is effectively based on the $r \to \infty$ limit. This means that $\theta_1 = \theta_2$ in the schematic below. In the calculation $\theta_1$ is computed over the range [0, $\pi$ / 2]. This is the same range for $\theta$ in the plot above.

For reference, here is the plotting script used to generate the figure above using the output of the simulation data (dipole_farfields.npz).

import math
from typing import Tuple

import matplotlib.pyplot as plt
import numpy as np

AZIMUTHAL_INDEX = 0

if __name__ == "__main__":
    data = np.load("dipole_farfields.npz")
    polar_rad = data["polar_rad"]
    poynting_flux = data["poynting_flux"][:, AZIMUTHAL_INDEX]

    fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(6, 6))
    ax.plot(polar_rad, poynting_flux / np.max(poynting_flux), 'b-', label="Meep")
    ax.plot(polar_rad, np.cos(polar_rad), 'r--', label="cos(θ)")
    ax.plot(polar_rad, np.square(np.cos(polar_rad)), 'g--', label="cos$^2$(θ)")
    ax.set_theta_direction(-1) # Theta increases in the clockwise direction                           
    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.legend()

    fig.savefig(
        "poynting_flux_jx.png",
	dpi=150,
        bbox_inches="tight",
    )

[stevengj]

Note that for the small-R calculation you could just take the far-field point (orange), subtract the source point (red) in Cartesian coordinates, then convert to spherical to get $\theta_2$.

[mochen4]

Denote the angle of the orange dot $\theta_3$, then you have $\theta_1 = \theta_2+ \theta_3 $, and by law of sines you have $\frac{sin \theta_2 }{R} = \frac{sin \theta_3}{L}$. Then you have a formula relating $\theta_1$ from $\theta_2$ :

$$\theta_1 = \theta_2 + arcsin(L sin\theta_2/R)$$

To get $$\theta_2$$ from $$\theta_1$$, the easiest approach might be using the law of cosines. Denote the distance between orange and red dot $S$, then by the law of cosines $S^2 = L^2 + R^2-2LRcos(\pi-\theta_1) =L^2 + R^2+2LRcos\theta_1$, then again by law of sines: $\frac{sin(\pi-\theta_1)}{S}=\frac{sin\theta_2}{R}$. Note that $sin(\pi-\theta_1)=sin\theta_1$, so

$$\theta_2 = arcsin \frac{R\sin\theta_1}{\sqrt{L^2 + R^2+2LRcos\theta_1}}$$

As a sanity check, when $L=R$, $\theta_1= 2\theta_2$; when $R &gt;&gt; L$, $\theta_1 \approx \theta_2$.

[oskooi]

Here is a plot of the analytic expression for $\theta_2(L, R, \theta_1)$ vs. $\theta_1$ for a fixed $L$ and different values of $R$.

To generate a range of [0, π/2] for $\theta_2$ (for comparing the Meep data to the expected result of a dipole) requires making $R$ larger. This is a problem because larger $R$ means a faster oscillating integrand. Unfortunately this means that it may not be practical to use small values for $R$ for this calculation.

[mochen4]

To generate a range of [0, π/2] for θ 2 (for comparing the Meep data to the expected result of a dipole) requires making R larger.

For $\theta_2$ to attain π/2, we must have R > L. Alternatively, for R < L, instead of comparing the full range of data, maybe we can just compare the data between [0, $\theta_2^{max}$], where $\theta_2^{max} = arcsin(\frac{R}{L})$.