From 1e77defd0c7e993caabdf68201b616d6ea6c661f Mon Sep 17 00:00:00 2001
From: Eliza Leung <87836029+elizaleung830@users.noreply.github.com>
Date: Fri, 15 Dec 2023 23:42:25 +1030
Subject: [PATCH 01/10] Add files via upload

---
 docs/photonics/examples/effective_area.py | 99 +++++++++++++++++++++++
 1 file changed, 99 insertions(+)
 create mode 100644 docs/photonics/examples/effective_area.py

diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py
new file mode 100644
index 00000000..101e1782
--- /dev/null
+++ b/docs/photonics/examples/effective_area.py
@@ -0,0 +1,99 @@
+# ---
+# 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"]
+import matplotlib.pyplot as plt
+import pandas as pd
+import numpy as np
+import shapely.affinity
+from matplotlib.ticker import MultipleLocator
+
+from shapely.ops import clip_by_rect
+from femwell.mesh import mesh_from_OrderedDict
+from collections import OrderedDict
+from skfem import Basis, ElementTriP0
+from skfem.io.meshio import from_meshio
+from femwell.maxwell.waveguide import compute_modes
+
+# %% [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 = 0.5
+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(50,700,10)] # x*1e-3 for x in range(50,700,10)
+neff_list = []
+aeff_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=clip_by_rect(nc.buffer(5.0, resolution=5), -np.inf, 0, np.inf, np.inf),
+    )
+
+    resolutions = dict(core={"resolution": 0.03, "distance": 0.1})
+
+    mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=10))
+
+    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=1, order=2)
+
+    for mode in modes:
+        print(f"Effective refractive index: {mode.n_eff:.4f}")
+        #print(f"Effective mode area: {mode.calculate_effective_area():.4f}")
+        neff_list.append(mode.n_eff)
+        #aeff_list.append(mode.calculate_effective_area())
+# %% [markdown]
+# Plot the result
+reference_neff_500nm = pd.read_csv("../reference_data/Rukhlenko/fig_1c_neff/h_500nm.csv",dtype=np.float64)
+ref_x,ref_y = np.split(reference_neff_500nm.values, 2, axis=1)
+
+fig, (ax1, ax2) = plt.subplots(2,1, figsize=(9,10))
+#ax1.scatter(w_list, aeff_list, c='b',label='stimulated 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)
+
+
+ax2.scatter(w_list, neff_list, c='r',label='stimulated neff')
+ax2.scatter(ref_x,ref_y, c='b',label='reference neff')
+ax2.set_title('neff at h = 500nm')
+ax2.yaxis.set_major_locator(MultipleLocator(0.4))
+ax2.yaxis.set_minor_locator(MultipleLocator(0.2))
+ax2.set_ylim(0, 2.8)
+
+plt.legend()
+plt.show()
\ No newline at end of file

From dfe594132fc0360da4c908ee4720fea3222975e6 Mon Sep 17 00:00:00 2001
From: Eliza Leung <87836029+elizaleung830@users.noreply.github.com>
Date: Fri, 15 Dec 2023 23:44:37 +1030
Subject: [PATCH 02/10] Add files via upload

---
 .../Rukhlenko/fig_1b_aeff/0.5 Eq1.csv         | 25 +++++++++++
 .../Rukhlenko/fig_1b_aeff/0.5 Eq2.csv         | 28 ++++++++++++
 .../Rukhlenko/fig_1b_aeff/0.5 Eq3.csv         | 31 +++++++++++++
 .../Rukhlenko/fig_1b_aeff/0.7 Eq1.csv         | 19 ++++++++
 .../Rukhlenko/fig_1b_aeff/0.7 Eq2.csv         | 28 ++++++++++++
 .../Rukhlenko/fig_1b_aeff/0.7 Eq3.csv         | 18 ++++++++
 .../Rukhlenko/fig_1c_neff/h_500nm.csv         | 44 +++++++++++++++++++
 .../Rukhlenko/fig_1c_neff/h_700nm.csv         | 42 ++++++++++++++++++
 8 files changed, 235 insertions(+)
 create mode 100644 docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.5 Eq1.csv
 create mode 100644 docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.5 Eq2.csv
 create mode 100644 docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.5 Eq3.csv
 create mode 100644 docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.7 Eq1.csv
 create mode 100644 docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.7 Eq2.csv
 create mode 100644 docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.7 Eq3.csv
 create mode 100644 docs/photonics/reference_data/Rukhlenko/fig_1c_neff/h_500nm.csv
 create mode 100644 docs/photonics/reference_data/Rukhlenko/fig_1c_neff/h_700nm.csv

