From c4e6e21de45b6bd6cfc9c376d09fdd55ac178a90 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Fri, 13 May 2022 21:05:33 +0000 Subject: [PATCH 1/3] tutorial and unit test for mode source using DiffractedPlanewave object --- doc/docs/Python_Tutorials/Eigenmode_Source.md | 12 + python/tests/test_diffracted_planewave.py | 217 +++++++++++++----- 2 files changed, 169 insertions(+), 60 deletions(-) diff --git a/doc/docs/Python_Tutorials/Eigenmode_Source.md b/doc/docs/Python_Tutorials/Eigenmode_Source.md index c54140087..a86cb139e 100644 --- a/doc/docs/Python_Tutorials/Eigenmode_Source.md +++ b/doc/docs/Python_Tutorials/Eigenmode_Source.md @@ -288,3 +288,15 @@ Shown below are the steady-state field profiles generated by the planewave for t
![](../images/eigenmode_planewave.png)
+ +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(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))] +``` diff --git a/python/tests/test_diffracted_planewave.py b/python/tests/test_diffracted_planewave.py index 962a3e3ea..a3da7d06b 100644 --- a/python/tests/test_diffracted_planewave.py +++ b/python/tests/test_diffracted_planewave.py @@ -7,39 +7,65 @@ class TestDiffractedPlanewave(unittest.TestCase): - def run_binary_grating_diffraction(self, gp, gh, gdc, theta, bands, orders): - resolution = 50 # pixels/um - - dpml = 1.0 # PML thickness - dsub = 3.0 # substrate thickness - dpad = 3.0 # length of padding between grating and PML - - sx = dpml+dsub+gh+dpad+dpml - sy = gp - - cell_size = mp.Vector3(sx,sy,0) - absorber_layers = [mp.Absorber(thickness=dpml,direction=mp.X)] - - wvl = 0.5 # center wavelength - fcen = 1/wvl # center frequency - df = 0.05*fcen # frequency width - - ng = 1.5 - glass = mp.Medium(index=ng) - + @classmethod + def setUp(cls): + cls.resolution = 50 # pixels/μm + + cls.dpml = 1.0 # PML thickness + cls.dsub = 3.0 # substrate thickness + cls.dpad = 3.0 # length of padding between grating and PML + cls.gp = 2.6 # grating period + cls.gh = 0.4 # grating height + cls.gdc = 0.5 # grating duty cycle + + cls.sx = cls.dpml+cls.dsub+cls.gh+cls.dpad+cls.dpml + cls.sy = cls.gp + + cls.cell_size = mp.Vector3(cls.sx,cls.sy,0) + cls.absorber_layers = [mp.PML(thickness=cls.dpml,direction=mp.X)] + + cls.wvl = 0.5 # center wavelength + cls.fcen = 1/cls.wvl # center frequency + cls.df = 0.05*cls.fcen # frequency width + + cls.ng = 1.5 + cls.glass = mp.Medium(index=cls.ng) + + cls.eig_parity = mp.ODD_Z + + cls.src_pt = mp.Vector3(-0.5*cls.sx+cls.dpml,0,0) + cls.tran_pt = mp.Vector3(0.5*cls.sx-cls.dpml,0,0) + + cls.geometry = [mp.Block(material=cls.glass, + size=mp.Vector3(cls.dpml+cls.dsub,mp.inf,mp.inf), + center=mp.Vector3(-0.5*cls.sx+0.5*(cls.dpml+cls.dsub),0,0)), + mp.Block(material=cls.glass, + size=mp.Vector3(cls.gh,cls.gdc*cls.gp,mp.inf), + center=mp.Vector3(-0.5*cls.sx+cls.dpml+cls.dsub+0.5*cls.gh,0,0))] + + cls.tol = 1e-8 # tolerance for mp.stop_dft_decayed termination criteria + + + def run_mode_decomposition(self,theta,bands,orders): + """Computes the transmittance of the diffraction orders + of a binary grating given an incident planewave at + angle theta. The transmittance is computed using + mode decomposition in two different ways: + (1) band number (bands) and (2) diffraction order (orders) + via a `DiffractedPlanewave` object. The test verifies + that these two methods produce equivalent results. + """ # rotation angle of incident planewave # counter clockwise (CCW) about Z axis, 0 degrees along +X axis theta_in = math.radians(theta) - eig_parity = mp.ODD_Z - # k (in source medium) with correct length (plane of incidence: XY) - k = mp.Vector3(fcen*ng).rotate(mp.Vector3(z=1), theta_in) + k = mp.Vector3(self.fcen*self.ng).rotate(mp.Vector3(z=1),theta_in) symmetries = [] if theta_in == 0: k = mp.Vector3() - eig_parity += mp.EVEN_Y + self.eig_parity += mp.EVEN_Y symmetries = [mp.Mirror(direction=mp.Y)] def pw_amp(k,x0): @@ -47,61 +73,52 @@ def _pw_amp(x): return cmath.exp(1j*2*math.pi*k.dot(x+x0)) return _pw_amp - src_pt = mp.Vector3(-0.5*sx+dpml,0,0) - sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df), + sources = [mp.Source(mp.GaussianSource(self.fcen,fwidth=self.df), component=mp.Ez, - center=src_pt, - size=mp.Vector3(0,sy,0), - amp_func=pw_amp(k,src_pt))] + center=self.src_pt, + size=mp.Vector3(0,self.sy,0), + amp_func=pw_amp(k,self.src_pt))] - sim = mp.Simulation(resolution=resolution, - cell_size=cell_size, - boundary_layers=absorber_layers, + sim = mp.Simulation(resolution=self.resolution, + cell_size=self.cell_size, + boundary_layers=self.absorber_layers, k_point=k, - default_material=glass, + default_material=self.glass, sources=sources, symmetries=symmetries) - tran_pt = mp.Vector3(0.5*sx-dpml,0,0) - tran_flux = sim.add_flux(fcen, + tran_flux = sim.add_flux(self.fcen, 0, 1, - mp.FluxRegion(center=tran_pt, - size=mp.Vector3(0,sy,0))) + mp.FluxRegion(center=self.tran_pt, + size=mp.Vector3(0,self.sy,0))) - sim.run(until_after_sources=10) + sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=self.tol)) input_flux = mp.get_fluxes(tran_flux) sim.reset_meep() - geometry = [mp.Block(material=glass, - size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), - center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub),0,0)), - mp.Block(material=glass, - size=mp.Vector3(gh,gdc*gp,mp.inf), - center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh,0,0))] - - sim = mp.Simulation(resolution=resolution, - cell_size=cell_size, - boundary_layers=absorber_layers, - geometry=geometry, + sim = mp.Simulation(resolution=self.resolution, + cell_size=self.cell_size, + boundary_layers=self.absorber_layers, + geometry=self.geometry, k_point=k, sources=sources, symmetries=symmetries) - tran_flux = sim.add_mode_monitor(fcen, + tran_flux = sim.add_mode_monitor(self.fcen, 0, 1, - mp.FluxRegion(center=tran_pt, - size=mp.Vector3(0,sy,0))) + mp.FluxRegion(center=self.tran_pt, + size=mp.Vector3(0,self.sy,0))) - sim.run(until_after_sources=100) + sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=self.tol)) for band,order in zip(bands,orders): res = sim.get_eigenmode_coefficients(tran_flux, [band], - eig_parity=eig_parity) + eig_parity=self.eig_parity) tran_ref = abs(res.alpha[0,0,0])**2/input_flux[0] if (theta_in == 0): tran_ref = 0.5*tran_ref @@ -121,17 +138,97 @@ def _pw_amp(x): vg_dp = res.vgrp[0] err = abs(tran_ref-tran_dp)/tran_ref - print("tran:, {:2d} (band), {:2d} (order), {:10.8f} (eigensolver), {:10.8f} (planewave), {:10.8f} (error)".format(band,order,tran_ref,tran_dp,err)) + print("info:, {:2d} (band number), {:2d} (diffraction order), " + "{:10.8f} [trans. (band number)], " + "{:10.8f} [trans. (diffracted order)], " + "{:10.8f} (error)".format(band, + order, + tran_ref, + tran_dp, + err)) self.assertAlmostEqual(vg_ref,vg_dp,places=5) self.assertAlmostEqual(tran_ref,tran_dp,places=5) - def test_diffracted_planewave(self): - self.run_binary_grating_diffraction(2.6,0.4,0.6,0,range(1,6),range(0,5)) - self.run_binary_grating_diffraction(2.6,0.4,0.6,13.4,range(1,6),[-2,-1,-3,0,-4]) - # self.run_binary_grating_diffraction(10.0,0.5,0.5,0,[2,4,6],[1,3,5]) - # self.run_binary_grating_diffraction(10.0,0.5,0.5,10.7,[1,4,8],[-6,-4,-2]) + def run_mode_source(self,m,diffpw): + """Computes the transmitted Poynting flux of a + binary grating given an incident planewave + specified by the diffraction order m in the + y direction. The incident planewave is defined + using a mode source with either a band number + or `DiffractedPlanewave` object. + """ + ky = m/self.gp + theta = 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) + + symmetries = [] + if theta == 0: + k = mp.Vector3() + self.eig_parity += mp.EVEN_Y + symmetries = [mp.Mirror(direction=mp.Y)] + + if diffpw: + # the *zeroth* diffraction order defines a planewave with a wavevector + # equal to the `k_point` of the `Simulation` object + sources = [mp.EigenModeSource(mp.GaussianSource(self.fcen,fwidth=self.df), + center=self.src_pt, + size=mp.Vector3(0,self.sy,0), + eig_band=mp.DiffractedPlanewave((0,0,0), + mp.Vector3(0,1,0), + 1, + 0))] + else: + sources = [mp.EigenModeSource(mp.GaussianSource(self.fcen,fwidth=self.df), + center=self.src_pt, + size=mp.Vector3(0,self.sy,0), + direction=mp.NO_DIRECTION, + eig_kpoint=k, + eig_band=1, + eig_parity=self.eig_parity)] + + sim = mp.Simulation(resolution=self.resolution, + cell_size=self.cell_size, + boundary_layers=self.absorber_layers, + k_point=k, + geometry=self.geometry, + sources=sources, + symmetries=symmetries) + + tran_flux = sim.add_flux(self.fcen, + 0, + 1, + mp.FluxRegion(center=self.tran_pt, + size=mp.Vector3(0,self.sy,0))) + + sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=self.tol)) + + tran = mp.get_fluxes(tran_flux)[0] + + sim.reset_meep() + + return tran + + + def test_mode_decomposition(self): + self.run_mode_decomposition(0,range(1,6),range(0,5)) + self.run_mode_decomposition(13.4,range(1,6),[-2,-1,-3,0,-4]) + + + def test_mode_source(self): + m = 3 + tran_bandnum = self.run_mode_source(m,False) + tran_diffpw = self.run_mode_source(m,True) + print("info:, {} (diffraction order), " + "{:.5f} [trans. (band number)]," + "{:.5f} [trans. (diffraction order)]".format(m, + tran_bandnum, + tran_diffpw)) + self.assertAlmostEqual(tran_bandnum,tran_diffpw,places=4) + if __name__ == '__main__': unittest.main() From f700429b5b085c691c8d0781cbeb3883f92aa5e0 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Fri, 13 May 2022 18:30:07 -0700 Subject: [PATCH 2/3] slightly increase tolerance for single-precision build --- doc/docs/Python_Tutorials/Eigenmode_Source.md | 14 +++++++----- python/examples/oblique-planewave.py | 22 ++++++++++++++----- python/tests/test_diffracted_planewave.py | 17 ++++++++------ 3 files changed, 35 insertions(+), 18 deletions(-) diff --git a/doc/docs/Python_Tutorials/Eigenmode_Source.md b/doc/docs/Python_Tutorials/Eigenmode_Source.md index a86cb139e..2824ec290 100644 --- a/doc/docs/Python_Tutorials/Eigenmode_Source.md +++ b/doc/docs/Python_Tutorials/Eigenmode_Source.md @@ -234,6 +234,8 @@ The simulation script is in [examples/oblique-planewave.py](https://github.com/N ```py import meep as mp import numpy as np +import matplotlib +matplotlib.use('agg') import matplotlib.pyplot as plt resolution = 50 # pixels/μm @@ -243,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) @@ -258,8 +260,7 @@ sources = [mp.EigenModeSource(src=mp.ContinuousSource(fsrc), 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, @@ -269,11 +270,12 @@ 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, + plot_boundaries_flag=False, output_plane=nonpml_vol) if mp.am_master(): @@ -289,7 +291,7 @@ Shown below are the steady-state field profiles generated by the planewave for t ![](../images/eigenmode_planewave.png) -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: +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 an $E_z$-polarized source (equivalent to the $\mathcal{S}$ polarization given the plane of incidence of $xy$), can also be defined using: ```py sources = [mp.EigenModeSource(src=mp.ContinuousSource(fsrc), diff --git a/python/examples/oblique-planewave.py b/python/examples/oblique-planewave.py index 4b3b971d4..eeedac345 100644 --- a/python/examples/oblique-planewave.py +++ b/python/examples/oblique-planewave.py @@ -1,7 +1,10 @@ import meep as mp import numpy as np +import matplotlib +matplotlib.use('agg') import matplotlib.pyplot as plt + resolution = 50 # pixels/μm cell_size = mp.Vector3(14,10,0) @@ -9,9 +12,9 @@ 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) @@ -24,8 +27,16 @@ 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)] + +## 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, @@ -35,11 +46,12 @@ 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, + plot_boundaries_flag=False, output_plane=nonpml_vol) if mp.am_master(): diff --git a/python/tests/test_diffracted_planewave.py b/python/tests/test_diffracted_planewave.py index a3da7d06b..a5aace85e 100644 --- a/python/tests/test_diffracted_planewave.py +++ b/python/tests/test_diffracted_planewave.py @@ -49,9 +49,9 @@ def setUp(cls): def run_mode_decomposition(self,theta,bands,orders): """Computes the transmittance of the diffraction orders of a binary grating given an incident planewave at - angle theta. The transmittance is computed using + angle `theta`. The transmittance is computed using mode decomposition in two different ways: - (1) band number (bands) and (2) diffraction order (orders) + (1) band number (`bands`) and (2) diffraction order (`orders`) via a `DiffractedPlanewave` object. The test verifies that these two methods produce equivalent results. """ @@ -154,10 +154,11 @@ def _pw_amp(x): def run_mode_source(self,m,diffpw): """Computes the transmitted Poynting flux of a binary grating given an incident planewave - specified by the diffraction order m in the + specified by the diffraction order `m` in the y direction. The incident planewave is defined using a mode source with either a band number - or `DiffractedPlanewave` object. + or `DiffractedPlanewave` object specified by + the boolean flag `diffpw`. """ ky = m/self.gp theta = math.asin(ky/(self.fcen*self.ng)) @@ -172,8 +173,8 @@ def run_mode_source(self,m,diffpw): symmetries = [mp.Mirror(direction=mp.Y)] if diffpw: - # the *zeroth* diffraction order defines a planewave with a wavevector - # equal to the `k_point` of the `Simulation` object + # the *zeroth* diffraction order defines a planewave with a + # wavevector equal to the `k_point` of the `Simulation` object sources = [mp.EigenModeSource(mp.GaussianSource(self.fcen,fwidth=self.df), center=self.src_pt, size=mp.Vector3(0,self.sy,0), @@ -227,7 +228,9 @@ def test_mode_source(self): "{:.5f} [trans. (diffraction order)]".format(m, tran_bandnum, tran_diffpw)) - self.assertAlmostEqual(tran_bandnum,tran_diffpw,places=4) + self.assertAlmostEqual(tran_bandnum, + tran_diffpw, + places=3 if mp.is_single_precision() else 4) if __name__ == '__main__': From a35f86774871732aaf1bb386767adcc2afea9f85 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Fri, 3 Jun 2022 08:37:58 -0700 Subject: [PATCH 3/3] fixes --- python/tests/test_diffracted_planewave.py | 253 +++++++++++----------- 1 file changed, 129 insertions(+), 124 deletions(-) diff --git a/python/tests/test_diffracted_planewave.py b/python/tests/test_diffracted_planewave.py index a5aace85e..a558baf65 100644 --- a/python/tests/test_diffracted_planewave.py +++ b/python/tests/test_diffracted_planewave.py @@ -9,150 +9,147 @@ class TestDiffractedPlanewave(unittest.TestCase): @classmethod def setUp(cls): - cls.resolution = 50 # pixels/μm + cls.resolution = 50 # pixels/μm - cls.dpml = 1.0 # PML thickness - cls.dsub = 3.0 # substrate thickness - cls.dpad = 3.0 # length of padding between grating and PML - cls.gp = 2.6 # grating period - cls.gh = 0.4 # grating height - cls.gdc = 0.5 # grating duty cycle + cls.dpml = 1.0 # PML thickness + cls.dsub = 3.0 # substrate thickness + cls.dpad = 3.0 # length of padding between grating and PML - cls.sx = cls.dpml+cls.dsub+cls.gh+cls.dpad+cls.dpml - cls.sy = cls.gp - - cls.cell_size = mp.Vector3(cls.sx,cls.sy,0) - cls.absorber_layers = [mp.PML(thickness=cls.dpml,direction=mp.X)] - - cls.wvl = 0.5 # center wavelength - cls.fcen = 1/cls.wvl # center frequency - cls.df = 0.05*cls.fcen # frequency width + cls.wvl = 0.5 # center wavelength + cls.fcen = 1/cls.wvl # center frequency cls.ng = 1.5 cls.glass = mp.Medium(index=cls.ng) - cls.eig_parity = mp.ODD_Z - - cls.src_pt = mp.Vector3(-0.5*cls.sx+cls.dpml,0,0) - cls.tran_pt = mp.Vector3(0.5*cls.sx-cls.dpml,0,0) - - cls.geometry = [mp.Block(material=cls.glass, - size=mp.Vector3(cls.dpml+cls.dsub,mp.inf,mp.inf), - center=mp.Vector3(-0.5*cls.sx+0.5*(cls.dpml+cls.dsub),0,0)), - mp.Block(material=cls.glass, - size=mp.Vector3(cls.gh,cls.gdc*cls.gp,mp.inf), - center=mp.Vector3(-0.5*cls.sx+cls.dpml+cls.dsub+0.5*cls.gh,0,0))] - - cls.tol = 1e-8 # tolerance for mp.stop_dft_decayed termination criteria + cls.pml_layers = [mp.PML(thickness=cls.dpml,direction=mp.X)] - def run_mode_decomposition(self,theta,bands,orders): - """Computes the transmittance of the diffraction orders - of a binary grating given an incident planewave at - angle `theta`. The transmittance is computed using - mode decomposition in two different ways: - (1) band number (`bands`) and (2) diffraction order (`orders`) - via a `DiffractedPlanewave` object. The test verifies - that these two methods produce equivalent results. + def run_binary_grating_diffraction(self,gp,gh,gdc,theta): + """Computes the mode coefficient of the transmitted orders of + a binary grating given an incident planewave using the two + approaches of band number and `DiffractedPlanewave` object in + `get_eigenmode_coefficients` and verifies that they produce + the same result. """ + 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) + k = mp.Vector3(self.fcen*self.ng).rotate(mp.Vector3(z=1), theta_in) - symmetries = [] - if theta_in == 0: + eig_parity = mp.ODD_Z + if theta == 0: k = mp.Vector3() - self.eig_parity += mp.EVEN_Y + eig_parity += mp.EVEN_Y symmetries = [mp.Mirror(direction=mp.Y)] + else: + symmetries = [] def pw_amp(k,x0): def _pw_amp(x): return cmath.exp(1j*2*math.pi*k.dot(x+x0)) return _pw_amp - sources = [mp.Source(mp.GaussianSource(self.fcen,fwidth=self.df), + src_pt = mp.Vector3(-0.5*sx+self.dpml,0,0) + sources = [mp.Source(mp.GaussianSource(self.fcen,fwidth=0.1*self.fcen), component=mp.Ez, - center=self.src_pt, - size=mp.Vector3(0,self.sy,0), - amp_func=pw_amp(k,self.src_pt))] - - sim = mp.Simulation(resolution=self.resolution, - cell_size=self.cell_size, - boundary_layers=self.absorber_layers, - k_point=k, - default_material=self.glass, - sources=sources, - symmetries=symmetries) - - tran_flux = sim.add_flux(self.fcen, - 0, - 1, - mp.FluxRegion(center=self.tran_pt, - size=mp.Vector3(0,self.sy,0))) + center=src_pt, + size=mp.Vector3(0,sy,0), + amp_func=pw_amp(k,src_pt))] - sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=self.tol)) - - input_flux = mp.get_fluxes(tran_flux) - - sim.reset_meep() + 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=self.cell_size, - boundary_layers=self.absorber_layers, - geometry=self.geometry, + cell_size=cell_size, + boundary_layers=self.pml_layers, + geometry=geometry, k_point=k, sources=sources, symmetries=symmetries) + tran_pt = mp.Vector3(0.5*sx-self.dpml,0,0) tran_flux = sim.add_mode_monitor(self.fcen, 0, 1, - mp.FluxRegion(center=self.tran_pt, - size=mp.Vector3(0,self.sy,0))) + 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)) - sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=self.tol)) + m_plus = int(np.floor((self.fcen-k.y)*gp)) + m_minus = int(np.ceil((-self.fcen-k.y)*gp)) - for band,order in zip(bands,orders): + if theta == 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] + ids = np.flip(np.argsort(kxs)) + orders = [ms[d] for d in ids] + + for band,order in enumerate(orders): res = sim.get_eigenmode_coefficients(tran_flux, - [band], - eig_parity=self.eig_parity) - tran_ref = abs(res.alpha[0,0,0])**2/input_flux[0] - if (theta_in == 0): + [band+1], + eig_parity=eig_parity) + tran_ref = abs(res.alpha[0,0,0])**2 + if theta == 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)) - if res is not None: - tran_dp = abs(res.alpha[0,0,0])**2/input_flux[0] - if ((theta_in == 0) and (order == 0)): - tran_dp = 0.5*tran_dp - else: - tran_dp = 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] err = abs(tran_ref-tran_dp)/tran_ref - print("info:, {:2d} (band number), {:2d} (diffraction order), " - "{:10.8f} [trans. (band number)], " - "{:10.8f} [trans. (diffracted order)], " - "{:10.8f} (error)".format(band, - order, - tran_ref, - tran_dp, - err)) + print("tran:, {:2d} (band), {:2d} (order), " + "{:10.8f} (band num.), {:10.8f} (diff. pw.), " + "{:10.8f} (error)".format(band,order,tran_ref,tran_dp,err)) + + self.assertAlmostEqual(vg_ref,vg_dp,places=4) + self.assertAlmostEqual(tran_ref,tran_dp,places=4) + if theta == 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) + else: + self.assertAlmostEqual(kdom_ref.x,kdom_dp.x,places=5) + self.assertAlmostEqual(kdom_ref.y,kdom_dp.y,places=5) + self.assertAlmostEqual(kdom_ref.z,kdom_dp.z,places=5) - self.assertAlmostEqual(vg_ref,vg_dp,places=5) - self.assertAlmostEqual(tran_ref,tran_dp,places=5) + print("PASSED.") - def run_mode_source(self,m,diffpw): - """Computes the transmitted Poynting flux of a + def test_diffracted_planewave(self): + self.run_binary_grating_diffraction(2.6,0.4,0.6,0) + self.run_binary_grating_diffraction(2.6,0.4,0.6,13.4) + + # 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,diffpw): + """Computes the transmitted flux of a binary grating given an incident planewave specified by the diffraction order `m` in the y direction. The incident planewave is defined @@ -160,52 +157,67 @@ def run_mode_source(self,m,diffpw): or `DiffractedPlanewave` object specified by the boolean flag `diffpw`. """ - ky = m/self.gp + sx = self.dpml+self.dsub+gh+self.dpad+self.dpml + sy = gp + cell_size = mp.Vector3(sx,sy,0) + + ky = m/gp theta = 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) - symmetries = [] + eig_parity = mp.ODD_Z if theta == 0: k = mp.Vector3() - self.eig_parity += mp.EVEN_Y + eig_parity += mp.EVEN_Y symmetries = [mp.Mirror(direction=mp.Y)] + else: + symmetries = [] + src_pt = mp.Vector3(-0.5*sx+self.dpml,0,0) if diffpw: # the *zeroth* diffraction order defines a planewave with a # wavevector equal to the `k_point` of the `Simulation` object - sources = [mp.EigenModeSource(mp.GaussianSource(self.fcen,fwidth=self.df), - center=self.src_pt, - size=mp.Vector3(0,self.sy,0), + sources = [mp.EigenModeSource(mp.GaussianSource(self.fcen,fwidth=0.1*self.fcen), + center=src_pt, + size=mp.Vector3(0,sy,0), eig_band=mp.DiffractedPlanewave((0,0,0), mp.Vector3(0,1,0), 1, 0))] else: - sources = [mp.EigenModeSource(mp.GaussianSource(self.fcen,fwidth=self.df), - center=self.src_pt, - size=mp.Vector3(0,self.sy,0), + sources = [mp.EigenModeSource(mp.GaussianSource(self.fcen,fwidth=0.1*self.fcen), + center=src_pt, + size=mp.Vector3(0,sy,0), direction=mp.NO_DIRECTION, eig_kpoint=k, eig_band=1, - eig_parity=self.eig_parity)] + eig_parity=eig_parity)] + + 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=self.cell_size, - boundary_layers=self.absorber_layers, + cell_size=cell_size, + boundary_layers=self.pml_layers, k_point=k, - geometry=self.geometry, + 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=self.tran_pt, - size=mp.Vector3(0,self.sy,0))) + mp.FluxRegion(center=tran_pt, + size=mp.Vector3(0,sy,0))) - sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=self.tol)) + sim.run(until_after_sources=mp.stop_when_fields_decayed(20,mp.Ez,src_pt,1e-6)) tran = mp.get_fluxes(tran_flux)[0] @@ -214,22 +226,15 @@ def run_mode_source(self,m,diffpw): return tran - def test_mode_decomposition(self): - self.run_mode_decomposition(0,range(1,6),range(0,5)) - self.run_mode_decomposition(13.4,range(1,6),[-2,-1,-3,0,-4]) - - def test_mode_source(self): - m = 3 - tran_bandnum = self.run_mode_source(m,False) - tran_diffpw = self.run_mode_source(m,True) - print("info:, {} (diffraction order), " - "{:.5f} [trans. (band number)]," - "{:.5f} [trans. (diffraction order)]".format(m, - tran_bandnum, - tran_diffpw)) - self.assertAlmostEqual(tran_bandnum, - tran_diffpw, + tran_bn = self.run_mode_source(1.5,0.5,0.3,2,False) + tran_dp = self.run_mode_source(1.5,0.5,0.3,2,True) + print("mode-source:, " + "{:.5f} (band number), " + "{:.5f} (diffraction order)".format(tran_bn, + tran_dp)) + self.assertAlmostEqual(tran_bn, + tran_dp, places=3 if mp.is_single_precision() else 4)