add unit test and apply phase fix only for real fields and beta

python/tests/test_binary_grating.py

@@ -1,42 +1,55 @@
 import unittest
+import parameterized
 import meep as mp
 import math
 import cmath
 import numpy as np
 
+
 class TestEigCoeffs(unittest.TestCase):
+  @classmethod
+  def setUpClass(cls):
+    cls.resolution = 30        # pixels/μm
 
-  def run_binary_grating_oblique(self, theta):
+    cls.dpml = 1.0             # PML thickness
+    cls.dsub = 1.0             # substrate thickness
+    cls.dpad = 1.0             # padding thickness between grating and PML
+    cls.gp = 6.0               # grating period
+    cls.gh = 0.5               # grating height
+    cls.gdc = 0.5              # grating duty cycle
 
-    resolution = 30        # pixels/um
+    cls.sx = cls.dpml+cls.dsub+cls.gh+cls.dpad+cls.dpml
+    cls.sy = cls.gp
 
-    dpml = 1.0             # PML thickness
-    dsub = 1.0             # substrate thickness
-    dpad = 1.0             # length of padding between grating and pml
-    gp = 6.0               # grating period
-    gh = 0.5               # grating height
-    gdc = 0.5              # grating duty cycle
+    cls.cell_size = mp.Vector3(cls.sx,cls.sy,0)
 
-    sx = dpml+dsub+gh+dpad+dpml
-    sy = gp
+    # replace anisotropic PML with isotropic Absorber to
+    # attenuate parallel-directed fields of oblique source
+    cls.abs_layers = [mp.Absorber(thickness=cls.dpml,direction=mp.X)]
 
-    cell_size = mp.Vector3(sx,sy,0)
+    wvl = 0.5              # center wavelength
+    cls.fcen = 1/wvl       # center frequency
+    cls.df = 0.05*cls.fcen # frequency width
 
-    # replace anisotropic PML with isotropic Absorber to attenuate parallel-directed fields of oblique source
-    abs_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.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))]
 
-    ng = 1.5
-    glass = mp.Medium(index=ng)
 
-    # rotation angle of incident planewave; CCW about Y axis, 0 degrees along +X axis
+  @parameterized.parameterized.expand([(0,), (10.7,)])
+  def test_binary_grating_oblique(self, theta):
+    # rotation angle of incident planewave
+    # counterclockwise (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(math.cos(theta_in),math.sin(theta_in),0).scale(fcen*ng)
+    k = mp.Vector3(self.fcen*self.ng).rotate(mp.Vector3(0,0,1), theta_in)
 
     symmetries = []
     eig_parity = mp.ODD_Z
@@ -50,68 +63,186 @@ 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.3*dsub,0,0)
-    sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
-                         component=mp.Ez,
+    src_pt = mp.Vector3(-0.5*self.sx+self.dpml,0,0)
+    sources = [mp.Source(mp.GaussianSource(self.fcen,fwidth=self.df),
+                         component=mp.Ez, # S polarization
                          center=src_pt,
-                         size=mp.Vector3(0,sy,0),
+                         size=mp.Vector3(0,self.sy,0),
                          amp_func=pw_amp(k,src_pt))]
 
-    sim = mp.Simulation(resolution=resolution,
-                        cell_size=cell_size,
-                        boundary_layers=abs_layers,
+    sim = mp.Simulation(resolution=self.resolution,
+                        cell_size=self.cell_size,
+                        boundary_layers=self.abs_layers,
                         k_point=k,
-                        default_material=glass,
+                        default_material=self.glass,
                         sources=sources,
                         symmetries=symmetries)
 
-    refl_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub,0,0)
-    refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))
+    refl_pt = mp.Vector3(-0.5*self.sx+self.dpml+0.5*self.dsub,0,0)
+    refl_flux = sim.add_flux(self.fcen,
+                             0,
+                             1,
+                             mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,self.sy,0)))
 
