Validating the near-to-far field transformation function (`greencyl`) in cylindrical coordinates

led_disc_cyl_n2f.py

import argparse
from typing import Tuple

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

parser = argparse.ArgumentParser()
parser.add_argument('-res',
                    type=float,
                    default=50.0,
                    help='resolution (default: 50 pixels/um)')
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()

resolution = args.res
dair = args.dair
dpml = args.dpml

L = 6.0          # length of cell in r (excluding PML)                                                            
n = 2.4          # refractive index of surrounding medium                                                         
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(dmat: float, h: float, rpos: float,
                    m: int) -> Tuple[np.ndarray, np.ndarray]:
    """Computes the radiated fields of a point dipole embedded within                                             
       a dielectric disc 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 as fraction of dmat.                                               
        rpos: position of source in r direction.                                                                  
        m: angular φ dependence of the fields exp(imφ).                                                           
    """
    sr = L + dpml
    sz = dmat + dair + dpml
    cell_size = mp.Vector3(sr, 0, sz)

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

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

    geometry = [
        mp.Block(
            material=mp.Medium(index=n),
            center=mp.Vector3(0.5, 0, -0.5 * sz + 0.5 * dmat),
            size=mp.Vector3(1.0, mp.inf, dmat),
        )
    ]

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

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

    dft_mon = sim.add_dft_fields(
        [mp.Er],
        fcen,
        0,
        1,
        center=mp.Vector3(0.5 * L, 0, 0.5 * sz - dpml),
        size=mp.Vector3(L, 0, 0),
    )

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

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

    er_dft = sim.get_dft_array(dft_mon, mp.Er, 0)

    er_ff = []
    rs = np.linspace(0, L, int(resolution * L) + 1)

    for r in rs:
        er_pt = sim.get_farfield(n2f_mon, mp.Vector3(r, 0, 0.5 * sz - dpml))
        er_ff.append(er_pt[0])

    er_ff = np.array(er_ff)

    fig, ax = plt.subplots(1, 2)
    ax[0].plot(rs, np.real(er_dft), 'bo-', label='DFT')
    ax[0].plot(rs, np.real(er_ff), 'ro-', label='N2F')
    ax[0].set_xlabel('$r$')
    ax[0].set_ylabel('real($E_r$)')
    ax[0].legend()
    ax[1].plot(rs, np.imag(er_dft), 'bo-', label='DFT')
    ax[1].plot(rs, np.imag(er_ff), 'ro-', label='N2F')
    ax[1].set_xlabel('$r$')
    ax[1].set_ylabel('imag($E_r$)')
    ax[1].legend()
    fig.subplots_adjust(wspace=0.3)
    fig.suptitle(f'm = {m}', y=0.95, verticalalignment='bottom')
    fig.savefig(f'er_dft_vs_n2f_m{m}_dair{dair}_res{resolution}.png',
		dpi=150, bbox_inches='tight')

    return er_dft, er_ff


if __name__ == "__main__":
    layer_thickness = 0.7 * wvl / n
    dipole_height = 0.5

    rpos = 0.2
    m = 2
    er_dft, er_ff = radiated_fields(
        layer_thickness,
        dipole_height,
        rpos,
        m,
    )

    rel_err_real = ((np.linalg.norm(np.real(er_dft) - np.real(er_ff))) /
                    np.linalg.norm(np.real(er_dft)))

    rel_err_imag = ((np.linalg.norm(np.imag(er_dft) - np.imag(er_ff))) /
                    np.linalg.norm(np.imag(er_dft)))

    print(f"norm:, {rel_err_real}, {rel_err_imag}")