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 0000000..733af7b
Binary files /dev/null and b/tests/data/PM1050NEG_aeff_slow.npz differ
diff --git a/tests/data/PM1050NEG_dispersion.npz b/tests/data/PM1050NEG_dispersion.npz
new file mode 100644
index 0000000..2c781b4
Binary files /dev/null and b/tests/data/PM1050NEG_dispersion.npz differ
diff --git a/tests/data/PM1050NEG_loss.npz b/tests/data/PM1050NEG_loss.npz
new file mode 100644
index 0000000..3f49b26
Binary files /dev/null and b/tests/data/PM1050NEG_loss.npz differ
diff --git a/tests/test_dispersion.py b/tests/test_dispersion.py
new file mode 100644
index 0000000..bc2593a
--- /dev/null
+++ b/tests/test_dispersion.py
@@ -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
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():