-    sim.run(until_after_sources=100)
+    sim.run(until_after_sources=mp.stop_when_dft_decayed())
 
     input_flux = mp.get_fluxes(refl_flux)
     input_flux_data = sim.get_flux_data(refl_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=abs_layers,
-                        geometry=geometry,
+    sim = mp.Simulation(resolution=self.resolution,
+                        cell_size=self.cell_size,
+                        boundary_layers=self.abs_layers,
+                        geometry=self.geometry,
                         k_point=k,
                         sources=sources,
                         symmetries=symmetries)
 
-    refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))
+    refl_flux = sim.add_flux(self.fcen,
+                             0,
+                             1,
+                             mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,self.sy,0)))
+
     sim.load_minus_flux_data(refl_flux,input_flux_data)
 
-    tran_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad,0,0)
-    tran_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,sy,0)))
+    tran_pt = mp.Vector3(0.5*self.sx-self.dpml-0.5*self.dpad,0,0)
+    tran_flux = sim.add_flux(self.fcen,
+                             0,
+                             1,
+                             mp.FluxRegion(center=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())
 
-    nm_r = np.floor((fcen*ng-k.y)*gp)-np.ceil((-fcen*ng-k.y)*gp) # number of reflected orders
+    # number of reflected orders
+    nm_r = np.floor((self.fcen*self.ng-k.y)*self.gp)-np.ceil((-self.fcen*self.ng-k.y)*self.gp)
     if theta_in == 0:
       nm_r = nm_r/2 # since eig_parity removes degeneracy in y-direction
     nm_r = int(nm_r)
 
-    res = sim.get_eigenmode_coefficients(refl_flux, range(1,nm_r+1), eig_parity=eig_parity)
+    res = sim.get_eigenmode_coefficients(refl_flux,
+                                         range(1,nm_r+1),
+                                         eig_parity=eig_parity)
     r_coeffs = res.alpha
 
     Rsum = 0
     for nm in range(nm_r):
       Rsum += abs(r_coeffs[nm,0,1])**2/input_flux[0]
 
-    nm_t = np.floor((fcen-k.y)*gp)-np.ceil((-fcen-k.y)*gp)       # number of transmitted orders
+    # number of transmitted orders
+    nm_t = np.floor((self.fcen-k.y)*self.gp)-np.ceil((-self.fcen-k.y)*self.gp)
     if theta_in == 0:
       nm_t = nm_t/2 # since eig_parity removes degeneracy in y-direction
     nm_t = int(nm_t)
 
-    res = sim.get_eigenmode_coefficients(tran_flux, range(1,nm_t+1), eig_parity=eig_parity)
+    res = sim.get_eigenmode_coefficients(tran_flux,
+                                         range(1,nm_t+1),
+                                         eig_parity=eig_parity)
+    t_coeffs = res.alpha
+
+    Tsum = 0
+    for nm in range(nm_t):
+      Tsum += abs(t_coeffs[nm,0,0])**2/input_flux[0]
+
+    r_flux = mp.get_fluxes(refl_flux)
+    t_flux = mp.get_fluxes(tran_flux)
+    Rflux = -r_flux[0]/input_flux[0]
+    Tflux =  t_flux[0]/input_flux[0]
+
+    self.assertAlmostEqual(Rsum,Rflux,places=2)
+    self.assertAlmostEqual(Tsum,Tflux,places=2)
+
+
+  @parameterized.parameterized.expand([(13.2,)])
+  def test_binary_grating_special_kz(self, theta):
+    # rotation angle of incident planewave
+    # counterclockwise (CCW) about Y axis, 0 degrees along +X axis
+    theta_in = math.radians(theta)
+
+    # k (in source medium) with correct length (plane of incidence: XZ)
+    k = mp.Vector3(self.fcen*self.ng).rotate(mp.Vector3(0,1,0), theta_in)
+
+    symmetries = [mp.Mirror(mp.Y)]
+    eig_parity = mp.EVEN_Y
+
+    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*self.sx+self.dpml,0,0)
+    sources = [mp.Source(mp.GaussianSource(self.fcen,fwidth=self.df),
+                         component=mp.Ez, # P polarization
+                         center=src_pt,
+                         size=mp.Vector3(0,self.sy,0),
+                         amp_func=pw_amp(k,src_pt))]
+
+    sim = mp.Simulation(resolution=self.resolution,
+                        cell_size=self.cell_size,
+                        boundary_layers=self.abs_layers,
+                        k_point=k,
+                        default_material=self.glass,
+                        sources=sources,
+                        symmetries=symmetries,
+                        kz_2d="3d")
+
+    refl_pt = mp.Vector3(-0.5*self.sx+self.dpml+0.5*self.dsub,0,0)
+    refl_flux = sim.add_flux(self.fcen,
+                             0,
+                             1,
+                             mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,self.sy,0)))
+
+    sim.run(until_after_sources=mp.stop_when_dft_decayed())
+
+    input_flux = mp.get_fluxes(refl_flux)
+    input_flux_data = sim.get_flux_data(refl_flux)
+
+    sim.reset_meep()
+
+    sim = mp.Simulation(resolution=self.resolution,
+                        cell_size=self.cell_size,
+                        boundary_layers=self.abs_layers,
+                        geometry=self.geometry,
+                        k_point=k,
+                        sources=sources,
+                        symmetries=symmetries,
+                        kz_2d="3d")
+
+    refl_flux = sim.add_flux(self.fcen,
+                             0,
+                             1,
+                             mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,self.sy,0)))
+
+    sim.load_minus_flux_data(refl_flux,input_flux_data)
+
+    tran_pt = mp.Vector3(0.5*self.sx-self.dpml-0.5*self.dpad,0,0)
+    tran_flux = sim.add_flux(self.fcen,
+                             0,
+                             1,
+                             mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,self.sy,0)))
+
+    sim.run(until_after_sources=mp.stop_when_dft_decayed())
+
+    # number of reflected orders
+    nm_r = np.floor((np.sqrt((self.fcen*self.ng)**2-k.z**2)-k.y)*self.gp)-np.ceil((-np.sqrt((self.fcen*self.ng)**2-k.z**2)-k.y)*self.gp)
+    nm_r = int(nm_r/2) # eig_parity removes degeneracy in y-direction
+
+    res = sim.get_eigenmode_coefficients(refl_flux,
+                                         range(1,nm_r+1),
+                                         eig_parity=eig_parity)
+    r_coeffs = res.alpha
+
+    Rsum = 0
+    for nm in range(nm_r):
+      Rsum += abs(r_coeffs[nm,0,1])**2/input_flux[0]
+
+    # number of transmitted orders
+    nm_t = np.floor((np.sqrt(self.fcen**2-k.z**2)-k.y)*self.gp)-np.ceil((-np.sqrt(self.fcen**2-k.z**2)-k.y)*self.gp)
+    nm_t = int(nm_t/2) # eig_parity removes degeneracy in y-direction
+
+    res = sim.get_eigenmode_coefficients(tran_flux,
+                                         range(1,nm_t+1),
+                                         eig_parity=eig_parity)
     t_coeffs = res.alpha
 
     Tsum = 0