diff --git a/docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.5 Eq1.csv b/docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.5 Eq1.csv
new file mode 100644
index 00000000..66ee323c
--- /dev/null
+++ b/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
diff --git a/docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.5 Eq2.csv b/docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.5 Eq2.csv
new file mode 100644
index 00000000..e4b7f658
--- /dev/null
+++ b/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
diff --git a/docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.5 Eq3.csv b/docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.5 Eq3.csv
new file mode 100644
index 00000000..c43de285
--- /dev/null
+++ b/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
diff --git a/docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.7 Eq1.csv b/docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.7 Eq1.csv
new file mode 100644
index 00000000..7bc44a5d
--- /dev/null
+++ b/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
diff --git a/docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.7 Eq2.csv b/docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.7 Eq2.csv
new file mode 100644
index 00000000..96bf20a1
--- /dev/null
+++ b/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
diff --git a/docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.7 Eq3.csv b/docs/photonics/reference_data/Rukhlenko/fig_1b_aeff/0.7 Eq3.csv
new file mode 100644
index 00000000..e8563903
--- /dev/null
+++ b/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
diff --git a/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/h_500nm.csv b/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/h_500nm.csv
new file mode 100644
index 00000000..26f8cfaa
--- /dev/null
+++ b/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
diff --git a/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/h_700nm.csv b/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/h_700nm.csv
new file mode 100644
index 00000000..6d4e5d2e
--- /dev/null
+++ b/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/h_700nm.csv
@@ -0,0 +1,42 @@
+12.07, 1.047
+29.44, 1.078
+46.81, 1.112
+64.19, 1.172
+81.56, 1.248
+98.93, 1.341
+116.3, 1.437
+133.7, 1.559
+151.0, 1.625
+185.8, 1.812
+203.2, 1.890
+220.5, 1.963
+237.9, 2.027
+255.3, 2.081
+272.6, 2.135
+290.0, 2.182
+307.4, 2.215
+324.8, 2.254
+342.1, 2.296
+359.5, 2.318
+376.9, 2.345
+394.2, 2.389
+411.6, 2.398
+429.0, 2.416
+446.4, 2.435
+463.7, 2.453
+481.1, 2.469
+498.5, 2.492
+515.9, 2.499
+533.2, 2.511
+550.6, 2.522
+568.0, 2.534
+585.3, 2.546
+602.7, 2.559
+620.1, 2.565
+637.5, 2.576
+654.8, 2.585
+672.2, 2.591
+689.6, 2.602
+699.7, 2.598
+171.6, 1.769
+163.6, 1.701

From 9fafd5ee040949b97f988a5c857a4cd7f1fae6a0 Mon Sep 17 00:00:00 2001
From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com>
Date: Fri, 15 Dec 2023 13:46:49 -0800
Subject: [PATCH 03/10] implement effective_area, adjust example

---
 docs/photonics/examples/effective_area.py | 80 +++++++++++++++--------
 femwell/maxwell/waveguide.py              | 34 ++++++++++
 2 files changed, 86 insertions(+), 28 deletions(-)

diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py
index 101e1782..8a4bd123 100644
--- a/docs/photonics/examples/effective_area.py
+++ b/docs/photonics/examples/effective_area.py
@@ -15,18 +15,19 @@
 # # 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 pandas as pd
 import numpy as np
+import pandas as pd
 import shapely.affinity
 from matplotlib.ticker import MultipleLocator
-
 from shapely.ops import clip_by_rect
-from femwell.mesh import mesh_from_OrderedDict
-from collections import OrderedDict
 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.
@@ -44,12 +45,14 @@
 n_air = 1
 n_nc = 1.6
 n_silica = 1.45
-w_list = [x for x in range(50,700,10)] # x*1e-3 for x in range(50,700,10)
+w_list = [x for x in range(50, 700, 10)]  # x*1e-3 for x in range(50,700,10)
 neff_list = []
 aeff_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)
+    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)
 
@@ -57,43 +60,64 @@
         core=nc,
         silicon=silicon,
         silica=silica,
-        air=clip_by_rect(nc.buffer(5.0, resolution=5), -np.inf, 0, np.inf, np.inf),
+        air=nc.buffer(5.0, resolution=5),
     )
 
-    resolutions = dict(core={"resolution": 0.03, "distance": 0.1})
+    resolutions = dict(
+        core={"resolution": 0.03, "distance": 0.1},
+        silicon={"resolution": 0.03, "distance": 0.1},
+        silica={"resolution": 0.03, "distance": 0.1},
+    )
 
-    mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=10))
+    mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=1))
 
     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=1, order=2)
+    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=4, order=1)
 
     for mode in modes:
-        print(f"Effective refractive index: {mode.n_eff:.4f}")
-        #print(f"Effective mode area: {mode.calculate_effective_area():.4f}")
-        neff_list.append(mode.n_eff)
-        #aeff_list.append(mode.calculate_effective_area())
+        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}")
+            neff_list.append(np.real(mode.n_eff))
+            aeff_list.append(mode.calculate_effective_area())
+            break
+    else:
+        print(f"no TM mode found for {width}")
 # %% [markdown]
 # Plot the result
-reference_neff_500nm = pd.read_csv("../reference_data/Rukhlenko/fig_1c_neff/h_500nm.csv",dtype=np.float64)
-ref_x,ref_y = np.split(reference_neff_500nm.values, 2, axis=1)
+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
+)
+ref_neff_x, ref_neff_y = np.split(reference_neff_500nm.values, 2, axis=1)
+ref_aeff_x, ref_aeff_y = np.split(reference_aeff_500nm.values, 2, axis=1)
 
-fig, (ax1, ax2) = plt.subplots(2,1, figsize=(9,10))
-#ax1.scatter(w_list, aeff_list, c='b',label='stimulated aeff')
-#ax1.set_title('aeff at h = 500nm')
+fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(9, 10))
+ax1.scatter(w_list, aeff_list, c="r", label="stimulated aeff")
+ax1.scatter(ref_aeff_x, ref_aeff_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()
 
-
-ax2.scatter(w_list, neff_list, c='r',label='stimulated neff')
-ax2.scatter(ref_x,ref_y, c='b',label='reference neff')
-ax2.set_title('neff at h = 500nm')
+ax2.scatter(w_list, neff_list, c="r", label="stimulated neff")
+ax2.scatter(ref_neff_x, ref_neff_y, c="b", label="reference neff")
+ax2.set_title("neff at h = 500nm")
 ax2.yaxis.set_major_locator(MultipleLocator(0.4))
 ax2.yaxis.set_minor_locator(MultipleLocator(0.2))
 ax2.set_ylim(0, 2.8)
+ax2.legend()
 
-plt.legend()
-plt.show()
\ No newline at end of file
+plt.show()
diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py
index b3854dc3..1e8ad0d9 100644
--- a/femwell/maxwell/waveguide.py
+++ b/femwell/maxwell/waveguide.py
@@ -112,6 +112,40 @@ def overlap(w):
             delta_epsilon=self.basis_epsilon_r.interpolate(delta_epsilon),
         )
 
+    def calculate_effective_area(self, field="xy"):
+        if field == "xy":
+
+            @Functional(dtype=complex)
+            def E2(w):
+                return dot(np.conj(w["E"][0]), w["E"][0])
+
+            @Functional(dtype=complex)
+            def E4(w):
+                return dot(np.conj(w["E"][0]), w["E"][0]) ** 2
+
+        else:
+            num = 0 if (field == "x") else 1
+
+            @Functional(dtype=complex)
+            def E2(w):
+                return np.conj(w["E"][0][num]) * w["E"][0][num]
+
+            @Functional(dtype=complex)
+            def E4(w):
+                return (np.conj(w["E"][0][num]) * w["E"][0][num]) ** 2
+
+        return np.real(
+            E2.assemble(
+                self.basis,
+                E=self.basis.interpolate(self.E),
+            )
+            ** 2
+            / E4.assemble(
+                self.basis,
+                E=self.basis.interpolate(self.E),
+            )
+        )
+
     def calculate_propagation_loss(self, distance):
         return -20 / np.log(10) * self.k0 * np.imag(self.n_eff) * distance
 

From de6671a468c7b7ad5be76f5260d613631abadf35 Mon Sep 17 00:00:00 2001
From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com>
Date: Fri, 15 Dec 2023 14:03:15 -0800
Subject: [PATCH 04/10] add effective_area example to docs

---
 docs/_toc.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/_toc.yml b/docs/_toc.yml
index e096cc38..a2dde925 100644
--- a/docs/_toc.yml
+++ b/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

From 52ab1f0148a634ce5280877096c54b966aa35998 Mon Sep 17 00:00:00 2001
From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com>
Date: Fri, 15 Dec 2023 14:03:52 -0800
Subject: [PATCH 05/10] plot calculated results with a line in-between

---
 docs/photonics/examples/effective_area.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py
index 8a4bd123..16191b4f 100644
--- a/docs/photonics/examples/effective_area.py
+++ b/docs/photonics/examples/effective_area.py
@@ -104,7 +104,7 @@
 ref_aeff_x, ref_aeff_y = np.split(reference_aeff_500nm.values, 2, axis=1)
 
 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(9, 10))
-ax1.scatter(w_list, aeff_list, c="r", label="stimulated aeff")
+ax1.plot(w_list, aeff_list, "-o", c="r", label="stimulated aeff")
 ax1.scatter(ref_aeff_x, ref_aeff_y, c="b", label="reference aeff")
 ax1.set_title("aeff at h = 500nm")
 ax1.yaxis.set_major_locator(MultipleLocator(0.05))
@@ -112,7 +112,7 @@
 ax1.set_ylim(0, 0.3)
 ax1.legend()
 
-ax2.scatter(w_list, neff_list, c="r", label="stimulated neff")
+ax2.plot(w_list, neff_list, "-o", c="r", label="stimulated neff")
 ax2.scatter(ref_neff_x, ref_neff_y, c="b", label="reference neff")
 ax2.set_title("neff at h = 500nm")
 ax2.yaxis.set_major_locator(MultipleLocator(0.4))

From 3843ad71fa632f3cf5c65fe838812b71f0117334 Mon Sep 17 00:00:00 2001
From: elizaleung830 <elizaleung830@gmail.com>
Date: Sat, 16 Dec 2023 21:29:55 +1030
Subject: [PATCH 06/10] add tm_fraction in reference data

---
 .../Rukhlenko/fig_1c_neff/tm_h_500nm.csv      | 111 ++++++++++++++++++
 .../Rukhlenko/fig_1c_neff/tm_h_700nm.csv      |  85 ++++++++++++++
 2 files changed, 196 insertions(+)
 create mode 100644 docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_500nm.csv
 create mode 100644 docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_700nm.csv

