Post

Ybranch

Ybranch

Y branch

splitter and combiner

Source Code


1
2
3
4
5
6
7
# importing necessary packages
import gplugins.modes as gmode
import numpy as np
import matplotlib.pyplot as plt
import meep as mp
import gdsfactory as gf
import gplugins.gmeep as gm
1
2
Using MPI version 4.1, 1 processes
2025-10-02 17:22:26.906 | INFO     | gplugins.gmeep:<module>:39 - Meep '1.31.0' installed at ['/home/ramprakash/anaconda3/envs/si_photo/lib/python3.13/site-packages/meep']

Y-branch from gdsfactory

Source Code


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
gf.clear_cache()
ybranch = gf.Component("ybranch")

inwg = gf.components.straight(length=5)
outwg = gf.components.straight(length=2)
splitter = gf.components.mmi1x2_with_sbend()

inwg_ref = ybranch.add_ref(inwg)
splitter_ref = ybranch.add_ref(splitter)
outwg_top_ref = ybranch.add_ref(outwg)
outwg_bot_ref = ybranch.add_ref(outwg)

# Connecting ports
inwg_ref.connect("o2", splitter_ref.ports["o1"])
outwg_top_ref.connect("o1", splitter_ref.ports["o2"])
outwg_bot_ref.connect("o1", splitter_ref.ports["o3"])

# adding port name
ybranch.add_port(name=f"o1", port=inwg_ref.ports["o1"])
ybranch.add_port(name=f"o2", port=outwg_top_ref.ports["o2"])
ybranch.add_port(name=f"o3", port=outwg_bot_ref.ports["o2"])
ybranch.auto_rename_ports()

# plotting
ybranch.draw_ports()
ybranch.plot()

png

Meep simulation

S-parameter calculation and steady-state fields

50-50 splitter

Input port 1, will not be exact 50-50 because of the losses at the split

Source Code


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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
from gdsfactory.technology import LayerLevel, LayerStack

mp.verbosity(0)

# Set up frequency points for simulation
npoints = 100
lcen = 1.55
dlam = 0.02
wl = np.linspace(lcen - dlam / 2, lcen + dlam / 2, npoints)
fcen = 1 / lcen
fwidth = 3 * dlam / lcen**2
fpoints = 1 / wl  # Center frequency
mode_parity = mp.EVEN_Y + mp.ODD_Z
dpml = 1
dpad = 1
resolution = 20

# Define materials
Si = mp.Medium(index=3.45)
SiO2 = mp.Medium(index=1.45)

cell_size = mp.Vector3(ybranch.xsize + 2 * dpml, ybranch.ysize + 2 * dpml + 2 * dpad, 0)

ybranch = gf.components.extend_ports(ybranch, port_names=["o1", "o2", "o3"], length=1)
ybranch = ybranch.copy()
ybranch.flatten()
ybranch.center = (0, 0)

geometry = gm.get_meep_geometry.get_meep_geometry_from_component(ybranch)

geometry = [
    mp.Prism(geom.vertices, geom.height, geom.axis, geom.center, material=Si)
    for geom in geometry
]

src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)
source = [
    mp.EigenModeSource(
        src=src,
        eig_band=1,
        eig_parity=mode_parity,
        size=mp.Vector3(0, 1),
        center=mp.Vector3(ybranch.ports["o1"].x + dpml + 1, ybranch["o1"].y),
    )
]
sim = mp.Simulation(
    resolution=resolution,
    cell_size=cell_size,
    boundary_layers=[
        mp.PML(dpml)
    ],  # the boundary layers to absorb fields that leave the simulation
    sources=source,  # The sources
    geometry=geometry,  # The geometry
    default_material=SiO2,
)
m1 = mp.Volume(
    center=mp.Vector3(ybranch.ports["o1"].x + dpml + 1 + 0.5, ybranch["o1"].y),
    size=mp.Vector3(0, 1),
)
m2 = mp.Volume(
    center=mp.Vector3(ybranch.ports["o2"].x - dpml - 1 - 0.5, ybranch["o2"].y),
    size=mp.Vector3(0, 1),
)
m3 = mp.Volume(
    center=mp.Vector3(ybranch.ports["o3"].x - dpml - 1 - 0.5, ybranch["o3"].y),
    size=mp.Vector3(0, 1),
)

mode_monitor_1 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=m1))
mode_monitor_2 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=m2))
mode_monitor_3 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=m3))
whole_dft = sim.add_dft_fields([mp.Ez], fcen, 0, 1, center=mp.Vector3(), size=cell_size)

sim.plot2D(labels=False)

sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=1e-4))

# Finds the S parameters
norm_mode_coeff = sim.get_eigenmode_coefficients(
    mode_monitor_1, [1], eig_parity=mode_parity
).alpha[0, :, 0]
port1_coeff = (
    sim.get_eigenmode_coefficients(mode_monitor_1, [1], eig_parity=mode_parity).alpha[
        0, :, 1
    ]
    / norm_mode_coeff
)

port2_coeff = (
    sim.get_eigenmode_coefficients(mode_monitor_2, [1], eig_parity=mode_parity).alpha[
        0, :, 0
    ]
    / norm_mode_coeff
)
port3_coeff = (
    sim.get_eigenmode_coefficients(mode_monitor_3, [1], eig_parity=mode_parity).alpha[
        0, :, 0
    ]
    / norm_mode_coeff
)

# Calculates the transmittance based off of the S parameters
port1_trans = abs(port1_coeff) ** 2
port2_trans = abs(port2_coeff) ** 2
port3_trans = abs(port3_coeff) ** 2
insertion_loss_dB = 10*np.log10(port2_trans)
fig = plt.figure(figsize=(15, 4))

ax_loss = fig.add_subplot(1,2,1)
ax_loss.set_title("Insertion Loss P2")
ax_loss.plot(wl,insertion_loss_dB)
ax_loss.set_xlabel(r"Wavelength (\mu m)")
ax_loss.set_ylabel("Loss (dB)")

# ax_trans1 = fig.add_subplot(1, 3, 1)
# ax_trans1.set_title("Transmission per Port")
# ax_trans1.plot(wl, port1_trans, label=r"s11^2")
# ax_trans1.plot(wl, port2_trans, label=r"s21^2")
# ax_trans1.plot(wl, port3_trans, label=r"s31^2")
# ax_trans1.set_xlabel(r"Wavelength (\mu m)")
# ax_trans1.set_ylabel(r"Transmission")
# ax_trans1.legend()

# ax_phase1 = fig.add_subplot(1, 3, 2)
# ax_phase1.set_title("Phase per Port")
# ax_phase1.plot(
#     wl, np.unwrap(np.angle(port1_coeff) * 180 / np.pi), label=r"phase of s11"
# )
# ax_phase1.plot(
#     wl, np.unwrap(np.angle(port2_coeff) * 180 / np.pi), label=r"phase of s21"
# )
# ax_phase1.plot(
#     wl, np.unwrap(np.angle(port3_coeff) * 180 / np.pi), label=r"phase of s31"
# )
# ax_phase1.set_xlabel(r"Wavelength (\mu m)")
# ax_phase1.set_ylabel("Angle (deg)")
# ax_phase1.legend()

# sim.reset_meep()

# sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=1e-4))

eps_data = sim.get_epsilon()  #
ez_data = sim.get_dft_array(
    whole_dft, mp.Ez, 0
)  # Values for the component of the E-field in the z direction (in/out of screen)

# Creates the plot
ax_field = fig.add_subplot(1, 2, 2)
ax_field.set_title("Steady State Fields")
ax_field.imshow(np.transpose(eps_data), interpolation="spline36", cmap="binary")
ax_field.imshow(
    np.flipud(np.transpose(np.real(ez_data))),
    interpolation="spline36",
    cmap="RdBu",
    alpha=0.9,
)
ax_field.axis("off")
plt.show()

png

png

Input in single guide

Splitter

Source Code


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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
from gdsfactory.technology import LayerLevel, LayerStack

mp.verbosity(0)

# Set up frequency points for simulation
npoints = 100
lcen = 1.55
dlam = 0.02
wl = np.linspace(lcen - dlam / 2, lcen + dlam / 2, npoints)
fcen = 1 / lcen
fwidth = 3 * dlam / lcen**2
fpoints = 1 / wl  # Center frequency
mode_parity = mp.EVEN_Y + mp.ODD_Z
dpml = 1
dpad = 1
resolution = 20

