From 5cf040cc08278b327027cd6f96b897208fccce9c Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Sun, 19 Nov 2023 19:24:11 -0800 Subject: [PATCH] add unit test for broadband sources at fixed angle (BFAST) --- python/tests/test_refl_angular.py | 214 +++++++++++++++++++++--------- 1 file changed, 150 insertions(+), 64 deletions(-) diff --git a/python/tests/test_refl_angular.py b/python/tests/test_refl_angular.py index 6e6bad0b6..ee14df981 100644 --- a/python/tests/test_refl_angular.py +++ b/python/tests/test_refl_angular.py @@ -1,4 +1,5 @@ import math +from typing import Tuple import unittest import numpy as np @@ -8,7 +9,7 @@ import meep as mp -class TestReflAngular(ApproxComparisonTestCase): +class TestReflectanceAngular(ApproxComparisonTestCase): @classmethod def setUpClass(cls): cls.resolution = 400 # pixels/μm @@ -18,35 +19,76 @@ def setUpClass(cls): 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.size_z = cls.dz + 2 * cls.dpml + + 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[np.ndarray, np.ndarray, np.ndarray]: + """Computes properties of the reflected planewave from an interface. + + 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 as 1D arrays. + """ + theta_rad = math.radians(theta_deg) + + if need_bfast: + bfast_k_bar = (np.sin(theta_rad), 0, 0) + + k = mp.Vector3() + + if theta_deg > 40.0: + # TODO(oskooi): determine whether a lower value is necessary. + Courant = 0.05 + else: + Courant = 0.5 + + else: + bfast_k_bar = (0, 0, 0) + + # 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) + ) + + Courant = 0.5 - dimensions = 1 if theta == 0 else 3 - cell_size = mp.Vector3(z=self.sz) + dimensions = 1 if theta_deg == 0 else 3 + cell_size = mp.Vector3(z=self.size_z) pml_layers = [mp.PML(self.dpml)] + # 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.dpml), ) ] @@ -58,25 +100,34 @@ 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.dpml + 0.25 * self.dz + 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-9 ) - 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 +140,93 @@ def refl_angular(self, theta): sources=sources, boundary_layers=pml_layers, k_point=k, - geometry=geometry, + 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) + 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) - sim.run(until_after_sources=termination_cond) + flux_monitor_flux = mp.get_fluxes(flux_monitor) + freqs = mp.get_flux_freqs(flux_monitor) - refl_flux = mp.get_fluxes(refl) - freqs = mp.get_flux_freqs(refl) + reflectance = -np.array(flux_monitor_flux) / np.array(empty_flux) - Rs = -np.array(refl_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) + ] - thetas = [math.asin(k.x / (self.n1 * freqs[i])) for i in range(self.nfreq)] - return freqs, thetas, Rs + return freqs, theta_in_rad, reflectance - @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) + @parameterized.parameterized.expand([(0, True), (20.6, 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) - # Fresnel reflectance for P polarization in medium n2 for - # an incident planewave in medium n1 at angle theta_in - Rfresnel = lambda theta_in: ( + # 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 + ) + + # 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) + if need_bfast: + tol = 0.1 + else: + tol = 0.005 if mp.is_single_precision() else 0.004 + + self.assertClose(reflectance_meep, reflectance_analytic, epsilon=tol) if __name__ == "__main__":