Tutorial for extraction efficiency of a collection of dipoles in a disc in cylindrical coordinates

[oskooi]

Closes #2690.

This script seems to be working and generates the radiation pattern below. The extraction efficiency computed from this radiation pattern is 0.9569.

[stevengj]

The fact that your vacuum radiation patterns are not independent of the location of the source looks wrong.

Not sure what the bug is — maybe you should have summed E and H over m before computing the E*xH power?

[oskooi]

Here are three radiation patterns scaled by $s(r) = 1/(2(r/r_0)^2$ for an $E_r$ dipole source at $\lambda = 0.45$ μm in vacuum at three different $r$ positions: $\approx$ 0, 0.5 μm, and 1.0 μm. Also included is the radiation pattern for a circularly polarized electric dipole (consisting of $Er + jE{\phi})$ obtained from antenna theory: $cos^2(\theta) + 1$.

There is good agreement in the radiation pattern between the theoretical result and the simulation for the dipole at $\approx$ 0. However, there are still noticeable differences in the radiation patterns for dipoles at $r > 0$ even after (1) increasing the resolution to 200 pixels/μm, (2) reducing the threshold flux used to truncate the Fourier-series expansion (i.e., including more terms), and (3) reducing the field tolerance to 1e-9 in the termination criteria meep.stop_when_fields_decayed. I also tried increasing the PML thickness and the size of the cell to no avail.

Schematic of the simulation cell.

[stevengj]

Would be good to overlay the analytical formula (for a circularly polarized source).

(Given the formula for the field of a dipole, add it to its rotation by 90° multiplied by i to get the field of a circularly polarized dipole. Then plug in the Poynting flux. Someone has probably worked it out already.)

[stevengj]

A good check is that summing the fields over m and then computing the Poynting flux ExH* should give the same answer as computing the Poynting flux and then summing over m.

[oskooi]

A good check is that summing the fields over m and then computing the Poynting flux ExH* should give the same answer as computing the Poynting flux and then summing over m.

The procedure for summing the fields first and then computing the Poynting flux is described in Tutorial/Nonaxisymmetric Dipoles Sources:

Note that in Step 2, the $|m| > 0$ terms in the expansion involve taking the real part of the fields. Also, the fields $E{tot}$ and $H{tot}$ and flux $P(\theta, \phi)$ are always computed at $\phi = 0$.

Using this approach, the radiation pattern of the dipoles at $r > 0$ is significantly different than the alternative approach (shown previously) of summing the Poynting flux of the individual $m$ calculations:

The same simulation parameters (resolution, field decay threshold for termination, and flux threshold for truncating the Fourier-series expansion) as before were used. For the $r = 0.5 \mu m$ case, there are $M + 1 = 14$ terms in the expansion, and for $r = 1.0 \mu m$ there are $M + 1 = 23$ terms.

[oskooi]

There is a subtlety in the calculation of the radiation pattern $P(\theta)$ based on summing the results for the Poynting flux $P{r,m}$ and $P{z,m}$ using the Fourier-series expansion over $m$. There are two ways to obtain $P(\theta)$ as shown below. The plots of the radiation pattern shown above were obtained using the method on the left. For comparison, I also tried the method on the right. The results were identical. This is expected because of the orthogonality in the flux for $P_r$, $P_z$, and $P(\theta)$ for different $m$.

[oskooi]

Would be good to overlay the analytical formula (for a circularly polarized source).

I have added the radiation pattern for a circularly polarized electric dipole obtained using antenna theory to the plot in the previous comment. There is good agreement between the theoretical result and the simulated radiation pattern for the dipole at $r_0 \approx 0$. This confirms that the simulation is set up correctly.

[oskooi]

edit (1/8/2024): disregard. Using accurate_fields_near_cyl_origin=True and smaller Courant had no effect on the results.

It seems that to obtain agreement in the radiation pattern in vacuum for dipoles at $r > 0$ it is necessary to set accurate_fields_near_cyl_origin to True when performing the Fourier-series expansion for terms with $|m| > 1$. Also, as described in #2644, the Courant factor needs to be made smaller with $m$.

These two modifications are:

    if abs(m) > 1:
        accurate_fields_near_cylorigin = True
        Courant = 1 / (abs(m) + 1)
    else:
        accurate_fields_near_cylorigin = False
        Courant = 0.5

With these changes, the agreement in the radiation pattern is noticeably improved compared to the original result.

[oskooi]

For debugging purposes, here is the script I used to generate the radiation patterns (shown above) for a dipole in vacuum at $r > 0$ using a Fourier-series expansion of the fields via $e^{im\phi}$.

"""Radiation pattern of a dipole in vacuum in cylindrical coordinates."""

import math

import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
import meep as mp
import numpy as np

RESOLUTION_UM = 100
WAVELENGTH_UM = 1.0
NUM_FARFIELD_PTS = 200
FARFIELD_RADIUS = 1e6 * WAVELENGTH_UM
DEBUG_OUTPUT = False

farfield_angles = np.linspace(0, 0.5 * math.pi, NUM_FARFIELD_PTS)

def plot_radiation_pattern_polar(radiation_pattern: np.ndarray):
    """Plots the radiation pattern in polar coordinates.                                                          

    The angles increase clockwise with zero at the pole (+z direction)                                            
    and π/2 at the equator (+r direction).                                                                        

    Args:                                                                                                         
        radiation_pattern: radial flux of the far fields in polar coordinates.                                    
    """
    radiation_pattern_theory = ((np.cos(farfield_angles)**2 + 1) *
                radiation_pattern[0] * 0.5)

    fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(6,6))
    ax.plot(
        farfield_angles,
        radiation_pattern,
        "b-",
        label="Meep"
    )
    ax.plot(
        farfield_angles,
        radiation_pattern_theory,
        "r-",
        label="theory"
    )
    ax.set_theta_direction(-1)
    ax.set_theta_offset(0.5 * math.pi)
    ax.set_thetalim(0, 0.5 * math.pi)
    ax.grid(True)
    ax.set_rlabel_position(22)
    ax.set_ylabel("radial flux (a.u.)")
    ax.set_title("radiation pattern in polar coordinates")
    ax.legend()

    if mp.am_master():
    fig.savefig(
            "vacuum_cyl_radpattern.png",
            dpi=150,
            bbox_inches="tight",
    )

