diff --git a/femwell/examples/Aeff.py b/femwell/examples/Aeff.py index 9d343aa1..d9c8096c 100644 --- a/femwell/examples/Aeff.py +++ b/femwell/examples/Aeff.py @@ -1,23 +1,19 @@ """Aeff analysis based on https://optics.ansys.com/hc/en-us/articles/15100783091731-Spontaneous-Four-wave-Mixing-SWFM-Microring-Resonator-Photon-Source.""" from collections import OrderedDict + import numpy as np +from scipy.constants import c, epsilon_0 from shapely.geometry import box from shapely.ops import clip_by_rect -from skfem import Basis, ElementTriP0 -from skfem import Functional +from skfem import Basis, ElementTriP0, Functional from skfem.helpers import dot from skfem.io.meshio import from_meshio + from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from scipy.constants import c, epsilon_0 -def calculate_sfwm_Aeff( - basis: Basis, - mode_p, - mode_s, - mode_i -) -> np.complex64: - + +def calculate_sfwm_Aeff(basis: Basis, mode_p, mode_s, mode_i) -> np.complex64: """ Calculates the overlap integral for SFWM process by considering the interactions between pump, signal, and idler modes in the xy plane. @@ -30,46 +26,43 @@ def calculate_sfwm_Aeff( np.complex64: The Aeff result for the SFWM process(1/overlap integral). """ - def normalize_mode(mode): + def normalization_factor_mode(mode): @Functional def E2(w): return dot(w["E"][0], np.conj(w["E"][0])) - E_xy, _ = mode.basis.interpolate(mode.E) # [0]=xy [1]=z - E_squared_integral = E2.assemble(mode.basis, E=E_xy) + E = mode.basis.interpolate(mode.E) # [0]=xy [1]=z + E_squared_integral = E2.assemble(mode.basis, E=E) normalization_factor = 1 / np.sqrt(E_squared_integral) return normalization_factor # Return the normalization factor instead of modifying the mode - # Calculate normalization factors for each mode - norm_factor_p = normalize_mode(mode_p) - norm_factor_s = normalize_mode(mode_s) - norm_factor_i = normalize_mode(mode_i) - # Apply normalization factors to the electric fields for the overlap calculation - E_p, _ = mode_p.basis.interpolate(mode_p.E * norm_factor_p) - E_s, _ = mode_s.basis.interpolate(mode_s.E * norm_factor_s) - E_i, _ = mode_i.basis.interpolate(mode_i.E * norm_factor_i) - @Functional(dtype=np.complex64) def sfwm_overlap(w): - overlap_SFWM = w["E_p"][0] * w["E_p"][0] * np.conj(w["E_s"][0]) * np.conj(w["E_i"][0]) - return overlap_SFWM + return dot(w["E_p"][0], w["E_p"][0]) * dot(np.conj(w["E_s"][0]), np.conj(w["E_i"][0])) - overlap_result = sfwm_overlap.assemble(basis, E_p=E_p, E_s=E_s, E_i=E_i) + overlap_result = sfwm_overlap.assemble( + basis, + E_p=mode_p.basis.interpolate(mode_p.E * normalization_factor_mode(mode_p)), + E_s=mode_s.basis.interpolate(mode_s.E * normalization_factor_mode(mode_s)), + E_i=mode_i.basis.interpolate(mode_i.E * normalization_factor_mode(mode_i)), + ) return 1 / overlap_result -#Dispersion relations of materials -#Core +# Dispersion relations of materials +# Core def n_X(wavelength): x = wavelength - return (1 - + 2.19244563 / (1 - (0.20746607 / x) ** 2) - + 0.13435116 / (1 - (0.3985835 / x) ** 2) - + 2.20997784 / (1 - (0.20747044 / x) ** 2) + return ( + 1 + + 2.19244563 / (1 - (0.20746607 / x) ** 2) + + 0.13435116 / (1 - (0.3985835 / x) ** 2) + + 2.20997784 / (1 - (0.20747044 / x) ** 2) ) ** 0.5 + # Box def n_silicon_dioxide(wavelength): x = wavelength @@ -80,27 +73,33 @@ def n_silicon_dioxide(wavelength): + 0.8974794 / (1 - (9.896161 / x) ** 2) ) ** 0.5 -Clad=1 + +Clad = 1 -core = box(0, 0, 0.5, 0.39) #500x390nm +core = box(0, 0, 1.05, 0.5) # 1050x500nm +# core = box(0, 0, .5, 0.39) # 500x390nm polygons = OrderedDict( core=core, box=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, -np.inf, np.inf, 0), clad=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, 0, np.inf, np.inf), ) -resolutions = {"core": {"resolution": 0.025, "distance": 2.}} +resolutions = {"core": {"resolution": 0.025, "distance": 2.0}} -mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.6)) +mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.2)) num_modes = 2 -#For SFWM we have energy conservation and momemtum(k) conservation for 2pumps and signal+idler +# For SFWM we have energy conservation and momemtum(k) conservation for 2pumps and signal+idler lambda_p0 = 1.4 lambda_s0 = 1.097 lambda_i0 = 1.686 +# lambda_p0 = 1.55 +# lambda_s0 = 1.55 +# lambda_i0 = 1.55 + basis0 = Basis(mesh, ElementTriP0()) epsilon_p = basis0.zeros() @@ -108,7 +107,9 @@ def n_silicon_dioxide(wavelength): epsilon_i = basis0.zeros() -for wavelength, epsilon in zip([lambda_p0, lambda_s0, lambda_i0], [epsilon_p, epsilon_s, epsilon_i]): +for wavelength, epsilon in zip( + [lambda_p0, lambda_s0, lambda_i0], [epsilon_p, epsilon_s, epsilon_i] +): for subdomain, n_func in { "core": n_X, "box": n_silicon_dioxide, @@ -122,6 +123,7 @@ def n_silicon_dioxide(wavelength): modes_s = compute_modes(basis0, epsilon_s, wavelength=lambda_s0, num_modes=num_modes, order=1) modes_i = compute_modes(basis0, epsilon_i, wavelength=lambda_i0, num_modes=num_modes, order=1) +# modes_p[0].show(modes_p[0].E.real) mode_p = max(modes_p, key=lambda mode: mode.te_fraction) mode_s = max(modes_s, key=lambda mode: mode.te_fraction) @@ -130,14 +132,15 @@ def n_silicon_dioxide(wavelength): A_eff = calculate_sfwm_Aeff(basis0, mode_p, mode_s, mode_i) print("Aeff in um2:", A_eff) -#Calculation for non-linear coef +# Calculation for non-linear coef chi_3 = 5e-21 # m^2/V^2 #7e-20? -lambda_p0_m = lambda_p0 * 1e-6 #to m +chi_3 = 7e-20 # ? +lambda_p0_m = lambda_p0 * 1e-6 # to m n_p0 = np.real(mode_p.n_eff) -A_eff_m2 = A_eff * 1e-12 #to m^2 +A_eff_m2 = A_eff * 1e-12 # to m^2 omega_p0 = 2 * np.pi * c / lambda_p0_m gamma = (3 * chi_3 * omega_p0) / (4 * epsilon_0 * c**2 * n_p0**2 * A_eff_m2) -print("gamma:",gamma) \ No newline at end of file +print("gamma:", gamma)