⚠️ 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.

@@            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     

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

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

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

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

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.

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

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",
    )

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.

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.

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.

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.

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",
    )

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

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

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.

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

There seems to be an inconsistency when calculating the radiation profile of an $E_r$ dipole in vacuum in cylindrical coordinates and comparing with the expected result of $\cos^2(\theta)$ (i.e. "donut") used to debug this PR which is based on a "3D" calculation. The radiation profile calculated in cylindrical coordinates for an $E_r$ (or $E_\phi$) dipole at $r = 0$ and $m=-1$ is equivalent to $(1 + \cos^2\theta)/2$ which is shown in the plot below. The same radiation profile is computed at any $r$ using a Fourier-series expansion. The expression $(1 + \cos^2\theta)/2$ is what we used in Tutorial/Extraction Efficiency of a Collection of Dipoles in a Disc to empirically determine the weighting factor $s(r_n)$ for a series of $N$ dipoles at position $r_n$ along a line in the $r$ direction.

Why is the radiation profile of a dipole different when computing the result in cylindrical vs. Cartesian coordinates: $(1 + \cos^2\theta)/2$ vs. $\cos^2(\theta)$ ? The radiation pattern should be independent of the coordinate system. See also the related discussion in #2656.

Simulation layout of a a dipole at $r = 0$ (red dot) in vacuum in cylindrical coordinates. The blue line denotes the near-field monitor used to enclose the emitter.