def radiation_pattern_flux(radiation_pattern: np.ndarray) -> float:
    """Computes the Poynting flux from the radiation pattern."""

    dtheta = 0.5 * math.pi / (NUM_FARFIELD_PTS - 1)
    dphi = 2 * math.pi
    flux_far = (np.sum(radiation_pattern * np.sin(farfield_angles)) *
                FARFIELD_RADIUS * FARFIELD_RADIUS * dtheta * dphi)

    # We need to multiply the result by two because this integral is only                                         
    # for θ in [0, 90°] whereas we need it to be [0, 180°].                                                       
    flux_far *= 2

    return flux_far

def radiation_pattern_from_farfields(
        sim: mp.Simulation, n2f_mon: mp.DftNear2Far
) -> np.ndarray:
    """Computes the radiation pattern from the far fields.                                                        

    Args:                                                                                                         
        sim: a `Simulation` object.                                                                               
        n2f_mon: a `DftNear2Far` object returned by `Simulation.add_near2far`.                                    

    Returns:                                                                                                      
        Array of radial Poynting flux, one for each point on the circumference                                    
        of a quarter circle with angular range of [0, π/2]. 0 rad is the +z                                       
        direction and π/2 is +r.                                                                                  
    """
    e_field = np.zeros((NUM_FARFIELD_PTS, 3), dtype=np.complex128)
    h_field = np.zeros((NUM_FARFIELD_PTS, 3), dtype=np.complex128)
    for n in range(NUM_FARFIELD_PTS):
        ff = sim.get_farfield(
            n2f_mon,
            mp.Vector3(
                FARFIELD_RADIUS * math.sin(farfield_angles[n]),
                0,
                FARFIELD_RADIUS * math.cos(farfield_angles[n])
            )
        )
        e_field[n, :] = [ff[j] for j in range(3)]
        h_field[n, :] = [ff[j + 3] for j in range(3)]

    x_flux = (np.real(np.conj(e_field[:, 1]) * h_field[:, 2] -
                      np.conj(e_field[:, 2]) * h_field[:, 1]))
    z_flux = (np.real(np.conj(e_field[:, 0]) * h_field[:, 1] -
                      np.conj(e_field[:, 1]) * h_field[:, 0]))
    r_flux = np.sqrt(np.square(x_flux) + np.square(z_flux))

    return r_flux

