Implement effective_area + docs example

docs/_toc.yml

@@ -37,6 +37,7 @@ parts:
     - file: photonics/examples/grating_coupler.py
     - file: photonics/examples/coupled_mode_theory.py
     - file: photonics/examples/depletion_waveguide.py
+    - file: photonics/examples/effective_area.py
     - file: photonics/examples/refinement.py
     - file: photonics/examples/propagation_loss.py
     - file: electronics/examples/capacitor.py

docs/photonics/examples/effective_area.py

@@ -0,0 +1,185 @@
+# ---
+# 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]
+# # Effective area and effective refractive index of a Si-NCs waveguide from
+# https://opg.optica.org/ol/abstract.cfm?uri=ol-37-12-2295#:~:text=A%20simple%20and%20physically%20meaningful,of%20Si%2DNC%20slot%20waveguides.
+# %% tags=["hide-input"]
+from collections import OrderedDict
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import shapely.affinity
+from matplotlib.ticker import MultipleLocator
+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
+
+# %% [markdown]
+# We describe the geometry using shapely.
+# In this case it's simple: we use a shapely.box for the waveguide according to the parameter provided in the paper
+# For the surrounding we buffer the core and clip it to the part above the waveguide for the air
+# Width is a sweep parameter from 50nm to 700nm
+
+# %%
+wavelength = 1.55
+capital_w = 1.4
+capital_h = 0.3
+h_list = [0.5, 0.7]
+t = 0.1
+n_silicon = 3.48
+n_air = 1
+n_nc = 1.6
+n_silica = 1.45
+w_list = [x for x in range(250, 700, 100)]
+neff_dict = dict()
+aeff_dict = dict()
+tm_dict = dict()
+for h in h_list:
+    neff_list = []
+    aeff_list = []
+    tm_list = []
+    for width in w_list:
+        width = width * 1e-3
+        nc = shapely.geometry.box(
+            -width / 2, capital_h + (h - t) / 2, +width / 2, capital_h + (h - t) / 2 + t
+        )
+        silica = shapely.geometry.box(-capital_w / 2, 0, +capital_w / 2, capital_h)
+        silicon = shapely.geometry.box(-width / 2, capital_h, +width / 2, capital_h + h)
+
+        polygons = OrderedDict(
+            core=nc,
+            silicon=silicon,
+            silica=silica,
+            air=nc.buffer(10.0, resolution=4),
+        )
+
+        resolutions = dict(
+            core={"resolution": 0.02, "distance": 0.3},
+            silicon={"resolution": 0.02, "distance": 0.1},
+            silica={"resolution": 0.04, "distance": 0.5},
+        )
+
+        mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=2))
+
+        basis0 = Basis(mesh, ElementTriP0())
+        epsilon = basis0.zeros()
+        for subdomain, n in {
+            "core": n_nc,
+            "silicon": n_silicon,
+            "air": n_air,
+            "silica": n_silica,
+        }.items():
+            epsilon[basis0.get_dofs(elements=subdomain)] = n**2
+        modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=3, order=1)
+
+        for mode in modes:
+            if mode.tm_fraction > 0.5:
+                # mode.show(np.real(mode.E))
+                print(f"Effective refractive index: {mode.n_eff:.4f}")
+                print(f"Effective mode area: {mode.calculate_effective_area(field='y'):.4f}")
+                print(f"Mode transversality: {mode.transversality}")
+                neff_list.append(np.real(mode.n_eff))
+                aeff_list.append(mode.calculate_effective_area())
+                tm_list.append(mode.transversality)
+                break
+        else:
+            print(f"no TM mode found for {width}")
+    neff_dict[str(h)] = neff_list
+    aeff_dict[str(h)] = aeff_list
+    tm_dict[str(h)] = tm_list
+# %% [markdown]
+# Plot the result
+reference_neff_500nm = pd.read_csv(
+    "../reference_data/Rukhlenko/fig_1c_neff/h_500nm.csv", dtype=np.float64
+)
+reference_aeff_500nm = pd.read_csv(
+    "../reference_data/Rukhlenko/fig_1b_aeff/0.5 Eq2.csv", dtype=np.float64
+)
+reference_tm_500nm = pd.read_csv(
+    "../reference_data/Rukhlenko/fig_1c_neff/tm_h_500nm.csv", dtype=np.float64
+)
+
+reference_neff_700nm = pd.read_csv(
+    "../reference_data/Rukhlenko/fig_1c_neff/h_700nm.csv", dtype=np.float64
+)
+reference_aeff_700nm = pd.read_csv(
+    "../reference_data/Rukhlenko/fig_1b_aeff/0.7 Eq2.csv", dtype=np.float64
+)
+reference_tm_700nm = pd.read_csv(
+    "../reference_data/Rukhlenko/fig_1c_neff/tm_h_700nm.csv", dtype=np.float64
+)
+
+ref_neff_500nm_x, ref_neff_500nm_y = np.split(reference_neff_500nm.values, 2, axis=1)
+ref_aeff_500nm_x, ref_aeff_500nm_y = np.split(reference_aeff_500nm.values, 2, axis=1)
+ref_tm_500nm_x, ref_tm_500nm_y = np.split(reference_tm_500nm.values, 2, axis=1)
+
+ref_neff_700nm_x, ref_neff_700nm_y = np.split(reference_neff_700nm.values, 2, axis=1)
+ref_aeff_700nm_x, ref_aeff_700nm_y = np.split(reference_aeff_700nm.values, 2, axis=1)
+ref_tm_700nm_x, ref_tm_700nm_y = np.split(reference_tm_700nm.values, 2, axis=1)
+
+fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(9, 20))
+
+ax1.plot(w_list, aeff_dict["0.5"], "-o", c="r", label="stimulated aeff")
+ax1.scatter(ref_aeff_500nm_x, ref_aeff_500nm_y, c="b", label="reference aeff")
+ax1.set_title("aeff at h = 500nm")
+ax1.yaxis.set_major_locator(MultipleLocator(0.05))
+ax1.yaxis.set_minor_locator(MultipleLocator(0.01))
+ax1.set_ylim(0, 0.3)
+ax1.legend()
+
+ax2b = ax2.twinx()
+ax2b.set_ylabel("mode transversality")
+ax2b.scatter(ref_tm_500nm_x, ref_tm_500nm_y, marker="v", c="b", label="reference transversality")
+ax2b.plot(w_list, tm_dict["0.5"], "-v", c="r", label="stimulated transversality")
+ax2b.set_ylim(0.775, 1)
+ax2b.legend()
+
+ax2.plot(w_list, neff_dict["0.5"], "-o", c="r", label="stimulated neff")
+ax2.scatter(ref_neff_500nm_x, ref_neff_500nm_y, c="b", label="reference neff")
+ax2.set_title("neff and mode transversality at h = 500nm")
+ax2.set_ylabel("neff")
+ax2.yaxis.set_major_locator(MultipleLocator(0.4))
+ax2.yaxis.set_minor_locator(MultipleLocator(0.2))
+ax2.set_ylim(0, 2.8)
+ax2.legend()
+
+ax3.plot(w_list, aeff_dict["0.7"], "-o", c="r", label="stimulated aeff")
+ax3.scatter(ref_aeff_700nm_x, ref_aeff_700nm_y, c="b", label="reference aeff")
+ax3.set_title("aeff at h = 700nm")
+ax3.yaxis.set_major_locator(MultipleLocator(0.05))
+ax3.yaxis.set_minor_locator(MultipleLocator(0.01))
+ax3.set_ylim(0, 0.3)
+ax3.legend()
+
+ax4b = ax4.twinx()
+ax4b.set_ylabel("mode transversality")
+ax4b.scatter(ref_tm_700nm_x, ref_tm_700nm_y, marker="v", c="b", label="reference transversality")
+ax4b.plot(w_list, tm_dict["0.7"], "-v", c="r", label="stimulated transversality")
+ax4b.set_ylim(0.775, 1)
+ax4b.legend()
+
+
+ax4.plot(w_list, neff_dict["0.7"], "-o", c="r", label="stimulated aeff")
+ax4.scatter(ref_neff_700nm_x, ref_neff_700nm_y, c="b", label="reference aeff")
+ax4.set_title("neff and mode transversality at h = 700nm")
+ax2.yaxis.set_major_locator(MultipleLocator(0.4))
+ax2.yaxis.set_minor_locator(MultipleLocator(0.2))
+ax4.set_ylim(0, 2.8)
+ax4.legend()
+
+plt.show()

docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.5 Eq1.csv

@@ -0,0 +1,25 @@
+141.1, 0.2904
+146.7, 0.2706
+150.2, 0.2518
+156.8, 0.2360
+164.0, 0.2224
+171.2, 0.2034
+177.0, 0.1903
+192.9, 0.1728
+208.8, 0.1618
+234.8, 0.1499
+266.9, 0.1448
+298.4, 0.1444
+330.2, 0.1468
+362.0, 0.1501
+393.8, 0.1550
+425.6, 0.1607
+457.4, 0.1664
+489.2, 0.1724
+521.0, 0.1784
+552.8, 0.1859
+584.6, 0.1928
+616.4, 0.1998
+648.2, 0.2075
+680.0, 0.2149
+698.8, 0.2195

docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.5 Eq2.csv

@@ -0,0 +1,28 @@
+120.9, 0.2868
+125.0, 0.2552
+130.8, 0.2277
+127.9, 0.2425
+135.1, 0.2087
+140.4, 0.1858
+148.1, 0.1657
+152.2, 0.1532
+162.6, 0.1392
+172.7, 0.1232
+177.6, 0.1145
+246.4, 0.08712
+278.2, 0.08491
+310.0, 0.08515
+341.8, 0.08789
+382.7, 0.08975
+411.2, 0.09183
+443.0, 0.09622
+474.8, 0.09965
+506.6, 0.1032
+538.4, 0.1071
+570.2, 0.1108
+602.0, 0.1149
+633.8, 0.1189
+665.6, 0.1229
+691.6, 0.1266
+224.8, 0.09111
+201.7, 0.09971

docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.5 Eq3.csv

@@ -0,0 +1,31 @@
+74.40, 0.2953
+80.18, 0.2480
+76.81, 0.2747
+89.34, 0.1981
+85.48, 0.2240
+99.94, 0.1615
+96.08, 0.1787
+107.6, 0.1357
+123.1, 0.1087
+140.7, 0.08927
+165.5, 0.07725
+197.3, 0.07098
+218.7, 0.06869
+246.4, 0.07154
+278.2, 0.07468
+309.3, 0.07899
+342.5, 0.08385
+373.6, 0.08967
+404.9, 0.09503
+437.2, 0.1019
+469.0, 0.1082
+500.8, 0.1142
+532.6, 0.1215
+564.4, 0.1277
+596.7, 0.1333
+628.0, 0.1420
+659.8, 0.1492
+688.7, 0.1557
+176.9, 0.07332
+150.5, 0.08264
+115.8, 0.1213

docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.7 Eq1.csv

@@ -0,0 +1,19 @@
+110.5, 0.2993
+119.9, 0.2753
+130.8, 0.2531
+143.8, 0.2339
+161.3, 0.2191
+191.5, 0.2067
+223.3, 0.2013
+255.1, 0.2016
+286.9, 0.2052
+318.7, 0.2114
+350.5, 0.2184
+382.3, 0.2268
+414.1, 0.2352
+445.9, 0.2454
+477.7, 0.2552
+509.5, 0.2652
+541.3, 0.2759
+570.2, 0.2857
+596.2, 0.2966

docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.7 Eq2.csv

@@ -0,0 +1,28 @@
+85.97, 0.2939
+91.95, 0.2394
+88.38, 0.2718
+90.30, 0.2538
+101.1, 0.1846
+96.08, 0.2160
+98.97, 0.2050
+110.5, 0.1520
+107.6, 0.1675
+122.5, 0.1293
+145.2, 0.1108
+174.1, 0.1008
+205.9, 0.1014
+237.7, 0.1067
+269.5, 0.1126
+301.3, 0.1207
+333.1, 0.1294
+364.9, 0.1378
+396.7, 0.1467
+428.5, 0.1555
+460.3, 0.1642
+492.1, 0.1730
+523.9, 0.1815
+555.7, 0.1903
+587.5, 0.1991
+619.3, 0.2073
+651.1, 0.2160
+682.9, 0.2246

docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.7 Eq3.csv

@@ -0,0 +1,18 @@
+61.39, 0.2521
+72.94, 0.1885
+82.84, 0.1492
+91.09, 0.1242
+102.6, 0.1113
+109.2, 0.1027
+119.1, 0.09769
+142.2, 0.09694
+152.1, 0.09836
+178.5, 0.1040
+201.7, 0.1133
+224.8, 0.1247
+251.2, 0.1368
+277.6, 0.1517
+300.7, 0.1653
+350.2, 0.1938
+399.7, 0.2237
+500.3, 0.2850

docs/photonics/reference_data/Rukhlenko/fig_1c_neff/h_500nm.csv

@@ -0,0 +1,44 @@
+5.245, 1.010
+23.65, 1.004
+41.02, 1.029
+58.39, 1.070
+93.14, 1.162
+110.5, 1.204
+127.9, 1.256
+145.3, 1.316
+162.6, 1.371
+178.5, 1.422
+214.7, 1.549
+232.1, 1.573
+249.5, 1.615
+266.9, 1.670
+283.9, 1.702
+301.6, 1.746
+318.1, 1.774
+336.3, 1.809
+352.8, 1.832
+369.6, 1.856
+388.5, 1.883
+405.0, 1.900
+422.3, 1.922
+440.3, 1.939
+458.2, 1.956
+476.2, 1.973
+493.6, 1.991
+510.1, 2.003
+528.3, 2.016
+545.1, 2.027
+561.9, 2.038
+578.1, 2.047
+596.9, 2.067
+615.2, 2.073
+628.8, 2.075
+649.0, 2.090
+666.1, 2.095
+683.8, 2.100
+696.8, 2.100
+201.4, 1.499
+187.7, 1.454
+275.9, 1.696
+79.94, 1.138
+68.48, 1.104