A test to verify the correctness of the near-to-far field transformation in cylindrical coordinates

freespace_reflector_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=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_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),
        ),
    )

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

    fig, ax = plt.subplots()
    sim.plot2D(ax=ax)
    if mp.am_master():
        fig.savefig(f'cyl_n2f_dair{args.dair}_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, len(er_dft)) # 0.1 * L, 0.9 * L                                                        

    for r in rs:
        er_pt = sim.get_farfield(n2f_mon, mp.Vector3(r, 0, 0.5 * sz - args.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.35)
    fig.suptitle(f'm = {m}', y=0.9, verticalalignment='bottom')
    fig.savefig(f'er_dft_vs_n2f_m{m}_dair{args.dair}_res{args.res}.png',
                dpi=150, bbox_inches='tight')

    return er_dft, er_ff


if __name__ == "__main__":
    h = 0.15
    rpos = 0.5
    er_dft, er_ff = radiated_fields(h, rpos, args.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}")