-
Notifications
You must be signed in to change notification settings - Fork 30
/
waveguide_modes.py
112 lines (89 loc) · 3.19 KB
/
waveguide_modes.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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
# ---
# jupyter:
# jupytext:
# formats: py:percent,md:myst
# text_representation:
# extension: .py
# format_name: percent
# format_version: '1.3'
# jupytext_version: 1.15.0
# kernelspec:
# display_name: Python 3
# name: python3
# ---
# %% [markdown]
# # Optical waveguide modes
# %% tags=["remove-stderr", "hide-input", "thebe-init"]
from collections import OrderedDict
import matplotlib.pyplot as plt
import numpy as np
import shapely
import shapely.affinity
from scipy.constants import epsilon_0, speed_of_light
from shapely.ops import clip_by_rect
from skfem import Basis, ElementTriP0
from skfem.io.meshio import from_meshio
from femwell.maxwell.waveguide import compute_modes
from femwell.mesh import mesh_from_OrderedDict
from femwell.visualization import plot_domains
# %% [markdown]
# We describe the geometry using shapely.
# In this case it's simple: we use a shapely.box for the waveguide.
# For the surrounding we buffer the core and clip it to the part below the waveguide for the box.
# The remaining buffer is used as the clad.
# For the core we set the resolution to 30nm and let it fall of over 500nm
# %%
wg_width = 2.5
wg_thickness = 0.3
core = shapely.geometry.box(-wg_width / 2, 0, +wg_width / 2, wg_thickness)
env = shapely.affinity.scale(core.buffer(5, resolution=8), xfact=0.5)
polygons = OrderedDict(
core=core,
box=clip_by_rect(env, -np.inf, -np.inf, np.inf, 0),
clad=clip_by_rect(env, -np.inf, 0, np.inf, np.inf),
)
resolutions = dict(core={"resolution": 0.03, "distance": 0.5})
mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=10))
mesh.draw().show()
# %% [markdown]
# Let's also plot the domains
# %%
plot_domains(mesh)
plt.show()
# %% [markdown]
# On this mesh, we define the epsilon. We do this by setting domainwise the epsilon to the squared refractive index.
# %%
basis0 = Basis(mesh, ElementTriP0())
epsilon = basis0.zeros()
for subdomain, n in {"core": 1.9963, "box": 1.444, "clad": 1}.items():
epsilon[basis0.get_dofs(elements=subdomain)] = n**2
basis0.plot(epsilon, colorbar=True).show()
# %% [markdown]
# And now we call `compute_modes` to calculate the modes of the waveguide we set up.
# As modes can have complex fields as soon as the epsilon gets complex, so we get a complex field for each mode.
# Here we show only the real part of the mode.
# %%
wavelength = 1.55
modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=2, order=2)
for mode in modes:
print(f"Effective refractive index: {mode.n_eff:.4f}")
mode.show("E", part="real", colorbar=True)
mode.show("E", part="imag", colorbar=True)
# %% [markdown]
# The intensity can be plotted directly from the mode object
# +
# %%
modes[0].show("I", colorbar=True)
# %% [markdown]
# -
# %% [markdown]
# Now, let's calculate with the modes:
# What percentage of the mode is within the core for the calculated modes?
# %%
powers_in_waveguide = []
confinement_factors_waveguide = []
for mode in modes:
powers_in_waveguide.append(mode.calculate_power(elements="core"))
confinement_factors_waveguide.append(mode.calculate_confinement_factor(elements="core"))
print(powers_in_waveguide)
print(confinement_factors_waveguide)