@@ -123,12 +254,11 @@ def _pw_amp(x):
     Rflux = -r_flux[0]/input_flux[0]
     Tflux =  t_flux[0]/input_flux[0]
 
+    print("refl:, {}, {}".format(Rsum,Rflux))
+    print("tran:, {}, {}".format(Tsum,Tflux))
     self.assertAlmostEqual(Rsum,Rflux,places=2)
     self.assertAlmostEqual(Tsum,Tflux,places=2)
 
-  def test_binary_grating(self):
-    self.run_binary_grating_oblique(0)
-    self.run_binary_grating_oblique(10.7)
 
 if __name__ == '__main__':
   unittest.main()

python/tests/test_special_kz.py

@@ -12,17 +12,18 @@ def refl_planar(self, theta, kz_2d):
         resolution = 100  # pixels/um
 
         dpml = 1.0
-        sx = 3+2*dpml
+        sx = 3.0 + 2*dpml
         sy = 1/resolution
         cell_size = mp.Vector3(sx,sy)
         pml_layers = [mp.PML(dpml,direction=mp.X)]
 
         fcen = 1.0
 
-        k_point = mp.Vector3(z=math.sin(theta)).scale(fcen)
+        # plane of incidence is XZ
+        k_point = mp.Vector3(1,0,0).rotate(mp.Vector3(0,1,0),theta).scale(fcen)
 
         sources = [mp.Source(mp.GaussianSource(fcen,fwidth=0.2*fcen),
-                             component=mp.Ez,
+                             component=mp.Ez, # P-polarization
                              center=mp.Vector3(-0.5*sx+dpml),
                              size=mp.Vector3(y=sy))]
 
