better sort_axis, extinction distance

src/scgenerator/physics/fiber.py

@@ -2,6 +2,7 @@ from typing import Any, Iterable, Literal, TypeVar
 
 import numpy as np
 from numpy.fft import fft, ifft
+from numpy import e
 from numpy.polynomial.chebyshev import Chebyshev, cheb2poly
 from scipy.interpolate import interp1d
 
@@ -792,7 +793,7 @@ def compute_capillary_loss(
     interpolation_range: tuple[float, float],
     he_mode: tuple[int, int],
 ) -> np.ndarray:
-    mask = l < interpolation_range[1]
+    mask = (l < interpolation_range[1]) & (l > 0)
     alpha = capillary_loss(l[mask], he_mode, core_radius)
     out = np.zeros_like(l)
     out[mask] = alpha
@@ -1099,14 +1100,12 @@ def effective_radius_HCARF(core_radius, t, f1, f2, wl_for_disp):
     return f1 * core_radius * (1 - f2 * wl_for_disp ** 2 / (core_radius * t))
 
 
-def capillary_loss(
+def capillary_loss(wl: np.ndarray, he_mode: tuple[int, int], core_radius: float) -> np.ndarray:
-    wl_for_disp: np.ndarray, he_mode: tuple[int, int], core_radius: float
-) -> np.ndarray:
     """computes the wavelenth dependent capillary loss according to Marcatili
 
     Parameters
     ----------
-    wl_for_disp : np.ndarray, shape (n, )
+    wl : np.ndarray, shape (n, )
         wavelength array
     he_mode : tuple[int, int]
         tuple of mode (n, m)
@@ -1118,9 +1117,10 @@ def capillary_loss(
     np.ndarray
         loss in 1/m
     """
-    alpha = np.zeros_like(wl_for_disp)
+    chi_silica = mat.sellmeier(wl, utils.load_material_dico("silica"))
-    mask = wl_for_disp > 0
-    chi_silica = mat.sellmeier(wl_for_disp[mask], utils.load_material_dico("silica"))
     nu_n = 0.5 * (chi_silica + 2) / np.sqrt(chi_silica)
-    alpha[mask] = nu_n * (u_nm(*he_mode) * wl_for_disp[mask] / pipi) ** 2 * core_radius ** -3
+    return nu_n * (u_nm(*he_mode) * wl / pipi) ** 2 * core_radius ** -3
-    return alpha
+
+
+def extinction_distance(loss: T, ratio=1 / e) -> T:
+    return np.log(ratio) / -loss

src/scgenerator/physics/pulse.py

@@ -275,7 +275,7 @@ def conform_pulse_params(
         return width, t0, peak_power, energy
     else:
         if soliton_num is None:
-            soliton_num = np.sqrt(peak_power * gamma * t0 ** 2 / abs(beta2))
+            soliton_num = np.sqrt(peak_power * gamma * L_D(t0, beta2))
         return width, t0, peak_power, energy, soliton_num
 
 

src/scgenerator/utils/parameter.py

@@ -1095,6 +1095,12 @@ class PlotRange:
     unit: Callable[[float], float] = Parameter(units.is_unit, converter=units.get_unit)
     conserved_quantity: bool = Parameter(boolean, default=True)
 
+    def __post_init__(self):
+        if self.left >= self.right:
+            raise ValueError(
+                f"left value {self.left!r} must be strictly smaller than right value {self.right!r}"
+            )
+
     def __str__(self):
         return f"{self.left:.1f}-{self.right:.1f} {self.unit.__name__}"
 
@@ -1110,51 +1116,40 @@ class PlotRange:
 def sort_axis(
     axis: np.ndarray, plt_range: PlotRange
 ) -> tuple[np.ndarray, np.ndarray, tuple[float, float]]:
-    """
+    """convert an array according to the given range
-    given an axis, returns this axis cropped according to the given range, converted and sorted
 
     Parameters
     ----------
-    axis : 1D array containing the original axis (usual the w or t array)
+    axis : np.ndarray, shape (n,)
-    plt_range : tupple (min, max, conversion_function) used to crop the axis
+        array
+    plt_range : PlotRange
+        range to crop in
 
     Returns
     -------
-    cropped : the axis cropped, converted and sorted
+    np.ndarray
-    indices : indices to use to slice and sort other array in the same fashion
+        new array converted to the desired unit and cropped in the given range
-    extent : tupple with min and max of cropped
+    np.ndarray
+        indices of the concerved values
+    tuple[float, float]
+        actual minimum and maximum of the new axis
 
     Example
     -------
-    w = np.append(np.linspace(0, -10, 20), np.linspace(0, 10, 20))
+    >> sort_axis([18.0, 19.0, 20.0, 13.0, 15.2], PlotRange(1400, 1900, "cm"))
-    t = np.linspace(-10, 10, 400)
+    ([1520.0, 1800.0, 1900.0], [4, 0, 1], (1520.0, 1900.0))
-    W, T = np.meshgrid(w, t)
-    y = np.exp(-W**2 - T**2)
-
-    # Define ranges
-    rw = (-4, 4, s)
-    rt = (-2, 6, s)
-
-    w, cw = sort_axis(w, rw)
-    t, ct = sort_axis(t, rt)
-
-    # slice y according to the given ranges
-    y = y[ct][:, cw]
     """
     if isinstance(plt_range, tuple):
         plt_range = PlotRange(*plt_range)
-    r = np.array((plt_range.left, plt_range.right), dtype="float")
 
-    indices = np.arange(len(axis))[
+    masked = np.ma.array(axis, mask=~np.isfinite(axis))
-        (axis <= np.max(plt_range.unit(r))) & (axis >= np.min(plt_range.unit(r)))
+    converted = plt_range.unit.inv(masked)
-    ]
+    converted[(converted < plt_range.left) | (converted > plt_range.right)] = np.ma.masked
-    cropped = axis[indices]
+    indices = np.arange(len(axis))[~converted.mask]
-    order = np.argsort(plt_range.unit.inv(cropped))
+    cropped = converted.compressed()
-    indices = indices[order]
+    order = cropped.argsort()
-    cropped = cropped[order]
-    out_ax = plt_range.unit.inv(cropped)
 
-    return out_ax, indices, (out_ax[0], out_ax[-1])
+    return cropped[order], indices[order], (cropped.min(), cropped.max())
 
 
 def get_arg_names(func: Callable) -> list[str]: