Add unit test for phase of mode coefficients

python/tests/test_mode_decomposition.py

@@ -1,10 +1,13 @@
 import cmath
+from enum import Enum
 import math
+import parameterized
 import unittest
 
+import meep as mp
 import numpy as np
 
-import meep as mp
+Polarization = Enum("Polarization", "S P")
 
 
 class TestModeDecomposition(unittest.TestCase):
@@ -182,14 +185,14 @@ def test_oblique_waveguide_backward_mode(self):
 
         flux = mp.get_fluxes(mode)[0]
         coeff = sim.get_eigenmode_coefficients(
-            mode, [1], direction=mp.NO_DIRECTION, kpoint_func=lambda f, n: kpoint
+            mode, [1], direction=mp.NO_DIRECTION, kpoint_func=lambda *not_used: kpoint
         ).alpha[0, 0, 0]
         flux_decimated = mp.get_fluxes(mode_decimated)[0]
         coeff_decimated = sim.get_eigenmode_coefficients(
             mode_decimated,
             [1],
             direction=mp.NO_DIRECTION,
-            kpoint_func=lambda f, n: kpoint,
+            kpoint_func=lambda *not_used: kpoint,
         ).alpha[0, 0, 0]
 
         print(f"oblique-waveguide-flux:, {-flux:.6f}, {abs(coeff) ** 2:.6f}")
@@ -665,6 +668,201 @@ def _pw_amp(x):
         self.assertAlmostEqual(Tsum, Tflux, places=3)
         self.assertAlmostEqual(Rsum + Tsum, 1.00, places=2)
 
+    @parameterized.parameterized.expand(
+        [
+            (Polarization.S, 54.3, 0.4),
+            (Polarization.P, 48.5, 1.2),
+        ]
+    )
+    def test_phase(self, pol: Polarization, theta: float, L: float):
+        """Unit test for phase of mode coefficients.
+
+        Verifies that the phase of a total internal reflected (TIR) mode of
+        a flat interface of two lossless materials given an incident planewave
+        at oblique incidence matches the Fresnel equations.
+
+        Args:
+            pol: polarization of the incident planewave (S or P).
+            theta: angle of the incident planewave (degrees).
+            L: position of the mode monitor relative to the flat interface.
+        """
+        resolution = 50.0
+
+        sx = 7.0
+        sy = 3.0
+        dpml = 2.0
+        cell_size = mp.Vector3(sx + 2 * dpml, sy, 0)
+        pml_layers = [mp.PML(dpml, direction=mp.X)]
+
+        n1 = 1.5
+        n2 = 1.0
+
+        # angle of incident planewave at center frequency
+        # 0° is +x; rotated CCW about z axis
+        theta = np.radians(theta)
+
+        fcen = 1.0  # center frequency
+        df = 0.1 * fcen
+
+        # k (in source medium) with correct length
+        # plane of incidence is xy
+        k = mp.Vector3(n1 * fcen, 0, 0).rotate(mp.Vector3(0, 0, 1), theta)
+
+        def pw_amp(k, x0):
+            def _pw_amp(x):
+                return cmath.exp(1j * 2 * math.pi * k.dot(x + x0))
+
+            return _pw_amp
+
+        src_pt = mp.Vector3(-0.5 * sx, 0, 0)
+
+        if pol.name == "S":
+            src_cmpt = mp.Ez
+            eig_parity = mp.ODD_Z
+        elif pol.name == "P":
+            src_cmpt = mp.Hz
+            eig_parity = mp.EVEN_Z
+        else:
+            raise ValueError("pol must be S or P, only.")
+
+        sources = [
+            mp.Source(
+                mp.GaussianSource(fcen, fwidth=df),
+                component=src_cmpt,
+                center=src_pt,
+                size=mp.Vector3(0, cell_size.y, 0),
+                amp_func=pw_amp(k, src_pt),
+            ),
+        ]
+
+        sim = mp.Simulation(
+            resolution=resolution,
+            cell_size=cell_size,
+            default_material=mp.Medium(index=n1),
+            boundary_layers=pml_layers,
+            k_point=k,
+            sources=sources,
+        )
+
+        # DFT monitor for incident fields
+        mode_mon = sim.add_mode_monitor(
+            fcen,
+            0,
+            1,
+            mp.FluxRegion(
+                center=mp.Vector3(-L, 0, 0),
+                size=mp.Vector3(0, cell_size.y, 0),
+            ),
+        )
+
+        sim.run(
+            until_after_sources=mp.stop_when_fields_decayed(
+                50,
+                src_cmpt,
+                mp.Vector3(-L, 0, 0),
+                1e-6,
+            ),
+        )
+
+        res = sim.get_eigenmode_coefficients(
+            mode_mon,
+            bands=[1],
+            eig_parity=eig_parity,
+            kpoint_func=lambda *not_used: k,
+            direction=mp.NO_DIRECTION,
+        )
+
+        input_mode_coeff = res.alpha[0, 0, 0]
+        input_flux_data = sim.get_flux_data(mode_mon)
+
+        sim.reset_meep()
+
+        geometry = [
+            mp.Block(
+                material=mp.Medium(index=n1),
+                center=mp.Vector3(-0.25 * (sx + 2 * dpml), 0, 0),
+                size=mp.Vector3(0.5 * (sx + 2 * dpml), mp.inf, mp.inf),
+            ),
+            mp.Block(
+                material=mp.Medium(index=n2),
+                center=mp.Vector3(0.25 * (sx + 2 * dpml), 0, 0),
+                size=mp.Vector3(0.5 * (sx + 2 * dpml), mp.inf, mp.inf),
+            ),
+        ]
+
+        sim = mp.Simulation(
+            resolution=resolution,
+            cell_size=cell_size,
+            boundary_layers=pml_layers,
+            k_point=k,
+            sources=sources,
+            geometry=geometry,
+        )
+
+        # DFT monitor for reflected fields
+        mode_mon = sim.add_mode_monitor(
+            fcen,
+            0,
+            1,
+            mp.FluxRegion(
+                center=mp.Vector3(-L, 0, 0),
+                size=mp.Vector3(0, cell_size.y, 0),
+            ),
+        )
+
+        sim.load_minus_flux_data(mode_mon, input_flux_data)
+
+        sim.run(
+            until_after_sources=mp.stop_when_fields_decayed(
+                50,
+                mp.Ez,
+                mp.Vector3(-L, 0, 0),
+                1e-6,
+            ),
+        )
+
+        res = sim.get_eigenmode_coefficients(
+            mode_mon,
+            bands=[1],
+            eig_parity=eig_parity,
+            kpoint_func=lambda *not_used: k,
+            direction=mp.NO_DIRECTION,
+        )
+
+        # mode coefficient of reflected planewave
+        refl_mode_coeff = res.alpha[0, 0, 1]
+
+        # reflection coefficient
+        refl_coeff = refl_mode_coeff / input_mode_coeff
+
+        # apply phase correction factor
+        refl_coeff /= cmath.exp(1j * k.x * 2 * math.pi * 2 * L)
+
+        # reflection coefficient (Fresnel equations)
+        if pol.name == "S":
+            refl_coeff_Fresnel = (
+                math.cos(theta) - ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5
+            ) / (math.cos(theta) + ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5)
+        else:
+            refl_coeff_Fresnel = (
+                -((n2 / n1) ** 2) * math.cos(theta)
+                + ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5
+            ) / (
+                (n2 / n1) ** 2 * math.cos(theta)
+                + ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5
+            )
+
+        print(
+            f"phase:, {pol.name}, {cmath.phase(refl_coeff)} (Meep), "
+            f"{cmath.phase(refl_coeff_Fresnel)} (Fresnel)"
+        )
+
+        self.assertAlmostEqual(
+            cmath.phase(refl_coeff),
+            cmath.phase(refl_coeff_Fresnel),
+            delta=0.04,
+        )
+
 
 if __name__ == "__main__":
     unittest.main()