diff --git a/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_500nm.csv b/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_500nm.csv
new file mode 100644
index 00000000..f09d4909
--- /dev/null
+++ b/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_500nm.csv
@@ -0,0 +1,111 @@
+6.766, 0.9465
+13.15, 0.9463
+19.54, 0.9444
+25.92, 0.9427
+32.84, 0.9405
+38.69, 0.9374
+43.91, 0.9347
+47.21, 0.9299
+52.66, 0.9279
+56.36, 0.9195
+62.64, 0.9165
+66.89, 0.9082
+70.26, 0.9042
+77.00, 0.8945
+82.86, 0.8918
+87.22, 0.8832
+93.35, 0.8805
+97.22, 0.8739
+103.6, 0.8725
+106.8, 0.8666
+113.2, 0.8661
+117.8, 0.8614
+123.9, 0.8596
+130.2, 0.8571
+143.0, 0.8519
+149.1, 0.8473
+162.1, 0.8437
+168.5, 0.8408
+174.9, 0.8388
+181.3, 0.8369
+187.7, 0.8342
+194.1, 0.8313
+200.5, 0.8301
+206.8, 0.8281
+213.2, 0.8257
+219.6, 0.8232
+226.0, 0.8223
+232.4, 0.8203
+238.8, 0.8184
+245.2, 0.8163
+251.5, 0.8156
+257.9, 0.8139
+264.3, 0.8125
+270.7, 0.8107
+277.1, 0.8102
+283.5, 0.8087
+289.8, 0.8078
+296.2, 0.8063
+302.6, 0.8059
+309.0, 0.8049
+315.4, 0.8042
+321.8, 0.8033
+328.2, 0.8026
+334.5, 0.8019
+340.9, 0.8012
+347.3, 0.8000
+353.7, 0.8000
+360.1, 0.7993
+366.5, 0.7987
+372.9, 0.7982
+379.2, 0.7978
+385.6, 0.7973
+392.0, 0.7968
+398.4, 0.7962
+404.8, 0.7961
+411.2, 0.7957
+417.6, 0.7955
+423.9, 0.7951
+430.3, 0.7948
+436.7, 0.7945
+443.1, 0.7942
+449.5, 0.7939
+455.9, 0.7937
+462.2, 0.7933
+468.6, 0.7931
+475.0, 0.7927
+481.4, 0.7927
+487.3, 0.7924
+494.2, 0.7922
+500.6, 0.7915
+506.6, 0.7918
+513.3, 0.7916
+519.7, 0.7915
+526.1, 0.7910
+532.5, 0.7910
+538.9, 0.7909
+545.3, 0.7905
+551.6, 0.7904
+558.0, 0.7904
+564.4, 0.7903
+570.8, 0.7899
+577.2, 0.7898
+583.6, 0.7898
+590.0, 0.7897
+596.3, 0.7895
+602.7, 0.7892
+609.1, 0.7892
+615.5, 0.7892
+621.9, 0.7892
+628.3, 0.7888
+634.7, 0.7887
+641.0, 0.7887
+647.4, 0.7887
+653.8, 0.7887
+660.2, 0.7884
+666.6, 0.7881
+673.0, 0.7881
+679.3, 0.7881
+685.7, 0.7881
+692.1, 0.7881
+698.5, 0.7882
diff --git a/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_700nm.csv b/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_700nm.csv
new file mode 100644
index 00000000..ccce1bb5
--- /dev/null
+++ b/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_700nm.csv
@@ -0,0 +1,85 @@
+0.7004, 0.9426
+7.831, 0.9424
+14.22, 0.9427
+57.10, 0.9272
+63.40, 0.9251
+69.55, 0.9211
+74.56, 0.9180
+77.11, 0.9138
+82.71, 0.9121
+86.58, 0.9054
+92.97, 0.9030
+97.22, 0.8956
+102.9, 0.8939
+106.8, 0.8862
+112.5, 0.8842
+116.6, 0.8764
+122.8, 0.8745
+126.9, 0.8683
+132.9, 0.8676
+136.5, 0.8609
+143.0, 0.8597
+147.6, 0.8554
+153.6, 0.8550
+160.0, 0.8517
+166.4, 0.8494
+172.8, 0.8467
+179.2, 0.8461
+185.6, 0.8445
+191.9, 0.8432
+198.6, 0.8415
+204.7, 0.8417
+242.5, 0.8387
+249.4, 0.8374
+255.8, 0.8371
+262.2, 0.8364
+268.6, 0.8360
+274.9, 0.8356
+281.3, 0.8354
+287.7, 0.8349
+294.1, 0.8343
+300.5, 0.8344
+306.9, 0.8341
+313.3, 0.8333
+319.6, 0.8331
+326.0, 0.8327
+332.4, 0.8325
+338.8, 0.8321
+345.2, 0.8317
+351.6, 0.8314
+358.0, 0.8310
+364.3, 0.8307
+370.7, 0.8304
+377.1, 0.8300
+383.5, 0.8297
+389.9, 0.8293
+396.3, 0.8289
+402.7, 0.8291
+409.0, 0.8284
+415.4, 0.8280
+421.8, 0.8279
+428.2, 0.8275
+434.6, 0.8273
+441.0, 0.8270
+447.4, 0.8267
+453.7, 0.8265
+460.1, 0.8261
+466.5, 0.8259
+472.9, 0.8256
+479.3, 0.8254
+485.7, 0.8251
+492.6, 0.8249
+498.4, 0.8244
+504.3, 0.8245
+510.1, 0.8247
+20.07, 0.9426
+26.67, 0.9398
+33.27, 0.9388
+41.52, 0.9360
+46.47, 0.9332
+53.07, 0.9308
+223.0, 0.8397
+229.6, 0.8397
+235.4, 0.8397
+208.2, 0.8415
+212.3, 0.8415

