dispersion rewrite

This commit is contained in:
Benoît Sierro
2023-09-26 09:52:44 +02:00
parent d65afaa456
commit 59eba7220c
11 changed files with 228 additions and 205 deletions

View File

@@ -461,11 +461,13 @@ default_rules: list[Rule] = [
), ),
Rule("w0", units.m_rads, ["wavelength"]), Rule("w0", units.m_rads, ["wavelength"]),
Rule("l", units.m_rads, ["w"]), Rule("l", units.m_rads, ["w"]),
Rule("w0_ind", math.argclosest, ["w_for_disp", "w0"]), Rule("w0_ind", math.argclosest, ["w", "w0"]),
Rule("w_num", len, ["w"]), Rule("w_num", len, ["w"]),
Rule("dw", lambda w: w[1] - w[0]), Rule("dw", lambda w: w[1] - w[0]),
Rule(["fft", "ifft"], utils.fft_functions, priorities=1), Rule(["fft", "ifft"], utils.fft_functions, priorities=1),
Rule("wavelength_window", lambda dt: (max(100e-9, 2 * units.c * dt), 8e-6)), Rule("wavelength_window", lambda dt: (max(100e-9, 2 * units.c * dt), 8e-6)),
Rule("wavelength_window", fiber.valid_wavelength_window),
Rule("dispersion_ind", fiber.dispersion_indices),
# Pulse # Pulse
Rule("field_0", pulse.finalize_pulse), Rule("field_0", pulse.finalize_pulse),
Rule(["input_time", "input_field"], pulse.load_custom_field), Rule(["input_time", "input_field"], pulse.load_custom_field),
@@ -490,9 +492,8 @@ default_rules: list[Rule] = [
Rule("soliton_length", pulse.soliton_length), Rule("soliton_length", pulse.soliton_length),
Rule("c_to_a_factor", lambda: 1, priorities=-1), Rule("c_to_a_factor", lambda: 1, priorities=-1),
# Fiber Dispersion # Fiber Dispersion
Rule("w_for_disp", units.m_rads, ["wl_for_disp"]),
Rule("gas_info", materials.Gas), Rule("gas_info", materials.Gas),
Rule("chi_gas", lambda gas_info, wl_for_disp: gas_info.sellmeier.chi(wl_for_disp)), Rule("chi_gas", lambda gas_info, l: gas_info.sellmeier.chi(l)),
Rule("n_gas_2", materials.n_gas_2), Rule("n_gas_2", materials.n_gas_2),
Rule("n_eff", fiber.n_eff_hasan, conditions=dict(model="hasan")), Rule("n_eff", fiber.n_eff_hasan, conditions=dict(model="hasan")),
Rule("n_eff", fiber.n_eff_marcatili, conditions=dict(model="marcatili")), Rule("n_eff", fiber.n_eff_marcatili, conditions=dict(model="marcatili")),
@@ -511,7 +512,7 @@ default_rules: list[Rule] = [
Rule("beta2_arr", fiber.beta2), Rule("beta2_arr", fiber.beta2),
Rule( Rule(
"zero_dispersion_wavelength", "zero_dispersion_wavelength",
lambda beta2_arr, wl_for_disp: wl_for_disp[math.argclosest(beta2_arr, 0)], lambda beta2_arr, l: l[math.argclosest(beta2_arr, 0)],
), ),
# Fiber nonlinearity # Fiber nonlinearity
Rule("effective_area", fiber.effective_area_pcf), Rule("effective_area", fiber.effective_area_pcf),
@@ -560,15 +561,11 @@ envelope_rules = default_rules + [
Rule("pre_field_0", pulse.initial_field_envelope, priorities=1), Rule("pre_field_0", pulse.initial_field_envelope, priorities=1),
Rule("c_to_a_factor", pulse.c_to_a_factor), Rule("c_to_a_factor", pulse.c_to_a_factor),
# Dispersion # Dispersion
Rule(["wl_for_disp", "dispersion_ind"], fiber.lambda_for_envelope_dispersion), Rule("beta2_coefficients", fiber.auto_dispersion_coefficients),
Rule("beta2_coefficients", fiber.dispersion_coefficients),
Rule("beta2_arr", fiber.dispersion_from_coefficients), Rule("beta2_arr", fiber.dispersion_from_coefficients),
Rule("beta2", lambda beta2_coefficients: beta2_coefficients[0]), Rule("beta2", lambda beta2_coefficients: beta2_coefficients[0]),
Rule( Rule("beta2", lambda beta2_arr, w0_ind: beta2_arr[w0_ind]),
["wl_for_disp", "beta2_arr", "wavelength_window"], Rule(["beta2_arr", "dispersion_ind"], fiber.load_custom_dispersion, priorities=[2, 2, 2]),
fiber.load_custom_dispersion,
priorities=[2, 2, 2],
),
# Operators # Operators
Rule("gamma_op", operators.variable_gamma, priorities=2), Rule("gamma_op", operators.variable_gamma, priorities=2),
Rule("gamma_op", operators.constant_quantity, ["gamma_arr"], priorities=1), Rule("gamma_op", operators.constant_quantity, ["gamma_arr"], priorities=1),
@@ -595,7 +592,6 @@ full_field_rules = default_rules + [
Rule("pre_field_0", pulse.initial_full_field), Rule("pre_field_0", pulse.initial_full_field),
Rule("spectrum_factor", pulse.spectrum_factor_fullfield, priorities=-1), Rule("spectrum_factor", pulse.spectrum_factor_fullfield, priorities=-1),
# Dispersion # Dispersion
Rule(["wl_for_disp", "dispersion_ind"], fiber.lambda_for_full_field_dispersion),
Rule("frame_velocity", fiber.frame_velocity), Rule("frame_velocity", fiber.frame_velocity),
Rule("beta2", lambda beta2_arr, w0_ind: beta2_arr[w0_ind]), Rule("beta2", lambda beta2_arr, w0_ind: beta2_arr[w0_ind]),
# Nonlinearity # Nonlinearity

View File

@@ -67,12 +67,12 @@ class ConstantGas:
pressure: float, pressure: float,
temperature: float, temperature: float,
ideal_gas: bool, ideal_gas: bool,
wl_for_disp: np.ndarray, l: np.ndarray,
): ):
self.gas = materials.Gas(gas_name) self.gas = materials.Gas(gas_name)
self.pressure_const = pressure self.pressure_const = pressure
self.number_density_const = self.gas.number_density(temperature, pressure, ideal_gas) self.number_density_const = self.gas.number_density(temperature, pressure, ideal_gas)
self.n_gas_2_const = self.gas.sellmeier.n_gas_2(wl_for_disp, temperature, pressure) self.n_gas_2_const = self.gas.sellmeier.n_gas_2(l, temperature, pressure)
self.n2_const = self.gas.n2(temperature, pressure) self.n2_const = self.gas.n2(temperature, pressure)
self.chi3_const = self.gas.chi3(temperature, pressure) self.chi3_const = self.gas.chi3(temperature, pressure)
@@ -102,7 +102,7 @@ class PressureGradientGas:
pressure_out: float, pressure_out: float,
temperature: float, temperature: float,
ideal_gas: bool, ideal_gas: bool,
wl_for_disp: np.ndarray, l: np.ndarray,
length: float, length: float,
): ):
self.p_in = pressure_in self.p_in = pressure_in
@@ -110,7 +110,7 @@ class PressureGradientGas:
self.gas = materials.Gas(gas_name) self.gas = materials.Gas(gas_name)
self.temperature = temperature self.temperature = temperature
self.ideal_gas = ideal_gas self.ideal_gas = ideal_gas
self.wl_for_disp = wl_for_disp self.l = l
self.z_max = length self.z_max = length
def _pressure(self, z: float) -> float: def _pressure(self, z: float) -> float:
@@ -120,7 +120,7 @@ class PressureGradientGas:
return self.gas.number_density(self.temperature, self._pressure(z), self.ideal_gas) return self.gas.number_density(self.temperature, self._pressure(z), self.ideal_gas)
def square_index(self, z: float) -> np.ndarray: def square_index(self, z: float) -> np.ndarray:
return self.gas.sellmeier.n_gas_2(self.wl_for_disp, self.temperature, self._pressure(z)) return self.gas.sellmeier.n_gas_2(self.l, self.temperature, self._pressure(z))
def n2(self, z: float) -> np.ndarray: def n2(self, z: float) -> np.ndarray:
return self.gas.n2(self.temperature, self._pressure(z)) return self.gas.n2(self.temperature, self._pressure(z))
@@ -135,19 +135,19 @@ class PressureGradientGas:
def marcatili_refractive_index( def marcatili_refractive_index(
square_index: VariableQuantity, core_radius: float, wl_for_disp: np.ndarray square_index: VariableQuantity, core_radius: float, l: np.ndarray
) -> VariableQuantity: ) -> VariableQuantity:
def operate(z: float) -> np.ndarray: def operate(z: float) -> np.ndarray:
return fiber.n_eff_marcatili(wl_for_disp, square_index(z), core_radius) return fiber.n_eff_marcatili(l, square_index(z), core_radius)
return operate return operate
def marcatili_adjusted_refractive_index( def marcatili_adjusted_refractive_index(
square_index: VariableQuantity, core_radius: float, wl_for_disp: np.ndarray square_index: VariableQuantity, core_radius: float, l: np.ndarray
) -> VariableQuantity: ) -> VariableQuantity:
def operate(z: float) -> np.ndarray: def operate(z: float) -> np.ndarray:
return fiber.n_eff_marcatili_adjusted(wl_for_disp, square_index(z), core_radius) return fiber.n_eff_marcatili_adjusted(l, square_index(z), core_radius)
return operate return operate
@@ -155,7 +155,7 @@ def marcatili_adjusted_refractive_index(
def vincetti_refractive_index( def vincetti_refractive_index(
square_index: VariableQuantity, square_index: VariableQuantity,
core_radius: float, core_radius: float,
wl_for_disp: np.ndarray, l: np.ndarray,
wavelength: float, wavelength: float,
wall_thickness: float, wall_thickness: float,
tube_radius: float, tube_radius: float,
@@ -165,7 +165,7 @@ def vincetti_refractive_index(
) -> VariableQuantity: ) -> VariableQuantity:
def operate(z: float) -> np.ndarray: def operate(z: float) -> np.ndarray:
return fiber.n_eff_vincetti( return fiber.n_eff_vincetti(
wl_for_disp, l,
wavelength, wavelength,
square_index(z), square_index(z),
wall_thickness, wall_thickness,
@@ -200,7 +200,7 @@ def constant_polynomial_dispersion(
def constant_direct_dispersion( def constant_direct_dispersion(
w_for_disp: np.ndarray, w: np.ndarray,
w0: np.ndarray, w0: np.ndarray,
t_num: int, t_num: int,
n_eff: np.ndarray, n_eff: np.ndarray,
@@ -211,8 +211,8 @@ def constant_direct_dispersion(
Direct dispersion for when the refractive index is known Direct dispersion for when the refractive index is known
""" """
disp_arr = np.zeros(t_num, dtype=complex) disp_arr = np.zeros(t_num, dtype=complex)
w0_ind = math.argclosest(w_for_disp, w0) w0_ind = math.argclosest(w, w0)
disp_arr[dispersion_ind] = fiber.fast_direct_dispersion(w_for_disp, w0, n_eff, w0_ind)[2:-2] disp_arr[dispersion_ind] = fiber.fast_direct_dispersion(w[dispersion_ind], w0, n_eff, w0_ind)
left_ind, *_, right_ind = np.nonzero(disp_arr[w_order])[0] left_ind, *_, right_ind = np.nonzero(disp_arr[w_order])[0]
disp_arr[w_order] = math._polynom_extrapolation_in_place( disp_arr[w_order] = math._polynom_extrapolation_in_place(
disp_arr[w_order], left_ind, right_ind, 1 disp_arr[w_order], left_ind, right_ind, 1
@@ -222,19 +222,19 @@ def constant_direct_dispersion(
def direct_dispersion( def direct_dispersion(
w_for_disp: np.ndarray, w: np.ndarray,
w0: np.ndarray, w0: np.ndarray,
t_num: int, t_num: int,
n_eff_op: SpecOperator, n_eff_op: SpecOperator,
dispersion_ind: np.ndarray, dispersion_ind: np.ndarray,
) -> VariableQuantity: ) -> VariableQuantity:
disp_arr = np.zeros(t_num, dtype=complex) disp_arr = np.zeros(t_num, dtype=complex)
w0_ind = math.argclosest(w_for_disp, w0) w0_ind = math.argclosest(w, w0)
def operate(z: float) -> np.ndarray: def operate(z: float) -> np.ndarray:
disp_arr[dispersion_ind] = fiber.fast_direct_dispersion( disp_arr[dispersion_ind] = fiber.fast_direct_dispersion(
w_for_disp, w0, n_eff_op(z), w0_ind w[dispersion_ind], w0, n_eff_op(z), w0_ind
)[2:-2] )
return disp_arr return disp_arr
return operate return operate
@@ -247,13 +247,13 @@ def direct_dispersion(
def constant_wave_vector( def constant_wave_vector(
n_eff: np.ndarray, n_eff: np.ndarray,
w_for_disp: np.ndarray, w: np.ndarray,
w_num: int, w_num: int,
dispersion_ind: np.ndarray, dispersion_ind: np.ndarray,
w_order: np.ndarray, w_order: np.ndarray,
): ):
beta_arr = np.zeros(w_num, dtype=float) beta_arr = np.zeros(w_num, dtype=float)
beta_arr[dispersion_ind] = fiber.beta(w_for_disp, n_eff)[2:-2] beta_arr[dispersion_ind] = fiber.beta(w[dispersion_ind], n_eff)
left_ind, *_, right_ind = np.nonzero(beta_arr[w_order])[0] left_ind, *_, right_ind = np.nonzero(beta_arr[w_order])[0]
beta_arr[w_order] = math._polynom_extrapolation_in_place( beta_arr[w_order] = math._polynom_extrapolation_in_place(
beta_arr[w_order], left_ind, right_ind, 1.0 beta_arr[w_order], left_ind, right_ind, 1.0
@@ -267,7 +267,7 @@ def constant_wave_vector(
################################################## ##################################################
def envelope_raman(hr_w: np.ndarra, raman_fraction: float) -> FieldOperator: def envelope_raman(hr_w: np.ndarray, raman_fraction: float) -> FieldOperator:
def operate(field: np.ndarray, z: float) -> np.ndarray: def operate(field: np.ndarray, z: float) -> np.ndarray:
return raman_fraction * np.fft.ifft(hr_w * np.fft.fft(math.abs2(field))) return raman_fraction * np.fft.ifft(hr_w * np.fft.fft(math.abs2(field)))

View File

@@ -646,7 +646,6 @@ class Parameters:
"t", "t",
"z_targets", "z_targets",
"l", "l",
"wl_for_disp",
"alpha", "alpha",
"gamma_arr", "gamma_arr",
"effective_area_arr", "effective_area_arr",

View File

@@ -1,12 +1,10 @@
import warnings import warnings
from io import BytesIO
from typing import Iterable, TypeVar from typing import Iterable, TypeVar
import numpy as np import numpy as np
from numpy import e from numpy import e
from numpy.fft import fft from numpy.fft import fft
from numpy.polynomial.chebyshev import Chebyshev, cheb2poly from numpy.polynomial.chebyshev import Chebyshev, cheb2poly
from scipy.interpolate import interp1d
from scgenerator import io from scgenerator import io
from scgenerator.io import DataFile from scgenerator.io import DataFile
@@ -19,56 +17,12 @@ pipi = 2 * pi
T = TypeVar("T") T = TypeVar("T")
def lambda_for_envelope_dispersion( def valid_wavelength_window(l: np.ndarray, dispersion_ind: np.ndarray) -> tuple[float, float]:
l: np.ndarray, wavelength_window: tuple[float, float] return l[dispersion_ind].min(), l[dispersion_ind].max()
) -> tuple[np.ndarray, np.ndarray]:
"""
Returns a wl vector for dispersion calculation in envelope mode
Returns
-------
np.ndarray
subset of l in the interpolation range with two extra values on each side
to accomodate for taking gradients
np.ndarray
indices of the original l where the values are valid
(i.e. without the two extra on each side)
"""
raw_indices = np.where((l >= wavelength_window[0]) & (l <= wavelength_window[1]))[0]
if l[raw_indices].min() > 1.01 * wavelength_window[0]:
raise ValueError(
f"lower range of {1e9*wavelength_window[0]:.1f}nm is not reached by the grid. "
f"Minimum of grid is {1e9*l[raw_indices].min():.1f}nm. Try a finer grid"
)
sorted_ind = raw_indices[np.argsort(l[raw_indices])]
return l[sorted_ind], sorted_ind[2:-2]
def lambda_for_full_field_dispersion( def dispersion_indices(l: np.ndarray, wavelength_window) -> np.ndarray:
l: np.ndarray, wavelength_window: tuple[float, float] return np.where((l >= wavelength_window[0]) & (l <= wavelength_window[1]))[0]
) -> tuple[np.ndarray, np.ndarray]:
"""
Returns a wl vector for dispersion calculation in full field mode
Returns
-------
np.ndarray
subset of l in the interpolation range with two extra values on each side
to accomodate for taking gradients
np.ndarray
indices of the original l where the values are valid
(i.e. without the two extra on each side)
"""
su = np.where((l >= wavelength_window[0]) & (l <= wavelength_window[1]))[0]
if l[su].min() > 1.01 * wavelength_window[0]:
raise ValueError(
f"lower range of {1e9*wavelength_window[0]:.1f}nm is not reached by the grid. "
"try a finer grid"
)
fu = np.concatenate((su[:2] - 2, su, su[-2:] + 2))
fu = np.where(fu < 0, 0, fu)
fu = np.where(fu >= len(l), len(l) - 1, fu)
return l[fu], su
def is_dynamic_dispersion(pressure=None): def is_dynamic_dispersion(pressure=None):
@@ -92,45 +46,23 @@ def is_dynamic_dispersion(pressure=None):
return out return out
def gvd_from_n_eff(n_eff: np.ndarray, wl_for_disp: np.ndarray): def beta2_to_D(beta2, l):
""" """returns the D parameter corresponding to beta2(l)"""
computes the dispersion parameter D from an effective index of refraction n_eff return -(pipi * c) / (l**2) * beta2
Since computing gradients/derivatives of discrete arrays is not well defined on the boundary,
it is advised to chop off the two values on either end of the returned array
Parameters
----------
n_eff : 1D array
a wl-dependent index of refraction
wl_for_disp : 1D array
the wavelength array (len must match n_eff)
Returns
-------
D : 1D array
wl-dependent dispersion parameter as function of wl_for_disp
"""
return -wl_for_disp / c * (np.gradient(np.gradient(n_eff, wl_for_disp), wl_for_disp))
def beta2_to_D(beta2, wl_for_disp): def D_to_beta2(D, l):
"""returns the D parameter corresponding to beta2(wl_for_disp)""" """returns the beta2 parameters corresponding to D(l)"""
return -(pipi * c) / (wl_for_disp**2) * beta2 return -(l**2) / (pipi * c) * D
def D_to_beta2(D, wl_for_disp): def plasma_dispersion(l, number_density, simple=False):
"""returns the beta2 parameters corresponding to D(wl_for_disp)"""
return -(wl_for_disp**2) / (pipi * c) * D
def plasma_dispersion(wl_for_disp, number_density, simple=False):
""" """
computes dispersion (beta2) for constant plasma computes dispersion (beta2) for constant plasma
Parameters Parameters
---------- ----------
wl_for_disp : array-like l : array-like
wavelengths over which to calculate the dispersion wavelengths over which to calculate the dispersion
number_density : number of ionized atoms /m^3 number_density : number of ionized atoms /m^3
@@ -141,7 +73,7 @@ def plasma_dispersion(wl_for_disp, number_density, simple=False):
""" """
e2_me_e0 = 3182.60735 # e^2 /(m_e * epsilon_0) e2_me_e0 = 3182.60735 # e^2 /(m_e * epsilon_0)
w = units.m_rads(wl_for_disp) w = units.m_rads(l)
if simple: if simple:
w_pl = number_density * e2_me_e0 w_pl = number_density * e2_me_e0
return -(w_pl**2) / (c * w**2) return -(w_pl**2) / (c * w**2)
@@ -150,16 +82,16 @@ def plasma_dispersion(wl_for_disp, number_density, simple=False):
return beta2_arr return beta2_arr
def n_eff_marcatili(wl_for_disp, n_gas_2, core_radius, he_mode=(1, 1)): def n_eff_marcatili(l, n_gas_2, core_radius, he_mode=(1, 1)):
""" """
computes the effective refractive index according to the Marcatili model of a capillary computes the effective refractive index according to the Marcatili model of a capillary
Parameters Parameters
---------- ----------
wl_for_disp : ndarray, shape (n, ) l : ndarray, shape (n, )
wavelengths array (m) wavelengths array (m)
n_gas_2 : ndarray, shape (n, ) n_gas_2 : ndarray, shape (n, )
square of the refractive index of the gas as function of wl_for_disp square of the refractive index of the gas as function of l
core_radius : float core_radius : float
inner radius of the capillary (m) inner radius of the capillary (m)
he_mode : tuple, shape (2, ), optional he_mode : tuple, shape (2, ), optional
@@ -175,20 +107,20 @@ def n_eff_marcatili(wl_for_disp, n_gas_2, core_radius, he_mode=(1, 1)):
""" """
u = u_nm(*he_mode) u = u_nm(*he_mode)
return np.sqrt(n_gas_2 - (wl_for_disp * u / (pipi * core_radius)) ** 2) return np.sqrt(n_gas_2 - (l * u / (pipi * core_radius)) ** 2)
def n_eff_marcatili_adjusted(wl_for_disp, n_gas_2, core_radius, he_mode=(1, 1), fit_parameters=()): def n_eff_marcatili_adjusted(l, n_gas_2, core_radius, he_mode=(1, 1), fit_parameters=()):
""" """
computes the effective refractive index according to the Marcatili model of a capillary computes the effective refractive index according to the Marcatili model of a capillary
but adjusted at longer wavelengths but adjusted at longer wavelengths
Parameters Parameters
---------- ----------
wl_for_disp : ndarray, shape (n, ) l : ndarray, shape (n, )
wavelengths array (m) wavelengths array (m)
n_gas_2 : ndarray, shape (n, ) n_gas_2 : ndarray, shape (n, )
refractive index of the gas as function of wl_for_disp refractive index of the gas as function of l
core_radius : float core_radius : float
inner radius of the capillary (m) inner radius of the capillary (m)
he_mode : tuple, shape (2, ), optional he_mode : tuple, shape (2, ), optional
@@ -207,9 +139,9 @@ def n_eff_marcatili_adjusted(wl_for_disp, n_gas_2, core_radius, he_mode=(1, 1),
""" """
u = u_nm(*he_mode) u = u_nm(*he_mode)
corrected_radius = effective_core_radius(wl_for_disp, core_radius, *fit_parameters) corrected_radius = effective_core_radius(l, core_radius, *fit_parameters)
return np.sqrt(n_gas_2 - (wl_for_disp * u / (pipi * corrected_radius)) ** 2) return np.sqrt(n_gas_2 - (l * u / (pipi * corrected_radius)) ** 2)
def effective_area_marcatili(core_radius: float) -> float: def effective_area_marcatili(core_radius: float) -> float:
@@ -238,14 +170,14 @@ def capillary_spacing_hasan(
def resonance_thickness( def resonance_thickness(
wl_for_disp: np.ndarray, order: int, n_gas_2: np.ndarray, core_radius: float l: np.ndarray, order: int, n_gas_2: np.ndarray, core_radius: float
) -> float: ) -> float:
""" """
computes at which wall thickness the specified wl is resonant computes at which wall thickness the specified wl is resonant
Parameters Parameters
---------- ----------
wl_for_disp : np.ndarray l : np.ndarray
in m in m
order : int order : int
0, 1, ... 0, 1, ...
@@ -259,20 +191,20 @@ def resonance_thickness(
float float
thickness in m thickness in m
""" """
n_si_2 = mat.n_gas_2(wl_for_disp, "silica", None, None) n_si_2 = mat.n_gas_2(l, "silica", None, None)
return ( return (
wl_for_disp l
/ (4 * np.sqrt(n_si_2)) / (4 * np.sqrt(n_si_2))
* (2 * order + 1) * (2 * order + 1)
* (1 - n_gas_2 / n_si_2 + wl_for_disp**2 / (4 * n_si_2 * core_radius**2)) ** -0.5 * (1 - n_gas_2 / n_si_2 + l**2 / (4 * n_si_2 * core_radius**2)) ** -0.5
) )
def resonance_strength( def resonance_strength(
wl_for_disp: np.ndarray, core_radius: float, capillary_thickness: float, order: int l: np.ndarray, core_radius: float, capillary_thickness: float, order: int
) -> float: ) -> float:
a = 1.83 + (2.3 * capillary_thickness / core_radius) a = 1.83 + (2.3 * capillary_thickness / core_radius)
n_si = np.sqrt(mat.n_gas_2(wl_for_disp, "silica", None, None)) n_si = np.sqrt(mat.n_gas_2(l, "silica", None, None))
return ( return (
capillary_thickness capillary_thickness
/ (n_si * core_radius) ** (2.303 * a / n_si) / (n_si * core_radius) ** (2.303 * a / n_si)
@@ -281,19 +213,19 @@ def resonance_strength(
def capillary_resonance_strengths( def capillary_resonance_strengths(
wl_for_disp: np.ndarray, l: np.ndarray,
core_radius: float, core_radius: float,
capillary_thickness: float, capillary_thickness: float,
capillary_resonance_max_order: int, capillary_resonance_max_order: int,
) -> list[float]: ) -> list[float]:
return [ return [
resonance_strength(wl_for_disp, core_radius, capillary_thickness, o) resonance_strength(l, core_radius, capillary_thickness, o)
for o in range(1, capillary_resonance_max_order + 1) for o in range(1, capillary_resonance_max_order + 1)
] ]
def n_eff_hasan( def n_eff_hasan(
wl_for_disp: np.ndarray, l: np.ndarray,
n_gas_2: np.ndarray, n_gas_2: np.ndarray,
core_radius: float, core_radius: float,
capillary_num: int, capillary_num: int,
@@ -308,10 +240,10 @@ def n_eff_hasan(
Parameters Parameters
---------- ----------
wl_for_disp : np.ndarray, shape(n,) l : np.ndarray, shape(n,)
wavelenghs array wavelenghs array
n_gas_2 : ndarray, shape (n, ) n_gas_2 : ndarray, shape (n, )
squared refractive index of the gas as a function of wl_for_disp squared refractive index of the gas as a function of l
core_radius : float core_radius : float
radius of the core (from cented to edge of a capillary) radius of the core (from cented to edge of a capillary)
capillary_num : int capillary_num : int
@@ -346,18 +278,18 @@ def n_eff_hasan(
if capillary_nested > 0: if capillary_nested > 0:
f2 += 0.0045 * np.exp(-4.1589 / (capillary_nested * Rg)) f2 += 0.0045 * np.exp(-4.1589 / (capillary_nested * Rg))
R_eff = f1 * core_radius * (1 - f2 * wl_for_disp**2 / (core_radius * capillary_thickness)) R_eff = f1 * core_radius * (1 - f2 * l**2 / (core_radius * capillary_thickness))
n_eff_2 = n_gas_2 - (u * wl_for_disp / (pipi * R_eff)) ** 2 n_eff_2 = n_gas_2 - (u * l / (pipi * R_eff)) ** 2
chi_sil = mat.Sellmeier.load("silica").chi(wl_for_disp) chi_sil = mat.Sellmeier.load("silica").chi(l)
with np.errstate(divide="ignore", invalid="ignore"): with np.errstate(divide="ignore", invalid="ignore"):
for m, strength in enumerate(capillary_resonance_strengths): for m, strength in enumerate(capillary_resonance_strengths):
n_eff_2 += ( n_eff_2 += (
strength strength
* wl_for_disp**2 * l**2
/ (alpha + wl_for_disp**2 - chi_sil * (2 * capillary_thickness / (m + 1)) ** 2) / (alpha + l**2 - chi_sil * (2 * capillary_thickness / (m + 1)) ** 2)
) )
return np.sqrt(n_eff_2) return np.sqrt(n_eff_2)
@@ -479,15 +411,15 @@ def effective_area_from_V(core_radius: float, V_eff: T) -> T:
return out return out
def beta(w_for_disp: np.ndarray, n_eff: np.ndarray) -> np.ndarray: def beta(w: np.ndarray, n_eff: np.ndarray) -> np.ndarray:
return n_eff * w_for_disp / c return n_eff * w / c
def beta1(w_for_disp: np.ndarray, n_eff: np.ndarray) -> np.ndarray: def beta1(w: np.ndarray, n_eff: np.ndarray) -> np.ndarray:
return np.gradient(beta(w_for_disp, n_eff), w_for_disp) return np.gradient(beta(w, n_eff), w[1] - w[0])
def beta2(w_for_disp: np.ndarray, n_eff: np.ndarray) -> np.ndarray: def beta2(w: np.ndarray, n_eff: np.ndarray) -> np.ndarray:
""" """
computes the dispersion parameter beta2 according to the effective refractive computes the dispersion parameter beta2 according to the effective refractive
index of the fiber and the frequency range index of the fiber and the frequency range
@@ -503,7 +435,8 @@ def beta2(w_for_disp: np.ndarray, n_eff: np.ndarray) -> np.ndarray:
------- -------
beta2 : ndarray, shape (n, ) beta2 : ndarray, shape (n, )
""" """
return np.gradient(np.gradient(beta(w_for_disp, n_eff), w_for_disp), w_for_disp) dw = w[1] - w[0]
return np.gradient(np.gradient(beta(w, n_eff), dw), dw)
def frame_velocity(beta1_arr: np.ndarray, w0_ind: int) -> float: def frame_velocity(beta1_arr: np.ndarray, w0_ind: int) -> float:
@@ -511,7 +444,7 @@ def frame_velocity(beta1_arr: np.ndarray, w0_ind: int) -> float:
def HCPCF_dispersion( def HCPCF_dispersion(
wl_for_disp, l,
gas_name: str, gas_name: str,
model="marcatili", model="marcatili",
pressure=None, pressure=None,
@@ -523,7 +456,7 @@ def HCPCF_dispersion(
Parameters Parameters
---------- ----------
wl_for_disp : ndarray, shape (n, ) l : ndarray, shape (n, )
wavelengths over which to calculate the dispersion wavelengths over which to calculate the dispersion
gas_name : str gas_name : str
name of the filling gas in lower case name of the filling gas in lower case
@@ -532,7 +465,7 @@ def HCPCF_dispersion(
model_params : tuple model_params : tuple
to be passed on to the function in charge of computing the effective index of the fiber. to be passed on to the function in charge of computing the effective index of the fiber.
Every n_eff_* function has a signature Every n_eff_* function has a signature
`n_eff_(wl_for_disp, n_gas_2, radius, *args)` and model_params corresponds to args `n_eff_(l, n_gas_2, radius, *args)` and model_params corresponds to args
temperature : float temperature : float
Temperature of the material Temperature of the material
pressure : float pressure : float
@@ -544,13 +477,13 @@ def HCPCF_dispersion(
beta2 as function of wavelength beta2 as function of wavelength
""" """
w = units.m_rads(wl_for_disp) w = units.m_rads(l)
n_gas_2 = mat.Sellmeier.load(gas_name).n_gas_2(wl_for_disp, temperature, pressure) n_gas_2 = mat.Sellmeier.load(gas_name).n_gas_2(l, temperature, pressure)
n_eff_func = dict( n_eff_func = dict(
marcatili=n_eff_marcatili, marcatili_adjusted=n_eff_marcatili_adjusted, hasan=n_eff_hasan marcatili=n_eff_marcatili, marcatili_adjusted=n_eff_marcatili_adjusted, hasan=n_eff_hasan
)[model] )[model]
n_eff = n_eff_func(wl_for_disp, n_gas_2, **model_params) n_eff = n_eff_func(l, n_gas_2, **model_params)
return beta2(w, n_eff) return beta2(w, n_eff)
@@ -558,23 +491,27 @@ def HCPCF_dispersion(
def gamma_parameter(n2: float, w0: float, effective_area: T) -> T: def gamma_parameter(n2: float, w0: float, effective_area: T) -> T:
if isinstance(effective_area, np.ndarray): if isinstance(effective_area, np.ndarray):
effective_area_term = np.sqrt(effective_area * effective_area[0]) effective_area_term = np.sqrt(effective_area * effective_area[0])
else: return np.divide(
effective_area_term = effective_area n2 * w0,
effective_area_term * c,
where=effective_area_term != 0,
out=effective_area_term,
)
return n2 * w0 / (effective_area_term * c) return n2 * w0 / (effective_area * c)
def constant_effective_area_arr(l: np.ndarray, effective_area: float) -> np.ndarray: def constant_effective_area_arr(l: np.ndarray, effective_area: float) -> np.ndarray:
return np.ones_like(l) * effective_area return np.ones_like(l) * effective_area
def n_eff_pcf(wl_for_disp: np.ndarray, pcf_pitch: float, pcf_pitch_ratio: float) -> np.ndarray: def n_eff_pcf(l: np.ndarray, pcf_pitch: float, pcf_pitch_ratio: float) -> np.ndarray:
""" """
semi-analytical computation of the dispersion profile of a triangular Index-guiding PCF semi-analytical computation of the dispersion profile of a triangular Index-guiding PCF
Parameters Parameters
---------- ----------
wl_for_disp : np.ndarray, shape (n,) l : np.ndarray, shape (n,)
wavelengths over which to calculate the dispersion wavelengths over which to calculate the dispersion
pcf_pitch : float pcf_pitch : float
distance between air holes in m distance between air holes in m
@@ -601,18 +538,18 @@ def n_eff_pcf(wl_for_disp: np.ndarray, pcf_pitch: float, pcf_pitch_ratio: float)
r_eff = pcf_pitch / np.sqrt(3) r_eff = pcf_pitch / np.sqrt(3)
pi2a = pipi * r_eff pi2a = pipi * r_eff
ratio_l = wl_for_disp / pcf_pitch ratio_l = l / pcf_pitch
A, B = saitoh_paramters(pcf_pitch_ratio) A, B = saitoh_paramters(pcf_pitch_ratio)
V = A[0] + A[1] / (1 + A[2] * np.exp(A[3] * ratio_l)) V = A[0] + A[1] / (1 + A[2] * np.exp(A[3] * ratio_l))
W = B[0] + B[1] / (1 + B[2] * np.exp(B[3] * ratio_l)) W = B[0] + B[1] / (1 + B[2] * np.exp(B[3] * ratio_l))
n_FSM2 = 1.45**2 - (wl_for_disp * V / (pi2a)) ** 2 n_FSM2 = 1.45**2 - (l * V / (pi2a)) ** 2
n_eff2 = (wl_for_disp * W / (pi2a)) ** 2 + n_FSM2 n_eff2 = (l * W / (pi2a)) ** 2 + n_FSM2
n_eff = np.sqrt(n_eff2) n_eff = np.sqrt(n_eff2)
chi_mat = mat.Sellmeier.load("silica").chi(wl_for_disp) chi_mat = mat.Sellmeier.load("silica").chi(l)
return n_eff + np.sqrt(chi_mat + 1) return n_eff + np.sqrt(chi_mat + 1)
@@ -675,14 +612,21 @@ def load_custom_effective_area(effective_area_file: DataFile, l: np.ndarray) ->
wl, effective_area = effective_area_file.load_arrays( wl, effective_area = effective_area_file.load_arrays(
("wavelength", "wl"), ("A_eff", "effective_area") ("wavelength", "wl"), ("A_eff", "effective_area")
) )
return interp1d(wl, effective_area, fill_value=1, bounds_error=False)(l) return np.interp(l, wl, effective_area, left=0, right=0)
def load_custom_dispersion(dispersion_file: DataFile) -> tuple[np.ndarray, np.ndarray]: def load_custom_dispersion(
wl_for_disp, D = dispersion_file.load_arrays(("wavelength", "wl"), "dispersion") l: np.ndarray, dispersion_file: DataFile
interp_range = (np.min(wl_for_disp), np.max(wl_for_disp)) ) -> tuple[np.ndarray, np.ndarray]:
beta2 = D_to_beta2(D, wl_for_disp) """loads a custom dispersion and returns beta2(l) where it is specified, as well as indices"""
return wl_for_disp, beta2, interp_range out = np.zeros_like(l)
wl_file, D = dispersion_file.load_arrays(("wavelength", "wl"), "dispersion")
wl_file = units.normalize_wavelengths(wl_file)
ind = np.where((l >= wl_file.min()) & (l <= wl_file.max()))[0]
out[ind] = np.interp(l[ind], wl_file, D_to_beta2(D, wl_file))
return out, ind
def load_custom_loss(l: np.ndarray, loss_file: DataFile) -> np.ndarray: def load_custom_loss(l: np.ndarray, loss_file: DataFile) -> np.ndarray:
@@ -702,13 +646,24 @@ def load_custom_loss(l: np.ndarray, loss_file: DataFile) -> np.ndarray:
loss in 1/m units loss in 1/m units
""" """
wl, loss = loss_file.load_arrays(("wavelength", "wl"), "loss") wl, loss = loss_file.load_arrays(("wavelength", "wl"), "loss")
return interp1d(wl, loss, fill_value=0, bounds_error=False)(l) wl = units.normalize_wavelengths(wl)
return np.interp(l, wl, loss, left=0, right=0)
def auto_dispersion_coefficients(
w_c: np.ndarray,
dispersion_ind: np.ndarray,
beta2_arr: np.ndarray,
interpolation_degree: int,
) -> np.ndarray:
return dispersion_coefficients(
w_c[dispersion_ind], beta2_arr[dispersion_ind], interpolation_degree
)
def dispersion_coefficients( def dispersion_coefficients(
w_for_disp: np.ndarray, w_c: np.ndarray,
beta2_arr: np.ndarray, beta2_arr: np.ndarray,
w0: float,
interpolation_degree: int, interpolation_degree: int,
): ):
""" """
@@ -716,10 +671,10 @@ def dispersion_coefficients(
Parameters Parameters
---------- ----------
wl_for_disp : 1D array w : np.ndarray, shape (n,)
wavelength angular frequencies array
beta2 : 1D array beta2 : np.ndarray, shape (n,)
beta2 as function of wl_for_disp beta2 as function of w
w0 : float w0 : float
pump angular frequency pump angular frequency
interpolation_degree : int interpolation_degree : int
@@ -734,10 +689,7 @@ def dispersion_coefficients(
# we get the beta2 Taylor coeffiecients by making a fit around w0 # we get the beta2 Taylor coeffiecients by making a fit around w0
if interpolation_degree < 2: if interpolation_degree < 2:
raise ValueError(f"interpolation_degree must be at least 2, got {interpolation_degree}") raise ValueError(f"interpolation_degree must be at least 2, got {interpolation_degree}")
w_c = w_for_disp - w0
w_c = w_c[2:-2]
beta2_arr = beta2_arr[2:-2]
fit = Chebyshev.fit(w_c, beta2_arr, interpolation_degree) fit = Chebyshev.fit(w_c, beta2_arr, interpolation_degree)
poly_coef = cheb2poly(fit.convert().coef) poly_coef = cheb2poly(fit.convert().coef)
beta2_coef = poly_coef * np.cumprod([1] + list(range(1, interpolation_degree + 1))) beta2_coef = poly_coef * np.cumprod([1] + list(range(1, interpolation_degree + 1)))
@@ -815,7 +767,7 @@ def delayed_raman_t(t: np.ndarray, raman_type: str) -> np.ndarray:
print("Not able to find the measured Raman response function. Going with agrawal model") print("Not able to find the measured Raman response function. Going with agrawal model")
return delayed_raman_t(t, raman_type="agrawal") return delayed_raman_t(t, raman_type="agrawal")
hr_arr = interp1d(t_stored, hr_arr_stored, bounds_error=False, fill_value=0)(t) hr_arr = np.interp(t, t_stored, hr_arr_stored, left=0, right=0)
else: else:
print("invalid raman response function, aborting") print("invalid raman response function, aborting")
quit() quit()
@@ -895,13 +847,13 @@ def fast_direct_dispersion(w: np.ndarray, w0: float, n_eff: np.ndarray, w0_ind:
return -1j * (beta_arr - beta1_arr[w0_ind] * (w - w0) - beta_arr[w0_ind]) return -1j * (beta_arr - beta1_arr[w0_ind] * (w - w0) - beta_arr[w0_ind])
def effective_core_radius(wl_for_disp, core_radius, s=0.08, h=200e-9): def effective_core_radius(l, core_radius, s=0.08, h=200e-9):
""" """
return the variable core radius according to Eq. S2.2 from Köttig2017 return the variable core radius according to Eq. S2.2 from Köttig2017
Parameters Parameters
---------- ----------
wl_for_disp : ndarray, shape (n, ) l : ndarray, shape (n, )
array of wl over which to calculate the effective core radius array of wl over which to calculate the effective core radius
core_radius : float core_radius : float
physical core radius in m physical core radius in m
@@ -914,12 +866,12 @@ def effective_core_radius(wl_for_disp, core_radius, s=0.08, h=200e-9):
------- -------
effective_core_radius : ndarray, shape (n, ) effective_core_radius : ndarray, shape (n, )
""" """
return core_radius / (1 + s * wl_for_disp**2 / (core_radius * h)) return core_radius / (1 + s * l**2 / (core_radius * h))
def effective_radius_HCARF(core_radius, t, f1, f2, wl_for_disp): def effective_radius_HCARF(core_radius, t, f1, f2, l):
"""eq. 3 in Hasan 2018""" """eq. 3 in Hasan 2018"""
return f1 * core_radius * (1 - f2 * wl_for_disp**2 / (core_radius * t)) return f1 * core_radius * (1 - f2 * l**2 / (core_radius * t))
def scalar_loss(alpha: float, w_num: int) -> np.ndarray: def scalar_loss(alpha: float, w_num: int) -> np.ndarray:
@@ -952,15 +904,15 @@ def capillary_loss(wl: np.ndarray, he_mode: tuple[int, int], core_radius: float)
def safe_capillary_loss( def safe_capillary_loss(
wl_for_disp: np.ndarray, l: np.ndarray,
dispersion_ind: np.ndarray, dispersion_ind: np.ndarray,
w_num: int, w_num: int,
core_radius: float, core_radius: float,
he_mode: tuple[int, int], he_mode: tuple[int, int],
) -> np.ndarray: ) -> np.ndarray:
alpha = capillary_loss(wl_for_disp, he_mode, core_radius) alpha = capillary_loss(l[dispersion_ind], he_mode, core_radius)
loss_arr = np.zeros(w_num) loss_arr = np.zeros(w_num)
loss_arr[dispersion_ind] = alpha[2:-2] loss_arr[dispersion_ind] = alpha
return loss_arr return loss_arr
@@ -1128,7 +1080,7 @@ def n_eff_correction_vincetti(
def n_eff_vincetti( def n_eff_vincetti(
wl_for_disp: np.ndarray, l: np.ndarray,
wavelength: float, wavelength: float,
n_gas_2: np.ndarray, n_gas_2: np.ndarray,
thickness: float, thickness: float,
@@ -1141,7 +1093,7 @@ def n_eff_vincetti(
""" """
Parameters Parameters
---------- ----------
wl_for_disp : ndarray l : ndarray
wavelength (m) array over which to compute the refractive index wavelength (m) array over which to compute the refractive index
wavelength : float wavelength : float
center wavelength / pump wavelength center wavelength / pump wavelength
@@ -1182,18 +1134,16 @@ def n_eff_vincetti(
""" """
if n_clad_2 is None: if n_clad_2 is None:
n_clad_2 = mat.Sellmeier.load("silica").n_gas_2(wl_for_disp) n_clad_2 = mat.Sellmeier.load("silica").n_gas_2(l)
n_clad_0 = np.sqrt(n_clad_2[argclosest(wl_for_disp, wavelength)]) n_clad_0 = np.sqrt(n_clad_2[argclosest(l, wavelength)])
f = normalized_frequency_vincetti(wl_for_disp, thickness, n_clad_2, n_gas_2) f = normalized_frequency_vincetti(l, thickness, n_clad_2, n_gas_2)
r_co_eff = effective_core_radius_vincetti(wl_for_disp, f, tube_radius, gap, n_tubes) r_co_eff = effective_core_radius_vincetti(l, f, tube_radius, gap, n_tubes)
r_co = core_radius_from_capillaries(tube_radius, gap, n_tubes) r_co = core_radius_from_capillaries(tube_radius, gap, n_tubes)
d_n_eff = n_eff_correction_vincetti( d_n_eff = n_eff_correction_vincetti(l, f, thickness / tube_radius, r_co, n_clad_0, n_terms)
wl_for_disp, f, thickness / tube_radius, r_co, n_clad_0, n_terms
)
n_gas = np.sqrt(n_gas_2) n_gas = np.sqrt(n_gas_2)
# eq. (21) in [1] # eq. (21) in [1]
return n_gas - 0.125 / n_gas * (u_nm(1, 1) * wl_for_disp / (np.pi * r_co_eff)) ** 2 + d_n_eff return n_gas - 0.125 / n_gas * (u_nm(1, 1) * l / (np.pi * r_co_eff)) ** 2 + d_n_eff

View File

@@ -321,9 +321,9 @@ class Gas:
return f"{self.__class__.__name__}({self.name!r})" return f"{self.__class__.__name__}({self.name!r})"
def n_gas_2(wl_for_disp: np.ndarray, gas_name: str, pressure: float, temperature: float): def n_gas_2(l: np.ndarray, gas_name: str, pressure: float, temperature: float):
"""Returns the sqare of the index of refraction of the specified gas""" """Returns the sqare of the index of refraction of the specified gas"""
return Sellmeier.load(gas_name).n_gas_2(wl_for_disp, temperature, pressure) return Sellmeier.load(gas_name).n_gas_2(l, temperature, pressure)
def pressure_from_gradient(ratio: float, p0: float, p1: float) -> float: def pressure_from_gradient(ratio: float, p0: float, p1: float) -> float:

View File

@@ -461,3 +461,12 @@ def sort_axis(
out_ax = plt_range.unit.inv(cropped) out_ax = plt_range.unit.inv(cropped)
return out_ax, indices, (out_ax[0], out_ax[-1]) return out_ax, indices, (out_ax[0], out_ax[-1])
def normalize_wavelengths(wl: np.ndarray) -> np.ndarray:
"""attempts to guess the units of wl and return the same array in base SI units (m)"""
if np.all((wl > 0.01) & (wl < 20)):
return wl * 1e-6
elif np.all(wl >= 20):
return wl * 1e-9
return wl

Binary file not shown.

Binary file not shown.

Binary file not shown.

75
tests/test_dispersion.py Normal file
View File

@@ -0,0 +1,75 @@
import pytest
import scgenerator as sc
PARAMS1 = sc.Parameters(
name="Dudley2006",
wavelength=835e-9,
width=50e-15,
peak_power=10e3,
repetition_rate=20e6,
shape="sech",
# fiber
length=15e-2,
raman_type="measured",
beta2_coefficients=[
-11.830e-27,
8.1038e-42,
-9.5205e-57,
2.0737e-71,
-5.3943e-86,
1.3486e-100,
-2.5495e-115,
3.0524e-130,
-1.7140e-145,
],
gamma=0.11,
# simulation
tolerated_error=1e-6,
wavelength_window=[400e-9, 1500e-9],
t_num=8192,
quantum_noise=False,
z_num=128,
)
PARAMS2 = sc.Parameters(
name="dual comb",
wavelength=1056e-9,
width=78e-15,
peak_power=13e3,
repetition_rate=20e6,
shape="gaussian",
# fiber
length=20e-2,
n2=2.2e-20,
raman_type="measured",
dispersion_file="tests/data/PM1050NEG_dispersion.npz",
loss_file="tests/data/PM1050NEG_loss.npz",
effective_area_file="tests/data/PM1050NEG_aeff_slow.npz",
# simulation
interpolation_degree=12,
tolerated_error=1e-8,
wavelength_window=[500e-9, 1500e-9],
t_num=8192,
quantum_noise=True,
z_num=128,
)
def test_poly_dispersion_grid_size():
params = PARAMS1.compile()
disp_ind = params.compute("dispersion_ind")
assert disp_ind is not None
assert all((params.l[disp_ind] < 1500e-9) & (params.l[disp_ind] > 400e-9))
assert len(params.compute("beta2_arr")) == params.t_num
def test_custom_dispersion_grid_size():
params = PARAMS2.compile()
params.compute("t")
disp_ind = params.compute("dispersion_ind")
assert disp_ind is not None
assert all((params.l[disp_ind] < 1500e-9) & (params.l[disp_ind] > 400e-9))
assert len(params.compute("beta2_arr")) == params.t_num

View File

@@ -14,12 +14,6 @@ def test_scaling():
assert dt == pytest.approx(period / (nt - 1)) assert dt == pytest.approx(period / (nt - 1))
def test_wl_dispersion():
t = sc.tspace(t_num=1 << 15, dt=3.8e-15)
w = sc.wspace(t)
wl = sc.units.nm_rads(w + sc.units.nm_rads(1546)) * 1e-9
wl_disp, ind_disp = sc.fiber.lambda_for_envelope_dispersion(wl, (950e-9, 4000e-9))
assert all(np.diff(wl_disp) > 0)
def test_iwspace(): def test_iwspace():