phase with planewave at oblique incidence

[Din-Y]

Hello! I am trying to get phase with planewave at oblique incidence. I use two approaches to get phase such as mode decomposition and dft monitor. I am not sure about correct results of mode decomposition in 3d model without symmetries so I decided to compare two approaches on simple model which consists of three homogeneous layers (Air – Si – Air) with the same thickness. Results of simulations are converged, but its different from analytic results. Could you tell me where I can be mistaken? Thanks.

import cmath
import math

import meep as mp
import numpy as np
import matplotlib.pyplot as plt

resolution = 50        # pixels/μm

dpml = 1.0             # PML thickness
dsub = 3.0             # substrate thickness
dpad = 3.0             # length of padding between grating and PML
gp = 3.0              # grating period
gh = 0.5               # grating height
gdc = 0.5              # grating duty cycle

sx = 10  # 8.5
sy = 3

cell_size = mp.Vector3(sx,sy,0)
pml_layers = [mp.PML(thickness=dpml,direction=mp.X)]

wvl_min = 0.9
wvl_max = 1.1
fmax = 1 / wvl_min
fmin = 1 / wvl_max
fcen = 0.5*(fmin+fmax)
df = fmax-fmin 
nfreq = 50

n_s = 3.42  # refrective index of silicon
n_a = 1     # refrective index of air

src_pt = mp.Vector3(-3,0,0)

def pw_amp(k,x0):
    def _pw_amp(x):
        return cmath.exp(1j*2*math.pi*k.dot(x+x0))
    return _pw_amp

# rotation angle of incident planewave
th = 45
theta_in = math.radians(th)

k = mp.Vector3(math.cos(theta_in),0,0).scale(fcen * n_a)

sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
                     component=mp.Ez,
                     center=src_pt,
                     size=mp.Vector3(0,sy,0),
                     amp_func=pw_amp(k,src_pt))]

sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    boundary_layers=pml_layers,
                    k_point=k,
                    default_material=mp.Medium(index=n_a),
                    sources=sources,
                    )

tran_pt = mp.Vector3(3,0,0)

tran_flux = sim.add_flux(fcen, df, nfreq, mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,sy,0)))

dft_mon = sim.add_dft_fields([mp.Ez],fcen, df, nfreq, center=tran_pt, size=mp.Vector3(0,sy,0),)

sim.run(until=200)

flux_freqs = mp.get_flux_freqs(tran_flux)
wl = (1/np.array(flux_freqs[:]))

dft_cmp = sim.get_dft_array(dft_mon, mp.Ez, 0)
centr_x = int(0.5 * len(dft_cmp))

my_list = list(idx for idx in range(centr_x * 2))
x = my_list

field_dft_input = []

for nf in range(nfreq):
    dft_cmp = sim.get_dft_array(dft_mon, mp.Ez, nf)
    for i in range(len(dft_cmp)):
        dft_cmp[i]= dft_cmp[i] / cmath.exp(1j * x[i] * (k.x + k.y))
    field_dft_input.append(dft_cmp[centr_x])

res_tran = sim.get_eigenmode_coefficients(tran_flux, [1], eig_parity=mp.ODD_Z + mp.EVEN_Y)
input_mode_coeff_tran = (res_tran.alpha[0, :, 0])

sim.reset_meep()

geometry = [mp.Block(material=mp.Medium(index = n_a),
                     size=mp.Vector3(2,mp.inf,mp.inf),
                     center=mp.Vector3(-2,0,0)),
           mp.Block(material=mp.Medium(index = n_s),
                     size=mp.Vector3(2,mp.inf,mp.inf),
                     center=mp.Vector3(0,0,0)),
           mp.Block(material=mp.Medium(index = n_a),
                     size=mp.Vector3(2,mp.inf,mp.inf),
                     center=mp.Vector3(2,0,0)),]

sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    boundary_layers=pml_layers,
                    geometry=geometry,
                    k_point=k,
                    sources=sources,
                    )

tran_flux = sim.add_flux(fcen, df, nfreq, mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,sy,0)))

dft_mon = sim.add_dft_fields([mp.Ez],fcen, df, nfreq, center=tran_pt, size=mp.Vector3(0,sy,0),)

sim.run(until=200)

sim.plot2D(output_plane=mp.Volume(center=mp.Vector3(0, 0), size=mp.Vector3(cell_size.x, cell_size.y, 0)))

field_dft = []

for nf in range(nfreq):
    dft_cmp1 = sim.get_dft_array(dft_mon, mp.Ez, nf)
    for i in range(len(dft_cmp1)):
        dft_cmp1[i]= dft_cmp1[i] / cmath.exp(1j * x[i] * (k.x + k.y))
    field_dft.append(dft_cmp1[centr_x])

ph_dft = np.angle(np.array(field_dft) / np.array(field_dft_input))

res_tran = sim.get_eigenmode_coefficients(tran_flux, [1], eig_parity=mp.ODD_Z + mp.EVEN_Y)
Ps = np.angle(res_tran.alpha[0,:,0] / input_mode_coeff_tran[:] )

# Analytic phase
th = 45
theta_rad = th*np.pi/180

th_s = math.asin((n_a / n_s) * math.cos(theta_rad))

ph_theor = []
for n in range(nfreq):
    ph_theor.append(np.angle(cmath.exp(1j * (flux_freqs[n] * math.cos(theta_rad) * 
                                             n_a * 4 * 2 * np.pi + 2 * np.pi * flux_freqs[n] * 
                                             math.cos(math.asin((n_a / n_s) * math.cos(theta_rad))) * n_s * 2))))

plt.figure()
plt.plot(wl, Ps, "bo-", label="Phase eiginemode_coefficients")
plt.plot(wl, ph_dft, color = "C1", marker = '*', label="Phase dft_monitor")
plt.plot(wl, ph_theor, color = "C3", marker = '*', label="Phase analytic")
plt.legend(bbox_to_anchor=( 1.05 , 1 ), loc='upper left', borderaxespad= 0 )
plt.show()