Add unit test for broadband sources at fixed angle in a periodic cell…

… (BFAST) (#2722)

* add unit test for broadband sources at fixed angle (BFAST)

* include index of source medium in bfast_k_bar

* set Courant factor using analytic expression involving bfast_k_bar

python/tests/test_refl_angular.py

@@ -1,4 +1,5 @@
 import math
+from typing import List, Tuple
 import unittest
 
 import numpy as np
@@ -8,45 +9,75 @@
 import meep as mp
 
 
-class TestReflAngular(ApproxComparisonTestCase):
+class TestReflectanceAngular(ApproxComparisonTestCase):
     @classmethod
     def setUpClass(cls):
-        cls.resolution = 400  # pixels/μm
+        cls.resolution = 200  # pixels/μm
 
         cls.n1 = 1.4  # refractive index of medium 1
         cls.n2 = 3.5  # refractive index of medium 2
 
-        cls.dpml = 1.0
-        cls.dz = 7.0
-        cls.sz = cls.dz + 2 * cls.dpml
-
-        cls.wvl_min = 0.4
-        cls.wvl_max = 0.8
-        cls.fmin = 1 / cls.wvl_max
-        cls.fmax = 1 / cls.wvl_min
-        cls.fcen = 0.5 * (cls.fmin + cls.fmax)
-        cls.df = cls.fmax - cls.fmin
-        cls.nfreq = 11
-
-    def refl_angular(self, theta):
-        theta_r = math.radians(theta)
-
-        # wavevector (in source medium); plane of incidence is XZ
-        k = (
-            mp.Vector3(0, 0, 1)
-            .rotate(mp.Vector3(0, 1, 0), theta_r)
-            .scale(self.n1 * self.fmin)
-        )
+        cls.t_pml = 1.0
+        cls.length_z = 7.0
+        cls.size_z = cls.length_z + 2 * cls.t_pml
+
+        cls.wavelength_min = 0.4
+        cls.wavelength_max = 0.8
+        cls.frequency_min = 1 / cls.wavelength_max
+        cls.frequency_max = 1 / cls.wavelength_min
+        cls.frequency_center = 0.5 * (cls.frequency_min + cls.frequency_max)
+        cls.frequency_width = cls.frequency_max - cls.frequency_min
+        cls.num_freq = 11
+
+    def reflectance_angular(
+        self, theta_deg: float, need_bfast: bool
+    ) -> Tuple[List, List, np.ndarray]:
+        """Computes properties of the incident and reflected planewave.
+
+        Args:
+          theta_deg: angle of incident planewave.
+          need_bfast: whether to use the same angle for the incident planewave
+            for all frequencies. If False, the incident angle is frequency
+            dependent.
+
+        Returns:
+          A 3-tuple comprising the frequencies of the incident planewave,
+          angles of the incident planewave, and the reflectance.
+        """
+        theta_rad = math.radians(theta_deg)
+
+        if need_bfast:
+            bfast_k_bar = (self.n1 * np.sin(theta_rad), 0, 0)
+
+            Courant = (1 - bfast_k_bar[0]) / 3**0.5
+
+            k = mp.Vector3()
+        else:
+            bfast_k_bar = (0, 0, 0)
+
+            Courant = 0.5
+
+            # Wavevector in source medium with refractive index n1.
+            # Plane of incidence is XZ. Rotation is counter clockwise about
+            # Y axis. A rotation angle of zero is the +Z axis.
+            k = (
+                mp.Vector3(0, 0, 1)
+                .rotate(mp.Vector3(0, 1, 0), theta_rad)
+                .scale(self.n1 * self.frequency_min)
+            )
 
-        dimensions = 1 if theta == 0 else 3
-        cell_size = mp.Vector3(z=self.sz)
-        pml_layers = [mp.PML(self.dpml)]
+        dimensions = 1 if theta_deg == 0 else 3
+        cell_size = mp.Vector3(z=self.size_z)
+        pml_layers = [mp.PML(self.t_pml)]
+
+        # P polarization.
+        source_component = mp.Ex
 
         sources = [
             mp.Source(
-                mp.GaussianSource(self.fcen, fwidth=self.df),
-                component=mp.Ex,  # P polarization
-                center=mp.Vector3(z=-0.5 * self.sz + self.dpml),
+                mp.GaussianSource(self.frequency_center, fwidth=self.frequency_width),
+                component=source_component,
+                center=mp.Vector3(z=-0.5 * self.size_z + self.t_pml),
             )
         ]
 
@@ -58,25 +89,31 @@ def refl_angular(self, theta):
             sources=sources,
             boundary_layers=pml_layers,
             k_point=k,
+            need_bfast=need_bfast,
+            bfast_k_bar=bfast_k_bar,
+            Courant=Courant,
         )
 
