Skip to content

Commit 95c4a9f

Browse files
authored
fixes and tweaks for tutorial on phase of mode coefficients (#2507)
1 parent cedb7ae commit 95c4a9f

File tree

2 files changed

+87
-64
lines changed

2 files changed

+87
-64
lines changed

doc/docs/Python_Tutorials/Mode_Decomposition.md

Lines changed: 50 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -146,17 +146,17 @@ In the reflected-flux calculation, we apply our usual trick of first performing
146146
Phase of a Total Internal Reflected Mode
147147
----------------------------------------
148148

149-
For certain applications, such as [physically based ray tracing](https://pbr-book.org/), the phase of the reflection/transmission coefficient of a surface interface (flat or textured) is an important parameter used for modeling coherent effects (i.e., constructive or destructive interference from a collection of intersecting rays scattered off the surface). This tutorial demonstrates how to use mode decomposition to compute the phase of the reflection coefficient of a [total internal reflected](https://en.wikipedia.org/wiki/Total_internal_reflection) (TIR) mode. The simulated results are verified using the [Fresnel equations](https://en.wikipedia.org/wiki/Total_internal_reflection#Phase_shifts).
149+
For certain applications, such as [physically based ray tracing](https://pbr-book.org/), the phase of the reflection/transmission coefficient of a surface interface (flat or textured) is an important parameter used for modeling coherent effects (i.e., constructive or destructive interference from a collection of intersecting rays scattered off the surface). This tutorial demonstrates how to use mode decomposition to compute the phase of the reflection coefficient of a [total internal reflected](https://en.wikipedia.org/wiki/Total_internal_reflection) (TIR) mode. The simulated results are validated using the [Fresnel equations](https://en.wikipedia.org/wiki/Total_internal_reflection#Phase_shifts).
150150

151151
The example involves a 2d simulation of a flat interface of two lossless media with $n_1=1.5$ and $n_2=1.0$. The [critical angle](https://en.wikipedia.org/wiki/Total_internal_reflection#Critical_angle) for this interface is $\theta_c = 41.8^\circ$. The source is an incident planewave with linear polarization ($S$ or $P$). The 2d cell is periodic in $y$ with PML in the $x$ direction. The cell size does not affect the results and is therefore arbitrary. A schematic of the simulation layout is shown below.
152152

153-
The key item to consider in these types of calculations is the *location of the mode monitor relative to the interface*. The mode monitor is a line extending the entire length of the cell in the $y$ direction. In order to compare the result with the Fresnel equations, the phase of the reflection coefficient must be computed *exactly* at the interface. However, it is problematic to measure the reflection coefficient exactly at the interface because the amplitude of the left-going wave drops discontinuously to zero across the interface. For this reason, the mode monitor must be positioned away from the interface (somewhere within the $n_1$ medium) and the measured phase must be corrected in post processing to account for this offset.
154-
155153
![](../images/refl_coeff_flat_interface.png#center)
156154

155+
The key item to consider in these types of calculations is the *location of the mode monitor relative to the interface*. The mode monitor is a line extending the entire length of the cell in the $y$ direction. In order to compare the result with the Fresnel equations, the phase of the reflection coefficient must be computed *exactly* at the interface. However, it is problematic to measure the reflection coefficient exactly at the interface because the amplitude of the left-going wave drops discontinuously to zero across the interface. For this reason, the mode monitor must be positioned away from the interface (somewhere within the $n_1$ medium) and the measured phase must be corrected in post processing to account for this offset.
156+
157157
The calculation of the reflection coefficient requires two separate runs: (1) a normalization run involving just the source medium ($n_1$) to obtain the incident fields, and (2) the main run including the interface whereby the incident fields from (1) are first subtracted from the monitor to obtain only the reflected fields. The mode coefficient in (2) divided by (1) is, by definition, the reflection coefficient.
158158

159-
The phase of the reflection coefficient needs to be corrected for the offset of the mode monitor relative to the interface — labeled $L$ in the schematic above — using the formula $\exp(i 2\pi k_x 2L)$, where $k_x$ is the $x$ component of the planewave's wavevector (`k` in the script). The $2\pi$ factor is necessary because `k` does *not* include this factor (per convention in Meep). The factor 2 in front of the $L$ is necessary because the phase needs to be corrected for the monitors in the normalization and main runs separately. The correction factor is just the phase accumulated as the wave propagates in the surface-normal direction for a distance $L$ in a medium with index $n_1$.
159+
The phase of the reflection coefficient needs to be corrected for the offset of the mode monitor relative to the interface — labeled $L$ in the schematic above — using the formula $\exp(i 2\pi k_x 2L)$, where $k_x$ is the $x$ component of the planewave's wavevector (`k` in the script). The $2\pi$ factor is necessary because `k` does *not* include this factor (per convention in Meep). The factor 2 in front of the $L$ is necessary because the phase needs to be corrected for the monitors in the normalization and main runs separately. The correction factor is just the phase accumulated as the wave propagates in the surface-normal direction for a distance $L$ in a medium with index $n_1$. (A transmitted mode would involve a correction factor for a medium with index $n_2$.)
160160

161161
With this setup, we measure the phase of the reflection coefficient for two different source configurations (polarization and angle) and compare the result with the Fresnel equations. The location of the mode monitor is also varied in the two test configurations. Results are shown in the two tables below. There is good agreement between the simulation and theory. (As additional validation, we note that the magnitude of the reflection coefficient of a TIR mode must be one which is indeed the case.)
162162

@@ -189,15 +189,18 @@ n2 = 1.0
189189

190190

191191
def refl_coeff_meep(pol: Polarization, theta: float, L: float) -> complex:
192-
"""Computes the reflection coefficient of a TIR mode using mode
193-
decomposition.
192+
"""Returns the complex reflection coefficient of a TIR mode computed
193+
using mode decomposition.
194194
195195
Args:
196196
pol: polarization of the incident planewave (S or P).
197197
theta: angle of the incident planewave (radians).
198-
L: position of the mode monitor relative to the flat interface. 0 is
199-
interface position.
198+
L: position of the mode monitor relative to the flat interface.
200199
"""
200+
if theta < math.asin(n2 / n1):
201+
raise ValueError(f"incident angle of {math.degrees(theta):.2f}° is "
202+
f"not a total internal reflected mode.")
203+
201204
resolution = 50.0
202205

203206
# cell size is arbitrary
@@ -207,18 +210,17 @@ def refl_coeff_meep(pol: Polarization, theta: float, L: float) -> complex:
207210
cell_size = mp.Vector3(sx + 2 * dpml, sy, 0)
208211
pml_layers = [mp.PML(dpml, direction=mp.X)]
209212

210-
fcen = 1.0 # center frequency
213+
fcen = 1.0 # source/monitor frequency
211214
df = 0.1 * fcen
212215

213216
# k (in source medium) with correct length
214217
# plane of incidence is xy
215218
k = mp.Vector3(n1 * fcen, 0, 0).rotate(mp.Vector3(0, 0, 1), theta)
216219

217-
# planewave amplitude function (for source)
220+
# planewave amplitude function (for line source)
218221
def pw_amp(k, x0):
219222
def _pw_amp(x):
220223
return cmath.exp(1j * 2 * math.pi * k.dot(x + x0))
221-
222224
return _pw_amp
223225

224226
src_pt = mp.Vector3(-0.5 * sx, 0, 0)
@@ -322,7 +324,7 @@ def refl_coeff_meep(pol: Polarization, theta: float, L: float) -> complex:
322324
sim.run(
323325
until_after_sources=mp.stop_when_fields_decayed(
324326
50,
325-
mp.Ez,
327+
src_cmpt,
326328
mp.Vector3(-L, 0, 0),
327329
1e-6,
328330
),
@@ -349,12 +351,12 @@ def refl_coeff_meep(pol: Polarization, theta: float, L: float) -> complex:
349351

350352

351353
def refl_coeff_Fresnel(pol: Polarization, theta: float) -> complex:
352-
"""Computes the reflection coefficient of a TIR mode using the Fresnel
353-
equations.
354+
"""Returns the complex reflection coefficient of a TIR mode computed
355+
using the Fresnel equations.
354356
355357
Args:
356358
pol: polarization of the incident planewave (S or P).
357-
theta: angle of the incident planewave (degrees).
359+
theta: angle of the incident planewave (radians).
358360
"""
359361
if pol.name == "S":
360362
refl_coeff = (
@@ -373,31 +375,42 @@ def refl_coeff_Fresnel(pol: Polarization, theta: float) -> complex:
373375

374376

375377
if __name__ == "__main__":
376-
thetas = [54.3, 48.5] # angle of incident planewave (degrees)
377-
Ls = [0.4, 1.2] # position of mode monitor relative to flat interface
378-
pols = [Polarization.S, Polarization.P] # polarization of incident planewave
378+
# angle of incident planewave (degrees)
379+
thetas = [54.3, 48.5]
380+
381+
# position of mode monitor relative to flat interface
382+
Ls = [0.4, 1.2]
383+
384+
# polarization of incident planewave
385+
pols = [Polarization.S, Polarization.P]
379386

380387
for pol, theta, L in zip(pols, thetas, Ls):
381388
theta_rad = np.radians(theta)
382-
rc_m = refl_coeff_meep(pol, theta_rad, L)
383-
rc_f = refl_coeff_Fresnel(pol, theta_rad)
384-
385-
rc_m_str = f"{rc_m.real:.5f}{rc_m.imag:+.5f}j"
386-
rc_f_str = f"{rc_f.real:.5f}{rc_f.imag:+.5f}j"
387-
print(f"refl-coeff:, {pol.name}, {theta}, {rc_m_str} (Meep), "
388-
f"{rc_f_str} (Fresnel)")
389-
390-
mag_m = abs(rc_m)
391-
mag_f = abs(rc_f)
392-
err_mag = abs(mag_m - mag_f) / mag_f
393-
print(f"magnitude:, {mag_m:.5f} (Meep), {mag_f:.5f} (Fresnel), "
394-
f"{err_mag:.5f} (error)")
395-
396-
phase_m = cmath.phase(rc_m)
397-
phase_f = cmath.phase(rc_f)
398-
err_phase = abs(phase_m - phase_f) / abs(phase_f)
399-
print(f"phase:, {phase_m:.5f} (Meep), {phase_f:.5f} (Fresnel), "
400-
f"{err_phase:.5f} (error)")
389+
R_meep = refl_coeff_meep(pol, theta_rad, L)
390+
R_fres = refl_coeff_Fresnel(pol, theta_rad)
391+
392+
complex_to_str = lambda cnum: f"{cnum.real:.5f}{cnum.imag:+.5f}j"
393+
print(
394+
f"refl-coeff:, {pol.name}, {theta}, "
395+
f"{complex_to_str(R_meep)} (Meep), "
396+
f"{complex_to_str(R_fres)} (Fresnel)"
397+
)
398+
399+
mag_meep = abs(R_meep)
400+
mag_fres = abs(R_fres)
401+
err_mag = abs(mag_meep - mag_fres) / mag_fres
402+
print(
403+
f"magnitude:, {mag_meep:.5f} (Meep), {mag_fres:.5f} (Fresnel), "
404+
f"{err_mag:.5f} (error)"
405+
)
406+
407+
phase_meep = cmath.phase(R_meep)
408+
phase_fres = cmath.phase(R_fres)
409+
err_phase = abs(phase_meep - phase_fres) / abs(phase_fres)
410+
print(
411+
f"phase:, {phase_meep:.5f} (Meep), {phase_fres:.5f} (Fresnel), "
412+
f"{err_phase:.5f} (error)"
413+
)
401414
```
402415

403416

python/examples/mode_coeff_phase.py

Lines changed: 37 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -22,15 +22,20 @@
2222

2323

2424
def refl_coeff_meep(pol: Polarization, theta: float, L: float) -> complex:
25-
"""Computes the reflection coefficient of a TIR mode using mode
26-
decomposition.
25+
"""Returns the complex reflection coefficient of a TIR mode computed
26+
using mode decomposition.
2727
2828
Args:
2929
pol: polarization of the incident planewave (S or P).
3030
theta: angle of the incident planewave (radians).
31-
L: position of the mode monitor relative to the flat interface. 0 is
32-
interface position.
31+
L: position of the mode monitor relative to the flat interface.
3332
"""
33+
if theta < math.asin(n2 / n1):
34+
raise ValueError(
35+
f"incident angle of {math.degrees(theta):.2f}° is "
36+
f"not a total internal reflected mode."
37+
)
38+
3439
resolution = 50.0
3540

3641
# cell size is arbitrary
@@ -40,14 +45,14 @@ def refl_coeff_meep(pol: Polarization, theta: float, L: float) -> complex:
4045
cell_size = mp.Vector3(sx + 2 * dpml, sy, 0)
4146
pml_layers = [mp.PML(dpml, direction=mp.X)]
4247

43-
fcen = 1.0 # center frequency
48+
fcen = 1.0 # source/monitor frequency
4449
df = 0.1 * fcen
4550

4651
# k (in source medium) with correct length
4752
# plane of incidence is xy
4853
k = mp.Vector3(n1 * fcen, 0, 0).rotate(mp.Vector3(0, 0, 1), theta)
4954

50-
# planewave amplitude function (for source)
55+
# planewave amplitude function (for line source)
5156
def pw_amp(k, x0):
5257
def _pw_amp(x):
5358
return cmath.exp(1j * 2 * math.pi * k.dot(x + x0))
@@ -155,7 +160,7 @@ def _pw_amp(x):
155160
sim.run(
156161
until_after_sources=mp.stop_when_fields_decayed(
157162
50,
158-
mp.Ez,
163+
src_cmpt,
159164
mp.Vector3(-L, 0, 0),
160165
1e-6,
161166
),
@@ -182,12 +187,12 @@ def _pw_amp(x):
182187

183188

184189
def refl_coeff_Fresnel(pol: Polarization, theta: float) -> complex:
185-
"""Computes the reflection coefficient of a TIR mode using the Fresnel
186-
equations.
190+
"""Returns the complex reflection coefficient of a TIR mode computed
191+
using the Fresnel equations.
187192
188193
Args:
189194
pol: polarization of the incident planewave (S or P).
190-
theta: angle of the incident planewave (degrees).
195+
theta: angle of the incident planewave (radians).
191196
"""
192197
if pol.name == "S":
193198
refl_coeff = (
@@ -206,34 +211,39 @@ def refl_coeff_Fresnel(pol: Polarization, theta: float) -> complex:
206211

207212

208213
if __name__ == "__main__":
209-
thetas = [54.3, 48.5] # angle of incident planewave (degrees)
210-
Ls = [0.4, 1.2] # position of mode monitor relative to flat interface
211-
pols = [Polarization.S, Polarization.P] # polarization of incident planewave
214+
# angle of incident planewave (degrees)
215+
thetas = [54.3, 48.5]
216+
217+
# position of mode monitor relative to flat interface
218+
Ls = [0.4, 1.2]
219+
220+
# polarization of incident planewave
221+
pols = [Polarization.S, Polarization.P]
212222

213223
for pol, theta, L in zip(pols, thetas, Ls):
214224
theta_rad = np.radians(theta)
215-
rc_m = refl_coeff_meep(pol, theta_rad, L)
216-
rc_f = refl_coeff_Fresnel(pol, theta_rad)
225+
R_meep = refl_coeff_meep(pol, theta_rad, L)
226+
R_fres = refl_coeff_Fresnel(pol, theta_rad)
217227

218-
rc_m_str = f"{rc_m.real:.5f}{rc_m.imag:+.5f}j"
219-
rc_f_str = f"{rc_f.real:.5f}{rc_f.imag:+.5f}j"
228+
complex_to_str = lambda cnum: f"{cnum.real:.5f}{cnum.imag:+.5f}j"
220229
print(
221-
f"refl-coeff:, {pol.name}, {theta}, {rc_m_str} (Meep), "
222-
f"{rc_f_str} (Fresnel)"
230+
f"refl-coeff:, {pol.name}, {theta}, "
231+
f"{complex_to_str(R_meep)} (Meep), "
232+
f"{complex_to_str(R_fres)} (Fresnel)"
223233
)
224234

225-
mag_m = abs(rc_m)
226-
mag_f = abs(rc_f)
227-
err_mag = abs(mag_m - mag_f) / mag_f
235+
mag_meep = abs(R_meep)
236+
mag_fres = abs(R_fres)
237+
err_mag = abs(mag_meep - mag_fres) / mag_fres
228238
print(
229-
f"magnitude:, {mag_m:.5f} (Meep), {mag_f:.5f} (Fresnel), "
239+
f"magnitude:, {mag_meep:.5f} (Meep), {mag_fres:.5f} (Fresnel), "
230240
f"{err_mag:.5f} (error)"
231241
)
232242

233-
phase_m = cmath.phase(rc_m)
234-
phase_f = cmath.phase(rc_f)
235-
err_phase = abs(phase_m - phase_f) / abs(phase_f)
243+
phase_meep = cmath.phase(R_meep)
244+
phase_fres = cmath.phase(R_fres)
245+
err_phase = abs(phase_meep - phase_fres) / abs(phase_fres)
236246
print(
237-
f"phase:, {phase_m:.5f} (Meep), {phase_f:.5f} (Fresnel), "
247+
f"phase:, {phase_meep:.5f} (Meep), {phase_fres:.5f} (Fresnel), "
238248
f"{err_phase:.5f} (error)"
239249
)

0 commit comments

Comments
 (0)