def point_dipole_in_vacuum(rpos_um: float, m: int) -> np.ndarray:
    pml_um = 1.0
    padding_um = 10.0
    size_r = padding_um + pml_um
    size_z = pml_um + padding_um + pml_um
    cell_size = mp.Vector3(size_r, 0, size_z)
    pml_layers = [mp.PML(thickness=pml_um)]

    frequency = 1 / WAVELENGTH_UM

    src_pt = mp.Vector3(rpos_um, 0, 0)
    sources = [
        mp.Source(
            src=mp.GaussianSource(
                frequency, fwidth=0.2 * frequency, cutoff=10.0
            ),
            component=mp.Er,
            center=src_pt,
        )
    ]

    sim = mp.Simulation(
        resolution=RESOLUTION_UM,
        cell_size=cell_size,
        dimensions=mp.CYLINDRICAL,
        sources=sources,
        boundary_layers=pml_layers,
        m=m,
        force_complex_fields=True,
    )

    flux_mon = sim.add_flux(
        frequency,
        0,
        1,
        mp.FluxRegion(
            center=mp.Vector3(
                0.5 * (size_r - pml_um), 0, 0.5 * size_z - pml_um
            ),
            size=mp.Vector3(size_r - pml_um, 0, 0)
        ),
        mp.FluxRegion(
            center=mp.Vector3(size_r - pml_um, 0, 0),
            size=mp.Vector3(0, 0, size_z - 2 * pml_um),
        ),
        mp.FluxRegion(
            center=mp.Vector3(
                0.5 * (size_r - pml_um), 0, -0.5 * size_z + pml_um
            ),
            size=mp.Vector3(size_r - pml_um, 0, 0),
            weight=-1.0,
        ),
    )

    n2f_mon = sim.add_near2far(
        frequency,
        0,
        1,
        mp.Near2FarRegion(
            center=mp.Vector3(
                0.5 * (size_r - pml_um), 0, 0.5 * size_z - pml_um
            ),
            size=mp.Vector3(size_r - pml_um, 0, 0)
        ),
        mp.Near2FarRegion(
            center=mp.Vector3(size_r - pml_um, 0, 0),
            size=mp.Vector3(0, 0, size_z - 2 * pml_um),
        ),
        mp.Near2FarRegion(
            center=mp.Vector3(
                0.5 * (size_r - pml_um), 0, -0.5 * size_z + pml_um
            ),
            size=mp.Vector3(size_r - pml_um, 0, 0),
            weight=-1.0,
        ),
    )

    if DEBUG_OUTPUT:
        fig, ax = plt.subplots()
        sim.plot2D(ax=ax)
        if mp.am_master():
            fig.savefig(
                "cyl_simulation_layout.png", dpi=150, bbox_inches="tight"
            )

    sim.run(
        until_after_sources=mp.stop_when_fields_decayed(
            25.0, mp.Er, src_pt, 1e-8
        )
    )

    flux_near = mp.get_fluxes(flux_mon)[0]

    radiation_pattern = radiation_pattern_from_farfields(sim, n2f_mon)
    flux_far = radiation_pattern_flux(radiation_pattern)

    err = abs(flux_far - flux_near) / flux_near

    print(
        f"flux-m:, {rpos_um}, {m}, "
        f"{flux_near:.6f} (near), {flux_far:.6f} (far), {err:.6f} (error)"
    )

    return radiation_pattern

if __name__ == "__main__":
    # An Er source at r = 0 must be slightly offset due to #2704.                                                 
    # Requires a single simulation with m = ±1.                                                                   
    rpos_um = 1.5 / RESOLUTION_UM
    m = 1
    radiation_pattern = point_dipole_in_vacuum(rpos_um, m)

    plot_radiation_pattern_polar(
        FARFIELD_RADIUS**2 * radiation_pattern
    )

    if mp.am_master():
        fname = f'vacuum_cyl_rpos{rpos_um}_res{RESOLUTION_UM}.npz'

        np.savez(
            fname,
            m=m,
            resolution=RESOLUTION_UM,
            wavelength_um=WAVELENGTH_UM,
            num_farfield_pts=NUM_FARFIELD_PTS,
            farfield_angles=farfield_angles,
            farfield_radius=FARFIELD_RADIUS,
            radiation_pattern=radiation_pattern
        )

    # An Er source at r > 0 requires Fourier-series expansion of exp(imϕ).                                        
    rpos_um = 0.19
    m = 0
    flux_max = 0
    flux_thresh = 1e-2
    radiation_patterns = np.zeros(NUM_FARFIELD_PTS)
    radiation_pattern_total = np.zeros(NUM_FARFIELD_PTS)
    while True:
        radiation_pattern = point_dipole_in_vacuum(rpos_um, m)
        radiation_pattern_total += radiation_pattern * (1 if m == 0 else 2)

        if m == 0:
            radiation_patterns = radiation_pattern[:, None]
        else:
            radiation_patterns = np.hstack(
                (radiation_patterns, radiation_pattern[:, None])
            )

        flux_far = radiation_pattern_flux(radiation_pattern)
        if flux_far > flux_max:
            flux_max = flux_far

        if m > 0 and (flux_far / flux_max) < flux_thresh:
            break
        else:
            m += 1

    plot_radiation_pattern_polar(
        FARFIELD_RADIUS**2 * radiation_pattern_total
    )

    if mp.am_master():
        fname = f'vacuum_cyl_rpos{rpos_um}_res{RESOLUTION_UM}.npz'

        np.savez(
            fname,
            num_m=m,
            resolution=RESOLUTION_UM,
            wavelength_um=WAVELENGTH_UM,
            num_farfield_pts=NUM_FARFIELD_PTS,
            farfield_angles=farfield_angles,
            farfield_radius=FARFIELD_RADIUS,
            radiation_patterns=radiation_patterns
        )

