diff --git a/python/tests/test_refl_angular.py b/python/tests/test_refl_angular.py index 6e6bad0b6..eb86f953e 100644 --- a/python/tests/test_refl_angular.py +++ b/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__":