In [1]:
# -*- coding: utf-8 -*-

import meep as mp
import math
import cmath
import numpy as np

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 = 10.0              # grating period
gh = 0.5               # grating height
gdc = 0.5              # grating duty cycle

sx = dpml+dsub+gh+dpad+dpml
sy = gp

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

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

ng = 1.5
glass = mp.Medium(index=ng)

use_cw_solver = False  # CW solver or time stepping?
tol = 1e-6             # CW solver tolerance
max_iters = 2000       # CW solver max iterations
L = 10                 # CW solver L

# rotation angle of incident planewave; counter clockwise (CCW) about Z axis, 0 degrees along +X axis
theta_in = math.radians(10.7)

# k (in source medium) with correct length (plane of incidence: XY)
k = mp.Vector3(fcen*ng).rotate(mp.Vector3(z=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

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*sx+dpml+0.3*dsub,0,0)
sources = [mp.Source(mp.ContinuousSource(fcen,fwidth=df) if use_cw_solver else 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=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)))

if use_cw_solver:
  sim.init_sim()
  sim.solve_cw(tol, max_iters, L)
else:
  sim.run(until_after_sources=100)

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=pml_layers,
                    geometry=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)))
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)))

if use_cw_solver:
  sim.init_sim()
  sim.solve_cw(tol, max_iters, L)
else:
  sim.run(until_after_sources=200)

nm_r = np.floor((fcen*ng-k.y)*gp)-np.ceil((-fcen*ng-k.y)*gp) # number of reflected orders
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)
r_coeffs = res.alpha

Rsum = 0
for nm in range(nm_r):
  r_kdom = res.kdom[nm]
  Rmode = abs(r_coeffs[nm,0,1])**2/input_flux[0]
  r_angle = np.sign(r_kdom.y)*math.acos(r_kdom.x/(ng*fcen))
  print("refl:, {}, {:.2f}, {:.8f}".format(nm,math.degrees(r_angle),Rmode))
  Rsum += Rmode

nm_t = np.floor((fcen-k.y)*gp)-np.ceil((-fcen-k.y)*gp)       # number of transmitted orders
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)
t_coeffs = res.alpha

Tsum = 0
for nm in range(nm_t):
  t_kdom = res.kdom[nm]
  Tmode = abs(t_coeffs[nm,0,0])**2/input_flux[0]
  t_angle = np.sign(t_kdom.y)*math.acos(t_kdom.x/fcen)
  print("tran:, {}, {:.2f}, {:.8f}".format(nm,math.degrees(t_angle),Tmode))
  Tsum += Tmode

