capillary loss

src/scgenerator/__init__.py

@@ -1,8 +1,15 @@
-from .initialize import ParamSequence, RecoveryParamSequence, ContinuationParamSequence
-from .io import Paths, load_toml
+from .initialize import (
+    ParamSequence,
+    RecoveryParamSequence,
+    ContinuationParamSequence,
+    Params,
+    Config,
+)
+from .io import Paths, load_toml, load_params
 from .math import abs2, argclosest, span
 from .physics import fiber, materials, pulse, simulate, units
 from .physics.simulate import RK4IP, new_simulation, resume_simulations
 from .plotting import plot_avg, plot_results_1D, plot_results_2D, plot_spectrogram
 from .spectra import Pulse
-from . import utils
+from .utils.parameter import BareParams, BareConfig
+from . import utils, io, initialize, math

src/scgenerator/defaults.py

@@ -24,7 +24,7 @@ default_parameters = dict(
     repeat=1,
     tolerated_error=1e-11,
     lower_wavelength_interp_limit=100e-9,
-    upper_wavelength_interp_limit=1900e-9,
+    upper_wavelength_interp_limit=2000e-9,
     interp_degree=8,
     ideal_gas=False,
     readjust_wavelength=False,

src/scgenerator/initialize.py

@@ -43,6 +43,13 @@ class Params(BareParams):
             )
             logger.info("added some quantum noise")
 
+        if self.step_size is not None:
+            self.error_ok = self.step_size
+            self.adapt_step_size = False
+        else:
+            self.error_ok = self.tolerated_error
+            self.adapt_step_size = True
+
         self.spec_0 = np.fft.fft(self.field_0)
 
     def __build_sim_grid(self):
@@ -112,6 +119,8 @@ class Params(BareParams):
         if "raman" in self.behaviors:
             self.hr_w = fiber.delayed_raman_w(self.t, self.dt, self.raman_type)
 
+        self.alpha = fiber.compute_loss(self)
+
     def __compute_custom_pulse(self):
         logger = get_logger(__name__)
 
@@ -129,15 +138,9 @@ class Params(BareParams):
             self.wavelength = pulse.correct_wavelength(self.wavelength, self.w_c, self.field_0)
             logger.info(f"moved wavelength from {1e9*old_wl:.2f} to {1e9*self.wavelength:.2f}")
             self.w_c, self.w0, self.w, self.w_power_fact = update_frequency_domain(
-                self.t, self.wavelength
+                self.t, self.wavelength, self.interp_degree
             )
 
-        if self.step_size is not None:
-            self.error_ok = self.step_size
-            self.adapt_step_size = False
-        else:
-            self.error_ok = self.tolerated_error
-            self.adapt_step_size = True
         return did_set_custom_pulse
 
 
@@ -191,6 +194,10 @@ class Config(BareConfig):
             else:
                 for param in hc_model_specific_parameters[self.model]:
                     self.get_fiber(param)
+        if self.contains("loss"):
+            if self.loss == "capillary":
+                for param in ["core_radius", "he_mode"]:
+                    self.get_fiber(param)
         for param in ["length", "input_transmission"]:
             self.get(param)
 
