radiated flux from a point dipole at r=0 in cylindrical coordinates

freespace_cyl_flux.py

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

def radiated_flux(res: float, dpml_r: float, dpml_z: float, sr: float,
                  sz: float, fcen: float, m: int, rpos: float) -> float:
    """Computes the radiated flux for an Er point source in vacuum."""

    cell_size = mp.Vector3(sr+dpml_r,0,sz+2*dpml_z)

    pml_layers = [
        mp.PML(dpml_r, direction=mp.R),
        mp.PML(dpml_z, direction=mp.Z),
    ]

    sources = [
	mp.Source(
            src=mp.GaussianSource(
		frequency=fcen,
                fwidth=0.1*fcen,
            ),
            center=mp.Vector3(rpos,0,0),
            component=mp.Er,
        ),
    ]

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

    flux_plus_z = sim.add_flux(
        fcen,
	0,
        1,
        mp.FluxRegion(
            center=mp.Vector3(0.5*sr,0,0.5*sz),
            size=mp.Vector3(sr,0,0)
        )
    )

    flux_plus_r = sim.add_flux(
        fcen,
        0,
        1,
        mp.FluxRegion(
            center=mp.Vector3(sr,0,0),
            size=mp.Vector3(0,0,sz)
        )
    )

    flux_minus_z = sim.add_flux(
        fcen,
        0,
        1,
        mp.FluxRegion(
            center=mp.Vector3(0.5*sr,0,-0.5*sz),
            size=mp.Vector3(sr,0,0),
            weight=-1.0
        )
    )

    sim.run(
        until_after_sources=mp.stop_when_fields_decayed(
            50,
            mp.Er,
            mp.Vector3(),
            1e-8,
        ),
    )

    flux_plusz = mp.get_fluxes(flux_plus_z)[0]
    flux_plusr = mp.get_fluxes(flux_plus_r)[0]
    flux_minusz = mp.get_fluxes(flux_minus_z)[0]
    flux_tot = flux_plusz + flux_plusr + flux_minusz

    print(f"flux0:, {sim.meep_time():5.1f}, {sr}, {sz}, {rpos}, {m}, "
          f"{flux_plusz:.10f}, {flux_plusr:.10f}, {flux_minusz:.10f}, "
          f"{flux_tot:.10f}")

    return flux_tot


if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument(
        '-res',
        type=float,
        default=25.,
        help='resolution (default: 25.0)',
    )
    parser.add_argument(
        '-dpml_r',
        type=float,
        default=1.,
        help='PML thickness in r direction (default: 1.0)',
    )
    parser.add_argument(
        '-dpml_z',
        type=float,
        default=1.,
        help='PML thickness in z direction (default: 1.0)',
    )
    parser.add_argument(
        '-sr',
        type=float,
        default=5.,
        help='size of cell in r (default: 5.0)',
    )
    parser.add_argument(
        '-sz',
        type=float,
        default=5.,
        help='size of cell in z (default: 5.0)',
    )
    parser.add_argument(
        '-fcen',
        type=float,
        default=1.,
        help='frequency of pulsed source (default: 1.0)',
    )
    args = parser.parse_args()

    # r = 0                                                                                
    rpos = 0
    m = +1
    ref_out_flux = radiated_flux(args.res, args.dpml_r, args.dpml_z, args.sr,
                                 args.sz, args.fcen, m, rpos)
    print(f"flux0-m:, {m}, {ref_out_flux}")