# Define materials
Si = mp.Medium(index=3.45)
SiO2 = mp.Medium(index=1.45)

cell_size = mp.Vector3(ybranch.xsize + 2 * dpml, ybranch.ysize + 2 * dpml + 2 * dpad, 0)

ybranch = gf.components.extend_ports(ybranch, port_names=["o1", "o2", "o3"], length=1)
ybranch = ybranch.copy()
ybranch.flatten()
ybranch.center = (0, 0)

geometry = gm.get_meep_geometry.get_meep_geometry_from_component(ybranch)

geometry = [
    mp.Prism(geom.vertices, geom.height, geom.axis, geom.center, material=Si)
    for geom in geometry
]

src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)
source = [
    mp.EigenModeSource(
        src=src,
        eig_band=1,
        eig_parity=mode_parity,
        eig_kpoint = mp.Vector3(-1,0,0),
        direction = mp.NO_DIRECTION,
        size=mp.Vector3(0, 1),
        center=mp.Vector3(ybranch.ports["o2"].x - dpml - 1, ybranch["o2"].y),
        amplitude=1
    ),
    # mp.EigenModeSource(
    #     src=src,
    #     eig_band=1,
    #     eig_parity=mode_parity,
    #     eig_kpoint = mp.Vector3(-1,0,0),
    #     direction = mp.NO_DIRECTION,
    #     size=mp.Vector3(0, 1),
    #     center=mp.Vector3(ybranch.ports["o3"].x - dpml - 1, ybranch["o3"].y),
    #     amplitude=1
    # ),
]
sim = mp.Simulation(
    resolution=resolution,
    cell_size=cell_size,
    boundary_layers=[mp.PML(dpml)],  # 
    sources=source,  # The sources
    geometry=geometry,  # The geometry
    default_material=SiO2,
)
m1 = mp.Volume(
    center=mp.Vector3(ybranch.ports["o1"].x + dpml + 1 + 0.5, ybranch["o1"].y),
    size=mp.Vector3(0, 1),
)
m2 = mp.Volume(
    center=mp.Vector3(ybranch.ports["o2"].x - dpml - 1 - 0.5, ybranch["o2"].y),
    size=mp.Vector3(0, 1),
)
m3 = mp.Volume(
    center=mp.Vector3(ybranch.ports["o3"].x - dpml - 1 - 0.5, ybranch["o3"].y),
    size=mp.Vector3(0, 1),
)

mode_monitor_1 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=m1))
mode_monitor_2 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=m2))
mode_monitor_3 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=m3))
whole_dft = sim.add_dft_fields([mp.Ez], fcen, 0, 1, center=mp.Vector3(), size=cell_size)

sim.plot2D(labels=False)

sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=1e-4))

# Finds the S parameters
norm_mode_coeff = sim.get_eigenmode_coefficients(
    mode_monitor_2, [1], eig_parity=mode_parity
).alpha[0, :, 1]  # Adding both the ports input values

port1_coeff = (
    sim.get_eigenmode_coefficients(mode_monitor_1, [1], eig_parity=mode_parity).alpha[
        0, :, 1
    ]
    / norm_mode_coeff
)

port2_coeff = (
    sim.get_eigenmode_coefficients(mode_monitor_2, [1], eig_parity=mode_parity).alpha[
        0, :, 0
    ]
    / norm_mode_coeff
)
port3_coeff = (
    sim.get_eigenmode_coefficients(mode_monitor_3, [1], eig_parity=mode_parity).alpha[
        0, :, 0
    ]
    / norm_mode_coeff
)

# Calculates the transmittance based off of the S parameters
port1_trans = abs(port1_coeff) ** 2
port2_trans = abs(port2_coeff) ** 2
port3_trans = abs(port3_coeff) ** 2
insertion_loss_dB = 10*np.log10(port1_trans)

fig = plt.figure(figsize=(15, 4))

ax_loss = fig.add_subplot(1,2,1)
ax_loss.set_title("Insertion Loss P1")
ax_loss.plot(wl,insertion_loss_dB)
ax_loss.set_xlabel(r"Wavelength (\mu m)")
ax_loss.set_ylabel("Loss (dB)")

# ax_trans1 = fig.add_subplot(1, 3, 1)
# ax_trans1.set_title("Transmission per Port")
# ax_trans1.plot(wl, port1_trans, label=r"s1(23)^2") # Not the proper way to represent 
# ax_trans1.plot(wl, port2_trans, label=r"s22^2")
# ax_trans1.plot(wl, port3_trans, label=r"s33^2")
# ax_trans1.set_xlabel(r"Wavelength (\mu m)")
# ax_trans1.set_ylabel(r"Transmission")
# ax_trans1.legend()

# ax_phase1 = fig.add_subplot(1, 3, 2)
# ax_phase1.set_title("Phase per Port")
# ax_phase1.plot(
#     wl, np.unwrap(np.angle(port1_coeff) * 180 / np.pi), label=r"phase of s1(23)"
# )
# ax_phase1.plot(
#     wl, np.unwrap(np.angle(port2_coeff) * 180 / np.pi), label=r"phase of s22"
# )
# ax_phase1.plot(
#     wl, np.unwrap(np.angle(port3_coeff) * 180 / np.pi), label=r"phase of s33"
# )
# ax_phase1.set_xlabel(r"Wavelength (\mu m)")
# ax_phase1.set_ylabel("Angle (deg)")
# ax_phase1.legend()

# sim.reset_meep()

# sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=1e-4))

eps_data = sim.get_epsilon()  #
ez_data = sim.get_dft_array(
    whole_dft, mp.Ez, 0
)  # Values for the component of the E-field in the z direction (in/out of screen)

# Creates the plot
ax_field = fig.add_subplot(1, 2, 2)
ax_field.set_title("Steady State Fields")
ax_field.imshow(np.transpose(eps_data), interpolation="spline36", cmap="binary")
ax_field.imshow(
    np.flipud(np.transpose(np.real(ez_data))),
    interpolation="spline36",
    cmap="RdBu",
    alpha=0.9,
)
ax_field.axis("off")
plt.show()

png

png

Combiner

Out-phase

distructive interference

Source Code


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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
from gdsfactory.technology import LayerLevel, LayerStack

mp.verbosity(0)

# Set up frequency points for simulation
npoints = 100
lcen = 1.55
dlam = 0.02
wl = np.linspace(lcen - dlam / 2, lcen + dlam / 2, npoints)
fcen = 1 / lcen
fwidth = 3 * dlam / lcen**2
fpoints = 1 / wl  # Center frequency
mode_parity = mp.EVEN_Y + mp.ODD_Z
dpml = 1
dpad = 1
resolution = 20

# Define materials
Si = mp.Medium(index=3.45)
SiO2 = mp.Medium(index=1.45)

cell_size = mp.Vector3(ybranch.xsize + 2 * dpml, ybranch.ysize + 2 * dpml + 2 * dpad, 0)

ybranch = gf.components.extend_ports(ybranch, port_names=["o1", "o2", "o3"], length=1)
ybranch = ybranch.copy()
ybranch.flatten()
ybranch.center = (0, 0)

geometry = gm.get_meep_geometry.get_meep_geometry_from_component(ybranch)

geometry = [
    mp.Prism(geom.vertices, geom.height, geom.axis, geom.center, material=Si)
    for geom in geometry
]