@@ -612,21 +619,54 @@ def build_sim_grid(
     length: float,
     z_num: int,
     wavelength: float,
+    deg: int,
     time_window: float = None,
     t_num: int = None,
     dt: float = None,
-):
+) -> tuple[
+    np.ndarray, np.ndarray, float, int, float, np.ndarray, float, np.ndarray, np.ndarray, np.ndarray
+]:
     """computes a bunch of values that relate to the simulation grid
 
     Parameters
     ----------
-    params : dict
-        flattened parameter dictionary
+    length : float
+        length of the fiber in m
+    z_num : int
+        number of spatial points
+    wavelength : float
+        pump wavelength in m
+    deg : int
+        dispersion interpolation degree
+    time_window : float, optional
+        total width of the temporal grid in s, by default None
+    t_num : int, optional
+        number of temporal grid points, by default None
+    dt : float, optional
+        spacing of the temporal grid in s, by default None
 
     Returns
     -------
-    dict
-        updated parameter dictionary
+    z_targets : np.ndarray, shape (z_num, )
+        spatial points in m
+    t : np.ndarray, shape (t_num, )
+        temporal points in s
+    time_window : float
+        total width of the temporal grid in s, by default None
+    t_num : int
+        number of temporal grid points, by default None
+    dt : float
+        spacing of the temporal grid in s, by default None
+    w_c : np.ndarray, shape (t_num, )
+        centered angular frequencies in rad/s where 0 is the pump frequency
+    w0 : float
+        pump angular frequency
+    w : np.ndarray, shape (t_num, )
+        actual angualr frequency grid in rad/s
+    w_power_fact : np.ndarray, shape (deg, t_num)
+        set of all the necessaray powers of w_c
+    l : np.ndarray, shape (t_num)
+        wavelengths in m
     """
     t = tspace(time_window, t_num, dt)
 
@@ -634,7 +674,7 @@ def build_sim_grid(
     dt = t[1] - t[0]
     t_num = len(t)
     z_targets = np.linspace(0, length, z_num)
-    w_c, w0, w, w_power_fact = update_frequency_domain(t, wavelength)
+    w_c, w0, w, w_power_fact = update_frequency_domain(t, wavelength, deg)
     l = units.To.m(w)
     return z_targets, t, time_window, t_num, dt, w_c, w0, w, w_power_fact, l
 
@@ -653,12 +693,18 @@ def build_sim_grid_in_place(params: BareParams):
         params.w_power_fact,
         params.l,
     ) = build_sim_grid(
-        params.length, params.z_num, params.wavelength, params.time_window, params.t_num, params.dt
+        params.length,
+        params.z_num,
+        params.wavelength,
+        params.interp_degree,
+        params.time_window,
+        params.t_num,
+        params.dt,
     )
 
 
 def update_frequency_domain(
-    t: np.ndarray, wavelength: float
+    t: np.ndarray, wavelength: float, deg: int
 ) -> Tuple[np.ndarray, float, np.ndarray, np.ndarray]:
     """updates the frequency grid
 
@@ -668,6 +714,8 @@ def update_frequency_domain(
         time array
     wavelength : float
         wavelength
+    deg : int
+        interpolation degree of the dispersion
 
     Returns
     -------
@@ -677,5 +725,5 @@ def update_frequency_domain(
     w_c = wspace(t)
     w0 = units.m(wavelength)
     w = w_c + w0
-    w_power_fact = np.array([power_fact(w_c, k) for k in range(2, 11)])
+    w_power_fact = np.array([power_fact(w_c, k) for k in range(2, deg + 3)])
     return w_c, w0, w, w_power_fact

src/scgenerator/physics/fiber.py

@@ -1,21 +1,25 @@
-from typing import Any, Dict, Iterable, List, Literal, Tuple, Union
+from typing import Any, Dict, Iterable, List, Literal, Optional, Tuple, Union
 
 import numpy as np
+from numpy.ma import core
 import toml
 from numpy.fft import fft, ifft
 from numpy.polynomial.chebyshev import Chebyshev, cheb2poly
+from numpy.polynomial.polynomial import Polynomial
 from scipy.interpolate import interp1d
 
 from ..logger import get_logger
 
 from .. import io
 from ..math import abs2, argclosest, power_fact, u_nm
-from ..utils.parameter import BareParams, hc_model_specific_parameters
+from ..utils.parameter import BareConfig, BareParams, hc_model_specific_parameters
 from ..utils.cache import np_cache
 from . import materials as mat
 from . import units
 from .units import c, pi
 
+pipi = 2 * pi
+
 
 def lambda_for_dispersion(left: float, right: float) -> np.ndarray:
     """Returns a wl vector for dispersion calculation
@@ -90,12 +94,12 @@ def dispersion_parameter(n_eff: np.ndarray, lambda_: np.ndarray):
 
 def beta2_to_D(beta2, λ):
     """returns the D parameter corresponding to beta2(λ)"""
-    return -(2 * pi * c) / (λ ** 2) * beta2
+    return -(pipi * c) / (λ ** 2) * beta2
 
 
 def D_to_beta2(D, λ):
     """returns the beta2 parameters corresponding to D(λ)"""
-    return -(λ ** 2) / (2 * pi * c) * D
+    return -(λ ** 2) / (pipi * c) * D
 
 
 def plasma_dispersion(lambda_, number_density, simple=False):
@@ -148,7 +152,7 @@ def n_eff_marcatili(lambda_, n_gas_2, core_radius, he_mode=(1, 1)):
     """
     u = u_nm(*he_mode)
 
-    return np.sqrt(n_gas_2 - (lambda_ * u / (2 * pi * core_radius)) ** 2)
+    return np.sqrt(n_gas_2 - (lambda_ * u / (pipi * core_radius)) ** 2)
 
 
 def n_eff_marcatili_adjusted(lambda_, n_gas_2, core_radius, he_mode=(1, 1), fit_parameters=()):
@@ -179,7 +183,7 @@ def n_eff_marcatili_adjusted(lambda_, n_gas_2, core_radius, he_mode=(1, 1), fit_
 
     corrected_radius = effective_core_radius(lambda_, core_radius, *fit_parameters)
 
-    return np.sqrt(n_gas_2 - (lambda_ * u / (2 * pi * corrected_radius)) ** 2)
+    return np.sqrt(n_gas_2 - (lambda_ * u / (pipi * corrected_radius)) ** 2)
 
 
 @np_cache
@@ -242,7 +246,7 @@ def n_eff_hasan(
 
     R_eff = f1 * core_radius * (1 - f2 * lambda_ ** 2 / (core_radius * capillary_thickness))
 
-    n_eff_2 = n_gas_2 - (u * lambda_ / (2 * pi * R_eff)) ** 2
+    n_eff_2 = n_gas_2 - (u * lambda_ / (pipi * R_eff)) ** 2
 
     chi_sil = mat.sellmeier(lambda_, io.load_material_dico("silica"))
 
@@ -593,7 +597,7 @@ def PCF_dispersion(lambda_, pitch, ratio_d, w0=None, n2=None, A_eff=None):
 
     n_co = 1.45
     a_eff = pitch / np.sqrt(3)
-    pi2a = 2 * pi * a_eff
+    pi2a = pipi * a_eff
 
     ratio_l = lambda_ / pitch
 
@@ -643,7 +647,7 @@ def PCF_dispersion(lambda_, pitch, ratio_d, w0=None, n2=None, A_eff=None):
         if A_eff is None:
             V_eff = pi2a / lambda_ * np.sqrt(n_co ** 2 - n_FSM2)
             w_eff = a_eff * (0.65 + 1.619 / V_eff ** 1.5 + 2.879 / V_eff ** 6)
-            A_eff = interp1d(lambda_, w_eff, kind="linear")(units.m.inv(w0)) ** 2 * pi
+            A_eff = interp1d(lambda_, w_eff, kind="linear")(units.m.inv(w0)) ** pipi
 
         if n2 is None:
             n2 = 2.6e-20
@@ -652,6 +656,16 @@ def PCF_dispersion(lambda_, pitch, ratio_d, w0=None, n2=None, A_eff=None):
         return beta2, gamma
 
 
+def compute_loss(params: BareParams) -> Optional[np.ndarray]:
+    if params.loss == "capillary":
+        mask = params.l < params.upper_wavelength_interp_limit
+        alpha = capillary_loss(params.l[mask], params.he_mode, params.core_radius)
+        out = np.zeros_like(params.l)
+        out[mask] = alpha
+        return out
+    return None
+
+
 def compute_dispersion(params: BareParams):
     """dispatch function depending on what type of fiber is used
 
@@ -695,7 +709,7 @@ def compute_dispersion(params: BareParams):
             )
 
         else:
-            material_dico = toml.loads(io.Paths.gets("gas"))[params.gas_name]
+            material_dico = io.load_material_dico(params.gas_name)
             if params.dynamic_dispersion:
                 return dynamic_HCPCF_dispersion(
                     lambda_,
@@ -788,18 +802,20 @@ def dispersion_coefficients(
     logger.debug(
         f"interpolating dispersion between {lambda_[r].min()*1e9:.1f}nm and {lambda_[r].max()*1e9:.1f}nm"
     )
-    # import matplotlib.pyplot as plt
 
     # we get the beta2 Taylor coeffiecients by making a fit around w0
     w_c = units.m(lambda_) - w0
-    # interp = interp1d(w_c[r], beta2[r])
-    # w_c = np.linspace(w_c)
-    # fig, ax = plt.subplots()
-    # ax.plot(w_c[r], beta2[r])
-    # fig.show()
+    interp = interp1d(w_c[r], beta2[r])
+    w_c = np.linspace(w_c[r].min(), w_c[r].max(), len(w_c[r]))
 
-    fit = Chebyshev.fit(w_c[r], beta2[r], deg)
-    beta2_coef = cheb2poly(fit.convert().coef) * np.cumprod([1] + list(range(1, deg + 1)))
+    # import matplotlib.pyplot as plt
+
+    # ax = plt.gca()
+    # ax.plot(w_c, interp(w_c) * 1e28)
+
+    fit = Chebyshev.fit(w_c, interp(w_c), deg)
+    poly_coef = cheb2poly(fit.convert().coef)
+    beta2_coef = poly_coef * np.cumprod([1] + list(range(1, deg + 1)))
 
     return beta2_coef
 
@@ -935,13 +951,13 @@ def create_non_linear_op(behaviors, w_c, w0, gamma, raman_type="stolen", f_r=0,
 
     if isinstance(gamma, (float, int)):
 
-        def N_func(spectrum: np.ndarray, r=0):
+        def N_func(spectrum: np.ndarray, r=0) -> np.ndarray:
             field = ifft(spectrum)
             return -1j * gamma * (1 + ss_part) * fft(field * (spm_part(field) + raman_part(field)))
 
     else:
 
-        def N_func(spectrum: np.ndarray, r=0):
+        def N_func(spectrum: np.ndarray, r=0) -> np.ndarray:
             field = ifft(spectrum)
             return (
                 -1j * gamma(r) * (1 + ss_part) * fft(field * (spm_part(field) + raman_part(field)))
@@ -950,7 +966,7 @@ def create_non_linear_op(behaviors, w_c, w0, gamma, raman_type="stolen", f_r=0,
     return N_func
 
 
-def fast_dispersion_op(w_c, beta_arr, power_fact_arr, where=slice(None)):
+def fast_dispersion_op(w_c, beta_arr, power_fact_arr, where=slice(None), alpha=None):
     """
     dispersive operator
 
@@ -975,8 +991,10 @@ def fast_dispersion_op(w_c, beta_arr, power_fact_arr, where=slice(None)):
 
     out = np.zeros_like(dispersion)
     out[where] = dispersion[where]
-
-    return -1j * out
+    if alpha is None:
+        return -1j * out
+    else:
+        return -1j * out - alpha / 2
 
 
 def _fast_disp_loop(dispersion: np.ndarray, beta_arr, power_fact_arr):
@@ -1046,6 +1064,31 @@ def effective_radius_HCARF(core_radius, t, f1, f2, lambda_):
     return f1 * core_radius * (1 - f2 * lambda_ ** 2 / (core_radius * t))
 
 
+def capillary_loss(lambda_: np.ndarray, he_mode: tuple[int, int], core_radius: float) -> np.ndarray:
+    """computes the wavelenth dependent capillary loss according to Marcatili
+
+    Parameters
+    ----------
+    lambda_ : np.ndarray, shape (n, )
+        wavelength array
+    he_mode : tuple[int, int]
+        tuple of mode (n, m)
+    core_radius : float
+        in m
+
+    Returns
+    -------
+    np.ndarray
+        loss in 1/m
+    """
+    alpha = np.zeros_like(lambda_)
+    mask = lambda_ > 0
+    chi_silica = mat.sellmeier(lambda_[mask], io.load_material_dico("silica"))
+    nu_n = 0.5 * (chi_silica + 2) / np.sqrt(chi_silica)
+    alpha[mask] = nu_n * (u_nm(*he_mode) * lambda_[mask] / pipi) ** 2 * core_radius ** -3
+    return alpha
+
+
 if __name__ == "__main__":
     w = np.linspace(0, 1, 4096)
     c = np.arange(8)

src/scgenerator/physics/pulse.py

@@ -282,14 +282,24 @@ def gauss_pulse(t, t0, P0, offset=0):
     return np.sqrt(P0) * np.exp(-(((t - offset) / t0) ** 2))
 
 
-def photon_number(spectrum, w, dw, gamma):
+def photon_number(spectrum, w, dw, gamma) -> float:
     return np.sum(1 / gamma * abs2(spectrum) / w * dw)
 
 
-def pulse_energy(spectrum, w, dw, _):
+def photon_number_with_loss(spectrum, w, dw, gamma, alpha, h) -> float:
+    spec2 = abs2(spectrum)
+    return np.sum(1 / gamma * spec2 / w * dw) - h * np.sum(alpha / gamma * spec2 / w * dw)
+
+
+def pulse_energy(spectrum, dw) -> float:
     return np.sum(abs2(spectrum) * dw)
 
 
+def pulse_energy_with_loss(spectrum, dw, alpha, h) -> float:
+    spec2 = abs2(spectrum)
+    return np.sum(spec2 * dw) - h * np.sum(alpha * spec2 * dw)
+
+
 def technical_noise(rms_noise, relative_factor=0.4):
     """
     To implement technical noise as described in Grenier2019, we need to know the

src/scgenerator/physics/simulate.py

@@ -61,8 +61,11 @@ class RK4IP:
         self.save_data = save_data
 
         self.w_c = params.w_c
+        self.w = params.w
+        self.dw = self.w[1] - self.w[0]
         self.w0 = params.w0
         self.w_power_fact = params.w_power_fact
+        self.alpha = params.alpha
         self.spec_0 = params.spec_0
         self.z_targets = params.z_targets
         self.z_final = params.length
@@ -85,19 +88,38 @@ class RK4IP:
         )
 
         if self.dynamic_dispersion:
-            self.disp = lambda r: fast_dispersion_op(self.w_c, self.beta(r), self.w_power_fact)
+            self.disp = lambda r: fast_dispersion_op(
+                self.w_c, self.beta(r), self.w_power_fact, alpha=self.alpha
+            )
         else:
-            self.disp = lambda r: fast_dispersion_op(self.w_c, self.beta, self.w_power_fact)
+            self.disp = lambda r: fast_dispersion_op(
+                self.w_c, self.beta, self.w_power_fact, alpha=self.alpha
+            )
 
         # Set up which quantity is conserved for adaptive step size
         if self.adapt_step_size:
-            if "raman" in self.behaviors:
-                self.conserved_quantity_func = pulse.photon_number
+            if "raman" in self.behaviors and self.alpha is not None:
+                self.logger.debug("Conserved quantity : photon number with loss")
+                self.conserved_quantity_func = lambda spectrum, h: pulse.photon_number_with_loss(
+                    spectrum, self.w, self.dw, self.gamma, self.alpha, h
+                )
+            elif "raman" in self.behaviors:
+                self.logger.debug("Conserved quantity : photon number without loss")
+                self.conserved_quantity_func = lambda spectrum, h: pulse.photon_number(
+                    spectrum, self.w, self.dw, self.gamma
+                )
+            elif self.alpha is not None:
+                self.logger.debug("Conserved quantity : energy with loss")
+                self.conserved_quantity_func = lambda spectrum, h: pulse.pulse_energy_with_loss(
+                    spectrum, self.dw, self.alpha, h
+                )
             else:
-                self.logger.debug("energy conserved")
-                self.conserved_quantity_func = pulse.pulse_energy
+                self.logger.debug("Conserved quantity : energy without loss")
+                self.conserved_quantity_func = lambda spectrum, h: pulse.pulse_energy(
+                    spectrum, self.dw
+                )
         else:
-            self.conserved_quantity_func = lambda a, b, c, d: 0
+            self.conserved_quantity_func = lambda spectrum, h: 0.0
 
     def _setup_sim_parameters(self):
         # making sure to keep only the z that we want
@@ -114,12 +136,7 @@ class RK4IP:
         self.current_spectrum = self.spec_0.copy()
         self.stored_spectra = self.starting_num * [None] + [self.current_spectrum.copy()]
         self.cons_qty = [
-            self.conserved_quantity_func(
-                self.current_spectrum,
-                self.w_c + self.w0,
-                self.d_w,
-                self.gamma,
-            ),
+            self.conserved_quantity_func(self.current_spectrum, 0),
             0,
         ]
         self.size_fac = 2 ** (1 / 5)
@@ -255,9 +272,7 @@ class RK4IP:
             new_spectrum = expD * (A_I + k1 / 6 + k2 / 3 + k3 / 3) + k4 / 6
 
             if self.adapt_step_size:
-                self.cons_qty[step] = self.conserved_quantity_func(
-                    new_spectrum, self.w_c + self.w0, self.d_w, self.gamma
-                )
+                self.cons_qty[step] = self.conserved_quantity_func(new_spectrum, h)
                 curr_p_change = np.abs(self.cons_qty[step - 1] - self.cons_qty[step])
                 cons_qty_change_ok = self.error_ok * self.cons_qty[step - 1]
 
@@ -275,6 +290,8 @@ class RK4IP:
                 else:
                     keep = True
                     h_next_step = h
+            else:
+                keep = True
         return h, h_next_step, new_spectrum
 
     def step_saved(self):

src/scgenerator/scripts/__init__.py

@@ -66,13 +66,17 @@ def plot_dispersion(config_path: Path, lim: tuple[float, float] = None):
     right.grid()
     all_labels = []
     already_plotted = set()
+    loss_ax = None
+    plt.sca(left)
     for style, lbl, params in plot_helper(config_path):
+        if params.alpha is not None and loss_ax is None:
+            loss_ax = right.twinx()
         if (bbb := tuple(params.beta)) not in already_plotted:
             already_plotted.add(bbb)
         else:
             continue
 
-        lbl = plot_1_dispersion(lim, left, right, style, lbl, params)
+        lbl = plot_1_dispersion(lim, left, right, style, lbl, params, loss_ax)
         all_labels.append(lbl)
     finish_plot(fig, right, all_labels, params)
 
@@ -84,16 +88,19 @@ def plot_init(
     lim_disp: tuple[float, float] = None,
 ):
     fig, ((tl, tr), (bl, br)) = plt.subplots(2, 2, figsize=(14, 10))
+    loss_ax = None
     tl.grid()
     tr.grid()
     all_labels = []
     already_plotted = set()
     for style, lbl, params in plot_helper(config_path):
+        if params.alpha is not None and loss_ax is None:
+            loss_ax = tr.twinx()
         if (fp := fingerprint(params)) not in already_plotted:
             already_plotted.add(fp)
         else:
             continue
-        lbl = plot_1_dispersion(lim_disp, tl, tr, style, lbl, params)
+        lbl = plot_1_dispersion(lim_disp, tl, tr, style, lbl, params, loss_ax)
         lbl = plot_1_init_spec_field(lim_field, lim_spec, bl, br, style, lbl, params)
         all_labels.append(lbl)
     finish_plot(fig, tr, all_labels, params)
@@ -132,6 +139,7 @@ def plot_1_dispersion(
     style: dict[str, Any],
     lbl: list[str],
     params: BareParams,
+    loss: plt.Axes = None,
 ):
     beta_arr = fiber.dispersion_from_coefficients(params.w_c, params.beta)
     wl = units.m.inv(params.w)
@@ -150,13 +158,21 @@ def plot_1_dispersion(
     m &= wl >= (lim[0] if lim[0] < 1 else lim[0] * 1e-9)
     m &= wl <= (lim[1] if lim[1] < 1 else lim[1] * 1e-9)
 
-    left.annotate(
+    info_str = (
         rf"$\lambda_{{\mathrm{{min}}}}={np.min(params.l[params.l>0])*1e9:.1f}$ nm"
-        f"lower interpolation limit : {params.interp_range[0]*1e9:.1f} nm",
-        (0, 1),
+        + f"\nlower interpolation limit : {params.interp_range[0]*1e9:.1f} nm\n"
+        + f"max time delay : {params.t.max()*1e12:.1f} ps"
+    )
+
+    left.annotate(
+        info_str,
+        xy=(1, 1),
+        xytext=(-12, -12),
         xycoords="axes fraction",
+        textcoords="offset points",
         va="top",
-        ha="left",
+        ha="right",
+        backgroundcolor=(1, 1, 1, 0.4),
     )
 
     m = np.argwhere(m)[:, 0]
@@ -170,11 +186,18 @@ def plot_1_dispersion(
     right.set_ylabel(units.D_ps_nm_km.label)
 
     # plot beta
-    left.plot(units.To.Prad_s(params.w[m]), units.beta2_fs_cm.inv(beta_arr[m]), label=" ", **style)
+    left.plot(units.To.nm(params.w[m]), units.beta2_fs_cm.inv(beta_arr[m]), label=" ", **style)
     left.set_ylabel(units.beta2_fs_cm.label)
 
-    left.set_xlabel(units.Prad_s.label)
+    left.set_xlabel(units.nm.label)
     right.set_xlabel("wavelength (nm)")
+
+    if params.alpha is not None and loss is not None:
+        loss.plot(1e9 * wl[m], params.alpha[m], c="r", ls="--")
+        loss.set_ylabel("loss (1/m)", color="r")
+        loss.set_yscale("log")
+        loss.tick_params(axis="y", labelcolor="r")
+
     return lbl
 
 

src/scgenerator/utils/parameter.py

@@ -320,6 +320,7 @@ class BareParams:
     input_transmission: float = Parameter(in_range_incl(0, 1))
     gamma: float = Parameter(non_negative(float, int))
     n2: float = Parameter(non_negative(float, int))
+    loss: str = Parameter(literal("capillary"))
     effective_mode_diameter: float = Parameter(positive(float, int))
     A_eff: float = Parameter(non_negative(float, int))
     pitch: float = Parameter(in_range_excl(0, 1e-3))
@@ -383,6 +384,7 @@ class BareParams:
     # computed
     field_0: np.ndarray = Parameter(type_checker(np.ndarray))
     spec_0: np.ndarray = Parameter(type_checker(np.ndarray))
+    alpha: np.ndarray = Parameter(type_checker(np.ndarray))
     w: np.ndarray = Parameter(type_checker(np.ndarray))
     l: np.ndarray = Parameter(type_checker(np.ndarray))
     w_c: np.ndarray = Parameter(type_checker(np.ndarray))
@@ -421,7 +423,17 @@ class BareParams:
         dico : dict
             dictionary
         """
-        forbiden_keys = ["w_c", "w_power_fact", "field_0", "spec_0", "w", "t", "z_targets", "l"]
+        forbiden_keys = [
+            "w_c",
+            "w_power_fact",
+            "field_0",
+            "spec_0",
+            "w",
+            "t",
+            "z_targets",
+            "l",
+            "alpha",
+        ]
         types = (np.ndarray, float, int, str, list, tuple, dict)
         out = {}
         for key, value in dico.items():