print("mode-coeff:, {:.6f}, {:.6f}, {:.6f}".format(Rsum,Tsum,Rsum+Tsum))

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("poynting-flux:, {:.6f}, {:.6f}, {:.6f}".format(Rflux,Tflux,Rflux+Tflux))
-----------
Initializing structure...
Meep: using complex fields.
Meep progress: 10.11/200.0 = 5.1% done in 4.0s, 75.2s to go
Meep progress: 19.91/200.0 = 10.0% done in 8.0s, 72.4s to go
Meep progress: 30.03/200.0 = 15.0% done in 12.0s, 67.9s to go
Meep progress: 40.33/200.0 = 20.2% done in 16.0s, 63.4s to go
Meep progress: 50.46/200.0 = 25.2% done in 20.0s, 59.3s to go
Meep progress: 60.7/200.0 = 30.4% done in 24.0s, 55.1s to go
Meep progress: 70.97/200.0 = 35.5% done in 28.0s, 50.9s to go
Meep progress: 81.13/200.0 = 40.6% done in 32.0s, 46.9s to go
Meep progress: 91.17/200.0 = 45.6% done in 36.0s, 43.0s to go
Meep progress: 101.35000000000001/200.0 = 50.7% done in 40.0s, 39.0s to go
Meep progress: 111.68/200.0 = 55.8% done in 44.0s, 34.8s to go
Meep progress: 121.94/200.0 = 61.0% done in 48.0s, 30.7s to go
Meep progress: 132.32/200.0 = 66.2% done in 52.0s, 26.6s to go
Meep progress: 142.6/200.0 = 71.3% done in 56.0s, 22.6s to go
Meep progress: 152.43/200.0 = 76.2% done in 60.0s, 18.7s to go
Meep progress: 162.74/200.0 = 81.4% done in 64.0s, 14.7s to go
Meep progress: 172.87/200.0 = 86.4% done in 68.0s, 10.7s to go
Meep progress: 183.22/200.0 = 91.6% done in 72.0s, 6.6s to go
Meep progress: 193.47/200.0 = 96.7% done in 76.0s, 2.6s to go
run 0 finished at t = 200.0 (20000 timesteps)
-----------
Initializing structure...
Meep: using complex fields.
Meep progress: 9.91/300.0 = 3.3% done in 4.0s, 117.2s to go
Meep progress: 19.97/300.0 = 6.7% done in 8.0s, 112.3s to go
Meep progress: 29.560000000000002/300.0 = 9.9% done in 12.0s, 109.9s to go
Meep progress: 39.53/300.0 = 13.2% done in 16.0s, 105.5s to go
Meep progress: 49.58/300.0 = 16.5% done in 20.0s, 101.1s to go
Meep progress: 59.56/300.0 = 19.9% done in 24.0s, 97.0s to go
Meep progress: 69.63/300.0 = 23.2% done in 28.0s, 92.7s to go
Meep progress: 79.63/300.0 = 26.5% done in 32.0s, 88.6s to go
Meep progress: 89.79/300.0 = 29.9% done in 36.0s, 84.3s to go
Meep progress: 99.65/300.0 = 33.2% done in 40.0s, 80.5s to go
Meep progress: 109.44/300.0 = 36.5% done in 44.0s, 76.7s to go
Meep progress: 119.5/300.0 = 39.8% done in 48.0s, 72.5s to go
Meep progress: 129.28/300.0 = 43.1% done in 52.0s, 68.7s to go
Meep progress: 138.88/300.0 = 46.3% done in 56.0s, 65.0s to go
Meep progress: 148.73/300.0 = 49.6% done in 60.0s, 61.1s to go
Meep progress: 158.18/300.0 = 52.7% done in 64.0s, 57.4s to go
Meep progress: 168.13/300.0 = 56.0% done in 68.0s, 53.4s to go
Meep progress: 178.14000000000001/300.0 = 59.4% done in 72.0s, 49.3s to go
Meep progress: 188.20000000000002/300.0 = 62.7% done in 76.0s, 45.2s to go
Meep progress: 198.17000000000002/300.0 = 66.1% done in 80.0s, 41.1s to go
Meep progress: 208.09/300.0 = 69.4% done in 84.0s, 37.1s to go
Meep progress: 218.17000000000002/300.0 = 72.7% done in 88.0s, 33.0s to go
Meep progress: 228.17000000000002/300.0 = 76.1% done in 92.1s, 29.0s to go
Meep progress: 238.24/300.0 = 79.4% done in 96.1s, 24.9s to go
Meep progress: 248.07/300.0 = 82.7% done in 100.1s, 20.9s to go
Meep progress: 258.12/300.0 = 86.0% done in 104.1s, 16.9s to go
Meep progress: 268.3/300.0 = 89.4% done in 108.1s, 12.8s to go
Meep progress: 278.44/300.0 = 92.8% done in 112.1s, 8.7s to go
Meep progress: 288.27/300.0 = 96.1% done in 116.1s, 4.7s to go
Meep progress: 298.33/300.0 = 99.4% done in 120.1s, 0.7s to go
run 0 finished at t = 300.0 (30000 timesteps)
refl:, 0, -0.82, 0.00005706
refl:, 1, 1.09, 0.00000404
refl:, 2, -2.73, 0.00000529
refl:, 3, 3.00, 0.00006039
refl:, 4, -4.65, 0.00005601
refl:, 5, 4.91, 0.00000845
refl:, 6, -6.57, 0.00000815
refl:, 7, 6.83, 0.00006645
refl:, 8, -8.49, 0.00005697
refl:, 9, 8.76, 0.00015744
refl:, 10, -10.43, 0.00001317
refl:, 11, 10.70, 0.04414669
refl:, 12, -12.38, 0.00005969
refl:, 13, 12.65, 0.00041534
refl:, 14, -14.34, 0.00001986
refl:, 15, 14.62, 0.00009002
refl:, 16, -16.32, 0.00006408
refl:, 17, 16.60, 0.00012194
refl:, 18, -18.32, 0.00003124
refl:, 19, 18.60, 0.00011016
refl:, 20, -20.34, 0.00006994
refl:, 21, 20.63, 0.00011303
refl:, 22, -22.40, 0.00004960
refl:, 23, 22.69, 0.00013856
refl:, 24, -24.48, 0.00007689
refl:, 25, 24.77, 0.00014294
refl:, 26, -26.59, 0.00007941
refl:, 27, 26.89, 0.00017844
refl:, 28, -28.75, 0.00008474
refl:, 29, 29.06, 0.00020622
refl:, 30, -30.95, 0.00012869
refl:, 31, 31.27, 0.00023543
refl:, 32, -33.21, 0.00009471
refl:, 33, 33.53, 0.00032429
refl:, 34, -35.52, 0.00021572
refl:, 35, 35.85, 0.00032455
refl:, 36, -37.90, 0.00011932
refl:, 37, 38.24, 0.00060212
refl:, 38, -40.37, 0.00042988
refl:, 39, 40.72, 0.00054810
refl:, 40, -42.92, 0.00031748
refl:, 41, 43.29, 0.00174005
refl:, 42, -45.59, 0.00064620
refl:, 43, 45.97, 0.00100453
refl:, 44, -48.39, 0.00034986
refl:, 45, 48.79, 0.00105225
refl:, 46, -51.35, 0.00032074
refl:, 47, 51.78, 0.00098224
refl:, 48, -54.52, 0.00044285
refl:, 49, 54.98, 0.00071074
refl:, 50, -57.96, 0.00023208
refl:, 51, 58.47, 0.00109322
refl:, 52, -61.76, 0.00035650
refl:, 53, 62.33, 0.00079058
refl:, 54, -66.11, 0.00026145
refl:, 55, 66.78, 0.00071195
refl:, 56, -71.38, 0.00002471
refl:, 57, 72.24, 0.00028160
refl:, 58, -78.81, 0.00001352
tran:, 0, -1.23, 0.00092964
tran:, 1, 1.63, 0.01528639
tran:, 2, -4.10, 0.00767314
tran:, 3, 4.50, 0.00094525
tran:, 4, -6.98, 0.00092643
tran:, 5, 7.38, 0.04285359
tran:, 6, -9.88, 0.00453034
tran:, 7, 10.28, 0.00097440
tran:, 8, -12.80, 0.00093868
tran:, 9, 13.21, 0.38618480
tran:, 10, -15.75, 0.00292543
tran:, 11, 16.17, 0.00214517
tran:, 12, -18.75, 0.00095546
tran:, 13, 19.18, 0.38260810
tran:, 14, -21.81, 0.00198520
tran:, 15, 22.24, 0.00107209
tran:, 16, -24.93, 0.00098420
tran:, 17, 25.37, 0.04148387
tran:, 18, -28.13, 0.00137337
tran:, 19, 28.59, 0.00113879
tran:, 20, -31.43, 0.00101836
tran:, 21, 31.90, 0.01425444
tran:, 22, -34.85, 0.00094045
tran:, 23, 35.35, 0.00121503
tran:, 24, -38.43, 0.00105001
tran:, 25, 38.94, 0.00671291
tran:, 26, -42.18, 0.00061482
tran:, 27, 42.73, 0.00129876
tran:, 28, -46.18, 0.00106158
tran:, 29, 46.76, 0.00353975
tran:, 30, -50.49, 0.00037006
tran:, 31, 51.12, 0.00137924
tran:, 32, -55.24, 0.00100766
tran:, 33, 55.94, 0.00185785
tran:, 34, -60.63, 0.00022601
tran:, 35, 61.46, 0.00134905
tran:, 36, -67.15, 0.00077375
tran:, 37, 68.20, 0.00106048
tran:, 38, -76.29, 0.00022316
mode-coeff:, 0.061048, 0.937868, 0.998915
poynting-flux:, 0.061102, 0.938344, 0.999447
In [ ]: