Open oskooi opened 2 years ago
Can't you already pass arbitrary keyword arguments corresponding to get_eigenmode_coefficients
to the EigenmodeCoefficient
object?
Unfortunately, it's not just a simple matter of passing a DiffractedPlanewave
object as the mode
parameter of the EigenmodeCoefficient
class constructor to enable this feature because the __call__
and place_adjoint_source
functions have been set up specifically to use a band number for defining an eigenmode.
Note that for the inverse design of diffractive structures, an alternative to using an objective function based on an EigenmodeCoefficient
object with a DiffractedPlanewave
argument is to use a Near2FarFields
object (which obviously already exists and has been debugged). These are functionally equivalent (i.e., for a given diffractive structure, you know exactly where to compute its diffracted orders in the far field). Also, the computational cost of each approach is similar.
Perhaps this means that it is not absolutely necessary to add this feature after all?
Here is an attempt to verify the claim that it is possible to compute the diffracted orders of a 2D grating (a cylinder with axis along z positioned on a substrate with a unit cell that is periodic in the xy plane with normally incident planewave along z) using the near-to-far field transformation. The test involves using the Poynting theorem to verify that the sum of the transmittance in z of all the transmitted orders computed using mode decomposition is equivalent to the sum of the radial Poynting flux of the transmitted orders computed using the near-to-far field transformation. This demonstration is similar to Tutorial/Near to Far-field Spectra/Radiation Pattern of an Antenna which verified that the total flux in the far fields is independent of the actual shape of the bounding surface used to compute the flux.
Unfortunately, there is currently a large difference in the results from the two approaches because of a likely bug.
import meep as mp
import math
import numpy as np
resolution = 25 # pixels/μm
nSi = 3.45
Si = mp.Medium(index=nSi)
nSiO2 = 1.45
SiO2 = mp.Medium(index=nSiO2)
wvl = 0.5 # wavelength
fcen = 1/wvl
dpml = 2.0 # PML thickness
dsub = 3.0 # substrate thickness
dair = 3.0 # air padding
hcyl = 0.5 # cylinder height
rcyl = 0.2 # cylinder radius
sx = 1.4
sy = 0.8
sz = dpml+dsub+hcyl+dair+dpml
cell_size = mp.Vector3(sx,sy,sz)
boundary_layers = [mp.PML(thickness=dpml,direction=mp.Z)]
# periodic boundary conditions
k_point = mp.Vector3()
P_pol = True
# plane of incidence is XZ
# P/TM polarization: Ex, S/TE polarization: Ey
src_cmpt = mp.Ex if P_pol else mp.Ey
sources = [mp.Source(src=mp.GaussianSource(fcen,fwidth=0.2*fcen),
size=mp.Vector3(sx,sy,0),
center=mp.Vector3(0,0,-0.5*sz+dpml),
component=src_cmpt)]
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
sources=sources,
default_material=SiO2,
boundary_layers=boundary_layers,
k_point=k_point)
tran_pt = mp.Vector3(0,0,0.5*sz-dpml)
tran_flux = sim.add_mode_monitor(fcen,
0,
1,
mp.ModeRegion(center=tran_pt,
size=mp.Vector3(sx,sy,0)))
stop_cond = mp.stop_when_fields_decayed(10,src_cmpt,tran_pt,1e-7)
sim.run(until_after_sources=stop_cond)
input_flux = mp.get_fluxes(tran_flux)
input_flux_data = sim.get_flux_data(tran_flux)
sim.reset_meep()
geometry = [mp.Block(size=mp.Vector3(mp.inf,mp.inf,dpml+dsub),
center=mp.Vector3(0,0,-0.5*sz+0.5*(dpml+dsub)),
material=SiO2),
mp.Cylinder(height=hcyl,
radius=rcyl,
center=mp.Vector3(0,0,-0.5*sz+dpml+dsub+0.5*hcyl),
material=Si)]
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
sources=sources,
geometry=geometry,
boundary_layers=boundary_layers,
k_point=k_point)
nearfield = sim.add_near2far(fcen,
0,
1,
mp.Near2FarRegion(center=tran_pt,
size=mp.Vector3(sx,sy,0)),
nperiods=100)
tran_flux = sim.add_mode_monitor(fcen,
0,
1,
mp.ModeRegion(center=tran_pt,
size=mp.Vector3(sx,sy,0)))
sim.run(until_after_sources=stop_cond)
# length of far-field points
r = 1000*wvl
# sum the radial Poynting flux from the far fields of all transmitted orders
Pr_sum = 0
# sum the Poynting flux in z direction for all transmitted orders
Tsum = 0
# number of transmitted modes/orders in air in x and y directions (upper bound)
nm_x = int(fcen*sx) + 1
nm_y = int(fcen*sy) + 1
for m_x in range(nm_x):
for m_y in range(nm_y):
for S_pol in [False,True]:
res = sim.get_eigenmode_coefficients(tran_flux,
mp.DiffractedPlanewave([m_x,m_y,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:, {}, {}, {}, {:.6f}".format("s" if S_pol else "p",m_x,m_y,Tmode))
if m_x == 0 and m_y == 0:
Tsum += Tmode
elif (m_x != 0 and m_y == 0) or (m_x == 0 and m_y != 0):
Tsum += 2*Tmode
else:
Tsum += 4*Tmode
# determine the far-field point by scaling the wavevector of the diffracted mode
kmode = mp.Vector3(m_x/sx,m_y/sy,(fcen**2-(m_x/sx)**2-(m_y/sy)**2)**0.5)
rpt = kmode.unit().scale(r)
ff = sim.get_farfield(nearfield,rpt)
Efar = [np.conj(ff[j]) for j in range(3)]
Hfar = [ff[j+3] for j in range(3)]
Px = np.real(Efar[1]*Hfar[2]-Efar[2]*Hfar[1]) # Ey*Hz - Ez*Hy
Py = np.real(Efar[2]*Hfar[0]-Efar[0]*Hfar[2]) # Ez*Hx - Ex*Hz
Pz = np.real(Efar[0]*Hfar[1]-Efar[1]*Hfar[0]) # Ex*Hy - Ey*Hx
Pr = np.sqrt(np.square(Px)+np.square(Py)+np.square(Pz))
print("n2f:, {}, {}, {:.6f}".format(m_x,m_y,Pr))
if m_x == 0 and m_y == 0:
Pr_sum += Pr
elif (m_x != 0 and m_y == 0) or (m_x == 0 and m_y != 0):
Pr_sum += 2*Pr
else:
Pr_sum += 4*Pr
t_flux = mp.get_fluxes(tran_flux)
Tflux = t_flux[0]/input_flux[0]
# integrate the radial Poynting flux over a hemisphere
radial_flux = Pr_sum*2*np.pi*r*r/(nm_x*nm_y)
print("tran:, {:.6f}, {:.6f}, {:.6f}".format(Tsum,Tflux,radial_flux))
Note that for the inverse design of diffractive structures, an alternative to using an objective function based on an EigenmodeCoefficient object with a DiffractedPlanewave argument is to use a Near2FarFields object (which obviously already exists and has been debugged). These are functionally equivalent (i.e., for a given diffractive structure, you know exactly where to compute its diffracted orders in the far field). Also, the computational cost of each approach is similar.
This won't work. The diffracted orders are planewaves, so they fill all space above the surface — you can't look at individual points in the far field and only see individual diffracted orders.
Currently, the
EigenmodeCoefficient
objective function of the adjoint solver can only be specified using an eigenmode/band number (an integer) via itsmode
member parameter. It would be useful to extend this feature to also support aDiffractedPlanewave
object similar to what already exists for theEigenmodeSource
(via itseig_band
property) andget_eigenmode_coefficients
(via itsmode
proeprty). This would make it easier to use the adjoint solver for the design of diffractive structures.