src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)
source = [
    mp.EigenModeSource(
        src=src,
        eig_band=1,
        eig_parity=mode_parity,
        eig_kpoint = mp.Vector3(-1,0,0),
        direction = mp.NO_DIRECTION,
        size=mp.Vector3(0, 1),
        center=mp.Vector3(ybranch.ports["o2"].x - dpml - 1, ybranch["o2"].y),
        amplitude=1
    ),
    mp.EigenModeSource(
        src=src,
        eig_band=1,
        eig_parity=mode_parity,
        eig_kpoint = mp.Vector3(-1,0,0),
        direction = mp.NO_DIRECTION,
        size=mp.Vector3(0, 1),
        center=mp.Vector3(ybranch.ports["o3"].x - dpml - 1, ybranch["o3"].y),
        amplitude=-1
    ),
]
sim = mp.Simulation(
    resolution=resolution,
    cell_size=cell_size,
    boundary_layers=[mp.PML(dpml)],  # 
    sources=source,  # The sources
    geometry=geometry,  # The geometry
    default_material=SiO2,
)
m1 = mp.Volume(
    center=mp.Vector3(ybranch.ports["o1"].x + dpml + 1 + 0.5, ybranch["o1"].y),
    size=mp.Vector3(0, 1),
)
m2 = mp.Volume(
    center=mp.Vector3(ybranch.ports["o2"].x - dpml - 1 - 0.5, ybranch["o2"].y),
    size=mp.Vector3(0, 1),
)
m3 = mp.Volume(
    center=mp.Vector3(ybranch.ports["o3"].x - dpml - 1 - 0.5, ybranch["o3"].y),
    size=mp.Vector3(0, 1),
)

mode_monitor_1 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=m1))
mode_monitor_2 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=m2))
mode_monitor_3 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=m3))
whole_dft = sim.add_dft_fields([mp.Ez], fcen, 0, 1, center=mp.Vector3(), size=cell_size)

sim.plot2D(labels=False)

sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=1e-4))

# Finds the S parameters
norm_mode_coeff = sim.get_eigenmode_coefficients(
    mode_monitor_2, [1], eig_parity=mode_parity
).alpha[0, :, 1] + sim.get_eigenmode_coefficients(
    mode_monitor_3, [1], eig_parity=mode_parity
).alpha[0, :, 1] # Adding both the ports input values

port1_coeff = (
    sim.get_eigenmode_coefficients(mode_monitor_1, [1], eig_parity=mode_parity).alpha[
        0, :, 1
    ]
    / norm_mode_coeff
)

port2_coeff = (
    sim.get_eigenmode_coefficients(mode_monitor_2, [1], eig_parity=mode_parity).alpha[
        0, :, 0
    ]
    / norm_mode_coeff
)
port3_coeff = (
    sim.get_eigenmode_coefficients(mode_monitor_3, [1], eig_parity=mode_parity).alpha[
        0, :, 0
    ]
    / norm_mode_coeff
)

# Calculates the transmittance based off of the S parameters
port1_trans = abs(port1_coeff) ** 2
port2_trans = abs(port2_coeff) ** 2
port3_trans = abs(port3_coeff) ** 2
insertion_loss_dB = 10*np.log10(port1_trans)
fig = plt.figure(figsize=(15, 4))

ax_loss = fig.add_subplot(1,2,1)
ax_loss.set_title("Insertion Loss P1")
ax_loss.plot(wl,insertion_loss_dB)
ax_loss.set_xlabel(r"Wavelength (\mu m)")
ax_loss.set_ylabel("Loss (dB)")

# ax_trans1 = fig.add_subplot(1, 3, 1)
# ax_trans1.set_title("Transmission per Port")
# ax_trans1.plot(wl, port1_trans, label=r"s1(23)^2") # Not the proper way to represent 
# ax_trans1.plot(wl, port2_trans, label=r"s22^2")
# ax_trans1.plot(wl, port3_trans, label=r"s33^2")
# ax_trans1.set_xlabel(r"Wavelength (\mu m)")
# ax_trans1.set_ylabel(r"Transmission")
# ax_trans1.legend()

# ax_phase1 = fig.add_subplot(1, 3, 2)
# ax_phase1.set_title("Phase per Port")
# ax_phase1.plot(
#     wl, np.unwrap(np.angle(port1_coeff) * 180 / np.pi), label=r"phase of s1(23)"
# )
# ax_phase1.plot(
#     wl, np.unwrap(np.angle(port2_coeff) * 180 / np.pi), label=r"phase of s22"
# )
# ax_phase1.plot(
#     wl, np.unwrap(np.angle(port3_coeff) * 180 / np.pi), label=r"phase of s33"
# )
# ax_phase1.set_xlabel(r"Wavelength (\mu m)")
# ax_phase1.set_ylabel("Angle (deg)")
# ax_phase1.legend()

# sim.reset_meep()

# sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=1e-4))

eps_data = sim.get_epsilon()  #
ez_data = sim.get_dft_array(
    whole_dft, mp.Ez, 0
)  # Values for the component of the E-field in the z direction (in/out of screen)

# Creates the plot
ax_field = fig.add_subplot(1, 2, 2)
ax_field.set_title("Steady State Fields")
ax_field.imshow(np.transpose(eps_data), interpolation="spline36", cmap="binary")
ax_field.imshow(
    np.flipud(np.transpose(np.real(ez_data))),
    interpolation="spline36",
    cmap="RdBu",
    alpha=0.9,
)
ax_field.axis("off")
plt.show()

png

png

In-phase

constructive interference

Source Code


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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
from gdsfactory.technology import LayerLevel, LayerStack

mp.verbosity(0)

# Set up frequency points for simulation
npoints = 100
lcen = 1.55
dlam = 0.02
wl = np.linspace(lcen - dlam / 2, lcen + dlam / 2, npoints)
fcen = 1 / lcen
fwidth = 3 * dlam / lcen**2
fpoints = 1 / wl  # Center frequency
mode_parity = mp.EVEN_Y + mp.ODD_Z
dpml = 1
dpad = 1
resolution = 20

# Define materials
Si = mp.Medium(index=3.45)
SiO2 = mp.Medium(index=1.45)

cell_size = mp.Vector3(ybranch.xsize + 2 * dpml, ybranch.ysize + 2 * dpml + 2 * dpad, 0)

ybranch = gf.components.extend_ports(ybranch, port_names=["o1", "o2", "o3"], length=1)
ybranch = ybranch.copy()
ybranch.flatten()
ybranch.center = (0, 0)

geometry = gm.get_meep_geometry.get_meep_geometry_from_component(ybranch)

geometry = [
    mp.Prism(geom.vertices, geom.height, geom.axis, geom.center, material=Si)
    for geom in geometry
]

src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)
source = [
    mp.EigenModeSource(
        src=src,
        eig_band=1,
        eig_parity=mode_parity,
        eig_kpoint = mp.Vector3(-1,0,0),
        direction = mp.NO_DIRECTION,
        size=mp.Vector3(0, 1),
        center=mp.Vector3(ybranch.ports["o2"].x - dpml - 1, ybranch["o2"].y),
        amplitude=1
    ),
    mp.EigenModeSource(
        src=src,
        eig_band=1,
        eig_parity=mode_parity,
        eig_kpoint = mp.Vector3(-1,0,0),
        direction = mp.NO_DIRECTION,
        size=mp.Vector3(0, 1),
        center=mp.Vector3(ybranch.ports["o3"].x - dpml - 1, ybranch["o3"].y),
        amplitude=1
    ),
]
sim = mp.Simulation(
    resolution=resolution,
    cell_size=cell_size,
    boundary_layers=[mp.PML(dpml)],  # 
    sources=source,  # The sources
    geometry=geometry,  # The geometry
    default_material=SiO2,
)
m1 = mp.Volume(
    center=mp.Vector3(ybranch.ports["o1"].x + dpml + 1 + 0.5, ybranch["o1"].y),
    size=mp.Vector3(0, 1),
)
m2 = mp.Volume(
    center=mp.Vector3(ybranch.ports["o2"].x - dpml - 1 - 0.5, ybranch["o2"].y),
    size=mp.Vector3(0, 1),
)
m3 = mp.Volume(
    center=mp.Vector3(ybranch.ports["o3"].x - dpml - 1 - 0.5, ybranch["o3"].y),
    size=mp.Vector3(0, 1),
)

mode_monitor_1 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=m1))
mode_monitor_2 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=m2))
mode_monitor_3 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=m3))
whole_dft = sim.add_dft_fields([mp.Ez], fcen, 0, 1, center=mp.Vector3(), size=cell_size)

sim.plot2D(labels=False)

sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=1e-4))

# Finds the S parameters
norm_mode_coeff = sim.get_eigenmode_coefficients(
    mode_monitor_2, [1], eig_parity=mode_parity
).alpha[0, :, 1] + sim.get_eigenmode_coefficients(
    mode_monitor_3, [1], eig_parity=mode_parity
).alpha[0, :, 1] # Adding both the ports input values