-        mon_pt = -0.5 * self.sz + self.dpml + 0.25 * self.dz
-        refl_fr = mp.FluxRegion(center=mp.Vector3(z=mon_pt))
-        refl = sim.add_flux(self.fcen, self.df, self.nfreq, refl_fr)
+        monitor_point = -0.5 * self.size_z + self.t_pml + 0.25 * self.length_z
+        monitor_region = mp.FluxRegion(center=mp.Vector3(z=monitor_point))
+        flux_monitor = sim.add_flux(
+            self.frequency_center, self.frequency_width, self.num_freq, monitor_region
+        )
 
-        termination_cond = mp.stop_when_fields_decayed(
-            50, mp.Ex, mp.Vector3(z=mon_pt), 1e-9
+        termination_criteria = mp.stop_when_fields_decayed(
+            50, source_component, mp.Vector3(z=monitor_point), 1e-6
         )
-        sim.run(until_after_sources=termination_cond)
+        sim.run(until_after_sources=termination_criteria)
+
+        empty_data = sim.get_flux_data(flux_monitor)
+        empty_flux = mp.get_fluxes(flux_monitor)
 
-        empty_data = sim.get_flux_data(refl)
-        empty_flux = mp.get_fluxes(refl)
         sim.reset_meep()
 
         geometry = [
             mp.Block(
-                size=mp.Vector3(mp.inf, mp.inf, 0.5 * self.sz),
-                center=mp.Vector3(z=0.25 * self.sz),
+                size=mp.Vector3(mp.inf, mp.inf, 0.5 * self.size_z),
+                center=mp.Vector3(z=0.25 * self.size_z),
                 material=mp.Medium(index=self.n2),
             )
         ]
@@ -89,58 +126,88 @@ def refl_angular(self, theta):
             sources=sources,
             boundary_layers=pml_layers,
             k_point=k,
+            need_bfast=need_bfast,
+            bfast_k_bar=bfast_k_bar,
+            Courant=Courant,
             geometry=geometry,
         )
 
-        refl = sim.add_flux(self.fcen, self.df, self.nfreq, refl_fr)
-        sim.load_minus_flux_data(refl, empty_data)
-
-        sim.run(until_after_sources=termination_cond)
-
-        refl_flux = mp.get_fluxes(refl)
-        freqs = mp.get_flux_freqs(refl)
-
-        Rs = -np.array(refl_flux) / np.array(empty_flux)
-
-        thetas = [math.asin(k.x / (self.n1 * freqs[i])) for i in range(self.nfreq)]
-        return freqs, thetas, Rs
-
-    @parameterized.parameterized.expand([(0,), (20.6,)])
-    def test_refl_angular(self, theta):
-        fmeep, tmeep, Rmeep = self.refl_angular(theta)
-
-        # angle of refracted planewave in medium n2 for an
-        # incident planewave in medium n1 at angle theta_in
-        theta_out = lambda theta_in: math.asin(self.n1 * math.sin(theta_in) / self.n2)
+        flux_monitor = sim.add_flux(
+            self.frequency_center, self.frequency_width, self.num_freq, monitor_region
+        )
+        sim.load_minus_flux_data(flux_monitor, empty_data)
+
+        sim.run(until_after_sources=termination_criteria)
+
+        flux_monitor_flux = mp.get_fluxes(flux_monitor)
+        freqs = mp.get_flux_freqs(flux_monitor)
+
+        reflectance = -np.array(flux_monitor_flux) / np.array(empty_flux)
+
+        if need_bfast:
+            theta_in_rad = [theta_rad] * self.num_freq
+        else:
+            # Returns the angle of the incident planewave in medium n1 based
+            # on its frequency given a fixed wavevector component in X.
+            theta_in_rad = [
+                math.asin(k.x / (self.n1 * freqs[i])) for i in range(self.num_freq)
+            ]
+
+        return freqs, theta_in_rad, reflectance
+
+    @parameterized.parameterized.expand([(0, False), (20.6, False), (35.7, True)])
+    def test_reflectance_angular(self, theta_deg: float, need_bfast: bool):
+        (
+            frequency_meep,
+            theta_in_rad_meep,
+            reflectance_meep,
+        ) = self.reflectance_angular(theta_deg, need_bfast)
+
+        # Returns angle of refracted planewave in medium n2 given
+        # an incident planewave in medium n1 at angle theta_in_rad.
+        theta_out = lambda theta_in_rad: math.asin(
+            self.n1 * math.sin(theta_in_rad) / self.n2
+        )
 
-        # Fresnel reflectance for P polarization in medium n2 for
-        # an incident planewave in medium n1 at angle theta_in
-        Rfresnel = lambda theta_in: (
+        # Returns Fresnel reflectance for P polarization in medium n2
+        # for an incident planewave in medium n1 at angle theta_in_rad.
+        reflectance_fresnel = lambda theta_in_rad: (
             math.fabs(
-                (self.n1 * math.cos(theta_out(theta_in)) - self.n2 * math.cos(theta_in))
+                (
+                    self.n1 * math.cos(theta_out(theta_in_rad))
+                    - self.n2 * math.cos(theta_in_rad)
+                )
                 / (
-                    self.n1 * math.cos(theta_out(theta_in))
-                    + self.n2 * math.cos(theta_in)
+                    self.n1 * math.cos(theta_out(theta_in_rad))
+                    + self.n2 * math.cos(theta_in_rad)
                 )
             )
             ** 2
         )
 
-        Ranalytic = np.empty((self.nfreq,))
+        reflectance_analytic = np.empty((self.num_freq,))
         print(
-            "refl:, wavelength (μm), incident angle (°), reflectance (Meep), reflectance (analytic), error"
+            "refl:, wavelength (μm), incident angle (°), reflectance (Meep), "
+            "reflectance (analytic), error"
         )
-        for i in range(self.nfreq):
-            Ranalytic[i] = Rfresnel(tmeep[i])
-            err = abs(Rmeep[i] - Ranalytic[i]) / Ranalytic[i]
+        for i in range(self.num_freq):
+            reflectance_analytic[i] = reflectance_fresnel(theta_in_rad_meep[i])
+            err = (
+                abs(reflectance_meep[i] - reflectance_analytic[i])
+                / reflectance_analytic[i]
+            )
             print(
                 "refl:, {:4.2f}, {:4.2f}, {:8.6f}, {:8.6f}, {:6.4f}".format(
-                    1 / fmeep[i], math.degrees(tmeep[i]), Rmeep[i], Ranalytic[i], err
+                    1 / frequency_meep[i],
+                    math.degrees(theta_in_rad_meep[i]),
+                    reflectance_meep[i],
+                    reflectance_analytic[i],
+                    err,
                 )
             )
 
-        tol = 0.005 if mp.is_single_precision() else 0.004
-        self.assertClose(Rmeep, Ranalytic, epsilon=tol)
+        tol = 0.03
+        self.assertClose(reflectance_meep, reflectance_analytic, epsilon=tol)
 
 
 if __name__ == "__main__":