-
Notifications
You must be signed in to change notification settings - Fork 587
/
mode-decomposition.py
90 lines (68 loc) · 3.25 KB
/
mode-decomposition.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
# -*- coding: utf-8 -*-
import meep as mp
resolution = 61 # pixels/μm
w1 = 1 # width of waveguide 1
w2 = 2 # width of waveguide 2
Lw = 10 # length of waveguide 1/2
dair = 3.0 # length of air region
dpml_x = 6.0 # length of PML in x direction
dpml_y = 2.0 # length of PML in y direction
sy = dpml_y+dair+w2+dair+dpml_y
Si = mp.Medium(epsilon=12.0)
boundary_layers = [mp.PML(dpml_x,direction=mp.X),
mp.PML(dpml_y,direction=mp.Y)]
# mode wavelength
lcen = 6.67
# mode frequency
fcen = 1/lcen
symmetries = [mp.Mirror(mp.Y)]
for m in range(5):
Lt = 2**m
sx = dpml_x+Lw+Lt+Lw+dpml_x
cell_size = mp.Vector3(sx,sy,0)
src_pt = mp.Vector3(-0.5*sx+dpml_x+0.2*Lw,0,0)
sources = [mp.EigenModeSource(src=mp.GaussianSource(fcen,fwidth=0.2*fcen),
component=mp.Ez,
center=src_pt,
size=mp.Vector3(0,sy-2*dpml_y,0),
eig_match_freq=True,
eig_parity=mp.ODD_Z+mp.EVEN_Y)]
# straight waveguide
vertices = [mp.Vector3(-0.5*sx-1,0.5*w1),
mp.Vector3(0.5*sx+1,0.5*w1),
mp.Vector3(0.5*sx+1,-0.5*w1),
mp.Vector3(-0.5*sx-1,-0.5*w1)]
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=boundary_layers,
geometry=[mp.Prism(vertices,height=mp.inf,material=Si)],
sources=sources,
symmetries=symmetries)
mon_pt = mp.Vector3(-0.5*sx+dpml_x+0.7*Lw,0,0)
flux = sim.add_flux(fcen,0,1,mp.FluxRegion(center=mon_pt,size=mp.Vector3(0,sy-2*dpml_y,0)))
sim.run(until_after_sources=mp.stop_when_fields_decayed(50,mp.Ez,mon_pt,1e-9))
incident_coeffs, vgrp, kpoints = sim.get_eigenmode_coefficients(flux,[1],eig_parity=mp.ODD_Z+mp.EVEN_Y)
incident_flux = mp.get_fluxes(flux)
incident_flux_data = sim.get_flux_data(flux)
sim.reset_meep()
# linear taper
vertices = [mp.Vector3(-0.5*sx-1,0.5*w1),
mp.Vector3(-0.5*Lt,0.5*w1),
mp.Vector3(0.5*Lt,0.5*w2),
mp.Vector3(0.5*sx+1,0.5*w2),
mp.Vector3(0.5*sx+1,-0.5*w2),
mp.Vector3(0.5*Lt,-0.5*w2),
mp.Vector3(-0.5*Lt,-0.5*w1),
mp.Vector3(-0.5*sx-1,-0.5*w1)]
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=boundary_layers,
geometry=[mp.Prism(vertices,height=mp.inf,material=Si)],
sources=sources,
symmetries=symmetries)
refl_flux = sim.add_flux(fcen,0,1,mp.FluxRegion(center=mon_pt,size=mp.Vector3(0,sy-2*dpml_y,0)))
sim.load_minus_flux_data(refl_flux,incident_flux_data)
sim.run(until_after_sources=mp.stop_when_fields_decayed(50,mp.Ez,mon_pt,1e-9))
coeffs, vgrp, kpoints = sim.get_eigenmode_coefficients(refl_flux,[1],eig_parity=mp.ODD_Z+mp.EVEN_Y)
taper_flux = mp.get_fluxes(refl_flux)
print("refl:, {}, {:.8f}, {:.8f}".format(Lt,abs(coeffs[0,0,1])**2/abs(incident_coeffs[0,0,0])**2,-taper_flux[0]/incident_flux[0]))