See https://github.com/NanoComp/meep/issues/2489

freespace_reflector_cyl_debug.py

import argparse
from typing import Tuple

import matplotlib
matplotlib.use("agg")
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import meep as mp
import numpy as np

parser = argparse.ArgumentParser()
parser.add_argument('-res',
                    type=int,
                    default=50,
                    help='resolution (default: 50 pixels/um)')
parser.add_argument('-m',
                    type=int,
                    default=1,
                    help='angular dependence of fields exp(imφ) (default: 1)')
parser.add_argument('-dair',
                    type=float,
                    default=1.0,
                    help='thickness of air padding (default: 1.0 um)')
parser.add_argument('-dpml',
                    type=float,
                    default=1.0,
                    help='PML thickness (default: 1.0 um)')
args = parser.parse_args()


L = 6.0          # length of cell in r (excluding PML)                          
wvl = 1.0        # wavelength (in vacuum)                                       
fcen = 1 / wvl   # center frequency of source/monitor                           

# source properties                                                             
df = 0.1 * fcen
src = mp.GaussianSource(fcen, fwidth=df)

def radiated_fields(h: float, rpos: float, m: int) -> Tuple[np.ndarray,
                                                            np.ndarray]:
    """Computes the radiated fields of a point dipole above a lossless          
       ground plane in cylindrical coordinates using two different methods:     
       (1) a DFT monitor and (2) a near-to-far field transformation.            
                                                                                
    Args:                                                                       
        dmat: thickness of dielectric layer.                                    
        h: height of dipole above ground plane.                                 
        rpos: position of dipole in r direction.                                
        m: angular φ dependence of the fields exp(imφ).                         
    """
    if h > args.dair:
        raise ValueError("dipole is positioned within z-PML.")

    sr = L + args.dpml
    sz = args.dair + args.dpml
    cell_size = mp.Vector3(sr, 0, sz)

    boundary_layers = [
        mp.PML(args.dpml, direction=mp.R),
        mp.PML(args.dpml, direction=mp.Z, side=mp.High),
    ]

    src_cmpt = mp.Er
    src_cmpt_str = mp.component_name(src_cmpt)
    src_pt = mp.Vector3(rpos, 0, -0.5 * sz + h)
    sources = [mp.Source(src=src, component=src_cmpt, center=src_pt)]

    sim = mp.Simulation(
	resolution=args.res,
        cell_size=cell_size,
        dimensions=mp.CYLINDRICAL,
        m=m,
        boundary_layers=boundary_layers,
	sources=sources,
    )

    n2f_mon = sim.add_near2far(
        fcen,
        0,
        1,
        mp.Near2FarRegion(
            center=mp.Vector3(
                0.25 * L, 0, 0.5 * sz - args.dpml - 0.5 * args.dair
            ),
            size=mp.Vector3(0.5 * L, 0, 0),
        ),
        mp.Near2FarRegion(
            center=mp.Vector3(0.5 * L, 0, -0.5 * sz + 0.25 * args.dair),
            size=mp.Vector3(0, 0, 0.5 * args.dair),
        ),
    )

    mon_cmpt = mp.Ez
    mon_cmpt_str = mp.component_name(mon_cmpt)
    dft_mon = sim.add_dft_fields(
        [mon_cmpt],
        fcen,
        0,
        1,
        center=mp.Vector3(0.5 * sr, 0, 0),
        size=mp.Vector3(sr, 0, sz),
    )

    fig, ax = plt.subplots()
    sim.plot2D(ax=ax)
    if mp.am_master():
        fig.savefig(f'cyl_dft_dair{args.dair}_dpml{args.dpml}_plot2D.png',
                    dpi=150, bbox_inches='tight')

    def field_snapshot(sim):
        fig, ax = plt.subplots()
        sim.plot2D(
            ax=ax,
            fields=mon_cmpt,
            field_parameters={
                'colorbar':True,
                'cmap':'inferno',
                'post_process':lambda x: np.log10(np.abs(x)),
            },
            colorbar_parameters={'label':"log scale"},
            plot_monitors_flag=False,
            plot_sources_flag=False,
        )
        ax.set_title(
            f"|{mon_cmpt_str}($r$, $z$, $t$)|, "
            f"$t$ = {sim.meep_time():.2f}, "
            f"$m$ = {m}"
        )
        fig.savefig(
            f'{mon_cmpt_str}_t{sim.meep_time():.2f}.png',
            dpi=150,
            bbox_inches='tight',
        )

        e_t = sim.get_array(
            mon_cmpt,
            center=mp.Vector3(0.5 * sr, 0, -0.5 * sz + h),
            size=mp.Vector3(sr),
            cmplx=True,
        )
        fig, ax = plt.subplots()
        rs = np.linspace(0, sr, len(e_t))
        ax.semilogy(rs, np.abs(e_t), 'bo-')
        ax.set_xlabel('$r$')
        ax.set_ylabel(f"|{mon_cmpt_str}($r, $z$, t$)|, "
                      f"$z$ = {src_pt.z:.2f}, $t$ = {sim.meep_time():.2f}")
        ax.set_title(f"$m$ = {m}")
        ax.set_xlim(0, sr)
        fig.savefig(
            f'{mon_cmpt_str}_t{sim.meep_time():.2f}_z{abs(src_pt.z):.2f}.png',
            dpi=150,
            bbox_inches='tight',
        )


    sim.run(
        mp.at_every(5.4581274, field_snapshot),
        until_after_sources=mp.stop_when_fields_decayed(
            50,
            src_cmpt,
            src_pt,
            1e-9,
        ),
    )

    e_dft = sim.get_dft_array(dft_mon, mon_cmpt, 0)
    print(f"{mon_cmpt_str}_dft:, {np.shape(e_dft)}, ",
          f"{np.amin(np.real(e_dft)):.2f}, {np.amax(np.real(e_dft)):.2f}")

    fig, ax = plt.subplots()
    extent = [0, sr, -0.5 * sz, 0.5 * sz]
    im = ax.imshow(
        np.flipud(np.log10(np.abs(e_dft))),
        cmap='inferno',
        extent=extent,
    )
    ax.set_xlabel('$r$')
    ax.set_ylabel('$z$')
    ax.set_title(f'|{mon_cmpt_str}_dft($r$, $z$)|, $m$ = {m}')
    ax_divider = make_axes_locatable(ax)
    cax = ax_divider.append_axes(
        pad="2%",
        size="5%",
        position="right",
    )
    fig.colorbar(im, cax=cax)
    fig.savefig(f'src_{src_cmpt_str}_mon_{mon_cmpt_str}_'
                f'dft_m{m}_dair{args.dair}_res{args.res}.png',
                dpi=150, bbox_inches='tight')


    rs = np.linspace(0, sr, np.shape(e_dft)[1])
    zs = np.linspace(-0.5 * sz, 0.5 * sz, np.shape(e_dft)[0])
    fig, ax = plt.subplots()
    z_idx = int(h * args.res)
    ax.semilogy(rs, np.abs(e_dft[z_idx, :]), 'bo-')
    ax.set_xlabel('$r$')
    ax.set_ylabel(f"|{mon_cmpt_str}_dft($r$, $z$)|, $z$ = {zs[z_idx]:.2f}")
    ax.set_title(f"$m$ = {m}")
    ax.set_xlim(0, sr)
    fig.savefig(f'src_{src_cmpt_str}_mon_{mon_cmpt_str}_'
                f'dft_m{m}_dair{args.dair}_res{args.res}_z.png',
                dpi=150, bbox_inches='tight')

    return e_dft


if __name__ == "__main__":
    h = 0.15
    rpos = 0.5
    e_dft = radiated_fields(h, rpos, args.m)