port1_coeff = (
    sim.get_eigenmode_coefficients(mode_monitor_1, [1], eig_parity=mode_parity).alpha[
        0, :, 1
    ]
    / norm_mode_coeff
)

port2_coeff = (
    sim.get_eigenmode_coefficients(mode_monitor_2, [1], eig_parity=mode_parity).alpha[
        0, :, 0
    ]
    / norm_mode_coeff
)
port3_coeff = (
    sim.get_eigenmode_coefficients(mode_monitor_3, [1], eig_parity=mode_parity).alpha[
        0, :, 0
    ]
    / norm_mode_coeff
)

# Calculates the transmittance based off of the S parameters
port1_trans = abs(port1_coeff) ** 2
port2_trans = abs(port2_coeff) ** 2
port3_trans = abs(port3_coeff) ** 2
insertion_loss_dB = 10*np.log10(port1_trans)

fig = plt.figure(figsize=(15, 4))

ax_loss = fig.add_subplot(1,2,1)
ax_loss.set_title("Insertion Loss P1")
ax_loss.plot(wl,insertion_loss_dB)
ax_loss.set_xlabel(r"Wavelength (\mu m)")
ax_loss.set_ylabel("Loss (dB)")

# ax_trans1 = fig.add_subplot(1, 3, 1)
# ax_trans1.set_title("Transmission per Port")
# ax_trans1.plot(wl, port1_trans, label=r"s1(23)^2") # Not the proper way to represent 
# ax_trans1.plot(wl, port2_trans, label=r"s22^2")
# ax_trans1.plot(wl, port3_trans, label=r"s33^2")
# ax_trans1.set_xlabel(r"Wavelength (\mu m)")
# ax_trans1.set_ylabel(r"Transmission")
# ax_trans1.legend()

# ax_phase1 = fig.add_subplot(1, 3, 2)
# ax_phase1.set_title("Phase per Port")
# ax_phase1.plot(
#     wl, np.unwrap(np.angle(port1_coeff) * 180 / np.pi), label=r"phase of s1(23)"
# )
# ax_phase1.plot(
#     wl, np.unwrap(np.angle(port2_coeff) * 180 / np.pi), label=r"phase of s22"
# )
# ax_phase1.plot(
#     wl, np.unwrap(np.angle(port3_coeff) * 180 / np.pi), label=r"phase of s33"
# )
# ax_phase1.set_xlabel(r"Wavelength (\mu m)")
# ax_phase1.set_ylabel("Angle (deg)")
# ax_phase1.legend()

# sim.reset_meep()

# sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=1e-4))

eps_data = sim.get_epsilon()  #
ez_data = sim.get_dft_array(
    whole_dft, mp.Ez, 0
)  # Values for the component of the E-field in the z direction (in/out of screen)

# Creates the plot
ax_field = fig.add_subplot(1, 2, 2)
ax_field.set_title("Steady State Fields")
ax_field.imshow(np.transpose(eps_data), interpolation="spline36", cmap="binary")
ax_field.imshow(
    np.flipud(np.transpose(np.real(ez_data))),
    interpolation="spline36",
    cmap="RdBu",
    alpha=0.9,
)
ax_field.axis("off")
plt.show()

png

png

Source Code


1
2
3
4
5
6
7
8
9
10
11
12
norm_mode_coeff = sim.get_eigenmode_coefficients(
    mode_monitor_3, [1], eig_parity=mode_parity
).alpha[0, :, 1] 
port1_coeff = (
    sim.get_eigenmode_coefficients(mode_monitor_1, [1], eig_parity=mode_parity).alpha[
        0, :, 1
    ]
    / norm_mode_coeff
)
port1_trans = abs(port1_coeff) ** 2
insertion_loss_dB = 10*np.log10(port1_trans)
plt.plot(wl, insertion_loss_dB)
1
[<matplotlib.lines.Line2D at 0x7b118fc26fd0>]

png

Source Code


1
plt.plot(wl, port1_trans)
1
[<matplotlib.lines.Line2D at 0x7b118fca1d10>]

png

This post is licensed under CC BY 4.0 by the author.