From 706de64ec3f286440ed40fb66d502cfed94177da Mon Sep 17 00:00:00 2001
From: elizaleung830 <elizaleung830@gmail.com>
Date: Sat, 16 Dec 2023 21:30:55 +1030
Subject: [PATCH 07/10] add graph for h=700nm and mode transversality

---
 docs/photonics/examples/effective_area.py | 174 +++++++++++++++-------
 1 file changed, 118 insertions(+), 56 deletions(-)

diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py
index 16191b4f..9d476602 100644
--- a/docs/photonics/examples/effective_area.py
+++ b/docs/photonics/examples/effective_area.py
@@ -39,59 +39,69 @@
 wavelength = 1.55
 capital_w = 1.4
 capital_h = 0.3
-h = 0.5
+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(50, 700, 10)]  # x*1e-3 for x in range(50,700,10)
-neff_list = []
-aeff_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(5.0, resolution=5),
-    )
-
-    resolutions = dict(
-        core={"resolution": 0.03, "distance": 0.1},
-        silicon={"resolution": 0.03, "distance": 0.1},
-        silica={"resolution": 0.03, "distance": 0.1},
-    )
-
-    mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=1))
-
-    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=4, 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}")
-            neff_list.append(np.real(mode.n_eff))
-            aeff_list.append(mode.calculate_effective_area())
-            break
-    else:
-        print(f"no TM mode found for {width}")
+w_list = [x for x in range(50, 700, 150)]  # x*1e-3 for x in range(50,700,10)
+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(5.0, resolution=5),
+        )
+
+        resolutions = dict(
+            core={"resolution": 0.03, "distance": 0.1},
+            silicon={"resolution": 0.03, "distance": 0.1},
+            silica={"resolution": 0.03, "distance": 0.1},
+        )
+
+        mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=1))
+
+        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=4, 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"tm fraction: {mode.tm_fraction}")
+                neff_list.append(np.real(mode.n_eff))
+                aeff_list.append(mode.calculate_effective_area())
+                tm_list.append(mode.tm_fraction)
+                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(
@@ -100,24 +110,76 @@
 reference_aeff_500nm = pd.read_csv(
     "../reference_data/Rukhlenko/fig_1b_aeff/0.5 Eq2.csv", dtype=np.float64
 )
-ref_neff_x, ref_neff_y = np.split(reference_neff_500nm.values, 2, axis=1)
-ref_aeff_x, ref_aeff_y = np.split(reference_aeff_500nm.values, 2, axis=1)
+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
+)
 
-fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(9, 10))
-ax1.plot(w_list, aeff_list, "-o", c="r", label="stimulated aeff")
-ax1.scatter(ref_aeff_x, ref_aeff_y, c="b", label="reference aeff")
+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()
 
-ax2.plot(w_list, neff_list, "-o", c="r", label="stimulated neff")
-ax2.scatter(ref_neff_x, ref_neff_y, c="b", label="reference neff")
-ax2.set_title("neff at h = 500nm")
+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()

From 74485fb5f3a0f9fa713e2818add14f033d1374b4 Mon Sep 17 00:00:00 2001
From: elizaleung830 <elizaleung830@gmail.com>
Date: Sat, 16 Dec 2023 21:32:02 +1030
Subject: [PATCH 08/10] fixing wrong range

---
 docs/photonics/examples/effective_area.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py
index 9d476602..f63c2e9b 100644
--- a/docs/photonics/examples/effective_area.py
+++ b/docs/photonics/examples/effective_area.py
@@ -45,7 +45,7 @@
 n_air = 1
 n_nc = 1.6
 n_silica = 1.45
-w_list = [x for x in range(50, 700, 150)]  # x*1e-3 for x in range(50,700,10)
+w_list = [x for x in range(50, 700, 10)]  # x*1e-3 for x in range(50,700,10)
 neff_dict = dict()
 aeff_dict = dict()
 tm_dict = dict()

From 7b355eb4b49c967db44af56b631c1b77ca378177 Mon Sep 17 00:00:00 2001
From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com>
Date: Mon, 18 Dec 2023 21:12:03 -0800
Subject: [PATCH 09/10] exchange datafiles

---
 .../Rukhlenko/fig_1c_neff/tm_h_500nm.csv      | 196 ++++++++----------
 .../Rukhlenko/fig_1c_neff/tm_h_700nm.csv      | 196 ++++++++++--------
 2 files changed, 196 insertions(+), 196 deletions(-)

