From 59eba7220cc57f0f65f56289b29ba31cc273b8dc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Sierro?= Date: Tue, 26 Sep 2023 09:52:44 +0200 Subject: [PATCH] dispersion rewrite --- src/scgenerator/evaluator.py | 20 +- src/scgenerator/operators.py | 42 ++-- src/scgenerator/parameter.py | 1 - src/scgenerator/physics/fiber.py | 276 +++++++++++---------------- src/scgenerator/physics/materials.py | 4 +- src/scgenerator/physics/units.py | 9 + tests/data/PM1050NEG_aeff_slow.npz | Bin 0 -> 2310 bytes tests/data/PM1050NEG_dispersion.npz | Bin 0 -> 16542 bytes tests/data/PM1050NEG_loss.npz | Bin 0 -> 12866 bytes tests/test_dispersion.py | 75 ++++++++ tests/test_grid.py | 6 - 11 files changed, 228 insertions(+), 205 deletions(-) create mode 100644 tests/data/PM1050NEG_aeff_slow.npz create mode 100644 tests/data/PM1050NEG_dispersion.npz create mode 100644 tests/data/PM1050NEG_loss.npz create mode 100644 tests/test_dispersion.py diff --git a/src/scgenerator/evaluator.py b/src/scgenerator/evaluator.py index 3349e34..26413b6 100644 --- a/src/scgenerator/evaluator.py +++ b/src/scgenerator/evaluator.py @@ -461,11 +461,13 @@ default_rules: list[Rule] = [ ), Rule("w0", units.m_rads, ["wavelength"]), 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("dw", lambda w: w[1] - w[0]), 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", fiber.valid_wavelength_window), + Rule("dispersion_ind", fiber.dispersion_indices), # Pulse Rule("field_0", pulse.finalize_pulse), Rule(["input_time", "input_field"], pulse.load_custom_field), @@ -490,9 +492,8 @@ default_rules: list[Rule] = [ Rule("soliton_length", pulse.soliton_length), Rule("c_to_a_factor", lambda: 1, priorities=-1), # Fiber Dispersion - Rule("w_for_disp", units.m_rads, ["wl_for_disp"]), 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_eff", fiber.n_eff_hasan, conditions=dict(model="hasan")), 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( "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 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("c_to_a_factor", pulse.c_to_a_factor), # Dispersion - Rule(["wl_for_disp", "dispersion_ind"], fiber.lambda_for_envelope_dispersion), - Rule("beta2_coefficients", fiber.dispersion_coefficients), + Rule("beta2_coefficients", fiber.auto_dispersion_coefficients), Rule("beta2_arr", fiber.dispersion_from_coefficients), Rule("beta2", lambda beta2_coefficients: beta2_coefficients[0]), - Rule( - ["wl_for_disp", "beta2_arr", "wavelength_window"], - fiber.load_custom_dispersion, - priorities=[2, 2, 2], - ), + Rule("beta2", lambda beta2_arr, w0_ind: beta2_arr[w0_ind]), + Rule(["beta2_arr", "dispersion_ind"], fiber.load_custom_dispersion, priorities=[2, 2, 2]), # Operators Rule("gamma_op", operators.variable_gamma, priorities=2), 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("spectrum_factor", pulse.spectrum_factor_fullfield, priorities=-1), # Dispersion - Rule(["wl_for_disp", "dispersion_ind"], fiber.lambda_for_full_field_dispersion), Rule("frame_velocity", fiber.frame_velocity), Rule("beta2", lambda beta2_arr, w0_ind: beta2_arr[w0_ind]), # Nonlinearity diff --git a/src/scgenerator/operators.py b/src/scgenerator/operators.py index 9b2ac4e..e5e0131 100644 --- a/src/scgenerator/operators.py +++ b/src/scgenerator/operators.py @@ -67,12 +67,12 @@ class ConstantGas: pressure: float, temperature: float, ideal_gas: bool, - wl_for_disp: np.ndarray, + l: np.ndarray, ): self.gas = materials.Gas(gas_name) self.pressure_const = pressure 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.chi3_const = self.gas.chi3(temperature, pressure) @@ -102,7 +102,7 @@ class PressureGradientGas: pressure_out: float, temperature: float, ideal_gas: bool, - wl_for_disp: np.ndarray, + l: np.ndarray, length: float, ): self.p_in = pressure_in @@ -110,7 +110,7 @@ class PressureGradientGas: self.gas = materials.Gas(gas_name) self.temperature = temperature self.ideal_gas = ideal_gas - self.wl_for_disp = wl_for_disp + self.l = l self.z_max = length 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) 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: return self.gas.n2(self.temperature, self._pressure(z)) @@ -135,19 +135,19 @@ class PressureGradientGas: 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: 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 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: 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 @@ -155,7 +155,7 @@ def marcatili_adjusted_refractive_index( def vincetti_refractive_index( square_index: VariableQuantity, core_radius: float, - wl_for_disp: np.ndarray, + l: np.ndarray, wavelength: float, wall_thickness: float, tube_radius: float, @@ -165,7 +165,7 @@ def vincetti_refractive_index( ) -> VariableQuantity: def operate(z: float) -> np.ndarray: return fiber.n_eff_vincetti( - wl_for_disp, + l, wavelength, square_index(z), wall_thickness, @@ -200,7 +200,7 @@ def constant_polynomial_dispersion( def constant_direct_dispersion( - w_for_disp: np.ndarray, + w: np.ndarray, w0: np.ndarray, t_num: int, n_eff: np.ndarray, @@ -211,8 +211,8 @@ def constant_direct_dispersion( Direct dispersion for when the refractive index is known """ disp_arr = np.zeros(t_num, dtype=complex) - w0_ind = math.argclosest(w_for_disp, w0) - disp_arr[dispersion_ind] = fiber.fast_direct_dispersion(w_for_disp, w0, n_eff, w0_ind)[2:-2] + w0_ind = math.argclosest(w, w0) + 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] disp_arr[w_order] = math._polynom_extrapolation_in_place( disp_arr[w_order], left_ind, right_ind, 1 @@ -222,19 +222,19 @@ def constant_direct_dispersion( def direct_dispersion( - w_for_disp: np.ndarray, + w: np.ndarray, w0: np.ndarray, t_num: int, n_eff_op: SpecOperator, dispersion_ind: np.ndarray, ) -> VariableQuantity: 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: disp_arr[dispersion_ind] = fiber.fast_direct_dispersion( - w_for_disp, w0, n_eff_op(z), w0_ind - )[2:-2] + w[dispersion_ind], w0, n_eff_op(z), w0_ind + ) return disp_arr return operate @@ -247,13 +247,13 @@ def direct_dispersion( def constant_wave_vector( n_eff: np.ndarray, - w_for_disp: np.ndarray, + w: np.ndarray, w_num: int, dispersion_ind: np.ndarray, w_order: np.ndarray, ): 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] beta_arr[w_order] = math._polynom_extrapolation_in_place( 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: return raman_fraction * np.fft.ifft(hr_w * np.fft.fft(math.abs2(field))) diff --git a/src/scgenerator/parameter.py b/src/scgenerator/parameter.py index 14fbc56..93c716d 100644 --- a/src/scgenerator/parameter.py +++ b/src/scgenerator/parameter.py @@ -646,7 +646,6 @@ class Parameters: "t", "z_targets", "l", - "wl_for_disp", "alpha", "gamma_arr", "effective_area_arr", diff --git a/src/scgenerator/physics/fiber.py b/src/scgenerator/physics/fiber.py index c31b2d4..0d10621 100644 --- a/src/scgenerator/physics/fiber.py +++ b/src/scgenerator/physics/fiber.py @@ -1,12 +1,10 @@ import warnings -from io import BytesIO from typing import Iterable, TypeVar import numpy as np from numpy import e from numpy.fft import fft from numpy.polynomial.chebyshev import Chebyshev, cheb2poly -from scipy.interpolate import interp1d from scgenerator import io from scgenerator.io import DataFile @@ -19,56 +17,12 @@ pipi = 2 * pi T = TypeVar("T") -def lambda_for_envelope_dispersion( - l: np.ndarray, wavelength_window: tuple[float, float] -) -> 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 valid_wavelength_window(l: np.ndarray, dispersion_ind: np.ndarray) -> tuple[float, float]: + return l[dispersion_ind].min(), l[dispersion_ind].max() -def lambda_for_full_field_dispersion( - l: np.ndarray, wavelength_window: tuple[float, float] -) -> 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 dispersion_indices(l: np.ndarray, wavelength_window) -> np.ndarray: + return np.where((l >= wavelength_window[0]) & (l <= wavelength_window[1]))[0] def is_dynamic_dispersion(pressure=None): @@ -92,45 +46,23 @@ def is_dynamic_dispersion(pressure=None): return out -def gvd_from_n_eff(n_eff: np.ndarray, wl_for_disp: np.ndarray): - """ - computes the dispersion parameter D from an effective index of refraction n_eff - 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, l): + """returns the D parameter corresponding to beta2(l)""" + return -(pipi * c) / (l**2) * beta2 -def beta2_to_D(beta2, wl_for_disp): - """returns the D parameter corresponding to beta2(wl_for_disp)""" - return -(pipi * c) / (wl_for_disp**2) * beta2 +def D_to_beta2(D, l): + """returns the beta2 parameters corresponding to D(l)""" + return -(l**2) / (pipi * c) * D -def D_to_beta2(D, wl_for_disp): - """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): +def plasma_dispersion(l, number_density, simple=False): """ computes dispersion (beta2) for constant plasma Parameters ---------- - wl_for_disp : array-like + l : array-like wavelengths over which to calculate the dispersion 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) - w = units.m_rads(wl_for_disp) + w = units.m_rads(l) if simple: w_pl = number_density * e2_me_e0 return -(w_pl**2) / (c * w**2) @@ -150,16 +82,16 @@ def plasma_dispersion(wl_for_disp, number_density, simple=False): 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 Parameters ---------- - wl_for_disp : ndarray, shape (n, ) + l : ndarray, shape (n, ) wavelengths array (m) 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 inner radius of the capillary (m) 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) - 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 but adjusted at longer wavelengths Parameters ---------- - wl_for_disp : ndarray, shape (n, ) + l : ndarray, shape (n, ) wavelengths array (m) 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 inner radius of the capillary (m) 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) - 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: @@ -238,14 +170,14 @@ def capillary_spacing_hasan( 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: """ computes at which wall thickness the specified wl is resonant Parameters ---------- - wl_for_disp : np.ndarray + l : np.ndarray in m order : int 0, 1, ... @@ -259,20 +191,20 @@ def resonance_thickness( float 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 ( - wl_for_disp + l / (4 * np.sqrt(n_si_2)) * (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( - 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: 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 ( capillary_thickness / (n_si * core_radius) ** (2.303 * a / n_si) @@ -281,19 +213,19 @@ def resonance_strength( def capillary_resonance_strengths( - wl_for_disp: np.ndarray, + l: np.ndarray, core_radius: float, capillary_thickness: float, capillary_resonance_max_order: int, ) -> list[float]: 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) ] def n_eff_hasan( - wl_for_disp: np.ndarray, + l: np.ndarray, n_gas_2: np.ndarray, core_radius: float, capillary_num: int, @@ -308,10 +240,10 @@ def n_eff_hasan( Parameters ---------- - wl_for_disp : np.ndarray, shape(n,) + l : np.ndarray, shape(n,) wavelenghs array 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 radius of the core (from cented to edge of a capillary) capillary_num : int @@ -346,18 +278,18 @@ def n_eff_hasan( if capillary_nested > 0: 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"): for m, strength in enumerate(capillary_resonance_strengths): n_eff_2 += ( strength - * wl_for_disp**2 - / (alpha + wl_for_disp**2 - chi_sil * (2 * capillary_thickness / (m + 1)) ** 2) + * l**2 + / (alpha + l**2 - chi_sil * (2 * capillary_thickness / (m + 1)) ** 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 -def beta(w_for_disp: np.ndarray, n_eff: np.ndarray) -> np.ndarray: - return n_eff * w_for_disp / c +def beta(w: np.ndarray, n_eff: np.ndarray) -> np.ndarray: + return n_eff * w / c -def beta1(w_for_disp: np.ndarray, n_eff: np.ndarray) -> np.ndarray: - return np.gradient(beta(w_for_disp, n_eff), w_for_disp) +def beta1(w: np.ndarray, n_eff: np.ndarray) -> np.ndarray: + 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 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, ) """ - 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: @@ -511,7 +444,7 @@ def frame_velocity(beta1_arr: np.ndarray, w0_ind: int) -> float: def HCPCF_dispersion( - wl_for_disp, + l, gas_name: str, model="marcatili", pressure=None, @@ -523,7 +456,7 @@ def HCPCF_dispersion( Parameters ---------- - wl_for_disp : ndarray, shape (n, ) + l : ndarray, shape (n, ) wavelengths over which to calculate the dispersion gas_name : str name of the filling gas in lower case @@ -532,7 +465,7 @@ def HCPCF_dispersion( model_params : tuple to be passed on to the function in charge of computing the effective index of the fiber. 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 of the material pressure : float @@ -544,13 +477,13 @@ def HCPCF_dispersion( beta2 as function of wavelength """ - w = units.m_rads(wl_for_disp) - n_gas_2 = mat.Sellmeier.load(gas_name).n_gas_2(wl_for_disp, temperature, pressure) + w = units.m_rads(l) + n_gas_2 = mat.Sellmeier.load(gas_name).n_gas_2(l, temperature, pressure) n_eff_func = dict( marcatili=n_eff_marcatili, marcatili_adjusted=n_eff_marcatili_adjusted, hasan=n_eff_hasan )[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) @@ -558,23 +491,27 @@ def HCPCF_dispersion( def gamma_parameter(n2: float, w0: float, effective_area: T) -> T: if isinstance(effective_area, np.ndarray): effective_area_term = np.sqrt(effective_area * effective_area[0]) - else: - effective_area_term = effective_area + return np.divide( + 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: 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 Parameters ---------- - wl_for_disp : np.ndarray, shape (n,) + l : np.ndarray, shape (n,) wavelengths over which to calculate the dispersion pcf_pitch : float 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) pi2a = pipi * r_eff - ratio_l = wl_for_disp / pcf_pitch + ratio_l = l / pcf_pitch A, B = saitoh_paramters(pcf_pitch_ratio) 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)) - n_FSM2 = 1.45**2 - (wl_for_disp * V / (pi2a)) ** 2 - n_eff2 = (wl_for_disp * W / (pi2a)) ** 2 + n_FSM2 + n_FSM2 = 1.45**2 - (l * V / (pi2a)) ** 2 + n_eff2 = (l * W / (pi2a)) ** 2 + n_FSM2 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) @@ -675,14 +612,21 @@ def load_custom_effective_area(effective_area_file: DataFile, l: np.ndarray) -> wl, effective_area = effective_area_file.load_arrays( ("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]: - wl_for_disp, D = dispersion_file.load_arrays(("wavelength", "wl"), "dispersion") - interp_range = (np.min(wl_for_disp), np.max(wl_for_disp)) - beta2 = D_to_beta2(D, wl_for_disp) - return wl_for_disp, beta2, interp_range +def load_custom_dispersion( + l: np.ndarray, dispersion_file: DataFile +) -> tuple[np.ndarray, np.ndarray]: + """loads a custom dispersion and returns beta2(l) where it is specified, as well as indices""" + 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: @@ -702,13 +646,24 @@ def load_custom_loss(l: np.ndarray, loss_file: DataFile) -> np.ndarray: loss in 1/m units """ 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( - w_for_disp: np.ndarray, + w_c: np.ndarray, beta2_arr: np.ndarray, - w0: float, interpolation_degree: int, ): """ @@ -716,10 +671,10 @@ def dispersion_coefficients( Parameters ---------- - wl_for_disp : 1D array - wavelength - beta2 : 1D array - beta2 as function of wl_for_disp + w : np.ndarray, shape (n,) + angular frequencies array + beta2 : np.ndarray, shape (n,) + beta2 as function of w w0 : float pump angular frequency interpolation_degree : int @@ -734,10 +689,7 @@ def dispersion_coefficients( # we get the beta2 Taylor coeffiecients by making a fit around w0 if interpolation_degree < 2: 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) poly_coef = cheb2poly(fit.convert().coef) 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") 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: print("invalid raman response function, aborting") 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]) -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 Parameters ---------- - wl_for_disp : ndarray, shape (n, ) + l : ndarray, shape (n, ) array of wl over which to calculate the effective core radius core_radius : float 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, ) """ - 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""" - 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: @@ -952,15 +904,15 @@ def capillary_loss(wl: np.ndarray, he_mode: tuple[int, int], core_radius: float) def safe_capillary_loss( - wl_for_disp: np.ndarray, + l: np.ndarray, dispersion_ind: np.ndarray, w_num: int, core_radius: float, he_mode: tuple[int, int], ) -> 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[dispersion_ind] = alpha[2:-2] + loss_arr[dispersion_ind] = alpha return loss_arr @@ -1128,7 +1080,7 @@ def n_eff_correction_vincetti( def n_eff_vincetti( - wl_for_disp: np.ndarray, + l: np.ndarray, wavelength: float, n_gas_2: np.ndarray, thickness: float, @@ -1141,7 +1093,7 @@ def n_eff_vincetti( """ Parameters ---------- - wl_for_disp : ndarray + l : ndarray wavelength (m) array over which to compute the refractive index wavelength : float center wavelength / pump wavelength @@ -1182,18 +1134,16 @@ def n_eff_vincetti( """ 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) - r_co_eff = effective_core_radius_vincetti(wl_for_disp, f, tube_radius, gap, n_tubes) + f = normalized_frequency_vincetti(l, thickness, n_clad_2, n_gas_2) + 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) - d_n_eff = n_eff_correction_vincetti( - wl_for_disp, f, thickness / tube_radius, r_co, n_clad_0, n_terms - ) + d_n_eff = n_eff_correction_vincetti(l, f, thickness / tube_radius, r_co, n_clad_0, n_terms) n_gas = np.sqrt(n_gas_2) # 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 diff --git a/src/scgenerator/physics/materials.py b/src/scgenerator/physics/materials.py index 2ebf4f9..1def4c7 100644 --- a/src/scgenerator/physics/materials.py +++ b/src/scgenerator/physics/materials.py @@ -321,9 +321,9 @@ class Gas: 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""" - 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: diff --git a/src/scgenerator/physics/units.py b/src/scgenerator/physics/units.py index 612c759..dba43cc 100644 --- a/src/scgenerator/physics/units.py +++ b/src/scgenerator/physics/units.py @@ -461,3 +461,12 @@ def sort_axis( out_ax = plt_range.unit.inv(cropped) 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 diff --git a/tests/data/PM1050NEG_aeff_slow.npz b/tests/data/PM1050NEG_aeff_slow.npz new file mode 100644 index 0000000000000000000000000000000000000000..733af7bcd4d013d8b221f5edc6061202bf9fcee9 GIT binary patch literal 2310 zcmd6pYcQN?7RTK$iMS<0GMhn2(C&hf73a*@tcakQkde$*h|nSvE7~=#BQA@ewRwdNf`&wmLx=}O6`)04TzzHQo_zhg;<1z$4N`6ORdVj?R~cvkxzPfotBo0 zkxDQMqC^Hp7&+(}IR+mxGSf2(z7!D^5fJ8oDI$onb$&D;G?KD;9vKo4PTAZWTU%S3 zk<9dx^#14CMqB8pel<{yJv^NHVQ&pCQg3PZN7SI6t(SMxcQyFg$^?zB!@^W?z)(yS z3x)CX@@IxvIO7=UbVRQfKal)8krY#lN54rItNpzeW7BrLc)70*InJe)Q?Ye8p{RSY zZ=?>DhZ@)%gL-txwfdPBSC1FTT^~h{)?=mmA5Q=202@cfhl-sN*tk*bJ!W)=jh4Ps z+f|Jk@a!In%ZtPYESOwK{eHXw`ij%#j1_C!)V=tW4$u5`!icmIxR=A?S~dT z=)CLM(!o{?RAiDKX0>AR-~q~wN3Hm6PHT0oEe}^b?5DD?^3YJdS3CSM56v#rkxYNb z$3G5*xab_^;~7(PiUyUB(Yed1a+Q31%Y%<>a(vJHGN0hmh7Z+c*T!7h2tIA7bT`MY zli5aaYQx``?4mM%Z6kOI@JnANH_cT*a1&s~9>ci{RRV&a0EZ`A4ZL0o2#!KrXWILJ z>?S053eo&*FJo_wkl-prKYiB+GOvUL-*(KOxufxmdpp6o9ededc#N{z3Eu6v95S8N zvD!{>??8Ri(EbW?2f@Dsy%%0cuGDo99&}(jFX?c|S_k2S2s@|>v1Bh1;e`nQX_GIo zWs3+uM5wi*{p`ba5#dQEx`?Yu(kD6zUpnz}Jc%;f&`EgHiH(ke>XA2{gg;{Ze#2o( z=q)BZ665x7ZfY|c#e`2{bWC&iPLYTSue$Kh&i84-C%Xv0y70Q9-;jHA7vWhKjt)MN zw~_89eCtL@&LyWkzTJd(-DpK#O_XrE3IDp$>%H7i){QoOmm8X!xj_jIrUXV^h@tog zL{I|$U$aA3c|3Bp>?)KQ4cJBZUxhQ=lRphEUxmL$Wjd_t=71aLa}^oS99R?!jDJba z0bXw?FNmK5Ix4OyeLv@b&&NvHa$30{xqo7P`gktbCm&hnrsM)cwoC0paV`jUI9SIl z06wf$TH4C=i! zE-oy9ESm?bYyAZf>2uBM@@fG*{qz&VSo1<~mLgld2`mJYxx#+6%0jqSX>_^$ej&`O zQ6)?2MKGAEYccFm1a;Gn-1gieNE%2xuSY7e&uNw;!_OYF}}PcQ&9}^ zPpkUwP8Wk7gJB(QQ3AckXz8N~C7{JRb$I)52`Ko{%~^X(p>^SFfd;h{^1eFn`utWY zFft_mKUW3FtS{bF=%ONpE5@Z=$ z4j?LU+4H0v!n{+415YtPXCp>|d7A<9gKB=Kk5xds6)Cgv+X@ITuVdA?RBry)V%)=# zO1R)~lQHGXg!Rdzf{`ao$Wr?13dUcD1O@NgBbrs9P1(>a<5$7btg6E#sTxjX`kvjN zTn*w>r7pYs)o{=8Ko7;T27Kr5K1<52f%IX=McsuOa1Gt|Nb55eO!#}YMT=SBr8`8b zFsOz7O`dFKVJ)~k&eSZHuLFZI{sZfzI;hggXkA&YgY0R}M{aTTu(SKvw#toqXufA* ztDMaS2lF&@;qC??XDgbkwKl*N&-COnR}PG4ua4tiwCZK zB*s`cA86;-=yOYacrzUIMiA2mqKmsFeKG>r=k5}AyhH#l&y!e*7DDJ|RPMAM6~da+ z*o=K_JG}0>y%=fK0bXv>GUl@HPaRwH#MY}y%pzM$iuBg^{Tbvx?`|c}yStiO?L?(W bACi;$FFEApCa zWS)3PaLI8+c%AeO@(uP2_g@}-ERvf`itDB3zpwqTc*}X0lZ!hy*GaC^x<0;PhmPx7 zXz8vyY_4mlrF%H!c=&OzV9$`_KED51zR4>n%y(gVn7`LC--TRn<%$(64fPGR&S?GL zK17oG7?~5#GhvFWI6vuKCiuubAT4Rkgm=2h)Q3Ma!P4>)+qEkb@|~kl_fRHCdk5E? zpU#9cC%kKwc(b50-mSh|Gz(nsIW%mN%K}sm&EM6sVA)fX*6<}+@K3|iPSNFAz}2nN zpKX=}L30wLYpt{3ixls4y?qvJSUfTB=aL1bYkK(S_hfOj{cRaS#S@h8EgHM1wtPaOwUWL*J0r zEO>KXD^Q>(3)cQ47kYmr3vvsEB3934f#k#4=sLb^IQ?!g)>}Lq+FRS=r{uF?`*>s0 z4UKGgDEuZxS2r8<)E}oK<7~LKf|0q)JRAPl7UlP@&xRwr?~^WU$%fBIZcCzTsCC$q#4Hu zLRaR%(ui%YLKZo2Gr`t9 zV?z$i-e2QswlxR*DF(i8Tyx-Kjb@4T(fTKI;GE20 zY;{Zy^ekzM-+eg;cCT$r>c5@?PaI#T#NEyTql4AyYFRmOCybs+E6M@B%Z2$4)Eo#& zzeoC2nFCFwH>jtc=D?=M7uk}pbAa;c6w3dU18SW?HCEqqAb#dx?I%tSj7aRN5A4r@ z1KL{}cqVe-m8n(po%tLv-LkTEg-|YJ?pxCNLMj(Tk16*bP|k&j3*w_=i*lhgm3#W? zvRrT~9-G%$l?&`gUHny+xv=>CchMaibK!d1d+9F6T$rACrg&~wE*uuER98HZ3-2|` zbjl9r!a9>|eY=ocD6~s8ZjQ(WS@#6fuyeT(6BJ`9awQkK&K|ePx|s{RZXLEaPtS$w z{Jl=^vU92%v znLK#Fq-9F0=7Df+L4Ki59-R1omt<|22R{d|QyW(2!B)PD?4wqBz*0Pg_%`K%wn0$M zUFSTw1_x_RcIUz5_MP>wJoCWEcXPu*|2%kiYF+chu{^N2x}r5HG7s`I7I!X*%>!wg zQvbt*JcxcSHtKpS4?4bbO?ThT1J}OM`SZDXfVjH&l}qzLU+%jonUM$o=)RY>M|r?y z`Ao6pc^(8hSE!$Omj?}Ar8;7bd9dL`mVVC9JSdI-$JnAP50vg*F?~Oj2N%oFSo%-r z!N8MHn|Ypm*xTT1pDLOU&$`^5jO6lR<(!jijaok3liJ|!u_PY^7n^ww8|6dTnq|IM z%<`edK_gJhIv=(ikO^hk=K~{DAi~)tAGG48qd9x>A@O#9Y_wNCjF+^=%O1%GZ`_bn zd^{iCe5g&aIh_yY?I^wRLO$e9Q8PoX=7XeoLB7C0`4FXbmvsMrK5$lDr>@S=hwb(k z*mY(3@NjPg@@D0OZg61D6wZfR=RIq0yv&C`DX#Up@ADy`(7pkE&4*7_md(3b^I`ov zAWzK4+`Uo{VKEh5GyP;O5rPjJ~giC&Ef^HXT|9J5BUP9vF+rK&?o@o z-OZvBx&@Gaq+UADxB&QLY82O*7eMGuw))5Q1<;&J(mArF0PNYB`rJDTfL51ce0y&J zX#9vbH9k}TSB9f4p9dDeD1WHU{_p}gsO)P$8eIT&hVD*@aRp$u-pN%vsQ|KFZQLtT z3qah@)YB!S08T~e`*su*!0&77fidI)*p?+7Dqm3mm300H(&GYH@^UKLwzdFne(Q^E z`cMF~{lDTvzZHNVPeYPWTLFBKuT9D5EdbC*>1JaEP_&Mk`DU&FtKZ&aK3c-kSq-OWZLbzRXptgT)A@DZtsE^xF2q8V2 z8q~HHLen3MW}0gu*vS~TI_xV1$`YN?k@ad9o0OJ?E$LV+!Fw z*vP!qz zJPsV3*!m2y!f76N}eEFt^sE0-H)5fA%`dJZd-L<^6@J$gg{k1!-8;W32v_gMF zOA%bVE;4$wqX;Im|4i`>7QrFr@ciA$B6#e#>~LRWa!MUN-HqD~6O)r!8Zhi(x)7 z#722{F$88DvL|~MLj!fUlf8d2*gW6r+H$NIO21mWpNK35rM^|3VzI?=flJRfC!rYn zWz_;LZWY5`U8&IbcZ=bf1z&`JZZWKIo{XL+6~jHR-q=(|F$kXc8E=G&;duO~q?+f& z&~o=>ipRTR*i!y5efUc;(4UYquly_q&4%23t*&B7?7B^24Hd)q>{Y7sbTN2Ioo91+ zO5lyo2^1|_0_JN1YGma~Aje^UZLwMjNF3N+Z?mKXqC)K&8jVVTb8&5Rh*=4^+%{?z zur7gzC5t-m+n0bYmhWHfQUb{zgh%W4l)xNkZpz!M1Oldp=BENm;FDN8|IOniV6D|8 zs&~2sN>;s<#urLJ!TyP&+tm^{zn7`rmr??KM~ihX-Y@gt?IqrqK=YWJQ{LASuovFyx~{bZsA|^kAG=FH zeZ?xzBO@hn#a7psd$t5d-Bbf_^O3;wh-9d-I0@>`@J2kBCxO|`iRk?rB*@C`i5*=^ z0x|ZF_(Wq8oT~enq-{=u-z_gvD%O+0dAKUwWeW)^`N^3b+exrQIX6FMF9~iK-X_T( zBEgLHRVpcv1b(jP*tX#$_;5HJHAR!aD$2hmG>!yC*Y?#4C6Pcb%cVXel>}$$8yn0r zNYL}by!lN53EaLJw)&7s@VH;AbB0X6d*VzB^hzP)PoyQRECmsnV4KftOX1`aZ~LGPrSQw#&53twDL6Vh zxTd?70^1YZSMDnXoiG#6mp-MCbZM!t=h0G_zM~R2exek7OC>_D#gsyQHBZFi%cWrX zX*{~>dMOlibjR+vT?#TYE%9BMrEo^#W74^zQs`RrB1Mr}3Oh}!(#tAK0dFqPw0l|# z%l75uH`kRyYDgL>>{BUlUr3~id@qHd)U)g?PAN1Ng(36)Qm}pGSMzSX6w2y7Y7fts zf^yro`Z>Weh?}r&NRcXoK@sz21LZRC&@gO$vZxGRtkmjsUseVtc5?lLCS`DckI?94 z%Q6rOnw`?vSO(!|2j>}%W$+`pjo)ck893xOini@917^i*>C=bHV9}duMVXK?xcZZ> zUKCLVlOu&X>(7&xE zJu=)J!)%?m4l}4wtA1qa!9bSZ5TdZ4&zR1ny)03gO`^5@e{}Mb=4?A;U?};`|>9WH1cK)i+Qj!>{8HNGCPPkQRO0_=ye~ zoGv9(#r4R*bK|P1yCE6acP_DWSCAn%=bYu>YBK1QMj^|!WcbPqx4CRZhNSA?nh!Q) z*jVdtud#^?(;vNS102Xu{(YYl!KFGs>R~1&da-J_{nlegh+4=}QF4- zJu(=i_VBx8lHq3tC;o0O8B&XWh;|f`!IAbg$)uDF+>bs=$56t86! zSCT=y>2dl&M25!SRqCWC3!nRMeNG13iEr;sOp8{Kd zY6b2#p%4y00p^HYXuk!8@Bj)d5)_Vzv!xI&K!Iy2bJ1#!6v786FljUtOWR2yoPYv{ zY&h`_`zeGMP~h#(uSvgrDTEtPz{3Ac%Bf%q;Rh7RJN-CaGJ--l0tKY5GcxnfQV36= zKy-FdzEwPha0Lo?^9tUzFRCbnLr@@anOpsV8VcbN6lkz=Xc&7Sqe! z6BHo%n6&D&Pza}>fYPa@ot5nr!Ye3nAyK7&`v8S-3kvi9br^ z!ZE1u?AiFdf-sfv3@WVn+Ra}kO(k4|3iooPxkgJV>}OQ>M}#KY;yDJtP6RLJ?f z&DH%pmGBcPNOak{4_=`XjzWc~*)^V*lc|KKP=O<9=&Ny;O1KIYTy(So897wKSEx{B zCKu{Rq7u$R1ziWBh~IQ7;Vo21-ai{1^@vKi3l-*$4aQ18rxN}`g@B7~@r7@xgu_tb zQ(9w^^%pAPF;uWFd7aYmgG#sz6^c=H`q54*;WJcF_(0F(8=?|ULxuC4!u-2aRKjbh z&^vjLWWqxu+=d2w#BNewiO>kYp}|wlIQBtV8sRuJSiUL}O{mfc&!Iv3rlU1Ui)n=G z(13rhSM3rb8sR%MICgYb{Xc5X|9E6@n}c z_zw-#tCn;s?V%A4M1w1DmHW#N(g+Wt!RRmX(M1U+e;n@FE&y8Mlb$B+>{sqJfz02Wg8G8sSGYIOX&SX@o1$pyEchekzkj_!144|j2 zb$6bj5k5tOv({@oIlOejsp!zX!@xIMj81qJ9o!CU2Fl9Q3Adud<4D=iVs$#j81qM z9rka#LDh|>6D~%Fm);lI=pvo)F*>Y1c?#{iMkkz%4jEU1YI^^n6JACKk^2X0FFc?V zZbpZb4bi&=}P``C-{zDI)@HaYGdUo+2 z8Kn~rM~8yr-$l9S=!C~HK<3hW>D&Aa!sQrn=FT%kV+jV~a}4MtRj5B#U=U8nfSuK) zI{P&lgx4_uf6US!)nyQF$AD!W{}?B(U=V)CfYj+LrrK*6gyS)QTO!7?!iGV39s`0F z9k+3DU=XgyfG?{L+js0>5WdF%+s%8OVmugx^D&^z!`W5dhe3EB1C&EHxRZhyg!?fd zF4oM`_5_3QKL!l`v&^^Y41;h$26z-{1cqK>5FW^Y=MQB1AeG~Nt&r(5Prx2hn00H$xj%B zBQk)w5vS|dF$hm&z@k00%<7K}!W9{CHLxImcQb?VMFvd7+$Ht5F$ia5z@g;p)VMwd z;f)M58aS&t@4#@z%Orf02_#0d zsAn*fa84#DzN(iVKglG#lL;4^YZR}YWfJbmg#H1x`r>#d;h#)!=PlK#y1^tIlnFHo zS^7KfFbNN3!ir`87QxVXH^C&Fl?kusr9%&Mu?TNv!5UeCh&dq^;jS#mSvnn^BF!TF zl?4(O{jmlrEW%+~5b4w!|3sTbcq|Lr4mKpY8?XqMWr0h0ZOY&(7U8ojsJe{OFRx<} zPRoL&cd3~gb}YhcS&&?opU-e&5pK(ZxyR`wM>iJXw=D4goJ9S7fJHbi3qE$nvZMT1 zgy*urdiEreKE@(kmj%U=fi;DvScLDgKwihQ*7`h)a9$RiH*>9TxWXd5mj%6B>>G|I zvk3QP!S4N*&3t!Rg#WVO=`rKhyE!bvfmyKpqE4p?iA8uY3+|*T_P?UD2p47nUy11G z!AC5@hgon8%}-4{V-Zfwf^Y9f=9AvC2rp*ACQb+c(l0E+jafjQ{3iP72aE7y7O0EW zOYiJt5su7)c+DEcoMc9Nhv%&Oe zim|;coA72fWSxsQZBbwhrZ4tNBAf7S zHk`fjD_$#wO*l6jx^q4!vF@=6?`Fd;)~gifJT~FpYkv_HsRrH;Nedv6~AB;F3yHvr6j6NJ)7`xHheXRWj8jm2`6X6M(dL(ZXP|-U`CgD`0=7MYFe9 z1>x!y@RGK?bxN*+@bwB<^+LPzhI$3z>=lsFq|mRsw1V*V3UK)!eE^FRlehi<`v80b z!jAbK4+{th{(rR(Kwfg~i`qynoSr=tQ*}%WYnm0Ck^Ho9?kn&9dk3_z{(jZ`YrC{? z>EOrxF^*c;;mA!#~YWObcHd=J@mQYT?Z3XWF%FEl|YGS&h!LrL$ns}}E zk~a?vH8GCRJU0A~Cf;~WSyu0YCa%qvW$r(wiMO%E-EZyJ#4ZnS7^^#I;?M6BKI~Yd ziT8F~JQ2TG6YuL0bd-_S#J{+B2JHW6;8wo4XYac;u%F6Vf5Ao#?61LNeX~Xb2WZF1 zt5P&@|C+OGmj@bnz>?>BQ-TH_w2t%jiO|49wnFRr4{6|G`$U;jE(>w1L>IS}1|HpZ z6DR5~L|3^~c?AtTzB?@_gJ&UnsMzWcsN*o_rA~t17UGT~otrPz@d-D>X%)IU9^P!u zHlbYYPc1}ShpKKLb=+&}_UxsLI__Wl_3jQUbsS*G$y=eXjssWt zM_0yu{p)JqU@cut4ryo+(aMR$@xb*vK zc=vCU>!+@$;g6r;j@t<}?DEnfXSJ6aeqHHSD&wSvx0HIb2F=y*vy9-!Zx^fK4are; zBq=qFFJAhbFsq7LzN)g1I#luMJ-X_5KdRz#hgG_fSQQ&uS+Cq(qKX%5d-y|H$|s#u}zTy=!IDo$>=TD!|m6-&Ol{b|)oRV-L+By(6z6<@r!R&Aq@ zDvmz$*8Sv&3Z4pSI<)hL3LbI^$+}dlf{)ssF7aWiVBf`3U8z|rxK%-IB>b8RZkiyk zC^)5pKQ!6}e?6pvx7BSw#oDHV?Mdy7E{lcrMkOT&moDTl4uw6JR>8D`YL2Y2eOASv>mIo>0 z-8p&5e0!Afy`4J@-EEa|EC2kn(&frnk5%HyqoRz1{r60{@h|Lm<&;nQpc4Mrn)j#w zn-U%_*>uF@g%Vz~e}wlkRSCN-O$n*JuY~V^Sud=4MF~I5{1xsUu7sNcVt?j&D&Z*; z9mg&QC9E|2lx1k9ggr}-Z9l(A37IqFM;wmB6l&!xNans-LKSiu{ zA@k1f3Po(~Xv43Yr-%>8w;nrkLlK{8iES=Dt%#pqH$KVZr-)nk)c>;Ep@_M)L!GWz zDdJ5Z)jPiHDPrGyXvZQsMND!E>^(8Bfa?^*_EdB$V5yhvA@K$UygSl2(fhFiKCwYy zB9Ekit463drqdKKcb@06Z5I@BE5dzn7Pa`;*3so}YQj+1QnkymTpC9`SXxiJD~$!Zt7-T1av8gDuA^|0?!Y3vZW&uTzM8gD)GLpE|=3Lg&$R+8(L!s)x_ z)XE#BaF_n=Z6430u=yX8SV5{3F0K69TJt~(Ye$^3jZKuoSasD0ms3*MqndM2(^m>l z_}m`o+b)GGhulbSKnmv`UX+Fmq;S%=p7Z5OQuve$Gboc^3d?vXi`sD{@hbV}>$YG? zoRg7I8*)by%PiZrCpAnGzq@F#;fIqXzWPJZP)|=1r;5ym-r|wOS*C-mKg|+&YQ?Pk zRh9(KoYma3E?ELkOc+eR3X;G*xBPbAvzNdv#yMrd+7kFx6OKPlZ4CRD^M3 z@}~CnlS23fC-0JegAitz^Xtmch47F^w2peR5RTd=8trgQ2=m&ee~NGx!rugJuDmc5 z!nE2BD@9Qu{N+Ht^OY_^?C>ERyQf_QCg{q{(1LA+v8J+b1a01k%>T58n-xZ+|bBPK%tZ@ueVz4MZ2W9u$58WHsaMza)cWlc3)v}2X3npJQ_tWFUjbhtYy%FZac6wXd zBYJo-&DQpg_8VS&)Oc6nKp8JS*%KrC;|4FTi>PiIJkE=?`FOH;=y0mYgPRT z;K7b(g#@{sc<_M(TB4g*@Zbdafn#( zzMom#m~U>@@M9b|zT1~L_bh-L$CQ~Y`{u-r%dLl~GA7*EDCyvBv88msmu9_PZm-nac|ySebryt(ji z>$vdxE6Em{wYjk2QC%*EAQyhUq0lzIX&yCgU1{BuH;=k^R$4fR&LiW$^Q-enDxtDA zPg$}bwMbC)prZS6uSddrnch8`2Ny(a4l{2WJHC*CE!VC)fa)~jtcLqHOdMQUU zm_dQrYYi)Ar_rZVO*6c2r%~fT@r7F%(O8|P+=^D;W~!Sa3$EK8I7Us-5RQqf@8=-Ue>Ga*C^ut+Xp#{3P^kRZ%i9S7c0Kk zojE>=4n99-a~J+4eQpH~tzy=YF1Qx$$TO`AT(j7Nw7% z*z=Y9*F}t=2g3=U-FA+kU3`iB%#|akkXhO1E;)iyg;kO^^bI4exx)64x?yB@({JKe z(JE>ehuI_oc%~=|{J`j=(VbW%eRsW&04aGyZt< z;`4?0Qo={AUhA+6MsnMY%LN0(c`|^WsI@)QT{nQz+r+nX3l5;4%f}z~ zedtF^B3tyD9`vI|ySNiy1N%{;jNy|;^L|v~b#tIwpdZzd$JORO^%3VuAL?+0wf13s z=y!pZ-dWo|)H1`P$dc_tn?B^Ub8~u8dErK;BeNI9apxW`j_*YtdwCUwJbF=v>7v3@ z%X?9)l=a#!?p~A}tXdM((1UDWPpAzR_aKA7Rr9lFd(gS7Mo-k-deG-m-MDpzJt%Ga zPrVS{9^^iD`*l`RH*#$v|J+aQMin-fKN}`?Bh8}ywKIph(Ma-e#&@f3^dkE539Q(S zJ`|j;S~lE;9RBvBb)hFbzl3y)x=?}14k>nA7rN^GRB?-U7b^OA>-r($ zANHVc;wl_>#VI}s=PRDX!)*06xq>+T+1d_s&==b z8V}dq=6!7_&_SZFWv~tHzs6q8J=TT}^kfCRoN7Z3W;6Tl|7kJzZiuWt=%^B~Og%yFDCK`geZjzyI$4bsYZdxVY6_VMPAJL2XCOG(#4STiDG%GPIcSym-qN189OQRwvHzJC4yu=PzID5wgN{u%EVr8HpiPUgMW9qWa=jwJ z($;QA2Zx+iHLYq#-@dut$lcVA+@8Dl1-rK+Sfyd$c(ffY&2078Ki`h*oCOqb{Ik#; z{nz)~(ZBu9?WpSq>FMm}c0}6V|75M>IWJ-huWByZ61W?LZ>WM}FikoNv1PPPFUX zooM66{7X^loha>ZA8;pHbm?SFuUjV~*(aHFhjpTyM4uPaH#^bYcNY6Z$(@KR_1*jB z?>dps-+tVM+-Hx`ym%M77w9NUGU`HkV}HH|I(MPAYAI<(NEcdiReSTF8(k><_N{?> zY8Tqven`;ha~Co=aqGIwWEa|GryYM@wHtZKQ_c-SH=-@g=<4+DMi!C1PmWyfMpI&U z#7>gBQU9;p=OQ1vk;gUi3kbX?4l07sa$!MIAoai#!Du>D4FnBAca`KCfl; zqP>?LMht)Uq9Xm^^F^Y4$XI;nZ26i#L>~-`VI1m1q4I8h8*lU>pVMit_dn`G3r(16 zTs?hg;BP-)KXMsab8okOKk86Rdo6adA2BZPI=VHlA6*=ll}Py1kMNHh*T3-%pqba= z(nh8O$gDCUImUMY{m_3t+L}6mF5ltaZU1rr`TgxL8$i|?Tg}`Q29bl}e=4Oxbai`f zm)n*>w5rF&_KWu*l6t;R{8#iK>dW;1(w#bp4y)z5ucZv4D09h_*L8zPZ>;E2Zu=l= zq->0g<{m=Yl9G$>s12dups#NkYlcw22yeHC>kv9tWXR7GGK7Zg-Lwi5hLHH>@{#2H zA><|Kzee!s5K`_Kd}8)v2nh*&MFTTKXhQFCWuekA;%q5mdYKL*QPQAx?v7!!Lc;ud zd&n?yJ9c(P@7ge;$}VjjEgnXatk?AFm&54D-msfj+lSE~h-~@AH-c{JoHbaZGlF!4 z;!a()89_U)Ri0UTa0IE+G%fkhE{sni>14eJBdE=K%TmXMaSA72d?&qm1htEQIIlK4 zg3il7_j{!}iqzFBolaVgq5)s=T$0BqS}L(}!-LbKC}o1?qI~ba$49}n6-On%kD`d3 zL95@*jS|<#2uMNi0*Pts5J?C+BsnY$;#be{>U7`E^CAY>AbnjI@7Bib!7^v)lW7{ zluse)JifG?Pg5vcWxep3=_$18#_7{;TGQw$KdaAp(=-}X)mCy2nntV-52wk=(&fn~3bQC^t*Yuam_=E~!&_>6XHlr*ocq~pv&8v1i=??E zw~sc@qSmizugm!6kkd{Vt}??pv@^uKcWCDvs_q)p9EzSp&-uQJn&!_OxmDzg2o?2-yzf{WB&&Ymu*}BGvp647EQZXdV8TB@@tuldi)2KJYBMD zbN3%~vHD)XUfFrH|B-A=(7Jize4R(dyH|RCPMSxWV_ysm+4JaDa?E7xuX)6GJGnK8 zp9_nz=5j4mxUg)4zn_pX7ginnsNZ1Ag>~etdvo?I)O);7L{g4&;kUyZ^)JS9;l-^l z?8DQzaO8W3!#<>i<<)O4?tILJkDv3OKit5D<=68@MfPxE_Wb?PL@sXZfz3l6$a7=I zZ9R4$mvQ4kiA&WSYi|7Wfq9zBE^hpOdC$>pN4W75bbsT;v)q_V(?8-{3O61uQ*o9m z=Ef@|Ip_Q^H(r#Hw)@^EZtOYll`zu9jURX#K604n#?5$qAWM!12MjjGs4nBd1ui+1 z)8ZKhh>VkIwSo#0_#c&HmxRieclL0~>*2whP0xCya`EC%zFysu3cUDWs%_*I172*AX63WQh8I77p*ixv zjTdLWRBo{d;>AW&GePw8ytqQYWyUCt7vEae*#3&di#0qdPuM=;#giWNgFn9T;*a~v zww~_e#Rsn4Xd30=!{$i|&u%F3;b$mlr^s?XynVp=_XS%%Y%sX_t)@F4W@#DhzZ1-d zZ!TK0*>u4jG411zjyx;m7}u%lPrR1|{dycKldWaHe^5FF$thdnW04j2}0o6+Czr$B*SFL*J|2 z<;Sjj*71dt`SHKk6MkH;^gclQ8$Xtayn8QlkRQ(rIfs?;3t$LP%I8-Vz{XRx(H1KN zu-26`t|4{;SXz0-vfRA_nBG{M%M1~~dH!ck_goafoAs2&=k5sLo(}3Usd52)y}+qQ zyG8(;_;vc4eigu%-2yKz=@-BT3Ig+Mcm?q}{wwNwN`g3k@#v2+LqYsqHsTM*Mi6sM zTB~Vpg1B$H(}nYaf_Us0+iu%AL0pbSNz*BU_=tq3On#9d*4a-zwH*oKmJT&zm@!1EQGx-%)9%p6~ZgML`g+XLfA#pc(0VV z5Z3OqAniOUg!87CXP-8f04IdU5MfbnjU8eW3#nS3*>3SxWh;9M8I=lytVp7d-_*lT(0_b_v>C^ zEE^&9YHohP^;RyA?UoV2Lk|z7O)U|@l@aG%+Ra7qYVXUXsCT`S58` zolp@hd_F*|?t%#JthJ`Nq>A9Ejen+Sg(6rKJ=~yJC4%p-Ih`2wM&$qV$Hoe>P}7hI zKKDe|s$sz;_3IvlevlW%7kK|{A?b?ZAC~*IFRc~Dy*arLthS2cY0+Mlv3;VrB1hFS zAy^duTNjFA7w=VmP7CfC=DzCP{cKU}a8C7G9ZeK3D&Xyysu9Kej2>G2{w#_)7p`{q zw~OKrYmHhN6Qa0!;LX=)Au%k&Z>>D9B!;QoFX%b?Vwh#T=r*r~7_JJBI(^7N3?J{R z;>p@8hQD3o7(F~9hR3}oX6aE2+nEW!E=>@_$pJmBUg-<@{Wnwc#S3|Txm0b17>+pP z9kTY77~VH7_qMQ64Ci;Ig@5c2!vjq!s(;4CaQAGSuYiC!zT^DWwO39YZ^|^>@Igl$ zyKi?n6=ouigLuNrGwRxsm?w^&m%L08YrVX&O>mVY&OTygv&K>qH>~^6u+df$hXk*oi*J?0 zRqrmE(sxMWii4GV&G$;;fJKcu>0XjJ#^(8ImcJynyLqpR5+aGQl*fjIlal!UKv%2! z8A)vP=b-NO3zFDnD>si3@Bh3bz<>VN|3AMH#C-z(6V7q}=lOr$AK-uA_CH+&|L1Lm dw79tbXHS94R(^s1?2%{T=hQ-!TK}Kye*y9U;-CNk literal 0 HcmV?d00001 diff --git a/tests/data/PM1050NEG_loss.npz b/tests/data/PM1050NEG_loss.npz new file mode 100644 index 0000000000000000000000000000000000000000..3f49b267bec149003350ac04ce5134193d90f546 GIT binary patch literal 12866 zcmeI33sj9;+xI&VkxCJzb|p#<$(iVXM-maGh$uxwiO!Lv6qQaoDMczB)_JY_Uh7_~ zgvcp{h>#*mDI!AT)Mt;a{fuXi&)CoRyzle9&-=c6uCdl{%(<@nnrmM7J?5I%e~xLh zRJPwFso#pO)Qrl^uHGL-QEHe}kZa)P?VESF`)@Vc;k{d0N?9srcM9UFG~x^4coeUa;Szs=v<{kFP#Z~k5% zH*40^iG~w(_Uio8Wzha&*UCWuJOo7ydG}s7AD+ugv;wLN5G(assC2Xlb|*&BOLB{# z2}^MVvj}3>vz=$iBJ7oV{`x~n5f+^o-ZiDX2uhg$*0Zt*_ZqgpE2t?#&Yq0-jSq^j zRpVl}dUFw`lyrWudr^c>4x@VF-xi^~d(o#WpNkNi=KEQ?9|LyNv%kz$U_kTs)vp1= z7-;c%|BbK707qF@>iHN30{F|M$4p|tbonkB=NSyhy*w$KVah;tY<0gH3kK51f9|ia zoB@}s`U4iNW?@_^8F(WXDt~GV1E&g@gWm06VE@9}!A3z0ST#!x*}9JbmHiVH z^P(82*IA>~kjOxOS%k7$1_Rrj`9p1T7?|<(fyyBU17EY`hFugeaM@()@NdNoB-S~p z%|6e7!>*VSzE>I089|Ta++v`OeyrYF$AG|IVU$i21H0eO(sX>z!2IM*TB%(O44RVo z%Z*PA)Ks6+9w^I1){a)41%sJz8=^8gM3sq&98*2%T_VrotHWBpc4~w=3kz4ZZi|o zkKY)5@?xSgRLgjJAQLHCOQ!4yWn$fF|EY|FOpIBVGp#9siH@GC>6+0_QtUkd+ z=q$Zihl`l7xVwCIC7%h!fZ#d(ikY~pT4;iEOq`%@&D~YO#O%9^1E`x!NCo%@cHC#; zs%lP<{u3sWsj6Kr?Myh@e%PJ)h6%mbqxanDVd6!?@(@KC7DSVR_bwj9!rp5Iq5D-> zu<*RK@3aOBii5s}b&qD@E@Rw&;|VNeTiQo>Ok=^lbzfxuToxuDU>|s7&ca8%dk04@ zW8qx6O!P{77GhkDV&c}buu{tL&}BCkG)^3dmGWYt*_4bk31ETU&=Bvxmj(Yl0~2`% zSTNO?ez+~3g?=R)l13k4;krX?^7?ER(z?YdM+#YRPHRfN$z#D_y5f;R5*9jZ=Ar^mUCT#2~d&_*ffDJQeZT2rq*^vKg#aUy=Mom^AH+dZ!M@{m0*EX@SvFpH*YzjoaF7%lCHG=J2Wzy11IAi$pnLk! zzzr)oXm=bm=%^zHRL_jTw>NUIC&N`?h$jc;v*L$X`f;Few^%WJ4+pjWPnAj|I5@7V zJoJ4m2V1DQDw9(=FtT-5-FA$F51q-wPUdq^mT-3XLkp3s>a<-Pvoo7|);j;l3*uUDnPYRXn&T z?TG1F=F3GyEdA+VFc+HqiPO*T=c2`W&5X~7xZrj~%$$+JMPMv{*3K+0i1GEanJ2jz za8=Hvna#yb_o;KWMBn?7bDrmL`g=bGF{ViuzxT6{M%6Vg-ZVWX`qy%CDqLYc@sNxC zIE53Y+TU?|?hrm!WZWvKAI`^!Szil>|H4OO?Ks8?eLfie_RJU~KD<@KSQlpU zF^6Kazk&~G+k2c@mV8up%J6(v@R5>W#AiG3v2K#1pk)IeW3C+_wYT!o;XzVsy}$P( z|4^8+i;qyoK>B(Z9~PF=#RFpaP;A{GnV-bR-N@MDJxBRCp(iP!bNSd>-gK&y$%nD4 z;%P&QkDlH+XEvSU<9zn!v&YN%IAoescCV5TyN1)}hTi65)SkBUOY8Y~sy^((fo49~ z#l*$37kv0TcwXvx%ZF+Ak;~J*@DXA4q&^`G=#J3#=G z?m)S~$pU0wQW%gjO#qp+ya9q)0`%spbX=P&z(yZNM>`?F6*Ih)TPVN;d*uMa@_Y1GfEbq->i|tmmU#dn;dWW)=UAe*P5#x$QEEq zL7iG&o&ZPucaA746yTG2|B($W0i0E{M)vRoC~q=WR}%`LPhD1@CJ`Vue7%NEsQ{fe zT^e3x0@&$9k2-Wg07<8oX3-S^G|PmV7pnvaOIV`yq*{QdPK{dMZVRwvQot{x?hC;E zGDv%Ng8)OXAtWwF7V48jNf_GC zGq^FFM4VoY!3%W~Z@S#a4g7_~*+@N$v7k5)5sBsH z!Inp-lAw~57mH_-;Ae+Tylp~akMX=EuR)@*_tugj<|GzYyIYO7B*Do3YGt{Ugu;%L zrEV)o_?Q|l3%4Utt5CKq+kwQK27Bw$btI0npIhJGKtg6u*z)(ANNllETcNU*M5Q`s z#pLZICb!SGS>{b5rTCuBHh&TyqP%Pm?jo_?L3ZWI5EAF~kFGoyMq*62k=?^c6493~ z+I@;9@iNuXendP88@Jc?)00SurX5^mn?^!iR&(`^qa^m$kgH>llX#MA>A=V(vB>A4 z!=(ZeEHl3~O-vGs%JOS_xg>lWv)5`;B<}I1JI)c4m>XK@xcU@{6V@*4{LYf-r`5gg za5;&s9kJ`Vmq}bJ)pe?@Brzpc;`HnWi8M!R=YF?I^cXff>)s>b)U#^?>PeJWD7vg` zB%z;?@3O0z#36SR*R*FOI%nQ=rCyM*liRYfx|4*s_VdP$wcIiif)6}+k??4L5ROpsW1q%M*!`vmx6zXhP?stY!Fx4^J`er1B z+|Jsqidq!p$~?DC(52v!AnmbuEQK3R=^mR6DNLO-$#efC3cY9BqE1Yquqh!l`t%G6 zSDc2$)XkwVaS|)0+myoLFT^3$1r**@)gCflL}6`~XRNgqg)<(~aURPl=$fR*MXjW8 zP<~Q;{%Q)(@0Z7)ccieKu_ob>GX*l}Rl?_u6h>G^CXRHc5Ui$gc!np1M=kurD|b>b zr!A6p`cWu~tWSyyqA=LrCz-j2g16p)l*{`l-03=&(i}lycKOs)sb~tvlPXfR;wVVF zI;WW&rm)%gU7ABGg{!?WNBlD=7*%VhCmo}ZoK2_mawv3rEzPLPr?AfSNk$ukLYYF) z(f(WtdJPJhqe%+U?7U2ZrqB^QJ8OLj1skgyS-a0r5UINzJ93_adV9|?;Ux+o#fisn zUZwCjYHaq)>l7>;O0x&mP+;oYo-nveVMurDiG>d+_+ARhaeYMLZmLRdXcGk!H%4w& zD~0T7$SZzML00xoUhOLi?ls#^c6Cvx%I!U=^nrr0Pg?%O&lFP4CKN2`rSMVde8Co3 zA)Fdl7e)*a;ykaTFlVq3V?*~Bolz1Z#(E^+Q=!rsD_84;tjD=uVT;`Nd7eXmxJ-2bT5Pt4m+%NNl zFzew~s?Qf9zrws~riBm#GwQ1BmI&eHzVn*5wGg*v_P-u)E5r=Btm~{*LS)t&-?*|? zh;R9qZ#;DpVx#~1YUzzaTruye{$-006I7#b&h-!?u}SOJnjJ#CC55*Fe1%vOzN98O zP>9nujWzsWA#`*CZeI%(;y~x1JMH_0Xg`~CM=nYT>x7xLda*)KPS{Jr!nAKU4ad$2QPJZ5u;ii;ljyGx?^Xzr6q)kq z?KL6Z+Zi@0-x9)6udH!Wtq^Cs>>pdz2{F3-`D6EoLPRBnJ&AlG#0wX-rrZ`GRv2@d z&bA97?493S|5AvN)%Th|z7ZlM+w1A@ZXq6f$+k@WB!q?O(Uuk8gkUO6Zrv^;fIUX01>>b ztyRUlMYyTita^K|2-7=u4SO9f!qHO2;X@9J@byss@bQO4aCJ0MvrG`-vf)iNw`38< z_iP>!enfw@#pQ=PC zv`N<9zJ8|^oa0o&H(+@Uqx6ee_Y>Bn#P&?Qws!zG*rvSDy_m>6FYBG>lAbr%LQ-Br80dddGoAcYVmTH|uCz z@S1I-=t9H5^oGp@HyW`DZnlfJ(s*6pW4n1f4SROt%KhFnii5}6o$#ljWmRf-dKZmw zbzA$o5E?CQt@hnvG?o;HtWrHdgBPW;+Bk-WiUVV{bv%tgeK>d|(Wrla+aW5A261WI zn*5{R$CKK-=KOIQ1KrZrKFX!BecA-a&jmDY%A9u`$)Yj6X7#!mJQ_!HJJzkFXnghF zzka8fhKt!qr?^uzE-CSxm}hAiHkvzME~gRCt8;F?OrvY>&J9vkG*(;pchRb*aY{4G z#pE^(?G9sChkG<4N-w+m*VA}*X#K{dMjFlO*B0@bX)Jc%Xj%1)2506+%eEIZl;si@ z_kT^pzgBO_=yx>g@=KNwA88cK4A`>%3ync?gWY#aiLs+L$Nfk@F>3N>ZWRs`W2XPL ztv40K$THu!?WM99y(%BK4H_=S#-?}=19dSfNIlPmT4GEHFY$EM6(iAR#rDv#V!YFC z*`8%6#+uGxui{B!oIb0(qjriIIthh4x@L%Rz-iu2B@;2)C)Vto2r;a`xO*>IAO=u+)OBg}u;1~I&Y z)B>h#664MyPQY?^F=neR2=w$6qd((oa`a9ywgsi66!?j8eUV{md5{=W)yh&E_lR+% z#XjxJJ~2Lto~Nltim@Rw?8wY$F)rAtrQ5}cp|8hD_dYDfp|1HE@u^~Tp1+sD$`HdY z$?NEqV`4~LWHX=Uh@okGG*h}jjIds#tY4U7geaWYH41_ zm+<5gF@{)mhx5;f;j13E|Jr#m?zZVhv|kbWMvUsV z>|Ce2VoWcdo)`Q;jEtzty!1z6d|B;sQq&}di+=aXTdiVTdLNtr>bV$(mvjpTcZw00 zDk&KER*W~BtP3qZh_Px~bK%C%VwA}2F51^C#xFIBjAODA?9a_-lnjud)!T%5cd!Ij zW;dB{l_cORZDA=7lR&lcGi%aF34(Y@Y%5I(8ul7++;t={vp&O#93w%#rX4qToCNY6 zZQQdHCGaW@<<%QYaO==8{>SMO%y4A!htH89)6h&X)l`CSA8Q3G=1Z`#!js%?A;Fb& zX)0!k1QXoTsX}WB5@${lUa*zmUH=Qh$Ezf8s9huax>kbH{8u6kX9=|ZBk5TiC5SZF z5ZiB&piM;}_VJKlS(Ak%VTS}HSubJxN-#X!r??_eg55R)N?L*?c%*&2L}s4^3p%Hs z(vFaz@NC7Yc~KG!N^mY+8!N$1r+1}+i4x51Ih`q1fSaVK-&LX<8^P&WHK}#PiT$Mn)=*eTl>k^Do3wpBXmIV7+6q+{GO3*CI zYYMBAU~%N^=Hm|~;Mm=0KJ`QbWj(j2_gW+f=<0d;u3dt<^NB4(UrAt^G`7{KOM+aN z($=Nj637|bKHK_90*~I-X9vDXaN|ZuTb@iYre&+Nm&p|)-HXxQFsK-x=izyeVliA4 z?mSmhEyl(AZ7-&cC`S45z_aJZ6=V60HP$95#)EGr4MMwO*!OSM|)&pZFR7 zQ9pflkt4TUDJt+I`Zz|k@=}z&AMww!&(ELnkNz*MATGGe`Vl|lpV^l@yq`Ru@+JEC zGe3RpeRh597d5iXz}1)NV_*A!7W@4EzcPPx|1$%WeTkp(KigklKbeB;xF=_Qh<}zX zYabO~@$w=5t9amjrMO7bhxm8lo(Q*t;Z@$mKT98z50AZqy@{XkXMTTX{~v8vx!e56 z+R&T$@5;^1IVI}^JBfc6p414YM=9Jy6Q#zN4e!( zf95&^FXDd+3qL%{A9Y|m@xP7FjEZyWn>>lX_*{#ic$v%LO~y2c-MmOtt?eY(yc^`bxOTYdUipWgOI9j{L}{Gj zfTR9|3`*JG;_cpubJFbHKj!z{#DF;tN9I}w6P77ed7r0+{8;aY{K#YW{>V#1 ziK7hF1Mg3V5^H@-VrFdEM;sGgzv>YjMqJx`ukuq~7%{Q(`8ADqVMMIaw;LuF;lzvN zp|>25gcAYYpG$+Dh5yhE_y5o@BYx{1 zP4R{WOwAC2kPU9qE*1}ij;cJtKt=%6M4&9$bR5%W{oN_yfP-;E8SYA2#$8(TO zs7PzhOt_p(ypvvX@xz1^Vv5n~d9O-Ri2GU}iB`*0LN7XcgXZq*gxV^)l=PG_TgS>9 zpDj1b%x0;y%p}>rKfCeksf=IWia$)R! 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 diff --git a/tests/test_grid.py b/tests/test_grid.py index 64305d6..ee46170 100644 --- a/tests/test_grid.py +++ b/tests/test_grid.py @@ -14,12 +14,6 @@ def test_scaling(): 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():