Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

add support for out-of-plane wavevector in 2d cell to fields::get_eigenmode #1968

Merged
merged 3 commits into from
Mar 17, 2022
Merged

add support for out-of-plane wavevector in 2d cell to fields::get_eigenmode #1968

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
4f79698
add support for out-of-plane wavevector in 2d cell to fields::get_eig…
oskooi Feb 27, 2022
56dd42e
add unit test and apply phase fix only for real fields and beta
oskooi Mar 1, 2022
826a21a
sum S- and P-polarized diffraction orders when kz is non-zero
oskooi Mar 10, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
256 changes: 206 additions & 50 deletions python/tests/test_binary_grating.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,47 +1,61 @@
# -*- coding: utf-8 -*-

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
if theta_in == 0:
k = mp.Vector3(0,0,0)
symmetries = [mp.Mirror(mp.Y)]
eig_parity += mp.EVEN_Y

Expand All @@ -50,68 +64,84 @@ 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
Expand All @@ -123,12 +153,138 @@ 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))
print("sum:, {}, {}".format(Rsum+Tsum,Rflux+Tflux))

self.assertAlmostEqual(Rsum,Rflux,places=2)
self.assertAlmostEqual(Tsum,Tflux,places=2)
self.assertAlmostEqual(Rsum+Tsum,1.00,places=2)


@parameterized.parameterized.expand([(13.2,"real/imag"),
(17.7,"complex"),
(21.2, "3d")])
def test_binary_grating_special_kz(self, theta, kz_2d):
# 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)]

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,
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=kz_2d)

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=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=kz_2d)

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=mp.stop_when_dft_decayed())

# 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 = int(nm_r/2)

Rsum = 0
for nm in range(nm_r):
for S_pol in [False, True]:
res = sim.get_eigenmode_coefficients(refl_flux,
mp.DiffractedPlanewave([0,nm,0],
mp.Vector3(1,0,0),
1 if S_pol else 0,
0 if S_pol else 1))
r_coeffs = res.alpha
Rmode = abs(r_coeffs[0,0,1])**2/input_flux[0]
print("refl-order:, {}, {}, {}".format("s" if S_pol else "p",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 = int(nm_t/2)

Tsum = 0
for nm in range(nm_t):
for S_pol in [False, True]:
res = sim.get_eigenmode_coefficients(tran_flux,
mp.DiffractedPlanewave([0,nm,0],
mp.Vector3(1,0,0),
1 if S_pol else 0,
0 if S_pol else 1))
t_coeffs = res.alpha
Tmode = abs(t_coeffs[0,0,0])**2/input_flux[0]
print("tran-order:, {}, {}, {}".format("s" if S_pol else "p",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)
self.assertAlmostEqual(Rsum+Tsum,1.00,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()
Loading