diff --git a/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_500nm.csv b/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_500nm.csv
index f09d4909..ccce1bb5 100644
--- a/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_500nm.csv
+++ b/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_500nm.csv
@@ -1,111 +1,85 @@
-6.766, 0.9465
-13.15, 0.9463
-19.54, 0.9444
-25.92, 0.9427
-32.84, 0.9405
-38.69, 0.9374
-43.91, 0.9347
-47.21, 0.9299
-52.66, 0.9279
-56.36, 0.9195
-62.64, 0.9165
-66.89, 0.9082
-70.26, 0.9042
-77.00, 0.8945
-82.86, 0.8918
-87.22, 0.8832
-93.35, 0.8805
-97.22, 0.8739
-103.6, 0.8725
-106.8, 0.8666
-113.2, 0.8661
-117.8, 0.8614
-123.9, 0.8596
-130.2, 0.8571
-143.0, 0.8519
-149.1, 0.8473
-162.1, 0.8437
-168.5, 0.8408
-174.9, 0.8388
-181.3, 0.8369
-187.7, 0.8342
-194.1, 0.8313
-200.5, 0.8301
-206.8, 0.8281
-213.2, 0.8257
-219.6, 0.8232
-226.0, 0.8223
-232.4, 0.8203
-238.8, 0.8184
-245.2, 0.8163
-251.5, 0.8156
-257.9, 0.8139
-264.3, 0.8125
-270.7, 0.8107
-277.1, 0.8102
-283.5, 0.8087
-289.8, 0.8078
-296.2, 0.8063
-302.6, 0.8059
-309.0, 0.8049
-315.4, 0.8042
-321.8, 0.8033
-328.2, 0.8026
-334.5, 0.8019
-340.9, 0.8012
-347.3, 0.8000
-353.7, 0.8000
-360.1, 0.7993
-366.5, 0.7987
-372.9, 0.7982
-379.2, 0.7978
-385.6, 0.7973
-392.0, 0.7968
-398.4, 0.7962
-404.8, 0.7961
-411.2, 0.7957
-417.6, 0.7955
-423.9, 0.7951
-430.3, 0.7948
-436.7, 0.7945
-443.1, 0.7942
-449.5, 0.7939
-455.9, 0.7937
-462.2, 0.7933
-468.6, 0.7931
-475.0, 0.7927
-481.4, 0.7927
-487.3, 0.7924
-494.2, 0.7922
-500.6, 0.7915
-506.6, 0.7918
-513.3, 0.7916
-519.7, 0.7915
-526.1, 0.7910
-532.5, 0.7910
-538.9, 0.7909
-545.3, 0.7905
-551.6, 0.7904
-558.0, 0.7904
-564.4, 0.7903
-570.8, 0.7899
-577.2, 0.7898
-583.6, 0.7898
-590.0, 0.7897
-596.3, 0.7895
-602.7, 0.7892
-609.1, 0.7892
-615.5, 0.7892
-621.9, 0.7892
-628.3, 0.7888
-634.7, 0.7887
-641.0, 0.7887
-647.4, 0.7887
-653.8, 0.7887
-660.2, 0.7884
-666.6, 0.7881
-673.0, 0.7881
-679.3, 0.7881
-685.7, 0.7881
-692.1, 0.7881
-698.5, 0.7882
+0.7004, 0.9426
+7.831, 0.9424
+14.22, 0.9427
+57.10, 0.9272
+63.40, 0.9251
+69.55, 0.9211
+74.56, 0.9180
+77.11, 0.9138
+82.71, 0.9121
+86.58, 0.9054
+92.97, 0.9030
+97.22, 0.8956
+102.9, 0.8939
+106.8, 0.8862
+112.5, 0.8842
+116.6, 0.8764
+122.8, 0.8745
+126.9, 0.8683
+132.9, 0.8676
+136.5, 0.8609
+143.0, 0.8597
+147.6, 0.8554
+153.6, 0.8550
+160.0, 0.8517
+166.4, 0.8494
+172.8, 0.8467
+179.2, 0.8461
+185.6, 0.8445
+191.9, 0.8432
+198.6, 0.8415
+204.7, 0.8417
+242.5, 0.8387
+249.4, 0.8374
+255.8, 0.8371
+262.2, 0.8364
+268.6, 0.8360
+274.9, 0.8356
+281.3, 0.8354
+287.7, 0.8349
+294.1, 0.8343
+300.5, 0.8344
+306.9, 0.8341
+313.3, 0.8333
+319.6, 0.8331
+326.0, 0.8327
+332.4, 0.8325
+338.8, 0.8321
+345.2, 0.8317
+351.6, 0.8314
+358.0, 0.8310
+364.3, 0.8307
+370.7, 0.8304
+377.1, 0.8300
+383.5, 0.8297
+389.9, 0.8293
+396.3, 0.8289
+402.7, 0.8291
+409.0, 0.8284
+415.4, 0.8280
+421.8, 0.8279
+428.2, 0.8275
+434.6, 0.8273
+441.0, 0.8270
+447.4, 0.8267
+453.7, 0.8265
+460.1, 0.8261
+466.5, 0.8259
+472.9, 0.8256
+479.3, 0.8254
+485.7, 0.8251
+492.6, 0.8249
+498.4, 0.8244
+504.3, 0.8245
+510.1, 0.8247
+20.07, 0.9426
+26.67, 0.9398
+33.27, 0.9388
+41.52, 0.9360
+46.47, 0.9332
+53.07, 0.9308
+223.0, 0.8397
+229.6, 0.8397
+235.4, 0.8397
+208.2, 0.8415
+212.3, 0.8415
diff --git a/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_700nm.csv b/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_700nm.csv
index ccce1bb5..f09d4909 100644
--- a/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_700nm.csv
+++ b/docs/photonics/reference_data/Rukhlenko/fig_1c_neff/tm_h_700nm.csv
@@ -1,85 +1,111 @@
-0.7004, 0.9426
-7.831, 0.9424
-14.22, 0.9427
-57.10, 0.9272
-63.40, 0.9251
-69.55, 0.9211
-74.56, 0.9180
-77.11, 0.9138
-82.71, 0.9121
-86.58, 0.9054
-92.97, 0.9030
-97.22, 0.8956
-102.9, 0.8939
-106.8, 0.8862
-112.5, 0.8842
-116.6, 0.8764
-122.8, 0.8745
-126.9, 0.8683
-132.9, 0.8676
-136.5, 0.8609
-143.0, 0.8597
-147.6, 0.8554
-153.6, 0.8550
-160.0, 0.8517
-166.4, 0.8494
-172.8, 0.8467
-179.2, 0.8461
-185.6, 0.8445
-191.9, 0.8432
-198.6, 0.8415
-204.7, 0.8417
-242.5, 0.8387
-249.4, 0.8374
-255.8, 0.8371
-262.2, 0.8364
-268.6, 0.8360
-274.9, 0.8356
-281.3, 0.8354
-287.7, 0.8349
-294.1, 0.8343
-300.5, 0.8344
-306.9, 0.8341
-313.3, 0.8333
-319.6, 0.8331
-326.0, 0.8327
-332.4, 0.8325
-338.8, 0.8321
-345.2, 0.8317
-351.6, 0.8314
-358.0, 0.8310
-364.3, 0.8307
-370.7, 0.8304
-377.1, 0.8300
-383.5, 0.8297
-389.9, 0.8293
-396.3, 0.8289
-402.7, 0.8291
-409.0, 0.8284
-415.4, 0.8280
-421.8, 0.8279
-428.2, 0.8275
-434.6, 0.8273
-441.0, 0.8270
-447.4, 0.8267
-453.7, 0.8265
-460.1, 0.8261
-466.5, 0.8259
-472.9, 0.8256
-479.3, 0.8254
-485.7, 0.8251
-492.6, 0.8249
-498.4, 0.8244
-504.3, 0.8245
-510.1, 0.8247
-20.07, 0.9426
-26.67, 0.9398
-33.27, 0.9388
-41.52, 0.9360
-46.47, 0.9332
-53.07, 0.9308
-223.0, 0.8397
-229.6, 0.8397
-235.4, 0.8397
-208.2, 0.8415
-212.3, 0.8415
+6.766, 0.9465
+13.15, 0.9463
+19.54, 0.9444
+25.92, 0.9427
+32.84, 0.9405
+38.69, 0.9374
+43.91, 0.9347
+47.21, 0.9299
+52.66, 0.9279
+56.36, 0.9195
+62.64, 0.9165
+66.89, 0.9082
+70.26, 0.9042
+77.00, 0.8945
+82.86, 0.8918
+87.22, 0.8832
+93.35, 0.8805
+97.22, 0.8739
+103.6, 0.8725
+106.8, 0.8666
+113.2, 0.8661
+117.8, 0.8614
+123.9, 0.8596
+130.2, 0.8571
+143.0, 0.8519
+149.1, 0.8473
+162.1, 0.8437
+168.5, 0.8408
+174.9, 0.8388
+181.3, 0.8369
+187.7, 0.8342
+194.1, 0.8313
+200.5, 0.8301
+206.8, 0.8281
+213.2, 0.8257
+219.6, 0.8232
+226.0, 0.8223
+232.4, 0.8203
+238.8, 0.8184
+245.2, 0.8163
+251.5, 0.8156
+257.9, 0.8139
+264.3, 0.8125
+270.7, 0.8107
+277.1, 0.8102
+283.5, 0.8087
+289.8, 0.8078
+296.2, 0.8063
+302.6, 0.8059
+309.0, 0.8049
+315.4, 0.8042
+321.8, 0.8033
+328.2, 0.8026
+334.5, 0.8019
+340.9, 0.8012
+347.3, 0.8000
+353.7, 0.8000
+360.1, 0.7993
+366.5, 0.7987
+372.9, 0.7982
+379.2, 0.7978
+385.6, 0.7973
+392.0, 0.7968
+398.4, 0.7962
+404.8, 0.7961
+411.2, 0.7957
+417.6, 0.7955
+423.9, 0.7951
+430.3, 0.7948
+436.7, 0.7945
+443.1, 0.7942
+449.5, 0.7939
+455.9, 0.7937
+462.2, 0.7933
+468.6, 0.7931
+475.0, 0.7927
+481.4, 0.7927
+487.3, 0.7924
+494.2, 0.7922
+500.6, 0.7915
+506.6, 0.7918
+513.3, 0.7916
+519.7, 0.7915
+526.1, 0.7910
+532.5, 0.7910
+538.9, 0.7909
+545.3, 0.7905
+551.6, 0.7904
+558.0, 0.7904
+564.4, 0.7903
+570.8, 0.7899
+577.2, 0.7898
+583.6, 0.7898
+590.0, 0.7897
+596.3, 0.7895
+602.7, 0.7892
+609.1, 0.7892
+615.5, 0.7892
+621.9, 0.7892
+628.3, 0.7888
+634.7, 0.7887
+641.0, 0.7887
+647.4, 0.7887
+653.8, 0.7887
+660.2, 0.7884
+666.6, 0.7881
+673.0, 0.7881
+679.3, 0.7881
+685.7, 0.7881
+692.1, 0.7881
+698.5, 0.7882

From 0395f4cc67cff2888d32aecc6920a99326f81e8d Mon Sep 17 00:00:00 2001
From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com>
Date: Mon, 18 Dec 2023 21:12:47 -0800
Subject: [PATCH 10/10] adjust definition of transversality

---
 docs/photonics/examples/effective_area.py | 44 +++++++++++------------
 femwell/maxwell/waveguide.py              | 26 ++++++++++++++
 2 files changed, 48 insertions(+), 22 deletions(-)

diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py
index f63c2e9b..386ce167 100644
--- a/docs/photonics/examples/effective_area.py
+++ b/docs/photonics/examples/effective_area.py
@@ -39,13 +39,13 @@
 wavelength = 1.55
 capital_w = 1.4
 capital_h = 0.3
-h_list = [0.5,0.7]
+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(50, 700, 10)]  # x*1e-3 for x in range(50,700,10)
+w_list = [x for x in range(250, 700, 100)]
 neff_dict = dict()
 aeff_dict = dict()
 tm_dict = dict()
@@ -65,16 +65,16 @@
             core=nc,
             silicon=silicon,
             silica=silica,
-            air=nc.buffer(5.0, resolution=5),
+            air=nc.buffer(10.0, resolution=4),
         )
 
         resolutions = dict(
-            core={"resolution": 0.03, "distance": 0.1},
-            silicon={"resolution": 0.03, "distance": 0.1},
-            silica={"resolution": 0.03, "distance": 0.1},
+            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=1))
+        mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=2))
 
         basis0 = Basis(mesh, ElementTriP0())
         epsilon = basis0.zeros()
