diff --git a/doc/docs/Python_Tutorials/Eigenmode_Source.md b/doc/docs/Python_Tutorials/Eigenmode_Source.md index c54140087..ea74dde29 100644 --- a/doc/docs/Python_Tutorials/Eigenmode_Source.md +++ b/doc/docs/Python_Tutorials/Eigenmode_Source.md @@ -243,9 +243,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 +258,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,16 +268,14 @@ 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, output_plane=nonpml_vol) - -if mp.am_master(): - plt.axis('off') - plt.savefig('pw.png',bbox_inches='tight') +plt.axis('off') +plt.show() ``` Note that the line source spans the *entire* length of the cell. (Planewave sources extending into the PML region must include `is_integrated=True`.) This example involves a continuous-wave (CW) time profile. For a pulse profile, the oblique planewave is incident at a given angle for only a *single* frequency component of the source as described in [Tutorial/Basics/Angular Reflectance Spectrum of a Planar Interface](../Python_Tutorials/Basics.md#angular-reflectance-spectrum-of-a-planar-interface). This is a fundamental feature of FDTD simulations and not of Meep per se. To simulate an incident planewave at multiple angles for a given frequency $\omega$, you will need to do separate simulations involving different values of $\vec{k}$ (`k_point`) since each set of $(\vec{k},\omega)$ specifying the Bloch-periodic boundaries and the frequency of the source will produce a different angle of the planewave. For more details, refer to Section 4.5 ("Efficient Frequency-Angle Coverage") in [Chapter 4](https://arxiv.org/abs/1301.5366) ("Electromagnetic Wave Source Conditions") of [Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707). @@ -288,3 +285,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..cdf8685e3 100644 --- a/python/examples/oblique-planewave.py +++ b/python/examples/oblique-planewave.py @@ -2,6 +2,7 @@ import numpy as np import matplotlib.pyplot as plt + resolution = 50 # pixels/μm cell_size = mp.Vector3(14,10,0) @@ -9,9 +10,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 +25,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,13 +44,11 @@ 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, output_plane=nonpml_vol) - -if mp.am_master(): - plt.axis('off') - plt.savefig('pw.png',bbox_inches='tight') +plt.axis('off') +plt.show() diff --git a/python/tests/test_diffracted_planewave.py b/python/tests/test_diffracted_planewave.py index 7c0d8b3c0..cdbbca1fe 100644 --- a/python/tests/test_diffracted_planewave.py +++ b/python/tests/test_diffracted_planewave.py @@ -5,16 +5,12 @@ import numpy as np -# Computes the mode coefficient of the transmitted orders of -# a binary grating given an incident planewave and verifies -# that the results are the same when using either a band number -# or `DiffractedPlanewave` object in `get_eigenmode_coefficients`. - class TestDiffractedPlanewave(unittest.TestCase): @classmethod def setUp(cls): - cls.resolution = 50 # pixels/um + 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 @@ -29,6 +25,12 @@ def setUp(cls): 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) @@ -88,7 +90,7 @@ def _pw_amp(x): m_minus = int(np.ceil((-self.fcen-k.y)*gp)) if theta == 0: - orders = range(m_plus) + 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) @@ -137,6 +139,7 @@ def _pw_amp(x): print("PASSED.") + 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) @@ -144,5 +147,94 @@ def test_diffracted_planewave(self): # 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] + + 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) + + if __name__ == '__main__': unittest.main()