From ea8b0dd94b812564199b470ead040b75a72aa112 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Thu, 24 Mar 2022 05:59:43 +0000 Subject: [PATCH] revamp test_refl_angular.py to validate results against theory rather than hard-coded values --- doc/docs/Python_Tutorials/Basics.md | 15 ++- python/examples/refl-angular.py | 22 ++-- python/tests/test_refl_angular.py | 195 ++++++++++++++-------------- 3 files changed, 123 insertions(+), 109 deletions(-) diff --git a/doc/docs/Python_Tutorials/Basics.md b/doc/docs/Python_Tutorials/Basics.md index cedab140f..8d2fab091 100644 --- a/doc/docs/Python_Tutorials/Basics.md +++ b/doc/docs/Python_Tutorials/Basics.md @@ -433,7 +433,7 @@ A 1d cell must be used since a higher-dimensional cell will introduce [artificia Creating an oblique planewave source typically requires specifying two parameters: (1) for periodic structures, the Bloch-periodic wavevector $\vec{k}$ via [`k_point`](../FAQ.md#how-does-k_point-define-the-phase-relation-between-adjacent-unit-cells), and (2) the source amplitude function `amp_func` for setting the $e^{i\vec{k} \cdot \vec{r}}$ spatial dependence ($\vec{r}$ is the position vector). Since we have a 1d cell and the source is at a single point, it is not necessary to specify the source amplitude (see this [2d example](https://github.com/NanoComp/meep/blob/master/python/examples/pw-source.py) for how this is done). The magnitude of the Bloch-periodic wavevector is specified according to the dispersion relation formula for a planewave in homogeneous media with index $n$: $\omega=c|\vec{k}|/n$. As the source in this example is incident from air, $|\vec{k}|$ is simply equal to the frequency $\omega$. Note that specifying $\vec{k}$ corresponds to a single frequency. Any broadband source is therefore incident at a specified angle for only a *single* frequency. This is described in more detail in Section 4.5 ("Efficient Frequency-Angle Coverage") in [Chapter 4](https://arxiv.org/abs/1301.5366) ("Electromagnetic Wave Source Conditions") of the book [Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707). In this example, $\omega$ is set to the minimum frequency of the pulse to produce propagating fields at all pulse frequencies. In general, any pulse frequencies which are *less* than any non-zero component of $\vec{k}$ will result in *evanescent* fields (which carry zero power to the far field, and computationally will yield exponentially small power). -In this example, the plane of incidence which contains $\vec{k}$ and the surface normal vector is chosen to be $xz$. The source angle $\theta$ is defined in degrees in the counterclockwise (CCW) direction around the $y$ axis with 0 degrees along the $+z$ axis. In Meep, a 1d cell is defined along the $z$ direction. When $\vec{k}$ is not set, only the $E_x$ and $H_y$ field components are permitted. A non-zero $\vec{k}$ results in a 3d simulation where all field components are included and are complex valued (note that the fields are real, by default). A current source with $E_x$ polarization lies within the plane of incidence and corresponds to the convention of $\mathcal{P}$-polarization. In order to model the $\mathcal{S}$-polarization, we must use an $E_y$ source. This example involves just the $\mathcal{P}$-polarization. +In this example, the plane of incidence which contains $\vec{k}$ and the surface normal vector is chosen to be $xz$. The source angle $\theta$ is defined in degrees in the counterclockwise (CCW) direction around the $y$ axis with 0 degrees along the $+z$ axis. In Meep, a 1d cell is defined along the $z$ direction. When $\vec{k}$ is not set, only the $E_x$ and $H_y$ field components are permitted. A non-zero $\vec{k}$ results in a 3d simulation where all field components are included and are complex valued (note that the fields are real, by default). A current source with $E_x$ polarization lies within the plane of incidence and corresponds to the convention of $\mathcal{P}$ polarization. In order to model the $\mathcal{S}$ polarization, we must use an $E_y$ source. This example involves just the $\mathcal{P}$ polarization. The simulation script is [examples/refl-angular.py](https://github.com/NanoComp/meep/blob/master/python/examples/refl-angular.py). The notebook is [examples/refl-angular.ipynb](https://nbviewer.jupyter.org/github/NanoComp/meep/blob/master/python/examples/refl-angular.ipynb). @@ -442,13 +442,12 @@ import meep as mp import argparse import math -def main(args): +def main(args): resolution = args.res dpml = 1.0 # PML thickness - sz = 10 # size of cell (without PMLs) - sz = 10 + 2*dpml + sz = 10 + 2*dpml # size of cell (without PMLs) cell_size = mp.Vector3(0,0,sz) pml_layers = [mp.PML(dpml)] @@ -472,7 +471,9 @@ def main(args): else: dimensions = 3 - sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df), component=mp.Ex, center=mp.Vector3(0,0,-0.5*sz+dpml))] + sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df), + component=mp.Ex, + center=mp.Vector3(0,0,-0.5*sz+dpml))] sim = mp.Simulation(cell_size=cell_size, boundary_layers=pml_layers, @@ -491,7 +492,9 @@ def main(args): sim.reset_meep() # add a block with n=3.5 for the air-dielectric interface - geometry = [mp.Block(mp.Vector3(mp.inf,mp.inf,0.5*sz), center=mp.Vector3(0,0,0.25*sz), material=mp.Medium(index=3.5))] + geometry = [mp.Block(size=mp.Vector3(mp.inf,mp.inf,0.5*sz), + center=mp.Vector3(0,0,0.25*sz), + material=mp.Medium(index=3.5))] sim = mp.Simulation(cell_size=cell_size, geometry=geometry, diff --git a/python/examples/refl-angular.py b/python/examples/refl-angular.py index c08bceea7..a51791cc9 100644 --- a/python/examples/refl-angular.py +++ b/python/examples/refl-angular.py @@ -2,12 +2,12 @@ import argparse import math + def main(args): resolution = args.res dpml = 1.0 # PML thickness - sz = 10 # size of computational cell (without PMLs) - sz = 10 + 2*dpml + sz = 10 + 2*dpml # size of computational cell (without PMLs) cell_size = mp.Vector3(0,0,sz) pml_layers = [mp.PML(dpml)] @@ -18,20 +18,22 @@ def main(args): fcen = 0.5*(fmin+fmax) # center frequency df = fmax-fmin # frequency width nfreq = 50 # number of frequency bins - + # rotation angle (in degrees) of source: CCW around Y axis, 0 degrees along +Z axis theta_r = math.radians(args.theta) # plane of incidence is XZ k = mp.Vector3(math.sin(theta_r),0,math.cos(theta_r)).scale(fmin) - + # if normal incidence, force number of dimensions to be 1 if theta_r == 0: dimensions = 1 else: dimensions = 3 - - sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df), component=mp.Ex, center=mp.Vector3(0,0,-0.5*sz+dpml))] + + sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df), + component=mp.Ex, + center=mp.Vector3(0,0,-0.5*sz+dpml))] sim = mp.Simulation(cell_size=cell_size, boundary_layers=pml_layers, @@ -42,7 +44,7 @@ def main(args): refl_fr = mp.FluxRegion(center=mp.Vector3(0,0,-0.25*sz)) refl = sim.add_flux(fcen, df, nfreq, refl_fr) - + sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ex, mp.Vector3(0,0,-0.5*sz+dpml), 1e-9)) empty_flux = mp.get_fluxes(refl) @@ -50,7 +52,9 @@ def main(args): sim.reset_meep() # add a block with n=3.5 for the air-dielectric interface - geometry = [mp.Block(mp.Vector3(mp.inf,mp.inf,0.5*sz), center=mp.Vector3(0,0,0.25*sz), material=mp.Medium(index=3.5))] + geometry = [mp.Block(size=mp.Vector3(mp.inf,mp.inf,0.5*sz), + center=mp.Vector3(0,0,0.25*sz), + material=mp.Medium(index=3.5))] sim = mp.Simulation(cell_size=cell_size, geometry=geometry, @@ -70,7 +74,7 @@ def main(args): for i in range(nfreq): print("refl:, {}, {}, {}, {}".format(k.x,1/freqs[i],math.degrees(math.asin(k.x/freqs[i])),-refl_flux[i]/empty_flux[i])) - + if __name__ == '__main__': parser = argparse.ArgumentParser() parser.add_argument('-res', type=int, default=200, help='resolution (default: 200 pixels/um)') diff --git a/python/tests/test_refl_angular.py b/python/tests/test_refl_angular.py index 209dd16eb..ba89eed29 100644 --- a/python/tests/test_refl_angular.py +++ b/python/tests/test_refl_angular.py @@ -1,127 +1,134 @@ -import meep as mp -from utils import ApproxComparisonTestCase -import math +# -*- coding: utf-8 -*- + import unittest +import parameterized import numpy as np +import math +import meep as mp +from utils import ApproxComparisonTestCase class TestReflAngular(ApproxComparisonTestCase): + @classmethod + def setUpClass(cls): + cls.resolution = 400 # pixels/μm - def test_refl_angular(self): - resolution = 100 + cls.n1 = 1.4 # refractive index of medium 1 + cls.n2 = 3.5 # refractive index of medium 2 - dpml = 1.0 - sz = 10 - sz = sz + 2 * dpml - cell_size = mp.Vector3(z=sz) - pml_layers = [mp.PML(dpml)] + cls.dpml = 1.0 + cls.dz = 7.0 + cls.sz = cls.dz+2*cls.dpml - wvl_min = 0.4 - wvl_max = 0.8 - fmin = 1 / wvl_max - fmax = 1 / wvl_min - fcen = 0.5 * (fmin + fmax) - df = fmax - fmin - nfreq = 50 + 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 - theta_r = math.radians(0) + def refl_angular(self, theta): + theta_r = math.radians(theta) - k = mp.Vector3(math.sin(theta_r), 0, math.cos(theta_r)).scale(fcen) + # 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) - dimensions = 1 + if theta == 0: + dimensions = 1 + else: + dimensions = 3 - sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ex, - center=mp.Vector3(z=-0.5 * sz + dpml))] + cell_size = mp.Vector3(z=self.sz) + pml_layers = [mp.PML(self.dpml)] - sim = mp.Simulation(cell_size=cell_size, - boundary_layers=pml_layers, - sources=sources, - k_point=k, + 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))] + + sim = mp.Simulation(resolution=self.resolution, + cell_size=cell_size, dimensions=dimensions, - resolution=resolution) + default_material=mp.Medium(index=self.n1), + sources=sources, + boundary_layers=pml_layers, + k_point=k) - refl_fr = mp.FluxRegion(center=mp.Vector3(z=-0.25 * sz)) - refl = sim.add_flux(fcen, df, nfreq, refl_fr, decimation_factor=1) + 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) - sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ex, mp.Vector3(z=-0.5 * sz + dpml), 1e-9)) + termination_cond = mp.stop_when_fields_decayed(50, + mp.Ex, + mp.Vector3(z=mon_pt), + 1e-9) + sim.run(until_after_sources=termination_cond) empty_data = sim.get_flux_data(refl) + empty_flux = mp.get_fluxes(refl) sim.reset_meep() - geometry = [mp.Block(mp.Vector3(mp.inf, mp.inf, 0.5 * sz), center=mp.Vector3(z=0.25 * sz), - material=mp.Medium(index=3.5))] + geometry = [mp.Block(size=mp.Vector3(mp.inf,mp.inf,0.5*self.sz), + center=mp.Vector3(z=0.25*self.sz), + material=mp.Medium(index=self.n2))] - sim = mp.Simulation(cell_size=cell_size, - geometry=geometry, - boundary_layers=pml_layers, + sim = mp.Simulation(resolution=self.resolution, + cell_size=cell_size, + dimensions=dimensions, + default_material=mp.Medium(index=self.n1), sources=sources, + boundary_layers=pml_layers, k_point=k, - dimensions=dimensions, - resolution=resolution) + geometry=geometry) - refl = sim.add_flux(fcen, df, nfreq, refl_fr, decimation_factor=1) + 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=mp.stop_when_fields_decayed(50, mp.Ex, mp.Vector3(z=-0.5 * sz + dpml), 1e-9)) + sim.run(until_after_sources=termination_cond) refl_flux = mp.get_fluxes(refl) freqs = mp.get_flux_freqs(refl) - expected = [ - (1.25, -1.123696883299492e-6), - (1.2755102040816326, -2.5749667658387866e-6), - (1.3010204081632653, -5.70480204599006e-6), - (1.3265306122448979, -1.2220464827582253e-5), - (1.3520408163265305, -2.531247480206961e-5), - (1.3775510204081631, -5.069850309492639e-5), - (1.4030612244897958, -9.819256552437341e-5), - (1.4285714285714284, -1.8390448659017395e-4), - (1.454081632653061, -3.330762066794769e-4), - (1.4795918367346936, -5.833650417163753e-4), - (1.5051020408163263, -9.8807834237052e-4), - (1.5306122448979589, -0.001618472171445976), - (1.5561224489795915, -0.0025638388059825985), - (1.5816326530612241, -0.003927863989816029), - (1.6071428571428568, -0.005819831283556752), - (1.6326530612244894, -0.008339881000982728), - (1.658163265306122, -0.011558769654206626), - (1.6836734693877546, -0.015494308354153143), - (1.7091836734693873, -0.02008850084337135), - (1.73469387755102, -0.025190871516857616), - (1.7602040816326525, -0.030553756123198477), - (1.7857142857142851, -0.03584404966066722), - (1.8112244897959178, -0.040672967700428275), - (1.8367346938775504, -0.04464118393086191), - (1.862244897959183, -0.047392712128477496), - (1.8877551020408156, -0.048667403362887635), - (1.9132653061224483, -0.048341494285878264), - (1.938775510204081, -0.04644739000778679), - (1.9642857142857135, -0.043168390293742316), - (1.9897959183673462, -0.0388094755730579), - (2.0153061224489788, -0.03375052221907117), - (2.0408163265306114, -0.02839209067703472), - (2.066326530612244, -0.023104245646230648), - (2.0918367346938767, -0.01818725699718267), - (2.1173469387755093, -0.013849270759480073), - (2.142857142857142, -0.010201733597436358), - (2.1683673469387745, -0.007269616609175294), - (2.193877551020407, -0.005011210495189995), - (2.21938775510204, -0.0033417192031464896), - (2.2448979591836724, -0.0021557351734376256), - (2.270408163265305, -0.0013453062176115673), - (2.2959183673469377, -8.121742663131631e-4), - (2.3214285714285703, -4.7433135191915683e-4), - (2.346938775510203, -2.6799188013374266e-4), - (2.3724489795918355, -1.464781343401766e-4), - (2.397959183673468, -7.745339273024636e-5), - (2.423469387755101, -3.9621374769542025e-5), - (2.4489795918367334, -1.9608458558430508e-5), - (2.474489795918366, -9.38818477949983e-6), - (2.4999999999999987, -4.3484671364929225e-6), - ] - - tol = 1e-7 if mp.is_single_precision() else 1e-8 - self.assertClose(expected, list(zip(freqs, refl_flux)), epsilon=tol) + Rs = -np.array(refl_flux)/np.array(empty_flux) + + thetas = [] + for i in range(self.nfreq): + thetas.append(math.asin(k.x/(self.n1*freqs[i]))) + + 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) + + # Fresnel reflectance for P polarization in medium n2 for + # an incident planewave in medium n1 at angle theta_in + Rfresnel = lambda theta_in: (math.fabs((self.n1*math.cos(theta_out(theta_in))-self.n2*math.cos(theta_in)) / + (self.n1*math.cos(theta_out(theta_in))+self.n2*math.cos(theta_in)))**2) + + Ranalytic = np.empty((self.nfreq,)) + print("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] + 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)) + + tol = 0.004 + self.assertClose(Rmeep, Ranalytic, epsilon=tol) if __name__ == '__main__':