Merge pull request #44 from HelgeGehring/38-object-oriented-mode-solver

object oriented mode solver

docs/benchmarks/mode_solver.py

@@ -30,8 +30,8 @@
 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.mode_solver import compute_modes
 
 core_width = 3
 core_thickness = 1
@@ -86,9 +86,9 @@
     for subdomain, n in {"core": 3.44, "box": 3.40, "clad": 1}.items():
         epsilon[basis0.get_dofs(elements=subdomain)] = n**2
 
-    lams, basis, xs = compute_modes(basis0, epsilon, wavelength=1.15, num_modes=1, order=2)
+    modes = compute_modes(basis0, epsilon, wavelength=1.15, num_modes=1, order=2)
 
-    neff_values_femwell.append(np.real(lams[0]))
+    neff_values_femwell.append(np.real(modes[0].n_eff))
 
 pd.DataFrame(
     {

docs/photonics/examples/bent_waveguide.py

@@ -19,20 +19,14 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
-import scipy.constants
 from shapely import box
 from shapely.ops import clip_by_rect
 from skfem import Basis, ElementDG, ElementTriP1
 from skfem.io.meshio import from_meshio
 from tqdm import tqdm
 
+from femwell.maxwell.waveguide import compute_modes
 from femwell.mesh import mesh_from_OrderedDict
-from femwell.mode_solver import (
-    calculate_hfield,
-    calculate_overlap,
-    compute_modes,
-    plot_mode,
-)
 
 # -
 
@@ -93,15 +87,9 @@
 # and for mode overlap calculations between straight and bent waveguides.
 
 # +
-lams_straight, basis_straight, xs_straight = compute_modes(
+modes_straight = compute_modes(
     basis0, epsilon, wavelength=wavelength, num_modes=1, order=2, radius=np.inf
 )
-H_straight = calculate_hfield(
-    basis_straight,
-    xs_straight[0],
-    2 * np.pi / wavelength * lams_straight[0],
-    omega=2 * np.pi / wavelength * scipy.constants.speed_of_light,
-)
 # -
 
 # Now we calculate the modes of bent waveguides with different radii.
@@ -112,9 +100,9 @@
 radiuss = np.linspace(40, 5, 21)
 radiuss_lams = []
 overlaps = []
-lam_guess = lams_straight[0]
+lam_guess = modes_straight[0].n_eff
 for radius in tqdm(radiuss):
-    lams, basis, xs = compute_modes(
+    modes = compute_modes(
         basis0,
         epsilon,
         wavelength=wavelength,
@@ -124,18 +112,10 @@
         n_guess=lam_guess,
         solver="scipy",
     )
-    lam_guess = lams[0]
-    H_bent = calculate_hfield(
-        basis_straight,
-        xs[0],
-        2 * np.pi / wavelength * lams[0],
-        omega=2 * np.pi / wavelength * scipy.constants.speed_of_light,
-    )
-    radiuss_lams.append(lams[0])
+    lam_guess = modes[0].n_eff
+    radiuss_lams.append(modes[0].n_eff)
 
-    overlaps.append(
-        calculate_overlap(basis_straight, xs[0], H_bent, basis_straight, xs_straight[0], H_straight)
-    )
+    overlaps.append(modes_straight[0].calculate_overlap(modes[0]))
 # -
 
 # And now we plot it!
@@ -164,7 +144,7 @@
 # Here we show only the real part of the mode.
 
 # + tags=["hide-input"]
-for i, lam in enumerate(lams):
-    print(f"Effective refractive index: {lam:.14f}")
-    plot_mode(basis, xs[i].real, colorbar=True, direction="x")
+for mode in modes:
+    print(f"Effective refractive index: {mode.n_eff:.14f}")
+    mode.plot(mode.E.real, colorbar=True, direction="x")
     plt.show()

docs/photonics/examples/coupled_mode_theory.py

@@ -28,6 +28,7 @@
 
 # %% tags=["hide-input"]
 from collections import OrderedDict
+from itertools import chain
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -37,14 +38,8 @@
 from skfem import Basis, ElementTriP0, Mesh
 from skfem.io import from_meshio
 
+from femwell.maxwell.waveguide import compute_modes
 from femwell.mesh import mesh_from_OrderedDict
-from femwell.mode_solver import (
-    calculate_coupling_coefficient,
-    calculate_hfield,
-    calculate_overlap,
-    compute_modes,
-    plot_mode,
-)
 
 # %% [markdown]
 # Let's set up the geometry!
@@ -119,17 +114,18 @@
 epsilon = basis0.zeros() + 1.444**2
 epsilon[basis0.get_dofs(elements=("core_1", "core_2"))] = 3.4777**2
 # basis0.plot(epsilon, colorbar=True).show()
-lams_both, basis, xs_both = compute_modes(
-    basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=2
+modes_both = compute_modes(basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=2)
+modes_both[0].show(modes_both[0].E.real, direction="x")
+modes_both[1].show(modes_both[1].E.real, direction="x")
+print(
+    "Refractive index of symmetric and assymetric mode:",
+    modes_both[0].n_eff,
+    ", ",
+    modes_both[1].n_eff,
 )
-plot_mode(basis, np.real(xs_both[0]), direction="x")
-plt.show()
-plot_mode(basis, np.real(xs_both[1]), direction="x")
-plt.show()
-print("Refractive index of symmetric and assymetric mode:", lams_both)
 # https://www.fiberoptics4sale.com/blogs/wave-optics/directional-couplers
 print(
-    f"Maximum power transfer after {np.pi / (2 * np.pi / wavelength * np.real(lams_both[0] - lams_both[1]))} um prop length"
+    f"Maximum power transfer after {np.pi / (2 * np.pi / wavelength * np.real(modes_both[0].n_eff - modes_both[1].n_eff))} um prop length"
 )
 
 # %% [markdown]
@@ -139,60 +135,44 @@
 epsilon = basis0.zeros() + 1.444**2
 epsilon[basis0.get_dofs(elements="core_1")] = 3.4777**2
 # basis0.plot(epsilon, colorbar=True).show()
-lams_1, basis, xs_1 = compute_modes(basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=1)
-print("Effective refractive index of the mode of the first waveguide", lams_1)
-plot_mode(basis, np.real(xs_1[0]), direction="x")
-plt.show()
+modes_1 = compute_modes(basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=1)
+print("Effective refractive index of the mode of the first waveguide", modes_1[0].n_eff)
+modes_1[0].show(modes_1[0].E.real, direction="x")
 
 epsilon_2 = basis0.zeros() + 1.444**2
 epsilon_2[basis0.get_dofs(elements="core_2")] = 3.4777**2
 # basis0.plot(epsilon_2, colorbar=True).show()
-lams_2, basis, xs_2 = compute_modes(basis0, epsilon_2, wavelength=wavelength, mu_r=1, num_modes=1)
-print("Effective refractive index of the mode of the second waveguide", lams_2)
-plot_mode(basis, np.real(xs_2[0]), direction="x")
-plt.show()
+modes_2 = compute_modes(basis0, epsilon_2, wavelength=wavelength, mu_r=1, num_modes=1)
+print("Effective refractive index of the mode of the second waveguide", modes_2[0].n_eff)
+modes_2[0].show(modes_2[0].E.real, direction="x")
 
 # %%
 length = 200
 ts = np.linspace(0, length, 1000)
 
 # %%
 epsilons = [epsilon, epsilon_2]
-modes = [(lam, x, 0) for lam, x in zip(lams_1, xs_1)] + [
-    (lam, x, 1) for lam, x in zip(lams_2, xs_2)
-]
-
-overlap_integrals = np.zeros((len(modes), len(modes)), dtype=complex)
-for i, (lam_i, E_i, epsilon_i) in enumerate(modes):
-    for j, (lam_j, E_j, epsilon_j) in enumerate(modes):
-        H_i = calculate_hfield(
-            basis,
-            E_i,
-            lam_i * (2 * np.pi / 1.55),
-            omega=2 * np.pi / wavelength * speed_of_light,
-        )
-        H_j = calculate_hfield(
-            basis,
-            E_j,
-            lam_j * (2 * np.pi / 1.55),
-            omega=2 * np.pi / wavelength * speed_of_light,
-        )
-        overlap_integrals[i, j] = calculate_overlap(basis, E_i, H_i, basis, E_j, H_j)
+
+num_modes = len(modes_1) + len(modes_2)
+overlap_integrals = np.zeros((num_modes, num_modes), dtype=complex)
+for i, mode_i in enumerate(chain(modes_1, modes_2)):
+    for j, mode_j in enumerate(chain(modes_1, modes_2)):
+        overlap_integrals[i, j] = mode_i.calculate_overlap(mode_j)
 
 print("overlap", overlap_integrals)
 # plt.imshow(np.abs(overlap_integrals))
 # plt.colorbar()
 # plt.show()
 
-coupling_coefficients = np.zeros((len(modes), len(modes)), dtype=complex)
-for i, (lam_i, E_i, epsilon_i) in enumerate(modes):
-    for j, (lam_j, E_j, epsilon_j) in enumerate(modes):
+coupling_coefficients = np.zeros((num_modes, num_modes), dtype=complex)
+for i, mode_i in enumerate(chain(modes_1, modes_2)):
+    for j, mode_j in enumerate(chain(modes_1, modes_2)):
         coupling_coefficients[i, j] = (
             k0
             * speed_of_light
             * epsilon_0
-            * calculate_coupling_coefficient(
-                basis0, epsilons[(epsilon_j + 1) % 2] - 1.444**2, basis, E_i, E_j
+            * mode_i.calculate_coupling_coefficient(
+                mode_j, epsilons[(j // len(modes_1) + 1) % 2] - 1.444**2
             )
             * 0.5
         )
@@ -225,8 +205,10 @@
 )
 print(kappas)
 
-delta = 0.5 * (np.real(lams_1[0]) * k0 + kappas[1, 1] - (np.real(lams_2[0]) * k0 + kappas[0, 0]))
-print(delta, np.real(lams_1[0]) * k0, kappas[1, 1])
+delta = 0.5 * (
+    np.real(modes_1[0].n_eff) * k0 + kappas[1, 1] - (np.real(modes_2[0].n_eff) * k0 + kappas[0, 0])
+)
+print(delta, np.real(modes_1[0].n_eff) * k0, kappas[1, 1])
 
 beta_c = (kappas[0, 1] * kappas[1, 0] + delta**2) ** 0.5
 
@@ -246,8 +228,11 @@
 # %%
 def fun(t, y):
     phase_matrix = [
-        [np.exp(2j * np.pi / wavelength * (lam_i - lam_j) * t) for lam_j, E_j, epsilon_j in modes]
-        for lam_i, E_i, epsilon_i in modes
+        [
+            np.exp(2j * np.pi / wavelength * (mode_i.n_eff - mode_j.n_eff) * t)
+            for mode_j in chain(modes_1, modes_2)
+        ]
+        for mode_i in chain(modes_1, modes_2)
     ]
     matrix = (
         np.linalg.inv(overlap_integrals * phase_matrix)
@@ -269,30 +254,15 @@ def fun(t, y):
 # ## two modes
 
 # %%
-R = []
-
-lam_i = lams_1[0]
-E_i = xs_1[0]
-
-for lam_j, E_j in zip(lams_both, xs_both):
-    H_i = calculate_hfield(
-        basis,
-        E_i,
-        lam_i * (2 * np.pi / 1.55),
-        omega=2 * np.pi / wavelength * speed_of_light,
-    )
-    H_j = calculate_hfield(
-        basis,
-        E_j,
-        lam_j * (2 * np.pi / 1.55),
-        omega=2 * np.pi / wavelength * speed_of_light,
-    )
-    R.append(np.abs(calculate_overlap(basis, E_i, H_i, basis, E_j, H_j) ** 2))
+R = [np.abs(modes_1[0].calculate_overlap(mode_j) ** 2) for mode_j in modes_both]
 print(R)
 P = (
     R[0] ** 2
     + R[1] ** 2
-    + 2 * R[0] * R[1] * np.cos(2 * np.pi / wavelength * (lams_both[0] - lams_both[1]) * ts)
+    + 2
+    * R[0]
+    * R[1]
+    * np.cos(2 * np.pi / wavelength * (modes_both[0].n_eff - modes_both[1].n_eff) * ts)
 )
 
 plt.plot(ts, P)

docs/photonics/examples/culomb.py

@@ -27,8 +27,8 @@
 from tqdm import tqdm
 
 from femwell.culomb import solve_coulomb
+from femwell.maxwell.waveguide import compute_modes
 from femwell.mesh import mesh_from_OrderedDict
-from femwell.mode_solver import compute_modes
 
 # -
 
@@ -108,8 +108,8 @@
     epsilon **= 2
     # basis_epsilon_r.plot(epsilon, colorbar=True).show()
 
-    neffs, basis_modes, modes = compute_modes(basis_epsilon_r, epsilon, 1.55, 1, 1, order=1)
-    voltages_neffs.append(neffs[0])
+    modes = compute_modes(basis_epsilon_r, epsilon, wavelength=1.55, order=1)
+    voltages_neffs.append(modes[0].n_eff)
 
     # from mode_solver import plot_mode
     # plot_mode(basis_modes, modes[0])

docs/photonics/examples/depletion_waveguide.py

@@ -23,8 +23,8 @@
 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.mode_solver import compute_modes
 from femwell.pn_analytical import index_pn_junction, k_to_alpha_dB
 
 # -
@@ -106,8 +106,8 @@ def define_index(mesh, V, xpn, NA, ND, wavelength):
 neff_vs_V = []
 for V in voltages:
     basis0, epsilon = define_index(mesh, V, xpn, NA, ND, wavelength)
-    lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=1, order=2)
-    neff_vs_V.append(lams)
+    modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=1, order=2)
+    neff_vs_V.append(modes[0].n_eff)
 # -
 
 # + tags=["hide-input"]

docs/photonics/examples/fiber_overlap.py

@@ -25,8 +25,8 @@
 from tqdm import tqdm
 
 from femwell.fiber import e_field_gaussian, overlap
+from femwell.maxwell.waveguide import compute_modes
 from femwell.mesh import mesh_from_OrderedDict
-from femwell.mode_solver import compute_modes, plot_mode
 
 # -
 
@@ -61,9 +61,9 @@
 # Thus, we know, that we chose the cladding thick enough if the field vanishes at the outer boundaries.
 
 # +
-lams, basis, xs = compute_modes(basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=1)
+modes = compute_modes(basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=1)
 
-fig, axs = plot_mode(basis, np.real(xs[0]), direction="x")
+fig, axs = modes[0].plot(modes[0].E.real, direction="x")
 plt.tight_layout()
 plt.show()
 # -
@@ -84,7 +84,7 @@
     )
 
     efficiency = overlap(
-        basis_fiber, basis.interpolate(xs[0])[0][1], basis_fiber.interpolate(x_fiber)
+        basis_fiber, modes[0].basis.interpolate(modes[0].E)[0][1], basis_fiber.interpolate(x_fiber)
     )
     efficiencies.append(efficiency)
 

docs/photonics/examples/leaky_waveguide.py

@@ -28,8 +28,8 @@
 from skfem import Basis, ElementDG, ElementTriP1
 from skfem.io.meshio import from_meshio
 
+from femwell.maxwell.waveguide import compute_modes
 from femwell.mesh import mesh_from_OrderedDict
-from femwell.mode_solver import compute_modes, plot_mode
 from femwell.visualization import plot_domains
 
 # %% [markdown]
@@ -105,13 +105,13 @@
 # %%
 wavelength = 1.55
 
-lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=1, order=1)
-for i, lam in enumerate(lams):
+modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=1, order=1)
+for mode in modes:
     print(
-        f"Effective refractive index: {lam:.12f}, "
-        f"Loss: {-20/np.log(10)*2*np.pi/wavelength*np.imag(lam):4f} / dB/um"
+        f"Effective refractive index: {mode.n_eff:.12f}, "
+        f"Loss: {mode.calculate_propagation_loss(distance=1):4f} / dB/um"
     )
-    plot_mode(basis, xs[i].real, colorbar=True, direction="x")
+    mode.plot(mode.E.real, colorbar=True, direction="x")
     plt.show()
 
 # %%