1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
|
"""Tutorial example for dipole current sources in cylindrical coordinates.
tutorial reference:
https://meep.readthedocs.io/en/latest/Python_Tutorials/Cylindrical_Coordinates/#nonaxisymmetric-dipole-sources
"""
from typing import Tuple
import meep as mp
import numpy as np
RESOLUTION_UM = 50
WAVELENGTH_UM = 1.0
N_SLAB = 2.4
SLAB_THICKNESS_UM = 0.7 * WAVELENGTH_UM / N_SLAB
def dipole_in_slab(zpos: float, rpos_um: float, m: int) -> Tuple[float, float]:
"""Computes the flux from a dipole in a slab.
Args:
zpos: position of dipole as a fraction of layer thickness.
rpos_um: position of source in radial direction.
m: angular φ dependence of the fields exp(imφ).
Returns:
A 2-tuple of the radiated and total flux.
"""
pml_um = 1.0 # thickness of PML
padding_um = 1.0 # thickness of air padding
r_um = 20.0 # length of cell in r
frequency = 1 / WAVELENGTH_UM # center frequency of source/monitor
# runtime termination criteria
flux_decay_threshold = 1e-4
size_r = r_um + pml_um
size_z = SLAB_THICKNESS_UM + padding_um + pml_um
cell_size = mp.Vector3(size_r, 0, size_z)
boundary_layers = [
mp.PML(pml_um, direction=mp.R),
mp.PML(pml_um, direction=mp.Z, side=mp.High),
]
src_pt = mp.Vector3(rpos_um, 0, -0.5 * size_z + zpos * SLAB_THICKNESS_UM)
sources = [
mp.Source(
src=mp.GaussianSource(frequency, fwidth=0.05 * frequency),
component=mp.Er,
center=src_pt,
),
]
geometry = [
mp.Block(
material=mp.Medium(index=N_SLAB),
center=mp.Vector3(0, 0, -0.5 * size_z + 0.5 * SLAB_THICKNESS_UM),
size=mp.Vector3(mp.inf, mp.inf, SLAB_THICKNESS_UM),
)
]
sim = mp.Simulation(
resolution=RESOLUTION_UM,
cell_size=cell_size,
dimensions=mp.CYLINDRICAL,
m=m,
boundary_layers=boundary_layers,
sources=sources,
geometry=geometry,
force_complex_fields=True,
)
flux_mon = sim.add_flux(
frequency,
0,
1,
mp.FluxRegion(
center=mp.Vector3(0.5 * r_um, 0, 0.5 * size_z - pml_um),
size=mp.Vector3(r_um, 0, 0),
),
mp.FluxRegion(
center=mp.Vector3(r_um, 0, 0.5 * size_z - pml_um - 0.5 * padding_um),
size=mp.Vector3(0, 0, padding_um),
),
)
sim.run(
mp.dft_ldos(frequency, 0, 1),
until_after_sources=mp.stop_when_dft_decayed(tol=flux_decay_threshold),
)
radiated_flux = mp.get_fluxes(flux_mon)[0]
# volume of the ring current source
delta_vol = 2 * np.pi * rpos_um / (RESOLUTION_UM**2)
# total flux from point source via LDOS
source_flux = -np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) * delta_vol
print(
f"flux-cyl:, {rpos_um:.2f}, {m:3d}, " f"{source_flux:.6f}, {radiated_flux:.6f}"
)
return radiated_flux, source_flux
if __name__ == "__main__":
dipole_height = 0.5
# An Er source at r = 0 needs to be slightly offset.
# https://github.com/NanoComp/meep/issues/2704
dipole_rpos_um = 1.5 / RESOLUTION_UM
# Er source at r = 0 requires a single simulation with m = ±1.
m = 1
radiated_flux, source_flux = dipole_in_slab(
dipole_height,
dipole_rpos_um,
m,
)
extraction_efficiency = radiated_flux / source_flux
print(f"exteff:, {dipole_rpos_um}, {extraction_efficiency:.6f}")
# Er source at r > 0 requires Fourier-series expansion of φ.
# Threshold flux to determine when to truncate expansion.
flux_decay_threshold = 1e-2
dipole_rpos_um = [3.5, 6.7, 9.5]
for rpos_um in dipole_rpos_um:
source_flux_total = 0
radiated_flux_total = 0
radiated_flux_max = 0
m = 0
while True:
radiated_flux, source_flux = dipole_in_slab(
dipole_height,
rpos_um,
m,
)
radiated_flux_total += radiated_flux * (1 if m == 0 else 2)
source_flux_total += source_flux * (1 if m == 0 else 2)
if radiated_flux > radiated_flux_max:
radiated_flux_max = radiated_flux
if m > 0 and (radiated_flux / radiated_flux_max) < flux_decay_threshold:
break
else:
m += 1
extraction_efficiency = radiated_flux_total / source_flux_total
print(f"exteff:, {rpos_um}, {extraction_efficiency:.6f}")
|
