radiated flux of a point dipole in 3d Cartesian coordinates

freespace_3d_flux.py

import argparse
import meep as mp
import numpy as np

parser = argparse.ArgumentParser()
parser.add_argument(
    '-res',
    type=float,
    help='resolution (pixels/μm)',
    default=100.0,
)
args = parser.parse_args()

resolution = args.res  # pixels/um                                                         

s = 5.0
dpml = 1.0
cell_size = mp.Vector3(s+2*dpml,s+2*dpml,s+2*dpml)

pml_layers = [mp.PML(dpml)]

n = 1.0  # refractive index of surrounding medium                                          

wvl = 1.0  # source wavelength                                                             
fcen = 1/wvl

src_cmpt = mp.Ex
sources = [mp.Source(src=mp.GaussianSource(fcen,fwidth=0.1*fcen),
                     center=mp.Vector3(),
                     component=src_cmpt)]

if src_cmpt == mp.Ex:
    symmetries = [mp.Mirror(mp.X,phase=-1),
                  mp.Mirror(mp.Y),
                  mp.Mirror(mp.Z)]
elif src_cmpt == mp.Ez:
    symmetries = [mp.Mirror(mp.X),
                  mp.Mirror(mp.Y),
                  mp.Mirror(mp.Z,phase=-1)]
else:
    symmetries = []

sim = mp.Simulation(cell_size=cell_size,
                    resolution=resolution,
                    sources=sources,
                    symmetries=symmetries,
                    boundary_layers=pml_layers,
                    default_material=mp.Medium(index=n))

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

sim.run(
    mp.dft_ldos(fcen,0,1),
    until_after_sources=mp.stop_when_fields_decayed(
        50,
        src_cmpt,
        mp.Vector3(),
        1e-8,
    ),
)

flux = mp.get_fluxes(flux_box)[0]
dV = 1/(resolution**3)
ldos = -np.real(sim.ldos_Fdata[0]*np.conj(sim.ldos_Jdata[0]))*dV
ratio = flux/ldos

print(f"flux:, {int(resolution):d}, {ldos:.6f}, {flux:.6f}, {ratio:.6f}")