Skip to content

Instantly share code, notes, and snippets.

@oskooi
Last active May 23, 2023 16:36
Show Gist options
  • Save oskooi/971f2bac2538dac19342b7e3d3e26a97 to your computer and use it in GitHub Desktop.
Save oskooi/971f2bac2538dac19342b7e3d3e26a97 to your computer and use it in GitHub Desktop.
A test to verify the correctness of the near-to-far field transformation in cylindrical coordinates
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}")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment