Skip to content

Instantly share code, notes, and snippets.

@oskooi
Last active April 26, 2023 04:26
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save oskooi/07bf0800ef2c8146c7602c10b710f476 to your computer and use it in GitHub Desktop.
Save oskooi/07bf0800ef2c8146c7602c10b710f476 to your computer and use it in GitHub Desktop.
Radiation pattern of a point dipole within an finite-size LED
import argparse
import math
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=25.,
help='resolution (default: 25 pixels/um)')
parser.add_argument('-L',
type=float,
default=5.0,
help='length of non-PML region (default: 5.0 um)')
parser.add_argument('-dpml',
type=float,
default=1.0,
help='PML thickness in r and z (default: 1.0 um)')
args = parser.parse_args()
resolution = args.res
L = args.L
dpml = args.dpml
dair = 1.0 # thickness of air padding
n = 2.4 # refractive index of surrounding medium
wvl = 1.0 # wavelength (in vacuum)
fcen = 1 / wvl # center frequency of source/monitor
def radiated_flux(dmat: float, h: float, rpos: float,
m: int) -> Tuple[float, float]:
"""Computes the radiated flux of a point dipole embedded
within a dielectric layer above a lossless ground plane in
cylindrical coordinates.
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=mp.GaussianSource(fcen, fwidth=0.1*fcen),
component=src_cmpt,
center=src_pt,
),
]
geometry = [
mp.Block(
material=mp.Medium(index=n),
center=mp.Vector3(0, 0, -0.5 * sz + 0.5 * dmat),
size=mp.Vector3(0.25 * L, 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,
)
# configuration 1
n2f_mon = sim.add_near2far(
fcen,
0,
1,
mp.Near2FarRegion(
center=mp.Vector3(0.5 * L, 0, 0.5 * sz - dpml),
size=mp.Vector3(L, 0, 0),
),
mp.Near2FarRegion(
center=mp.Vector3(L, 0, 0.5 * sz - dpml - 0.5 * (dmat + dair)),
size=mp.Vector3(0, 0, dmat + dair),
),
)
flux_air_mon = sim.add_flux(
fcen,
0,
1,
mp.FluxRegion(
center=mp.Vector3(0.5 * L, 0, 0.5 * sz - dpml),
size=mp.Vector3(L, 0, 0),
),
mp.FluxRegion(
center=mp.Vector3(L, 0, 0.5 * sz - dpml - 0.5 * (dmat + dair)),
size=mp.Vector3(0, 0, dmat + dair),
),
)
# configuration 2
# 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 * dai\
r)),
# size=mp.Vector3(0, 0, dmat + 0.5 * dair),
# ),
# )
# flux_air_mon = sim.add_flux(
# fcen,
# 0,
# 1,
# mp.FluxRegion(
# center=mp.Vector3(0.25 * L, 0, 0.5 * sz - dpml - 0.5 * dair),
# size=mp.Vector3(0.5 * L, 0, 0),
# ),
# mp.FluxRegion(
# center=mp.Vector3(0.5 * L, 0, -0.5 * sz + 0.5 * (dmat + 0.5 * dai\
r)),
# size=mp.Vector3(0, 0, dmat + 0.5 * dair),
# ),
# )
sim.run(
mp.dft_ldos(fcen, 0, 1),
until_after_sources=mp.stop_when_fields_decayed(
50.,
src_cmpt,
src_pt,
1e-8,
),
)
fig, ax = plt.subplots()
sim.plot2D(ax=ax)
if mp.am_master():
fig.savefig('led_cyl_plot2D.png', dpi=150, bbox_inches='tight')
flux_air = mp.get_fluxes(flux_air_mon)[0]
npts = 100 # number of points in [0,pi/2] range of angles
angles = np.linspace(0, 0.5*math.pi, npts)
# radius of quarter-circle flux monitor
r = 1000 * wvl
E = np.zeros((npts,3), dtype=np.complex128)
H = np.zeros((npts,3), dtype=np.complex128)
for i in range(npts):
ff = sim.get_farfield(n2f_mon,
mp.Vector3(r*math.cos(angles[i]),
0,
r*math.sin(angles[i])))
E[i,:] = [np.conj(ff[j]) for j in range(3)]
H[i,:] = [ff[j+3] for j in range(3)]
Pr = np.real(E[:,1]*H[:,2]-E[:,2]*H[:,1])
Pz = np.real(E[:,0]*H[:,1]-E[:,1]*H[:,0])
Prz = np.sqrt(np.square(Pr)+np.square(Pz))
fig, ax = plt.subplots()
ax.plot(np.degrees(angles), Prz / np.amax(Prz), 'b-')
ax.set_xlabel('angle (degrees)')
ax.set_ylabel('radial flux (normalized by maximum flux)')
ax.set_xlim(0, 90)
ax.set_ylim(0, 1)
ax.set_title('radiation pattern of $E_r$ point source at $r = 0$')
fig.savefig(f'led_cyl_radpattern_h{h}.png',dpi=150,bbox_inches='tight')
fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
ax.plot(
angles,
Prz / np.amax(Prz),
'k-',
)
ax.set_rmax(1)
ax.set_rticks([0,0.5,1])
ax.set_thetalim(0, 0.5*math.pi)
ax.grid(True)
ax.set_rlabel_position(22)
ax.set_title('radiation pattern of $E_r$ point source at $r = 0$')
ax.legend()
fig.savefig(f'led_cyl_radpattern_h{h}_polar.png',dpi=150,bbox_inches='tight')
if rpos == 0:
dV = np.pi / (resolution**3)
else:
dV = 2 * np.pi * rpos / (resolution**2)
flux_src = -np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) * dV
return flux_air, flux_src
if __name__ == "__main__":
layer_thickness = 0.7 * wvl / n
dipole_height = 0.5
# r = 0
rpos = 0
m = +1
flux_air, flux_src = radiated_flux(
layer_thickness,
dipole_height,
rpos,
m,
)
ext_eff = flux_air / flux_src
print(f"exteff:, {rpos}, {ext_eff:.6f}")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment