tutorial and unit test for mode sources using DiffractedPlanewave object

doc/docs/Python_Tutorials/Eigenmode_Source.md

@@ -245,9 +245,9 @@ cell_size = mp.Vector3(14,10,0)
 pml_layers = [mp.PML(thickness=2,direction=mp.X)]
 
 # rotation angle (in degrees) of planewave, counter clockwise (CCW) around z-axis
-rot_angle = np.radians(0)
+rot_angle = np.radians(23.2)
 
-fsrc = 1.0 # frequency of planewave (wavelength = 1/fsrc)
+fsrc = 1.0 # frequency (wavelength = 1/fsrc)
 
 n = 1.5  # refractive index of homogeneous material
 default_material = mp.Medium(index=n)
@@ -259,9 +259,7 @@ sources = [mp.EigenModeSource(src=mp.ContinuousSource(fsrc),
                               size=mp.Vector3(y=10),
                               direction=mp.AUTOMATIC if rot_angle == 0 else mp.NO_DIRECTION,
                               eig_kpoint=k_point,
-                              eig_band=1,
-                              eig_parity=mp.EVEN_Y+mp.ODD_Z if rot_angle == 0 else mp.ODD_Z,
-                              eig_match_freq=True)]
+                              eig_parity=mp.EVEN_Y+mp.ODD_Z if rot_angle == 0 else mp.ODD_Z)]
 
 sim = mp.Simulation(cell_size=cell_size,
                     resolution=resolution,
@@ -271,20 +269,31 @@ sim = mp.Simulation(cell_size=cell_size,
                     default_material=default_material,
                     symmetries=[mp.Mirror(mp.Y)] if rot_angle == 0 else [])
 
-sim.run(until=100)
+sim.run(until=20.)
 
 nonpml_vol = mp.Volume(center=mp.Vector3(), size=mp.Vector3(10,10,0))
 
 sim.plot2D(fields=mp.Ez,
            output_plane=nonpml_vol)
 
-if mp.am_master():
-    plt.axis('off')
-    plt.savefig('pw.png',bbox_inches='tight')
+plt.axis('off')
+plt.savefig('pw.png',bbox_inches='tight')
 ```
 
 Note that the line source spans the *entire* length of the cell. (Planewave sources extending into the PML region must include `is_integrated=True`.) This example involves a continuous-wave (CW) time profile. For a pulse profile, the oblique planewave is incident at a given angle for only a *single* frequency component of the source as described in [Tutorial/Basics/Angular Reflectance Spectrum of a Planar Interface](../Python_Tutorials/Basics.md#angular-reflectance-spectrum-of-a-planar-interface). This is a fundamental feature of FDTD simulations and not of Meep per se. To simulate an incident planewave at multiple angles for a given frequency $\omega$, you will need to do separate simulations involving different values of $\vec{k}$ (`k_point`) since each set of $(\vec{k},\omega)$ specifying the Bloch-periodic boundaries and the frequency of the source will produce a different angle of the planewave. For more details, refer to Section 4.5 ("Efficient Frequency-Angle Coverage") in [Chapter 4](https://arxiv.org/abs/1301.5366) ("Electromagnetic Wave Source Conditions") of [Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707).
 
 Shown below are the steady-state field profiles generated by the planewave for the three rotation angles. Residues of the backward-propagating waves due to the discretization are slightly visible.
 
 ![](../images/eigenmode_planewave.png#center)
+
+An alternative method to launching a planewave in homogeneous media which provides more control over specifying its polarization (particularly in 3D) is to use a `DiffractedPlanewave` object for the `eig_band` property of the `EigenModeSource` instead of a band number and parity as in the previous example which used 1 and `mp.ODD_Z`, respectively. A planewave with wavevector $\vec{k}$ can be defined using the *zeroth* diffraction order combined with the `k_point` of the `Simulation` object set to $\vec{k}$. As a demonstration, the previous example which involved a $E_z$-polarized source (which is equivalent to the $\mathcal{S}$ polarization given the plane of incidence of $xy$) source, can also be defined using:
+
+```py
+sources = [mp.EigenModeSource(mp.ContinuousSource(fsrc),
+                              center=mp.Vector3(),
+                              size=mp.Vector3(y=10),
+                              eig_band=mp.DiffractedPlanewave((0,0,0),
+                                                              mp.Vector3(0,1,0),
+                                                              1,
+                                                              0))]
+```

python/examples/oblique-planewave.py

@@ -3,16 +3,17 @@
 
 import meep as mp
 
+
 resolution = 50  # pixels/μm
 
 cell_size = mp.Vector3(14, 10, 0)
 
 pml_layers = [mp.PML(thickness=2, direction=mp.X)]
 
 # rotation angle (in degrees) of planewave, counter clockwise (CCW) around z-axis
-rot_angle = np.radians(0)
+rot_angle = np.radians(23.2)
 
-fsrc = 1.0  # frequency of planewave (wavelength = 1/fsrc)
+fsrc = 1.0  # frequency (wavelength = 1/fsrc)
 
 n = 1.5  # refractive index of homogeneous material
 default_material = mp.Medium(index=n)
@@ -26,12 +27,20 @@
         size=mp.Vector3(y=10),
         direction=mp.AUTOMATIC if rot_angle == 0 else mp.NO_DIRECTION,
         eig_kpoint=k_point,
-        eig_band=1,
         eig_parity=mp.EVEN_Y + mp.ODD_Z if rot_angle == 0 else mp.ODD_Z,
-        eig_match_freq=True,
     )
 ]
 
+## equivalent definition
+# sources = [
+#     mp.EigenModeSource(
+#         src=mp.ContinuousSource(fsrc),
+#         center=mp.Vector3(),
+#         size=mp.Vector3(y=10),
+#         eig_band=mp.DiffractedPlanewave((0, 0, 0), mp.Vector3(0, 1, 0), 1, 0),
+#     )
+# ]
+
 sim = mp.Simulation(
     cell_size=cell_size,
     resolution=resolution,
@@ -42,12 +51,11 @@
     symmetries=[mp.Mirror(mp.Y)] if rot_angle == 0 else [],
 )
 
-sim.run(until=100)
+sim.run(until=20.0)
 
-nonpml_vol = mp.Volume(center=mp.Vector3(), size=mp.Vector3(10, 10, 0))
+nonpml_vol = mp.Volume(center=mp.Vector3(), size=mp.Vector3(10.0, 10.0, 0))
 
 sim.plot2D(fields=mp.Ez, output_plane=nonpml_vol)
 
-if mp.am_master():
-    plt.axis("off")
-    plt.savefig("pw.png", bbox_inches="tight")
+plt.axis("off")
+plt.savefig("pw.png", bbox_inches="tight")

python/tests/test_diffracted_planewave.py

@@ -6,11 +6,6 @@
 
 import meep as mp
 
-# Computes the mode coefficient of the transmitted orders of
-# a binary grating given an incident planewave and verifies
-# that the results are the same when using either a band number
-# or `DiffractedPlanewave` object in `get_eigenmode_coefficients`.
-
 
 class TestDiffractedPlanewave(unittest.TestCase):
     @classmethod
@@ -28,24 +23,33 @@ def setUp(cls):
 
         cls.pml_layers = [mp.PML(thickness=cls.dpml, direction=mp.X)]
 
-    def run_binary_grating_diffraction(self, gp, gh, gdc, theta):
+    def run_binary_grating_diffraction(self, gp, gh, gdc, theta_pw):
+        """Computes the mode coefficient of the transmitted orders of
+           a binary grating given an incident planewave and verifies
+           that the results are the same when using either a band number
+           or `DiffractedPlanewave` object in `get_eigenmode_coefficients`.
+
+        Args:
+          gp: grating periodicity (μm).
+          gh: grating height (μm).
+          gdc: grating duty cycle (dimensionless).
+          theta: angle (in degrees) of incident planewave rotated counter
+                 clockwise (CCW) about z axis. 0 degrees along +x axis.
+        """
         sx = self.dpml + self.dsub + gh + self.dpad + self.dpml
         sy = gp
         cell_size = mp.Vector3(sx, sy, 0)
 
-        # rotation angle of incident planewave
-        # counter clockwise (CCW) about Z axis, 0 degrees along +X axis
-        theta_in = math.radians(theta)
-
-        # k (in source medium) with correct length (plane of incidence: XY)
-        k = mp.Vector3(self.fcen * self.ng).rotate(mp.Vector3(z=1), theta_in)
-
-        eig_parity = mp.ODD_Z
-        if theta == 0:
+        eig_parity = mp.ODD_Z  # S polarization
+        if theta_pw == 0:
             k = mp.Vector3()
             eig_parity += mp.EVEN_Y
             symmetries = [mp.Mirror(direction=mp.Y)]
         else:
+            # k (in source medium) with correct length (plane of incidence: xy)
+            k = mp.Vector3(self.fcen * self.ng).rotate(
+                mp.Vector3(z=1), math.radians(theta_pw)
+            )
             symmetries = []
 
         def pw_amp(k, x0):
@@ -100,14 +104,14 @@ def _pw_amp(x):
         m_plus = int(np.floor((self.fcen - k.y) * gp))
         m_minus = int(np.ceil((-self.fcen - k.y) * gp))
 
-        if theta == 0:
-            orders = range(m_plus)
+        if theta_pw == 0:
+            orders = range(m_plus + 1)
         else:
-            # ordering of the modes computed by MPB is according to *decreasing*
-            # values of kx (i.e., closest to propagation direction of 0° or +x)
             ms = range(m_minus, m_plus + 1)
             kx = lambda m: np.power(self.fcen, 2) - np.power(k.y + m / gp, 2)
             kxs = [kx(m) for m in ms]
+            # Ordering of the modes computed by MPB is according to *decreasing*
+            # values of kx (i.e., closest to propagation direction of 0° or +x)
             ids = np.flip(np.argsort(kxs))
             orders = [ms[d] for d in ids]
 
@@ -116,18 +120,20 @@ def _pw_amp(x):
                 tran_flux, [band + 1], eig_parity=eig_parity
             )
             tran_ref = abs(res.alpha[0, 0, 0]) ** 2
-            if theta == 0:
+            # For a planewave at normal incidence, the ±m diffracted orders
+            # contain identical power which means that to obtain the power in
+            # each order requires halving the total power. This applies
+            # to all orders except m=0.
+            if theta_pw == 0 and band != 0:
                 tran_ref = 0.5 * tran_ref
             vg_ref = res.vgrp[0]
             kdom_ref = res.kdom[0]
 
             res = sim.get_eigenmode_coefficients(
                 tran_flux,
-                mp.DiffractedPlanewave((0, order, 0), mp.Vector3(0, 1, 0), 1, 0),
+                mp.DiffractedPlanewave([0, order, 0], mp.Vector3(0, 1, 0), 1, 0),
             )
             tran_dp = abs(res.alpha[0, 0, 0]) ** 2
-            if (theta == 0) and (order == 0):
-                tran_dp = 0.5 * tran_dp
             vg_dp = res.vgrp[0]
             kdom_dp = res.kdom[0]
 
@@ -140,7 +146,7 @@ def _pw_amp(x):
 
             self.assertAlmostEqual(vg_ref, vg_dp, places=4)
             self.assertAlmostEqual(tran_ref, tran_dp, places=4)
-            if theta == 0:
+            if theta_pw == 0:
                 self.assertAlmostEqual(abs(kdom_ref.x), kdom_dp.x, places=5)
                 self.assertAlmostEqual(abs(kdom_ref.y), kdom_dp.y, places=5)
                 self.assertAlmostEqual(abs(kdom_ref.z), kdom_dp.z, places=5)
@@ -158,6 +164,120 @@ def test_diffracted_planewave(self):
         # self.run_binary_grating_diffraction(10.0,0.5,0.5,0)
         # self.run_binary_grating_diffraction(10.0,0.5,0.5,10.7)
 
+    def run_mode_source(self, gp, gh, gdc, m, use_diffpw):
+        """Computes the flux of a transmitted order of a binary grating
+           given an incident planewave. The incident planewave is defined
+           using a mode source specified using a band number and parity or
+           `DiffractedPlanewave` object.
+
+        Args:
+          gp: grating periodicity (μm).
+          gh: grating height (μm).
+          gdc: grating duty cycle (dimensionless).
+          m: diffraction order in y direction (integer).
+          use_diffpw: use a `DiffractedPlanewave` (True) or band number
+                      and parity (False).
+        """
+        sx = self.dpml + self.dsub + gh + self.dpad + self.dpml
+        sy = gp
+        cell_size = mp.Vector3(sx, sy, 0)
+
+        eig_parity = mp.ODD_Z  # S polarization
+        if m == 0:
+            k = mp.Vector3()
+            eig_parity += mp.EVEN_Y
+            symmetries = [mp.Mirror(direction=mp.Y)]
+        else:
+            ky = m / gp
+            theta_pw = math.asin(ky / (self.fcen * self.ng))
+            # k (in source medium) with correct length (plane of incidence: xy)
+            k = mp.Vector3(self.fcen * self.ng).rotate(mp.Vector3(z=1), theta_pw)
+            symmetries = []
+
+        if use_diffpw:
+            # The *zeroth* diffraction order specifies a planewave with a
+            # wavevector equal to the `k_point` of the `Simulation` object.
+            eigsrc_args = {
+                "eig_band": mp.DiffractedPlanewave(
+                    [0, 0, 0],
+                    mp.Vector3(0, 1, 0),
+                    1,
+                    0,
+                )
+            }
+        else:
+            eigsrc_args = {
+                "eig_band": 1,
+                "eig_parity": eig_parity,
+                "direction": mp.NO_DIRECTION,
+                "eig_kpoint": k,
+            }
+
+        src_pt = mp.Vector3(-0.5 * sx + self.dpml, 0, 0)
+        sources = [
+            mp.EigenModeSource(
+                mp.GaussianSource(self.fcen, fwidth=0.1 * self.fcen),
+                center=src_pt,
+                size=mp.Vector3(0, sy, 0),
+                **eigsrc_args,
+            )
+        ]
+
+        geometry = [
+            mp.Block(
+                material=self.glass,
+                size=mp.Vector3(self.dpml + self.dsub, mp.inf, mp.inf),
+                center=mp.Vector3(-0.5 * sx + 0.5 * (self.dpml + self.dsub), 0, 0),
+            ),
+            mp.Block(
+                material=self.glass,
+                size=mp.Vector3(gh, gdc * gp, mp.inf),
+                center=mp.Vector3(-0.5 * sx + self.dpml + self.dsub + 0.5 * gh, 0, 0),
+            ),
+        ]
+
+        sim = mp.Simulation(
+            resolution=self.resolution,
+            cell_size=cell_size,
+            boundary_layers=self.pml_layers,
+            k_point=k,
+            geometry=geometry,
+            sources=sources,
+            symmetries=symmetries,
+        )
+
+        tran_pt = mp.Vector3(0.5 * sx - self.dpml, 0, 0)
+        tran_flux = sim.add_flux(
+            self.fcen, 0, 1, mp.FluxRegion(center=tran_pt, size=mp.Vector3(0, sy, 0))
+        )
+
+        sim.run(
+            until_after_sources=mp.stop_when_fields_decayed(20, mp.Ez, src_pt, 1e-6)
+        )
+
+        tran = mp.get_fluxes(tran_flux)[0]
+
+        # force garbage collection
+        sim.reset_meep()
+
+        return tran
+
+    def test_mode_source(self):
+        gp = 1.5
+        gh = 0.5
+        gdc = 0.3
+        m = 2
+        tran_bn = self.run_mode_source(gp, gh, gdc, m, False)
+        tran_dp = self.run_mode_source(gp, gh, gdc, m, True)
+        print(
+            f"mode-source:, "
+            f"{tran_bn:.5f} (band number), "
+            f"{tran_dp:.5f} (diffraction order)"
+        )
+        self.assertAlmostEqual(
+            tran_bn, tran_dp, places=3 if mp.is_single_precision() else 4
+        )
+
 
 if __name__ == "__main__":
     unittest.main()