[oskooi]

I figured out why the radiation pattern for an $E_r$ point-dipole current source in vacuum computed at various positions for $r > 0$ is not equivalent to the analytic result. It has to do with the way we are combining the results from multiple independent simulations using the Fourier-series expansion of the fields with $e^{im\phi}$.

Specifically, in Tutorial/Nonaxisymmetric Dipole Sources, the Fourier-series expansion is described as:

This approach produces incorrect results as shown above in the plots of the radiation pattern at various source positions. In order to produce the correct result which matches the radiation pattern from antenna theory, the Fourier-series expansion must be changed to:

Note: the $m=0$ term has twice the weight of the $m \neq 0$ terms. This is the exact opposite of what we have been using.

This requires changing a single line in the script above:

radiation_pattern_total += radiation_pattern * (1 if m == 0 else 2)

to

radiation_pattern_total += radiation_pattern * (2 if m == 0 else 1)

With this change, the radiation pattern at e.g. $r = 0.35$ μm matches the analytic result as expected.

[oskooi]

With the correction described in the previous comment to the way the Fourier-series expansion for $e^{im\phi}$ is computed, it is now possible to demonstrate that the radiation pattern of a dipole in vacuum at different positions with $r > 0$ is equivalent. This involves using the scaling factor $s(r) = (r_0 / r)^2$ for a dipole at position $r$ where $r_0$ is the dipole location at the origin (= 1.5 / resolution).

The simulation script as well as the description in the text of the tutorial have been updated with these changes.

[stevengj]

I think the issue is that, as an optimization, Meep only simulates the real parts of the field for m=0 (where the real and imaginary parts are decoupled), which effectively halves the current source for the m=0. To compensate, you could either use force_complex_fields=True for the m=0 simulation or multiply the m=0 power by a factor of 4. (This corresponds to 2x the empirical formula you wrote above.)

[stevengj]

For the emission at 90°, you might want to make sure that greencyl is doing the angular integral accurately. e.g. force it to use more angular points by setting N0 = 128 at https://github.com/NanoComp/meep/blob/4f67cb69592e7899877cabd3dc15dce27ee33f95/src/near2far.cpp#L291C5-L291C20

[stevengj]

If you want to validate this, an ideal tool here would be SCUFF, which should be very efficient for this type of problem.

[stevengj]

Note that a dielectric disk will have very high Q resonance, especially at higher frequencies. This might cause the fields to decay very slowly if you have a broadband source. If these aren't the frequencies of interest, then you might want to either (a) not excite the resonances by using a narrow-band source or (b) use stop_when_dft_decayed (or whatever it's called) to only look at the frequency components of interest when stopping the simulation. Or both.

[oskooi]

For the emission at 90°, you might want to make sure that greencyl is doing the angular integral accurately. e.g. force it to use more angular points by setting N0 = 128

Increasing N0 from 4 (currently) to 128 considerably reduces the discretization artifacts in the radiation pattern. The figure below shows these artifacts near $\theta \approx 90^\circ$ for N0 = 4 (blue line) but not for N0 = 128 (red) which is much smoother and matches the theoretical result at all angles.

Based on these results, it might therefore be useful to make N0 a user parameter in the routines get_farfield and get_farfields rather than a hard-coded internal value.

[oskooi]

I think the issue is that, as an optimization, Meep only simulates the real parts of the field for m=0 (where the real and imaginary parts are decoupled), which effectively halves the current source for the m=0. To compensate, you could either use force_complex_fields=True for the m=0 simulation or multiply the m=0 power by a factor of 4. (This corresponds to 2x the empirical formula you wrote above.)

As suggested, setting force_complex_fields=True in the Simulation constructor for the m = 0 run fixes the issue. With this change, using the existing expression for the Fourier-series expansion involving a weight of 1 for the results of the m = 0 run and 2 for the m ≠ 0 runs produces identical radiation patterns for dipoles at $r > 0$ in vacuum. Obtaining this agreement requires using the scaling factor $s(r) = 0.5(r_0/r)^2$.

[codecov-commenter]

Codecov Report

All modified and coverable lines are covered by tests :white_check_mark:

Comparison is base (4f67cb6) 74.06% compared to head (098ff7a) 73.78%. Report is 1 commits behind head on master.

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

Additional details and impacted files

```diff @@ Coverage Diff @@ ## master #2726 +/- ## ========================================== - Coverage 74.06% 73.78% -0.29% ========================================== Files 18 18 Lines 5399 5421 +22 ========================================== + Hits 3999 4000 +1 - Misses 1400 1421 +21 ``` [see 1 file with indirect coverage changes](https://app.codecov.io/gh/NanoComp/meep/pull/2726/indirect-changes?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=None)