unit test for mode decomposition with out-of-plane wavevector

test_binary_grating.py

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 = 50        # pixels/μm                                                                                                                   

    cls.dpml = 1.0             # PML thickness                                                                                                               
    cls.dsub = 1.0             # substrate thickness                                                                                                         
    cls.dpad = 1.0             # length of padding between grating and PML                                                                                   
    cls.gp = 0.5               # grating period                                                                                                              
    cls.gh = 0.5               # 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)

    # 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)]

    wvl = 0.5              # center wavelength                                                                                                               
    cls.fcen = 1/wvl       # center frequency                                                                                                                
    cls.df = 0.05*cls.fcen # frequency width                                                                                                                 

    cls.ng = 1.5
    cls.glass = mp.Medium(index=cls.ng)
    
     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))]

  @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)
    print("k_point:, ({}, {}, {})".format(k.x,k.y,k.z))

    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="complex")

    refl_pt = mp.Vector3(-0.5*self.sx+self.dpml+0.5*self.dsub,0,0)
    refl_flux = sim.add_mode_monitor(self.fcen,
                                     0,
                                     1,
                                     mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,self.sy,0)))

    sim.run(until_after_sources=10)

    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="complex")
                        
    refl_flux = sim.add_mode_monitor(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_mode_monitor(self.fcen,
                                     0,
                                     1,
                                     mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,self.sy,0)))

    sim.run(until_after_sources=500)

    # number of reflected orders                                                                                                                             
    nm_r = np.ceil((np.sqrt((self.fcen*self.ng)**2-k.z**2)-k.y)*self.gp)-np.floor((-np.sqrt((self.fcen*self.ng)**2-k.z**2)-k.y)*self.gp)
    nm_r = nm_r/2 # since eig_parity removes degeneracy in y-direction                                                                                       
    nm_r = int(nm_r)
    print("number of reflected orders: {}".format(nm_r))

    Rsum = 0
    for nm in range(nm_r):
      res = sim.get_eigenmode_coefficients(refl_flux,
                                           mp.DiffractedPlanewave([0,nm,0],
                                                                  mp.Vector3(1,0,0),
                                                                  0,
                                                                  1))
      r_coeffs = res.alpha
      Rmode = abs(r_coeffs[0,0,1])**2/input_flux[0]
      print("refl-order:, {}, {}".format(nm,Rmode))
      Rsum += Rmode if nm == 0 else 2*Rmode

    # number of transmitted orders                                                                                                                           
    nm_t = np.ceil((np.sqrt(self.fcen**2-k.z**2)-k.y)*self.gp)-np.floor((-np.sqrt(self.fcen**2-k.z**2)-k.y)*self.gp)
    nm_t = nm_t/2 # since eig_parity removes degeneracy in y-direction                                                                                       
    nm_t = int(nm_t)
    print("number of transmitted orders: {}".format(nm_t))

    Tsum = 0
    for nm in range(nm_t):
      res = sim.get_eigenmode_coefficients(tran_flux,
                                           mp.DiffractedPlanewave([0,nm,0],
                                                                  mp.Vector3(1,0,0),
                                                                  0,
                                                                  1))
      t_coeffs = res.alpha
      Tmode = abs(t_coeffs[0,0,0])**2/input_flux[0]
      print("tran-order:, {}, {}".format(nm,Tmode))
      Tsum += Tmode if nm == 0 else 2*Tmode

    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]

    print("refl:, {}, {}".format(Rsum,Rflux))
    print("tran:, {}, {}".format(Tsum,Tflux))
    print("sum:,  {}, {}".format(Rsum+Tsum,Rflux+Tflux))
    self.assertAlmostEqual(Rsum,Rflux,places=2)
    self.assertAlmostEqual(Tsum,Tflux,places=2)


if __name__ == '__main__':
  unittest.main()