@@ -84,22 +84,22 @@
             "air": n_air,
             "silica": n_silica,
         }.items():
-            epsilon[basis0.get_dofs(elements=subdomain)] = n ** 2
-        modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=4, order=1)
+            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"tm fraction: {mode.tm_fraction}")
+                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.tm_fraction)
+                tm_list.append(mode.transversality)
                 break
         else:
             print(f"no TM mode found for {width}")
-    neff_dict[str(h)]=neff_list
+    neff_dict[str(h)] = neff_list
     aeff_dict[str(h)] = aeff_list
     tm_dict[str(h)] = tm_list
 # %% [markdown]
@@ -134,7 +134,7 @@
 
 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.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))
@@ -144,12 +144,12 @@
 
 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.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.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")
@@ -158,7 +158,7 @@
 ax2.set_ylim(0, 2.8)
 ax2.legend()
 
-ax3.plot(w_list, aeff_dict['0.7'], "-o", c="r", label="stimulated aeff")
+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))
@@ -168,13 +168,13 @@
 
 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.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.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))
diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py
index 1e8ad0d9..68bffc52 100644
--- a/femwell/maxwell/waveguide.py
+++ b/femwell/maxwell/waveguide.py
@@ -92,6 +92,32 @@ def tm_fraction(self):
 
         return 1 - self.te_fraction
 
+    @cached_property
+    def transversality(self):
+        """TE-fraction of the mode"""
+
+        @Functional
+        def ex(w):
+            return np.abs(w.E[0][0]) ** 2
+
+        @Functional
+        def ey(w):
+            return np.abs(w.E[0][1]) ** 2
+
+        @Functional
+        def ez(w):
+            return np.abs(
+                w.E[0][0] * np.conj(w.E[0][0])
+                + w.E[0][1] * np.conj(w.E[0][1])
+                + w.E[1] * np.conj(w.E[1])
+            )
+
+        ex_sum = ex.assemble(self.basis, E=self.basis.interpolate(self.E))
+        ey_sum = ey.assemble(self.basis, E=self.basis.interpolate(self.E))
+        ez_sum = ez.assemble(self.basis, E=self.basis.interpolate(self.E))
+
+        return (ex_sum + ey_sum) / (ez_sum)
+
     def __repr__(self) -> str:
         return f"{self.__class__.__name__}(k: {self.k}, n_eff:{self.n_eff})"