fdtd-waveguide-simulation/mode_solver.py at main · vortex2jm/fdtd-waveguide-simulation · GitHub

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
import meep as mp
import numpy as np
import matplotlib.pyplot as plt
from meep import mpb

mp.verbosity(0)

# -------- Material parameters --------
n_core = 2.5   # Silicon
n_clad = 1.44  # (SiO2)
core = mp.Medium(index=n_core)
cladding = mp.Medium(index=n_clad)

# -------- Simulation geometry --------
resolution = 100  # pixels/μm
geometry_lattice = mp.Lattice(size=mp.Vector3(3, 2))  # simulation size (microns)

w = 0.7   # width
h = 0.25  # height

geometry = [
    mp.Block(center=mp.Vector3(), size=mp.Vector3(mp.inf, mp.inf), material=cladding),
    mp.Block(center=mp.Vector3(), size=mp.Vector3(w, h), material=core)
]

# -------- Optical parameters --------
wl0 = 1.55         # wave length
freq = 1.0 / wl0   # freq
num_modes = 1      # amount of modes to calculate


# -------- Solver config --------
ms = mpb.ModeSolver(
    geometry_lattice=geometry_lattice,
    geometry=geometry,
    resolution=resolution,
    num_bands=num_modes,
)

ms.init_params(mp.NO_PARITY, True)

# -------- Permittivity view --------
eps = ms.get_epsilon()
x = np.arange(eps.shape[0]) / resolution
y = np.arange(eps.shape[1]) / resolution

plt.contourf(x, y, eps.transpose(), cmap='binary')
plt.colorbar(label="Permittivity")
plt.xlabel("x (μm)")
plt.ylabel("y (μm)")
plt.title("Distribuição da Permissividade (ε)")
plt.savefig("Permittivity.png")
plt.close()

# -------- Mode calc --------
n_mode_guess = 0.5 * (n_core + n_clad)
k_guess = freq * n_mode_guess
k_min = freq * n_clad
k_max = freq * n_core
tol = 1e-5
neff = []
cmap = 'viridis'

for mode_num in range(1, num_modes + 1):
    # Finds the mode
    k_mpb = ms.find_k(mp.NO_PARITY, freq, mode_num, mode_num, mp.Vector3(0, 0, 1), tol, k_guess, k_min, k_max)
    neff_val = k_mpb[0] / freq
    neff.append(neff_val)

    # Fields
    E = ms.get_efield(which_band=mode_num)
    H = ms.get_hfield(which_band=mode_num)
    P = ms.get_poynting(which_band=mode_num)

    Ex = E[:,:,0,0]
    Ey = E[:,:,0,1]
    Ez = E[:,:,0,2]
    Pz = 0.5 * np.real(P[:, :, 0, 2])

    plt.contourf(x, y, np.abs(Ex.transpose()), 202, cmap=cmap)
    plt.colorbar()
    plt.xlabel("x (μm)")
    plt.ylabel("y (μm)")
    plt.title(f"Ex do modo {mode_num}")
    plt.savefig(f"Ex_mode_{mode_num}.png")
    plt.close()

    plt.contourf(x, y, np.abs(Ey.transpose()), 202, cmap=cmap)
    plt.colorbar()
    plt.xlabel("x (μm)")
    plt.ylabel("y (μm)")
    plt.title(f"Ey do modo {mode_num}")
    plt.savefig(f"Ey_mode_{mode_num}.png")
    plt.close()

    plt.contourf(x, y, np.abs(Ez.transpose()), 202, cmap=cmap)
    plt.colorbar()
    plt.xlabel("x (μm)")
    plt.ylabel("y (μm)")
    plt.title(f"Ez do modo {mode_num}")
    plt.savefig(f"Ez_mode_{mode_num}.png")
    plt.close()

    plt.contourf(x, y, Pz.transpose(), 202, cmap=cmap)
    plt.colorbar()
    plt.xlabel("x (μm)")
    plt.ylabel("y (μm)")
    plt.title(f"Poynting z do modo {mode_num}")
    plt.savefig(f"Poynting_mode_{mode_num}.png")
    plt.close()