@@ -79,26 +80,28 @@ def test_special_kz(self):
         Rfresnel = lambda theta_in: math.fabs((n1*math.cos(theta_out(theta_in))-n2*math.cos(theta_in))
                                               / (n1*math.cos(theta_out(theta_in))+n2*math.cos(theta_in)))**2
 
+        # angle of incident planewave; clockwise (CW) about Y axis, 0 degrees along +X axis
         theta = math.radians(23)
 
         start = time()
-        Rmeep_complex = self.refl_planar(theta, 'complex')
+        Rmeep_complex = self.refl_planar(theta, "complex")
         t_complex = time() - start
 
         start = time()
-        Rmeep_real_imag = self.refl_planar(theta, 'real/imag')
+        Rmeep_real_imag = self.refl_planar(theta, "real/imag")
         t_real_imag = time() - start
 
         Rfres = Rfresnel(theta)
 
         self.assertAlmostEqual(Rmeep_complex,Rfres,places=2)
         self.assertAlmostEqual(Rmeep_real_imag,Rfres,places=2)
 
-        # the real/imag algorithm should be faster, but on CI machines performance is too variable for this to reliably hold
+        # the real/imag algorithm should be faster, but on CI machines performance is too variable
+        # for this to reliably hold
         # self.assertLess(t_real_imag,t_complex) 
 
 
-    @parameterized.parameterized.expand(["complex","real/imag"])
+    @parameterized.parameterized.expand([("complex",),("real/imag",)])
     def test_eigsrc_kz(self, kz_2d):
         resolution = 30 # pixels/um
 

src/mpb.cpp

@@ -292,11 +292,11 @@ static int nextpow2357(int n) {
   }
 }
 
-void special_kz_phasefix(eigenmode_data *edata, bool neg) {
+void special_kz_phasefix(eigenmode_data *edata, bool phase_flip) {
   size_t n = edata->n[0] * edata->n[1] * edata->n[2];
   complex<mpb_real> *E = (complex<mpb_real> *)edata->fft_data_E;
   complex<mpb_real> *H = (complex<mpb_real> *)edata->fft_data_H;
-  complex<mpb_real> im(0, neg ? -1 : 1);
+  complex<mpb_real> im(0, phase_flip ? -1 : 1);
   for (size_t i = 0; i < n; ++i) {
     E[3*i + 2] *= im; // Ez
     H[3*i + 0] *= im; // Hx
@@ -821,7 +821,7 @@ void fields::add_eigenmode_source(component c0, const src_time &src, direction d
                                       NULL, NULL, dp);
   finished_working();
 
-  if (beta != 0) special_kz_phasefix(global_eigenmode_data, true);
+  if (is_real && beta != 0) special_kz_phasefix(global_eigenmode_data, true /* phase_flip */);
 
   if (global_eigenmode_data == NULL) meep::abort("MPB could not find the eigenmode");
 
@@ -937,7 +937,7 @@ void fields::get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, in
         continue;
       }
 
-      if (beta != 0) special_kz_phasefix((eigenmode_data *) mode_data, false);
+      if (is_real && beta != 0) special_kz_phasefix((eigenmode_data *) mode_data, false /* phase_flip */);
 
       double vg = get_group_velocity(mode_data);
       vec kfound = get_k(mode_data);