diff --git a/doc/docs/Python_Tutorials/Eigenmode_Source.md b/doc/docs/Python_Tutorials/Eigenmode_Source.md index c54140087..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(): @@ -288,3 +290,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 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), + 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/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 962a3e3ea..a558baf65 100644 --- a/python/tests/test_diffracted_planewave.py +++ b/python/tests/test_diffracted_planewave.py @@ -7,131 +7,236 @@ class TestDiffractedPlanewave(unittest.TestCase): - def run_binary_grating_diffraction(self, gp, gh, gdc, theta, bands, orders): - resolution = 50 # pixels/um + @classmethod + def setUp(cls): + cls.resolution = 50 # pixels/μm - dpml = 1.0 # PML thickness - dsub = 3.0 # substrate thickness - dpad = 3.0 # length of padding between grating and PML + cls.dpml = 1.0 # PML thickness + cls.dsub = 3.0 # substrate thickness + cls.dpad = 3.0 # length of padding between grating and PML - sx = dpml+dsub+gh+dpad+dpml - sy = gp + cls.wvl = 0.5 # center wavelength + cls.fcen = 1/cls.wvl # center frequency - cell_size = mp.Vector3(sx,sy,0) - absorber_layers = [mp.Absorber(thickness=dpml,direction=mp.X)] + cls.ng = 1.5 + cls.glass = mp.Medium(index=cls.ng) - wvl = 0.5 # center wavelength - fcen = 1/wvl # center frequency - df = 0.05*fcen # frequency width + cls.pml_layers = [mp.PML(thickness=cls.dpml,direction=mp.X)] - ng = 1.5 - glass = mp.Medium(index=ng) + + 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) - 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: + eig_parity = mp.ODD_Z + if theta == 0: k = mp.Vector3() 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 - src_pt = mp.Vector3(-0.5*sx+dpml,0,0) - sources = [mp.Source(mp.GaussianSource(fcen,fwidth=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=src_pt, size=mp.Vector3(0,sy,0), amp_func=pw_amp(k,src_pt))] - sim = mp.Simulation(resolution=resolution, - cell_size=cell_size, - boundary_layers=absorber_layers, - k_point=k, - default_material=glass, - sources=sources, - symmetries=symmetries) - - tran_pt = mp.Vector3(0.5*sx-dpml,0,0) - tran_flux = sim.add_flux(fcen, - 0, - 1, - mp.FluxRegion(center=tran_pt, - size=mp.Vector3(0,sy,0))) - - sim.run(until_after_sources=10) - - 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, + 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+dpml+dsub+0.5*gh,0,0))] + center=mp.Vector3(-0.5*sx+self.dpml+self.dsub+0.5*gh,0,0))] - sim = mp.Simulation(resolution=resolution, + sim = mp.Simulation(resolution=self.resolution, cell_size=cell_size, - boundary_layers=absorber_layers, + boundary_layers=self.pml_layers, geometry=geometry, k_point=k, sources=sources, symmetries=symmetries) - tran_flux = sim.add_mode_monitor(fcen, + tran_pt = mp.Vector3(0.5*sx-self.dpml,0,0) + tran_flux = sim.add_mode_monitor(self.fcen, 0, 1, mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,sy,0))) - sim.run(until_after_sources=100) + sim.run(until_after_sources=mp.stop_when_fields_decayed(20,mp.Ez,src_pt,1e-6)) + + 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], + [band+1], eig_parity=eig_parity) - tran_ref = abs(res.alpha[0,0,0])**2/input_flux[0] - if (theta_in == 0): + 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("tran:, {:2d} (band), {:2d} (order), {:10.8f} (eigensolver), {:10.8f} (planewave), {: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) + + print("PASSED.") - 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(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 + using a mode source with either a band number + or `DiffractedPlanewave` object specified by + the boolean flag `diffpw`. + """ + 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) + + eig_parity = mp.ODD_Z + if theta == 0: + k = mp.Vector3() + 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=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=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=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=cell_size, + boundary_layers=self.pml_layers, + k_point=k, + 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=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)) + + tran = mp.get_fluxes(tran_flux)[0] + + sim.reset_meep() + + return tran + + + def test_mode_source(self): + 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) - # 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]) if __name__ == '__main__': unittest.main()