Repository: radis/radis (language: python)

Issue:
Intensity units consistent with the wavespace unit
Turning an insightful comment from @AllanHOlesenBW  into a Todo/Request. 

---

**in short: I think it would be beneficial if the default wunit and Iunit in a spectrum would correspond to each other so a new user can jump directly to work with the units he is given.**

---

> Context :
> 
> My main expectation is that the unit describing the wavenum/wavelength/frequency corresponds to the dSomething part of the unit for the radiance/absorbance/whatever property. If this expectation is not met, I cannot create an integral of the property over the wavelength range. And the area in a plot will not represent anything meaningful.
> 
> The first place where this expectation is not met is rather obvious. It is not directly related to this issue #215, but I find it an important part of the context:
> 
> **If I create a spectrum with calc_spectrum, I get a wunit of cm-1 and a Iunit of mW/cm2/sr/nm.**
> 
> To me, this is meaningless. The '/nm' describes how many mW/cm2/sr I will get over a given wavelength range in the spectrum. But this spectrum is identified by wavenums in cm-1, not by wavelengths in nm.
> 
> Not only does this annoy me visually - it also means that I cannot calculate an integral using the the two units.
> 
> So I have to either rescale the wunit to nm or convert the Iunit to mW/cm2/sr/cm-1.
> 
> This was not intuitive to me when I started using the software. I spent a lot of time exploring my options for creating the spectrums with corresponding units so I didn't have to convert them afterwards. And when I did not succeed in that, I started learning how I could force the software to plot and/or export arrays with (to me) meaningful unit combinations. I will not say that this effort was wasted, but it  stalled me for several days and prevented me from using the software for something productive.
> 
> The non-corresponding choice of default units also creates extra work in the software for example in Spectrum.get_power() where the software has to go through an extra conversion step to create the integral.

_Originally posted by @AllanHOlesenBW in https://github.com/radis/radis/issues/215#issuecomment-1112985821_







----

### Implementation  notes

- [ ] default RADIS calculation unit is defined in https://github.com/radis/radis/blob/develop/radis/lbl/base.py#L107

Make sure radiance_noslit() is properly calculated if default unit changes  : 
- [ ] equilibirium : (https://github.com/radis/radis/blob/develop/radis/lbl/factory.py#L820) --> @erwanp  : calls calc_radiance() which won't work with /cm-1 
- [ ] non_eq_spectrum() converts to proper unit here : https://github.com/radis/radis/blob/develop/radis/lbl/factory.py#L1671
- [ ] GPU code uses same routine as equilibium https://github.com/radis/radis/blob/a0c9bf3d9444b05ad4196959398eddf1879f8d40/radis/lbl/factory.py#L1156
 - [ ] GPU_interactive may need a deeper look  


Discussion / hints:

Fix patch:
diff --git a/docs/_templates/buttons.html b/docs/_templates/buttons.html
index bfff07068..bce722e94 100644
--- a/docs/_templates/buttons.html
+++ b/docs/_templates/buttons.html
@@ -10,6 +10,12 @@
 
 
 <p>
+</a>
+<a href="https://radis.github.io/radis-lab/">
+    <img class="circle-icon"
+         src="{{ pathto('radis_ico.png', 1) }}"
+         alt="Radis-Website">
+</a>
 <a href="https://github.com/radis/radis">
     <img class="circle-icon"
          src="{{ pathto('_static/assets/general/button_github.svg', 1) }}"
diff --git a/docs/conf.py b/docs/conf.py
index 1710f9346..c5d5ccdb4 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -263,6 +263,7 @@ def setup(app):
     "fixed_sidebar": True,
     "analytics_id": "UA-113616205-1",
     "link": "#7A306C",
+    "extra_nav_links": {"RADIS Website": "https://radis.github.io/"},
 }
 
 # Add any paths that contain custom themes here, relative to this directory.
@@ -436,7 +437,7 @@ def setup(app):
         "RADIS Documentation",
         author,
         "RADIS",
-        "A fast line-by-line code for high resolution infrared molecular spectra",
+        "A fast line-by-line code for high resolution infrared molecular spectra: https://radis.github.io/",
         "Miscellaneous",
     ),
 ]
diff --git a/docs/spectrum/spectrum.rst b/docs/spectrum/spectrum.rst
index ad3325ab4..6bf1098cb 100644
--- a/docs/spectrum/spectrum.rst
+++ b/docs/spectrum/spectrum.rst
@@ -338,6 +338,11 @@ method to retrieve the quantity un the units you want::
     _, R = s.get('radiance_noslit', wunit='nm', Iunit='W/cm2/sr/nm',
                  medium='air')
 
+Use with `return_units` to get dimensioned Astropy Quantities ::
+
+    w, R  = s.get('radiance_noslit', return_units=True)
+    # w, R are astropy quantities
+
 See :ref:`spectral arrays <label_spectral _arrays>` for the list
 of spectral arrays.
 
diff --git a/radis/lbl/bands.py b/radis/lbl/bands.py
index bcc969d8a..0df41e63b 100644
--- a/radis/lbl/bands.py
+++ b/radis/lbl/bands.py
@@ -65,11 +65,12 @@
     from .loader import KNOWN_DBFORMAT, KNOWN_LVLFORMAT
 except ImportError:  # if ran from here
     from radis.lbl.loader import KNOWN_DBFORMAT, KNOWN_LVLFORMAT
+
 from radis.misc.basics import all_in, is_float
 from radis.misc.progress_bar import ProgressBar
 from radis.misc.warning import reset_warnings
 from radis.phys.constants import k_b
-from radis.phys.units import convert_emi2nm, convert_rad2nm
+from radis.phys.units import convert_universal
 from radis.phys.units_astropy import convert_and_strip_units
 from radis.spectrum.equations import calc_radiance
 from radis.spectrum.spectrum import Spectrum
@@ -321,6 +322,7 @@ def eq_bands(
                 Tgas,
                 unit=self.units["radiance_noslit"],
             )
+            assert self.units["abscoeff"] == "cm-1"
 
             # ----------------------------- Export:
 
@@ -354,34 +356,30 @@ def add_attr(attr):
                 add_attr("band_htrn")
                 add_attr("viblvl_l")
                 add_attr("viblvl_u")
+
+            quantities = {
+                "wavenumber": wavenumber,
+                "abscoeff": abscoeff,
+                "absorbance": absorbance,
+                "emissivity_noslit": emissivity_noslit,
+                "transmittance_noslit": transmittance_noslit,
+                "radiance_noslit": radiance_noslit,
+            }
+            conditions["default_output_unit"] = self.input_wunit
+
             s = Spectrum(
-                quantities={
-                    "abscoeff": (wavenumber, abscoeff),
-                    "absorbance": (wavenumber, absorbance),
-                    "emissivity_noslit": (wavenumber, emissivity_noslit),
-                    "transmittance_noslit": (wavenumber, transmittance_noslit),
-                    # (mW/cm2/sr/nm)
-                    "radiance_noslit": (wavenumber, radiance_noslit),
-                },
+                quantities=quantities,
                 conditions=conditions,
                 populations=populations,
                 lines=lines,
                 units=self.units,
                 cond_units=self.cond_units,
-                wunit=self.params.waveunit,  # cm-1
                 name=band,
                 # dont check input (much faster, and Spectrum
-                check_wavespace=False,
                 # is freshly baken so probably in a good format
+                check_wavespace=False,
             )
 
-            #            # update database if asked so
-            #            if self.autoupdatedatabase:
-            #                self.SpecDatabase.add(s)
-            #                                                     # Tvib=Trot=Tgas... but this way names in a database
-            #                                                     # generated with eq_spectrum are consistent with names
-            #                                                     # in one generated with non_eq_spectrum
-
             s_bands[band] = s
 
             pb.update(i)  # progress bar
@@ -664,17 +662,27 @@ def non_eq_bands(
                 radiance_noslit[b] = emisscoeff[b] * path_length
             else:
                 # Note that for k -> 0,
-                radiance_noslit = emisscoeff * path_length  # (mW/sr/cm2/cm-1)
-
-            # Convert `radiance_noslit` from (mW/sr/cm2/cm-1) to (mW/sr/cm2/nm)
-            radiance_noslit = convert_rad2nm(
-                radiance_noslit, wavenumber, "mW/sr/cm2/cm-1", "mW/sr/cm2/nm"
+                radiance_noslit = emisscoeff * path_length  # (mW/cm2/sr/cm-1)
+
+            # Convert `radiance_noslit` from (mW/sr/cm2/cm-1) to output unit
+            radiance_noslit = convert_universal(
+                radiance_noslit,
+                from_unit="mW/cm2/sr/cm-1",
+                to_unit=self.units["radiance_noslit"],
+                wavenum=wavenumber,
+                per_nm_is_like="mW/cm2/sr/nm",
+                per_cm_is_like="mW/cm2/sr/cm-1",
             )
-            # Convert 'emisscoeff' from (mW/sr/cm3/cm-1) to (mW/sr/cm3/nm)
-            emisscoeff = convert_emi2nm(
-                emisscoeff, wavenumber, "mW/sr/cm3/cm-1", "mW/sr/cm3/nm"
+            # Convert 'emisscoeff' from (mW/sr/cm3/cm-1) to output unit
+            emisscoeff = convert_universal(
+                emisscoeff,
+                from_unit="mW/cm3/sr/cm-1",
+                to_unit=self.units["emisscoeff"],
+                wavenum=wavenumber,
+                per_nm_is_like="mW/cm3/sr/nm",
+                per_cm_is_like="mW/cm3/sr/cm-1",
             )
-            # Note: emissivity not defined under non equilibrium
+            assert self.units["abscoeff"] == "cm-1"
 
             # ----------------------------- Export:
 
@@ -710,32 +718,30 @@ def add_attr(attr):
             add_attr("band_htrn")
             add_attr("viblvl_l")
             add_attr("viblvl_u")
+
+            quantities = {
+                "wavenumber": wavenumber,
+                "abscoeff": abscoeff,
+                "absorbance": absorbance,
+                "emisscoeff": emisscoeff,
+                "transmittance_noslit": transmittance_noslit,
+                "radiance_noslit": radiance_noslit,
+            }
+            conditions["default_output_unit"] = self.input_wunit
+
             s = Spectrum(
-                quantities={
-                    "abscoeff": (wavenumber, abscoeff),
-                    "absorbance": (wavenumber, absorbance),
-                    # (mW/cm3/sr/nm)
-                    "emisscoeff": (wavenumber, emisscoeff),
-                    "transmittance_noslit": (wavenumber, transmittance_noslit),
-                    # (mW/cm2/sr/nm)
-                    "radiance_noslit": (wavenumber, radiance_noslit),
-                },
+                quantities=quantities,
                 conditions=conditions,
                 populations=populations,
                 lines=lines,
                 units=self.units,
                 cond_units=self.cond_units,
-                wunit=self.params.waveunit,  # cm-1
                 name=band,
                 # dont check input (much faster, and Spectrum
-                check_wavespace=False,
                 # is freshly baken so probably in a good format
+                check_wavespace=False,
             )
 
-            #            # update database if asked so
-            #            if self.autoupdatedatabase:
-            #                self.SpecDatabase.add(s, add_info=['Tvib', 'Trot'], add_date='%Y%m%d')
-
             s_bands[band] = s
 
             pb.update(i)  # progress bar
diff --git a/radis/lbl/base.py b/radis/lbl/base.py
index 07e7147b8..b43d58b66 100644
--- a/radis/lbl/base.py
+++ b/radis/lbl/base.py
@@ -96,26 +96,29 @@
 
 class BaseFactory(DatabankLoader):
 
-    # Output units
     units = {
+        "waverange": "cm-1",  # should be "cm-1", "nm" [assumes in air], "nm_vac" [in vacuum]. Note that Radis caluclations will still happen in cm-1, units are converted at the export only.
         "absorbance": "",
         "abscoeff": "cm-1",
         "abscoeff_continuum": "cm-1",
-        # TODO: deal with case where 'cm-1' is given as input for a Spectrum
-        # (write a cast_unit of some kind)
-        # different in Specair (mw/cm2/sr) because slit
-        "radiance": "mW/cm2/sr/nm",
+        # # different in Specair (mw/cm2/sr) because slit
+        # "radiance": "mW/cm2/sr/cm-1",
         # function is not normalized to conserve energy
-        "radiance_noslit": "mW/cm2/sr/nm",  # it's actually a spectral radiance
-        "emisscoeff": "mW/cm3/sr/nm",
-        "emisscoeff_continuum": "mW/cm3/sr/nm",
-        "emissivity": "",
+        "radiance_noslit": "mW/cm2/sr/cm-1",  # it's actually a spectral radiance
+        "emisscoeff": "mW/cm3/sr/cm-1",
+        "emisscoeff_continuum": "mW/cm3/sr/cm-1",
+        # "emissivity": "",
         "emissivity_noslit": "",
-        "transmittance": "",
+        # "transmittance": "",
         "transmittance_noslit": "",
     }
+    """
+    Default output units
+    ... may be changed at the initialisation of the SpectrumFactory, for instance
+    ... if user gives wavelength units we want to return radiance in
+    ... "mW/cm2/sr/nm" units for consistency
+    """
 
-    # Calculation Conditions units
     cond_units = {
         "wavenum_min": "cm-1",
         "wavenum_max": "cm-1",
@@ -138,6 +141,9 @@ class BaseFactory(DatabankLoader):
         # The later is never stored in Factory, but exported in Spectrum at the end of the calculation
         "calculation_time": "s",
     }
+    """
+    Calculation Conditions units
+    """
 
     def __init__(self):
         """
@@ -149,6 +155,10 @@ def __init__(self):
 
         super(BaseFactory, self).__init__()  # initialize parent class
 
+        # Make units specific to this BaseFactory instance :
+        self.units = BaseFactory.units.copy()
+        self.cond_units = BaseFactory.cond_units.copy()
+
         # Define variable names
         # ... Note: defaults values are overwritten by SpectrumFactory input
         # ... values here are just to help autocompletion tools
@@ -3484,7 +3494,7 @@ def plot_populations(self, what="vib", isotope=None, nfig=None):
         fix_style()
 
 
-def get_waverange(
+def get_wavenumber_range(
     wmin=None,
     wmax=None,
     wunit=Default("cm-1"),
@@ -3493,6 +3503,7 @@ def get_waverange(
     wavelength_min=None,
     wavelength_max=None,
     medium="air",
+    return_input_wunit=False,
 ):
     """Returns wavenumber based on whatever input was given: either ν_min,
     ν_max directly, or λ_min, λ_max  in the given propagation ``medium``.
@@ -3513,11 +3524,15 @@ def get_waverange(
         wavelengths in given ``medium``
     Returns
     -------
-    wavenum_min, wavenum_max,: float
+    wavenum_min, wavenum_max: float
         wavenumbers
+    input_wunit: 'nm', 'nm_vac', 'cm-1'
+        in which wavespace was the input given before conversion (used to keep
+        default plot/get consistent with input units)
     """
 
     # Checking consistency of all input variables
+    assert medium in ["air", "vacuum"]
 
     w_present = wmin is not None and wmax is not None
     wavenum_present = wavenum_min is not None and wavenum_max is not None
@@ -3562,6 +3577,7 @@ def get_waverange(
             )
             assert wavelength_min.unit.is_equivalent(u.m)
             assert wavelength_max.unit.is_equivalent(u.m)
+        input_wunit = "wavelength"
 
     if wavenum_min is not None or wavenum_max is not None:
         assert wavenum_min is not None and wavenum_max is not None
@@ -3574,13 +3590,16 @@ def get_waverange(
             )
             assert wavenum_min.unit.is_equivalent(1 / u.cm)
             assert wavenum_max.unit.is_equivalent(1 / u.cm)
+        input_wunit = "wavenumber"
 
     if isinstance(wmin, u.Quantity) or isinstance(wmax, u.Quantity):
         assert wmin is not None and wmax is not None
         assert isinstance(wmin, u.Quantity) and isinstance(wmax, u.Quantity)
-        assert wmin.unit.is_equivalent(u.m) or wmin.unit.is_equivalent(1 / u.m)
-        assert wmax.unit.is_equivalent(u.m) or wmax.unit.is_equivalent(1 / u.m)
-        assert wmin.unit.is_equivalent(wmax.unit)
+        if wmin.unit.is_equivalent(u.m):
+            assert wmax.unit.is_equivalent(u.m)
+        else:
+            assert wmin.unit.is_equivalent(1 / u.m)
+            assert wmax.unit.is_equivalent(1 / u.m)
         if not isinstance(wunit, Default):
             if not wmin.unit.is_equivalent(u.Unit(wunit)):
                 raise ValueError("Conflicting units passed for wmin/wmax and wunit")
@@ -3591,20 +3610,26 @@ def get_waverange(
             if wmin.unit.is_equivalent(u.m):
                 wavelength_min = wmin
                 wavelength_max = wmax
+                input_wunit = "wavelength"
             else:
                 wavenum_min = wmin
                 wavenum_max = wmax
+                input_wunit = "wavenumber"
         else:
             if isinstance(wunit, Default):
                 wavenum_min = wmin * u.Unit(wunit.value)
                 wavenum_max = wmax * u.Unit(wunit.value)
+                assert wunit.value == "cm-1"
+                input_wunit = "wavenumber"
             else:
                 if u.Unit(wunit).is_equivalent(u.m):
                     wavelength_min = wmin * u.Unit(wunit)
                     wavelength_max = wmax * u.Unit(wunit)
+                    input_wunit = "wavelength"
                 else:
                     wavenum_min = wmin * u.Unit(wunit)
                     wavenum_max = wmax * u.Unit(wunit)
+                    input_wunit = "wavenumber"
 
     # We now have wavenum_min/max, or wavelength_min/max defined. Let's convert these to cm-1 (warning: propagating medium is required if we start from wavelengths!)
     if wavenum_min is not None or wavenum_max is not None:
@@ -3627,7 +3652,18 @@ def get_waverange(
             wavenum_min = nm2cm(wavelength_max)
             wavenum_max = nm2cm(wavelength_min)
 
-    return wavenum_min, wavenum_max
+    if return_input_wunit:
+        # get cm-1, nm, nm_vac
+        if input_wunit == "wavenumber":
+            input_wunit = "cm-1"
+        else:
+            assert input_wunit == "wavelength"
+            input_wunit = "nm" if medium == "air" else "nm_vac"
+
+        assert input_wunit in ["cm-1", "nm", "nm_vac"]
+        return wavenum_min, wavenum_max, input_wunit
+    else:
+        return wavenum_min, wavenum_max
 
 
 def linestrength_from_Einstein(
diff --git a/radis/lbl/calc.py b/radis/lbl/calc.py
index ee81ba170..64379ff32 100644
--- a/radis/lbl/calc.py
+++ b/radis/lbl/calc.py
@@ -23,11 +23,11 @@
 # from radis.db.classes import M
 
 try:  # Proper import
-    from .base import get_waverange
+    from .base import get_wavenumber_range
     from .factory import SpectrumFactory
 except ImportError:  # if ran from here
     from radis.lbl.factory import SpectrumFactory
-    from radis.lbl.base import get_waverange
+    from radis.lbl.base import get_wavenumber_range
 
 from radis.misc.basics import all_in
 from radis.misc.utils import Default
@@ -361,7 +361,7 @@ def calc_spectrum(
     # ... wavelengths / wavenumbers
 
     # Get wavenumber, based on whatever was given as input.
-    wavenum_min, wavenum_max = get_waverange(
+    wavenum_min, wavenum_max, input_wunit = get_wavenumber_range(
         wmin,
         wmax,
         wunit,
@@ -370,6 +370,7 @@ def calc_spectrum(
         kwargs.pop("wavelength_min") if "wavelength_min" in kwargs else None,
         kwargs.pop("wavelength_max") if "wavelength_max" in kwargs else None,
         medium,
+        return_input_wunit=True,
     )
 
     # Deal with Multi-molecule mode:
@@ -510,6 +511,7 @@ def _check_molecules_are_consistent(
         generated_spectrum = _calc_spectrum_one_molecule(
             wavenum_min=wavenum_min,
             wavenum_max=wavenum_max,
+            input_wunit=input_wunit,
             Tgas=Tgas,
             Tvib=Tvib,
             Trot=Trot,
@@ -563,6 +565,7 @@ def _check_molecules_are_consistent(
 def _calc_spectrum_one_molecule(
     wavenum_min,
     wavenum_max,
+    input_wunit,
     Tgas,
     Tvib,
     Trot,
@@ -589,7 +592,17 @@ def _calc_spectrum_one_molecule(
     return_factory=False,
     **kwargs,
 ) -> Spectrum:
-    """See :py:func:`~radis.lbl.calc.calc_spectrum`"""
+    """See :py:func:`~radis.lbl.calc.calc_spectrum`
+
+    Parameters
+    ----------
+    input_wunit: 'nm', 'nm_vac', 'cm-1'
+        in which wavespace was the input given before conversion (used to keep
+        default plot/get consistent with input units)
+    """
+
+    # Initialize Factory
+    # ------------------
 
     # Check inputs
 
@@ -653,6 +666,13 @@ def _is_at_equilibrium():
         export_lines=export_lines,
         **kwargs,
     )
+    # Have consistent output units
+    sf.input_wunit = input_wunit
+
+    # Load databank
+    # -------------
+
+    # Get databank
     if (
         databank
         in [
@@ -807,6 +827,8 @@ def _is_at_equilibrium():
     if overpopulation is not None or overpopulation != {}:
         sf.misc.export_rovib_fraction = True  # required to compute Partition fucntions with overpopulation being taken into account
 
+    # Calculate Spectrum
+    # ------------------
     # Use the standard eq_spectrum / non_eq_spectrum functions
     if _equilibrium:
         if mode == "cpu":
diff --git a/radis/lbl/factory.py b/radis/lbl/factory.py
index 1d66486c0..c27283061 100644
--- a/radis/lbl/factory.py
+++ b/radis/lbl/factory.py
@@ -91,16 +91,16 @@
 
 try:  # Proper import
     from .bands import BandFactory
-    from .base import get_waverange
+    from .base import get_wavenumber_range
 except ImportError:  # if ran from here
     from radis.lbl.bands import BandFactory
-    from radis.lbl.base import get_waverange
+    from radis.lbl.base import get_wavenumber_range
 
 from radis.misc.basics import flatten, is_float, is_range, list_if_float, round_off
 from radis.misc.utils import Default
 from radis.phys.constants import k_b
 from radis.phys.convert import conv2
-from radis.phys.units import convert_emi2nm, convert_rad2nm
+from radis.phys.units import convert_universal
 from radis.phys.units_astropy import convert_and_strip_units
 from radis.spectrum.equations import calc_radiance
 from radis.spectrum.spectrum import Spectrum
@@ -459,7 +459,7 @@ def __init__(
         # --------------------
 
         # Get wavenumber, based on whatever was given as input.
-        wavenum_min, wavenum_max = get_waverange(
+        wavenum_min, wavenum_max, input_wunit = get_wavenumber_range(
             wmin,
             wmax,
             wunit,
@@ -468,7 +468,12 @@ def __init__(
             wavelength_min,
             wavelength_max,
             medium,
+            return_input_wunit=True,
         )
+        # ... Make default Spectrum's output unit consistent with the input waverange
+        # ... see https://github.com/radis/radis/issues/456
+        # (note: may be overwritten by user after Factory creation)
+        self.input_wunit = input_wunit
 
         # Storing inital value of wstep if wstep != "auto"
         self._wstep = wstep
@@ -820,6 +825,7 @@ def eq_spectrum(
         radiance_noslit = calc_radiance(
             wavenumber, emissivity_noslit, Tgas, unit=self.units["radiance_noslit"]
         )
+        assert self.units["abscoeff"] == "cm-1"
 
         self.profiler.stop(
             "calc_other_spectral_quan", "Calculated other spectral quantities"
@@ -873,17 +879,16 @@ def eq_spectrum(
 
         # Spectral quantities
         quantities = {
-            "abscoeff": (wavenumber, abscoeff),
-            "absorbance": (wavenumber, absorbance),
-            "emissivity_noslit": (wavenumber, emissivity_noslit),
-            "transmittance_noslit": (wavenumber, transmittance_noslit),
-            # (mW/cm2/sr/nm)
-            "radiance_noslit": (wavenumber, radiance_noslit),
+            "wavenumber": wavenumber,
+            "abscoeff": abscoeff,
+            "absorbance": absorbance,
+            "emissivity_noslit": emissivity_noslit,
+            "transmittance_noslit": transmittance_noslit,
+            "radiance_noslit": radiance_noslit,
         }
         if I_continuum is not None and self._export_continuum:
-            quantities.update(
-                {"abscoeff_continuum": (wavenumber, I_continuum * density)}
-            )
+            quantities.update({"abscoeff_continuum": I_continuum * density})
+        conditions["default_output_unit"] = self.input_wunit
 
         # Store results in Spectrum class
         s = Spectrum(
@@ -893,13 +898,13 @@ def eq_spectrum(
             lines=lines,
             units=self.units,
             cond_units=self.cond_units,
-            wunit=self.params.waveunit,  # cm-1
             # dont check input (much faster, and Spectrum
-            check_wavespace=False,
             # is freshly baken so probably in a good format
+            check_wavespace=False,
             name=name,
             references=dict(self.reftracker),
         )
+        # OPTION 2.  Change a posteriori using a Spectrum.method. More universal. Can it be slower?
 
         # update database if asked so
         if self.autoupdatedatabase:
@@ -1156,6 +1161,7 @@ def eq_spectrum_gpu(
         radiance_noslit = calc_radiance(
             wavenumber, emissivity_noslit, Tgas, unit=self.units["radiance_noslit"]
         )
+        assert self.units["abscoeff"] == "cm-1"
 
         self.profiler.stop(
             "calc_other_spectral_quan", "Calculated other spectral quantities"
@@ -1203,15 +1209,16 @@ def eq_spectrum_gpu(
 
         # Spectral quantities
         quantities = {
-            "abscoeff": (wavenumber, abscoeff),
-            "absorbance": (wavenumber, absorbance),
-            "emissivity": (wavenumber, emissivity),
-            "emissivity_noslit": (wavenumber, emissivity_noslit),
-            "transmittance_noslit": (wavenumber, transmittance_noslit),
-            # (mW/cm2/sr/nm)
-            "radiance_noslit": (wavenumber, radiance_noslit),
-            "transmittance": (wavenumber, transmittance),
+            "wavenumber": wavenumber,
+            "abscoeff": abscoeff,
+            "absorbance": absorbance,
+            "emissivity": emissivity,
+            "emissivity_noslit": emissivity_noslit,
+            "transmittance_noslit": transmittance_noslit,
+            "radiance_noslit": radiance_noslit,
+            "transmittance": transmittance,
         }
+        conditions["default_output_unit"] = self.input_wunit
 
         # Store results in Spectrum class
         s = Spectrum(
@@ -1220,7 +1227,6 @@ def eq_spectrum_gpu(
             conditions=conditions,
             lines=lines,
             cond_units=self.cond_units,
-            wunit=self.params.waveunit,  # cm-1
             # dont check input (much faster, and Spectrum
             check_wavespace=False,
             # is freshly baken so probably in a good format
@@ -1645,9 +1651,6 @@ def non_eq_spectrum(
         abscoeff = abscoeff_v * density  # cm-1
         emisscoeff = emisscoeff_v * density  # mW/sr/cm3/cm-1
 
-        # ... # TODO: if the code is extended to multi-species, then density has to be added
-        # ... before lineshape broadening (as it would not be constant for all species)
-
         # get absorbance (technically it's the optical depth `tau`,
         #                absorbance `A` being `A = tau/ln(10)` )
 
@@ -1666,16 +1669,27 @@ def non_eq_spectrum(
             radiance_noslit[b] = emisscoeff[b] * path_length
         else:
             # Note that for k -> 0,
-            radiance_noslit = emisscoeff * path_length  # (mW/sr/cm2/cm-1)
-
-        # Convert `radiance_noslit` from (mW/sr/cm2/cm-1) to (mW/sr/cm2/nm)
-        radiance_noslit = convert_rad2nm(
-            radiance_noslit, wavenumber, "mW/sr/cm2/cm-1", "mW/sr/cm2/nm"
+            radiance_noslit = emisscoeff * path_length  # (mW/cm2/sr/cm-1)
+
+        # Convert `radiance_noslit` from (mW/sr/cm2/cm-1) to output unit
+        radiance_noslit = convert_universal(
+            radiance_noslit,
+            from_unit="mW/cm2/sr/cm-1",
+            to_unit=self.units["radiance_noslit"],
+            wavenum=wavenumber,
+            per_nm_is_like="mW/cm2/sr/nm",
+            per_cm_is_like="mW/cm2/sr/cm-1",
         )
-        # Convert 'emisscoeff' from (mW/sr/cm3/cm-1) to (mW/sr/cm3/nm)
-        emisscoeff = convert_emi2nm(
-            emisscoeff, wavenumber, "mW/sr/cm3/cm-1", "mW/sr/cm3/nm"
+        # Convert 'emisscoeff' from (mW/sr/cm3/cm-1) to output unit
+        emisscoeff = convert_universal(
+            emisscoeff,
+            from_unit="mW/cm3/sr/cm-1",
+            to_unit=self.units["emisscoeff"],
+            wavenum=wavenumber,
+            per_nm_is_like="mW/cm3/sr/nm",
+            per_cm_is_like="mW/cm3/sr/cm-1",
         )
+        assert self.units["abscoeff"] == "cm-1"
 
         self.profiler.stop(
             "calc_other_spectral_quan", "Calculated other spectral quantities"
@@ -1729,21 +1743,21 @@ def non_eq_spectrum(
 
         # Spectral quantities
         quantities = {
-            "abscoeff": (wavenumber, abscoeff),
-            "absorbance": (wavenumber, absorbance),
-            # (mW/cm3/sr/nm)
-            "emisscoeff": (wavenumber, emisscoeff),
-            "transmittance_noslit": (wavenumber, transmittance_noslit),
-            # (mW/cm2/sr/nm)
-            "radiance_noslit": (wavenumber, radiance_noslit),
+            "wavenumber": wavenumber,
+            "abscoeff": abscoeff,
+            "absorbance": absorbance,
+            "emisscoeff": emisscoeff,
+            "transmittance_noslit": transmittance_noslit,
+            "radiance_noslit": radiance_noslit,
         }
         if k_continuum is not None and self._export_continuum:
             quantities.update(
                 {
-                    "abscoeff_continuum": (wavenumber, k_continuum * density),
-                    "emisscoeff_continuum": (wavenumber, j_continuum * density),
+                    "abscoeff_continuum": k_continuum * density,
+                    "emisscoeff_continuum": j_continuum * density,
                 }
             )
+        conditions["default_output_unit"] = self.input_wunit
 
         # Store results in Spectrum class
         s = Spectrum(
@@ -1753,7 +1767,6 @@ def non_eq_spectrum(
             lines=lines,
             units=self.units,
             cond_units=self.cond_units,
-            wunit=self.params.waveunit,  # cm-1
             # dont check input (much faster, and Spectrum
             check_wavespace=False,
             # is freshly baken so probably in a good format
diff --git a/radis/lbl/overp.py b/radis/lbl/overp.py
index 508b67bc2..258b1807a 100644
--- a/radis/lbl/overp.py
+++ b/radis/lbl/overp.py
@@ -659,22 +659,17 @@ def rescale_updown_levels(
 
     Parameters
     ----------
-
     new_nu: float
         new upper state mole fraction
-
     old_nu: float
         current upper state mole fraction
-
     new_nl: float
         new lower state mole fraction
-
     old_nl: float
         current lower state mole fraction
 
     Other Parameters
     ----------------
-
     force: boolean
         if False, won't allow rescaling to 0 (not to loose information).
         Default ``False``
@@ -778,6 +773,7 @@ def rescale_updown_levels(
             units,
             extra,
             true_path_length,
+            assume_equilibrium=False,
         )
 
     if "emisscoeff" in recompute:
diff --git a/radis/los/slabs.py b/radis/los/slabs.py
index 31227d596..24c1a2091 100644
--- a/radis/los/slabs.py
+++ b/radis/los/slabs.py
@@ -137,7 +137,6 @@ def SerialSlabs(*slabs, **kwargs) -> Spectrum:
     # Check inputs, get defaults
     resample_wavespace = kwargs.pop("resample", "never")  # default 'never'
     out_of_bounds = kwargs.pop("out", "nan")  # default 'nan'
-    verbose = kwargs.pop("verbose", False)  # type: bool
     modify_inputs = kwargs.pop("modify_inputs", False)  # type: bool
     if len(kwargs) > 0:
         raise ValueError("Unexpected input: {0}".format(list(kwargs.keys())))
@@ -198,57 +197,50 @@ def SerialSlabs(*slabs, **kwargs) -> Spectrum:
         # Get all data
         # -------------
 
-        I, In, T, Tn = None, None, None, None
-
         # To make it easier, the radiative transfer equation is solved with 'radiance_noslit' and
         # 'transmittance_noslit' only. Here we first try to get these quantities:
 
+        # ... if not given in the spectra, they are recomputed automatically
+        # ... starting from RADIS 0.12 https://github.com/radis/radis/pull/413
+
         # ... get sn quantities
         try:
-            sn.update("transmittance_noslit", verbose=verbose)
-        except ValueError:
-            pass
-        else:
             Tn = sn.get(
                 "transmittance_noslit",
                 wunit=waveunit,
                 Iunit=unitsn["transmittance_noslit"],
                 copy=False,
             )[1]
-        try:
-            sn.update("radiance_noslit", verbose=verbose)
         except ValueError:
-            pass
-        else:
+            Tn = None
+        try:
             In = sn.get(
                 "radiance_noslit",
                 wunit=waveunit,
                 Iunit=unitsn["radiance_noslit"],
                 copy=False,
             )[1]
+        except ValueError:
+            In = None
         # ... get s quantities
         try:
-            s.update("transmittance_noslit", verbose=verbose)
-        except ValueError:
-            pass
-        else:
             T = s.get(
                 "transmittance_noslit",
                 wunit=waveunit,
                 Iunit=unitsn["transmittance_noslit"],
                 copy=False,
             )[1]
-        try:
-            s.update("radiance_noslit", verbose=verbose)
         except ValueError:
-            pass
-        else:
+            T = None
+        try:
             I = s.get(
                 "radiance_noslit",
                 wunit=waveunit,
                 Iunit=unitsn["radiance_noslit"],
                 copy=False,
             )[1]
+        except ValueError:
+            I = None
 
         # Solve radiative transfer equation
         # ---------------------------------
@@ -281,10 +273,14 @@ def SerialSlabs(*slabs, **kwargs) -> Spectrum:
                 conditions[cond] = s.conditions[cond] + sn.conditions[cond]
                 if cond in s.cond_units and cond in sn.cond_units:
                     assert s.cond_units == sn.cond_units
+        # ... for the parmaeters below, if N/A, remove the value
+        for cond in ["default_output_unit"]:
+            if cond in conditions and conditions[cond] == "N/A":
+                conditions.pop(cond)
         cond_units = intersect(s.cond_units, sn.cond_units)
 
-        # Update references
-        # -----------------
+        # Update bibliography references
+        # ------------------------------
         # (just add everything)
         references = {**s.references, **sn.references}
 
@@ -425,7 +421,7 @@ def same_wavespace(wl):
         # defined on the same range)
         slabsk = [s for s in slabs if k in s.get_vars()]  # note that these are
         # references to the actual Spectrum copy
-        wl = [s.get(k, wunit=waveunit, copy=False)[0] for s in slabsk]
+        wl = [s.get(k, wunit=waveunit, Iunit=s.units[k], copy=False)[0] for s in slabsk]
         if not same_wavespace(wl):
             # resample slabs if allowed
             if resample_wavespace == "never":
@@ -666,6 +662,10 @@ def MergeSlabs(*slabs, **kwargs) -> Spectrum:
             conditions = intersect(conditions, s.conditions)
             cond_units = intersect(cond_units, s.cond_units)
             # units = intersect(units0, s.units)  # we're actually using [slabs0].units insteads
+        # ... for the parmaeters below, if N/A, remove the value
+        for cond in ["default_output_unit"]:
+            if cond in conditions and conditions[cond] == "N/A":
+                conditions.pop(cond)
         # ... Add extensive parameters
         for cond in ["molecule"]:  # list of all
             if in_all(cond, [s.conditions for s in slabs]):
@@ -696,8 +696,8 @@ def MergeSlabs(*slabs, **kwargs) -> Spectrum:
         #  "mole_fractions":{'MergeSlabs':dict},  # make a dict, same for mole fractions?
         #  }
 
-        # Update references
-        # -----------------
+        # Update bibliography references
+        # ------------------------------
         # (just add everything)
         references = slabs[0].references
         for s in slabs[1:]:
@@ -714,19 +714,10 @@ def MergeSlabs(*slabs, **kwargs) -> Spectrum:
             recompute.append("absorbance")
             recompute.append("transmittance_noslit")
 
-        # To make it easier, we start from abscoeff and emisscoeff of all slabs
-        # Let's recompute them all
-        # TODO: if that changes the initial Spectra, maybe we should just work on copies
-        for s in slabs:
-            if "abscoeff" in recompute and not "abscoeff" in list(s._q.keys()):
-                s.update("abscoeff", verbose=False)
-                # that may crash if Spectrum doesnt have the correct inputs.
-                # let update() handle that
-            if "emisscoeff" in recompute and not "emisscoeff" in list(s._q.keys()):
-                s.update("emisscoeff", verbose=False)
-                # same
-
         # %% Calculate total emisscoeff and abscoeff
+        # To make it easier, we start from abscoeff and emisscoeff of all slabs
+        # ... if not given in the spectra, they are recomputed automatically
+        # ... starting from RADIS 0.12 https://github.com/radis/radis/pull/413
         added = {}
 
         # ... absorption coefficient (cm-1)
@@ -736,7 +727,9 @@ def MergeSlabs(*slabs, **kwargs) -> Spectrum:
                 printdbg("... merge: calculating abscoeff k=sum(k_i)")
             abscoeff_eq = np.sum(
                 [
-                    s.get("abscoeff", wunit=waveunit, Iunit=units0["abscoeff"])[1]
+                    s.take("abscoeff").get(
+                        "abscoeff", wunit=waveunit, Iunit=units0["abscoeff"]
+                    )[1]
                     for s in slabs
                 ],
                 axis=0,
@@ -750,7 +743,9 @@ def MergeSlabs(*slabs, **kwargs) -> Spectrum:
                 printdbg("... merge: calculating emisscoeff j=sum(j_i)")
             emisscoeff_eq = np.sum(
                 [
-                    s.get("emisscoeff", wunit=waveunit, Iunit=units0["emisscoeff"])[1]
+                    s.take("emisscoeff").get(
+                        "emisscoeff", wunit=waveunit, Iunit=units0["emisscoeff"]
+                    )[1]
                     for s in slabs
                 ],
                 axis=0,
diff --git a/radis/misc/arrays.py b/radis/misc/arrays.py
index 6cd9e4b7d..04a1cfdee 100644
--- a/radis/misc/arrays.py
+++ b/radis/misc/arrays.py
@@ -293,25 +293,56 @@ def centered_diff(w):
     return (hstack((dw, dw[-1])) + hstack((dw[0], dw))) / 2
 
 
-def evenly_distributed(w, tolerance=1e-5):
+def evenly_distributed(w, atolerance=1e-5):
     """Make sure array `w` is evenly distributed.
 
     Parameters
     ----------
-
     w : numpy array
         array to test
-    tolerance: float
+    atolerance: float
         absolute tolerance
 
     Returns
     -------
-
     out: bool
         ``True`` or ``False`` if ``w`` is evenly distributed.
+
+    See Also
+    --------
+    :py:func:`~radis.misc.arrays.evenly_distributed_fast`
     """
     mean_step = np.diff(w).mean()
-    return (np.abs((np.diff(w) - mean_step)) > tolerance).sum() == 0
+    return (np.abs((np.diff(w) - mean_step)) > atolerance).sum() == 0
+
+
+# TODO : rewrite so it works in one loop only ; Jit it
+
+
+def evenly_distributed_fast(w, rtolerance=1e-5):
+    """Make sure array `w` is evenly distributed, by looking only at the first
+    and last steps. i.e::
+
+        w[-1] - w[-2] == w[1] - w[0]
+
+    Parameters
+    ----------
+    w : numpy array
+        array to test
+    rtolerance: float
+        relative tolerance
+
+    Returns
+    -------
+    out: bool
+        ``True`` or ``False`` if ``w`` is evenly distributed.
+
+    See Also
+    --------
+    :py:func:`~radis.misc.arrays.evenly_distributed`
+    """
+
+    return np.isclose(w[-1] - w[-2], w[1] - w[0], rtol=rtolerance)
 
 
 def anynan(a):
diff --git a/radis/misc/signal.py b/radis/misc/signal.py
index 860a1b33f..9f09db311 100644
--- a/radis/misc/signal.py
+++ b/radis/misc/signal.py
@@ -109,7 +109,7 @@ def resample(
         reverse = True
     else:
         raise ValueError(
-            "Resampling requires wavespace to be sorted. It is not! Use .sort()? "
+            "Resampling requires wavespace to be sorted. It is not! If working with a Spectrum object, use `s.sort()`? "
         )
 
     if reverse:
@@ -249,39 +249,30 @@ def resample_even(
 
     Parameters
     ----------
-
     xspace: array
         space on which vector was generated
-
     vector: array
         quantity to resample
-
     resfactor: float
         increase of resolution. If 1, output vector has the same number of
         points as the input vector. Default 2.
-
     k: int
         order of spline interpolation. 3: cubic, 1: linear. Default 1.
-
     ext: 'error', 'extrapolate', 0
         Controls the value returned for elements of xspace_new not in the interval
         defined by xspace. If 'error', raise a ValueError. If 'extrapolate', well,
         extrapolate. If '0' or 0, then fill with 0. Default 'error'.
-
     energy_threshold: float
         if energy conservation (integrals) is above this threshold, raise an
         error
-
     print_conservation: boolean
         if True, prints energy conservation
 
 
     Returns
     -------
-
     xspace_new: array
         evenly spaced mapping of xspace (same min, same max)
-
     vector_new: array
         resampled vector on evenly spaced array. Number of element is conserved.
 
@@ -291,12 +282,14 @@ def resample_even(
 
     Examples
     --------
-
     Resample a :class:`~radis.spectrum.spectrum.Spectrum` radiance
     on an evenly spaced wavenumber space::
 
         w_nm, I_nm = s.get('radiance')
         w_cm, I_cm = resample_even(nm2cm(w_nm), I_nm)
+
+    You can also use :py:meth:`~radis.spectrum.spectrum.Spectrum.resample_even`
+    directly
     """
 
     # Get new evenly space array
diff --git a/radis/phys/blackbody.py b/radis/phys/blackbody.py
index b79e8f2fb..589305775 100644
--- a/radis/phys/blackbody.py
+++ b/radis/phys/blackbody.py
@@ -31,7 +31,7 @@
 from radis.phys.units import conv2
 
 
-def planck(lmbda, T, eps=1, unit="mW/sr/cm2/nm"):
+def planck(lmbda, T, eps=1, unit="mW/cm2/sr/nm"):
     r"""Planck function for blackbody radiation.
 
     .. math::
@@ -64,13 +64,13 @@ def planck(lmbda, T, eps=1, unit="mW/sr/cm2/nm"):
     )  # S.I  (W.sr-1.m-3)
     iplanck *= 1e-10  # W.sr-1.m-3 >>> mW.sr-1.cm-2.nm-1
 
-    if Q_(unit) != Q_("mW/sr/cm2/nm"):
-        iplanck = conv2(iplanck, "mW/sr/cm2/nm", unit)
+    if Q_(unit) != Q_("mW/cm2/sr/nm"):
+        iplanck = conv2(iplanck, "mW/cm2/sr/nm", unit)
 
     return iplanck
 
 
-def planck_wn(wavenum, T, eps=1, unit="mW/sr/cm2/cm-1"):
+def planck_wn(wavenum, T, eps=1, unit="mW/cm2/sr/cm-1"):
     r"""Planck function for blackbody radiation, wavenumber version.
 
     .. math::
@@ -109,8 +109,8 @@ def planck_wn(wavenum, T, eps=1, unit="mW/sr/cm2/cm-1"):
     # iplanck in erg/s/sr/cm2/cm-1
     iplanck *= 1e-4  # erg/s/sr/cm2/cm-1 > mW/sr/cm^2/cm-1
 
-    if Q_(unit) != Q_("mW/sr/cm2/cm-1"):
-        iplanck = conv2(iplanck, "mW/sr/cm2/cm-1", unit)
+    if Q_(unit) != Q_("mW/cm2/sr/cm-1"):
+        iplanck = conv2(iplanck, "mW/cm2/sr/cm-1", unit)
 
     return iplanck
 
diff --git a/radis/phys/units.py b/radis/phys/units.py
index 1bb76aa1f..33982aca4 100644
--- a/radis/phys/units.py
+++ b/radis/phys/units.py
@@ -70,22 +70,26 @@ def conv2(quantity, fromunit, tounit):
     want to let the users choose another output unit
     """
 
-    import astropy.units as u
+    if fromunit == tounit:
+        return quantity
+    else:
 
-    try:
-        a = quantity * Unit(fromunit)
-        a = a.to(Unit(tounit))
+        import astropy.units as u
+
+        try:
+            a = quantity * Unit(fromunit)
+            a = a.to(Unit(tounit))
 
-    except u.UnitConversionError as err:
-        suggestion = ""
-        if fromunit == "" or tounit == "":
-            suggestion += " Is one of the two arrays normalized ?"
-        raise TypeError(
-            f"Cannot convert quantity `{fromunit}` to the specified unit `{tounit}`. Please check the dimensions."
-            + suggestion
-        ) from err
+        except u.UnitConversionError as err:
+            suggestion = ""
+            if fromunit == "" or tounit == "":
+                suggestion += " Is one of the two arrays normalized ?"
+            raise TypeError(
+                f"Cannot convert quantity `{fromunit}` to the specified unit `{tounit}`. Please check the dimensions."
+                + suggestion
+            ) from err
 
-    return a.value
+        return a.value
 
 
 def is_homogeneous(unit1, unit2):
@@ -97,13 +101,17 @@ def is_homogeneous(unit1, unit2):
         units
     """
 
-    import astropy.units as u
-
-    try:
-        1 * Unit(unit1) + 1 * Unit(unit2)
+    if unit1 == unit2:
         return True
-    except u.UnitConversionError:
-        return False
+    else:
+
+        import astropy.units as u
+
+        try:
+            1 * Unit(unit1) + 1 * Unit(unit2)
+            return True
+        except u.UnitConversionError:
+            return False
 
 
 # %% Emission density conversion functions (with changes from ~1/nm in ~1/cm-1)
@@ -329,9 +337,9 @@ def convert_universal(
     I,
     from_unit,
     to_unit,
-    spec=None,
-    per_nm_is_like="mW/sr/cm2/nm",
-    per_cm_is_like="mW/sr/cm2/cm-1",
+    wavenum,
+    per_nm_is_like="mW/cm2/sr/nm",
+    per_cm_is_like="mW/cm2/sr/cm-1",
 ):
     """Return variable var in whatever unit, and converts to to_unit Also deal
     with cases where var is in ~1/nm (per_nm_is_like) or ~1/cm-1
@@ -343,57 +351,52 @@ def convert_universal(
         variable to get. Usually 'radiance' or 'radiance_noslit'
     to_unit: str
         unit to convert variable to
-
-    Other Parameters
-    ----------------
-    spec: :class:`~radis.spectrum.spectrum.Spectrum` object
-        needed to get wavenumber in case we need to do a change of variable
-        within the integral
-
-    Notes
-    -----
-    wavenumber is needed in case we convert from ~1/nm to ~1/cm-1 (requires
-    a change of variable in the integral)
+    wavenum: array (cm-1)
+        wavenumber is needed in case we convert from ~1/nm to ~1/cm-1 (requires
+        a change of variable in the integral)
     """
-    import astropy.units as u
 
-    Iunit0 = from_unit
-    Iunit = to_unit
-    try:
-        if is_homogeneous(Iunit0, per_nm_is_like) and is_homogeneous(
-            Iunit, per_cm_is_like
-        ):
-            w_cm = spec.get_wavenumber()
-            I = convert_rad2cm(
-                I,
-                w_cm,
-                Iunit0,
-                Iunit,
-                per_nm_is_like=per_nm_is_like,
-                per_cm_is_like=per_cm_is_like,
-            )
-            # note that there may still be a new DimensionalityError
-            # raise here if the input was non-sense.
-        elif is_homogeneous(Iunit0, per_cm_is_like) and is_homogeneous(
-            Iunit, per_nm_is_like
-        ):
-            w_cm = spec.get_wavenumber()
-            I = convert_rad2nm(
-                I,
-                w_cm,
-                Iunit0,
-                Iunit,
-                per_nm_is_like=per_nm_is_like,
-                per_cm_is_like=per_cm_is_like,
-            )
-            # note that there may still be a new DimensionalityError
-            # raise here if the input was non-sense.
-        else:  # general case: just convert
-            I = conv2(I, Iunit0, Iunit)
-    except u.UnitConversionError:
-        raise TypeError(Iunit0, Iunit)
-
-    return I
+    if from_unit == to_unit:
+        return I
+    else:
+
+        import astropy.units as u
+
+        Iunit0 = from_unit
+        Iunit = to_unit
+        try:
+            if is_homogeneous(Iunit0, per_nm_is_like) and is_homogeneous(
+                Iunit, per_cm_is_like
+            ):
+                I = convert_rad2cm(
+                    I,
+                    wavenum,
+                    Iunit0,
+                    Iunit,
+                    per_nm_is_like=per_nm_is_like,
+                    per_cm_is_like=per_cm_is_like,
+                )
+                # note that there may still be a new DimensionalityError
+                # raise here if the input was non-sense.
+            elif is_homogeneous(Iunit0, per_cm_is_like) and is_homogeneous(
+                Iunit, per_nm_is_like
+            ):
+                I = convert_rad2nm(
+                    I,
+                    wavenum,
+                    Iunit0,
+                    Iunit,
+                    per_nm_is_like=per_nm_is_like,
+                    per_cm_is_like=per_cm_is_like,
+                )
+                # note that there may still be a new DimensionalityError
+                # raise here if the input was non-sense.
+            else:  # general case: just convert
+                I = conv2(I, Iunit0, Iunit)
+        except u.UnitConversionError:
+            raise TypeError(Iunit0, Iunit)
+
+        return I
 
 
 # %% Test
diff --git a/radis/spectrum/compare.py b/radis/spectrum/compare.py
index fbfd0ffaf..490f22c4f 100644
--- a/radis/spectrum/compare.py
+++ b/radis/spectrum/compare.py
@@ -270,6 +270,7 @@ def get_residual(
     normalize=False,
     normalize_how="max",
     wunit="default",
+    Iunit="default",
 ) -> float:
     # type: (Spectrum, Spectrum, str, bool, int) -> np.array, np.array
     """Returns L2 norm of ``s1`` and ``s2``
@@ -295,7 +296,7 @@ def get_residual(
     s1, s2: :class:`~radis.spectrum.spectrum.Spectrum` objects
         if not on the same range, ``s2`` is resampled on ``s1``.
     var: str
-        spectral quantity
+        spectral array
     norm: 'L2', 'L1'
         which norm to use
 
@@ -319,11 +320,12 @@ def get_residual(
         noisy experimental spectra. ``'area'`` will normalize the integral to 1.
         ``'mean'`` will normalize by the mean amplitude value
     wunit: str
-        used if normalized is a range
+        used if normalized is a range. If ``'default'``, use first spectrum unit.
+    Iunit: str
+        if ``'default'``, use first spectrum unit
 
     Notes
     -----
-
     0 values for I1 yield nans except if I2 = I1 = 0
 
     when s1 and s2 dont have the size wavespace range, they are automatically
@@ -352,15 +354,17 @@ def get_residual(
             wrange = normalize
         else:
             wrange = ()
-        var, wunit, Iunit = get_default_units(s1, s2, var=var, wunit=wunit)
+        var, wunit, Iunit = get_default_units(s1, s2, var=var, wunit=wunit, Iunit=Iunit)
         s1 = s1.take(var).normalize(
             wrange=wrange, normalize_how=normalize_how, wunit=wunit
         )
         s2 = s2.take(var).normalize(
             wrange=wrange, normalize_how=normalize_how, wunit=wunit
         )
+        Iunit = s1.units[var]
+        assert Iunit == s2.units[var]
 
-    var, wunit, Iunit = get_default_units(s1, s2, var=var, wunit=wunit)
+    var, wunit, Iunit = get_default_units(s1, s2, var=var, wunit=wunit, Iunit=Iunit)
     wdiff, dI = get_diff(
         s1,
         s2,
@@ -391,7 +395,14 @@ def get_residual(
     return output
 
 
-def get_residual_integral(s1: Spectrum, s2: Spectrum, var, ignore_nan=False) -> float:
+def get_residual_integral(
+    s1: Spectrum,
+    s2: Spectrum,
+    var,
+    ignore_nan=False,
+    wunit="default",
+    Iunit="default",
+) -> float:
     # type: (Spectrum, Spectrum, str, bool) -> float
     """Returns integral of the difference between two spectra s1 and s2,
     relatively to the integral of spectrum s1.
@@ -423,6 +434,11 @@ def get_residual_integral(s1: Spectrum, s2: Spectrum, var, ignore_nan=False) ->
         if ``True``, ignore nan in the difference between s1 and s2 (ex: out of bound)
         when calculating residual. Default ``False``. Note: ``get_residual_integral``
         will still fail if there are nan in initial Spectrum.
+    wunit: str
+        If ``'default'``, use first spectrum unit.
+    Iunit: str
+        if ``'default'``, use first spectrum unit
+
 
     Notes
     -----
@@ -446,7 +462,7 @@ def get_residual_integral(s1: Spectrum, s2: Spectrum, var, ignore_nan=False) ->
     :func:`~radis.spectrum.compare.plot_diff`,
     :meth:`~radis.spectrum.spectrum.compare_with`
     """
-    var, wunit, Iunit = get_default_units(s1, s2, var=var)
+    var, wunit, Iunit = get_default_units(s1, s2, var=var, wunit=wunit, Iunit=Iunit)
 
     # Get data
     # ----
@@ -516,17 +532,21 @@ def get_default_units(
         s1.get(var)
     if var not in s2.get_vars():
         s2.get(var)
+    if wunit == "default":
+        wunit = s1.get_waveunit()
+    wunit = cast_waveunit(wunit)
     if Iunit == "default":
         try:
-            Iunit = s1.units[var]
+            s1.units[var]
         except KeyError:  # unit not defined in dictionary
             raise KeyError(
                 "Iunit not defined in spectrum for variable {0}. ".format(var)
                 + "Cant use default unit. Specify unit in s.units['{0}'].".format(var)
             )
-    if wunit == "default":
-        wunit = s1.get_waveunit()
-    wunit = cast_waveunit(wunit)
+        # get default (note that defaut "Iunit" may depend on waveunit used. Let s.get() handle this)
+        _, _, _, Iunit = s1.get(
+            var, wunit=wunit, Iunit="default", copy=False, return_units="as_str"
+        )
 
     return var, wunit, Iunit
 
@@ -1187,7 +1207,7 @@ def _compare_dataframes(df1, df2, name):
             for error message
         """
 
-        from pandas.util.testing import assert_frame_equal
+        from pandas.testing import assert_frame_equal
 
         try:
             assert_frame_equal(
@@ -1254,8 +1274,8 @@ def _display_difference(q, q0):
         # Compare spectral arrays
         # -----------
         for k in vars:
-            w, q = first.get(k, wunit=wunit)
-            w0, q0 = other.get(k, wunit=wunit)
+            w, q = first.get(k, wunit=wunit, Iunit=first.units[k])
+            w0, q0 = other.get(k, wunit=wunit, Iunit=first.units[k])
             if len(w) != len(w0):
                 print(
                     "Wavespaces have different length (for {0}: {1} vs {2})".format(
@@ -1305,8 +1325,8 @@ def _display_difference(q, q0):
         # Compare spectral variables
         # -----------
         for k in vars:
-            w, q = first.get(k, wunit=wunit)
-            w0, q0 = other.get(k, wunit=wunit)
+            w, q = first.get(k, wunit=wunit, Iunit=first.units[k])
+            w0, q0 = other.get(k, wunit=wunit, Iunit=first.units[k])
             if len(w) != len(w0):
                 print(
                     "Wavespaces have different length (for {0}: {1} vs {2})".format(
diff --git a/radis/spectrum/equations.py b/radis/spectrum/equations.py
index 4f6f319f7..7f2a8f238 100644
--- a/radis/spectrum/equations.py
+++ b/radis/spectrum/equations.py
@@ -10,8 +10,10 @@
 in post processing)
 """
 
-from radis.phys.blackbody import planck
+from radis.phys.blackbody import planck, planck_wn
+from radis.phys.constants import k_b
 from radis.phys.convert import cm2nm
+from radis.phys.units import is_homogeneous
 
 # ------- Get radiance from emissivity
 
@@ -34,6 +36,27 @@ def calc_radiance(wavenumber, emissivity, Tgas, unit="mW/sr/cm2/nm"):
     :py:func:`~radis.phys.blackbody.planck`
     """
 
-    radiance = emissivity * planck(cm2nm(wavenumber), Tgas, unit=unit)
+    if is_homogeneous(unit, "mW/sr/cm2/nm"):
+        radiance = emissivity * planck(cm2nm(wavenumber), Tgas, unit=unit)
+    elif is_homogeneous(unit, "mW/sr/cm2/cm-1"):
+        radiance = emissivity * planck_wn(wavenumber, Tgas, unit=unit)
 
     return radiance
+
+
+def abscoeff2xsection(abscoeff_cm1, Tgas_K, pressure_Pa, mole_fraction):
+    """Convert Absorption coefficient (cm-1) to Cross-section  (cm2)
+
+    Parameters
+    ----------
+    abscoeff_cm : cm-1
+    Tgas_K : K
+    pressure_Pa : Pa
+    mole_fraction : 0-1
+
+    Returns
+    -------
+    xsection: cm2
+    """
+
+    return abscoeff_cm1 * (k_b * Tgas_K / mole_fraction / pressure_Pa) * 1e6  # cm2
diff --git a/radis/spectrum/operations.py b/radis/spectrum/operations.py
index d961bac52..c91b98905 100644
--- a/radis/spectrum/operations.py
+++ b/radis/spectrum/operations.py
@@ -113,6 +113,9 @@ def Transmittance(s: Spectrum) -> Spectrum:
     :py:meth:`~radis.spectrum.Spectrum.take`
 
     """
+    warn(
+        "Transmittance(Spectrum) will be deprecated in future RADIS versions. Use Spectrum.take('transmittance')"
+    )
 
     return s.copy(copy_lines=True, quantity="transmittance")
 
@@ -152,6 +155,9 @@ def Transmittance_noslit(s: Spectrum) -> Spectrum:
     :py:meth:`~radis.spectrum.Spectrum.take`
 
     """
+    warn(
+        "Transmittance_noslit(Spectrum) will be deprecated in future RADIS versions. Use Spectrum.take('transmittance_noslit')"
+    )
 
     return s.copy(copy_lines=True, quantity="transmittance_noslit")
 
@@ -191,6 +197,9 @@ def Radiance(s: Spectrum) -> Spectrum:
     :py:meth:`~radis.spectrum.Spectrum.take`
 
     """
+    warn(
+        "Radiance(Spectrum) will be deprecated in future RADIS versions. Use Spectrum.take('radiance')"
+    )
 
     return s.copy(copy_lines=True, quantity="radiance")
 
@@ -230,6 +239,9 @@ def Radiance_noslit(s: Spectrum) -> Spectrum:
     :py:meth:`~radis.spectrum.Spectrum.take`
 
     """
+    warn(
+        "Radiance_noslit(Spectrum) will be deprecated in future RADIS versions. Use Spectrum.take('radiance_noslit')"
+    )
 
     return s.copy(copy_lines=True, quantity="radiance_noslit")
 
@@ -243,20 +255,13 @@ def PerfectAbsorber(s: Spectrum) -> Spectrum:
     Useful to get contribution of different slabs in line-of-sight
     calculations (see example).
 
-    .. note:
-
-        formerly named "Transmittance", but "Transmittance(s)" wouldnt
-        return the Transmittance exactly
-
     Parameters
     ----------
-
     s: Spectrum
         :class:`~radis.spectrum.spectrum.Spectrum` object
 
     Returns
     -------
-
     s_tr: Spectrum
         :class:`~radis.spectrum.spectrum.Spectrum` object, with only the ``transmittance``,
         ``absorbance`` and/or ``abscoeff`` part of ``s``, where ``radiance_noslit`` ,
@@ -264,7 +269,6 @@ def PerfectAbsorber(s: Spectrum) -> Spectrum:
 
     Examples
     --------
-
     Let's say you have a total line of sight::
 
         s_los = s1 > s2 > s3
@@ -279,7 +283,6 @@ def PerfectAbsorber(s: Spectrum) -> Spectrum:
         (s1 > PerfectAbsorber(s2>s3)).plot('radiance_noslit')
 
     See more examples in :ref:`Line-of-Sight module <label_los_index>`
-
     """
 
     s_tr = s.copy()
@@ -308,10 +311,8 @@ def crop(s: Spectrum, wmin=None, wmax=None, wunit=None, inplace=False) -> Spectr
 
     Parameters
     ----------
-
     s: Spectrum object
         object to crop
-
     wmin, wmax: float, or None
         boundaries of spectral range (in ``wunit``)
 
@@ -322,28 +323,23 @@ def crop(s: Spectrum, wmin=None, wmax=None, wunit=None, inplace=False) -> Spectr
 
     Other Parameters
     ----------------
-
     inplace: bool
         if ``True``, modifiy ``s`` directly. Else, returns a copy.
 
     Returns
     -------
-
     s_crop: Spectrum
         a cropped Spectrum.
         if using ``inplace``, then ``s_crop`` and ``s`` are still the same object
 
     Examples
     --------
-
     ::
 
         crop(s, 420, 480, 'nm', 'air')
-
     Or in ``cm-1``::
 
         crop(s, 2000, 2300, 'cm-1')
-
     """
 
     # Check inputs
@@ -492,7 +488,7 @@ def multiply(s, coef, unit=None, var=None, inplace=False):
         s = s.copy(quantity=var)
 
     # Multiply inplace       ( @dev: we have copied already if needed )
-    w, I = s.get(var, wunit=s.get_waveunit(), copy=False)
+    w, I = s.get(var, wunit=s.get_waveunit(), Iunit=s.units[var], copy=False)
     I *= coef  # @dev: updates the Spectrum directly because of copy=False
 
     # Convert Spectrum unit
@@ -555,7 +551,7 @@ def add_constant(s, cst, unit=None, var=None, inplace=False):
         s = s.copy(quantity=var)
 
     # Add inplace       ( @dev: we have copied already if needed )
-    w, I = s.get(var, wunit=s.get_waveunit(), copy=False)
+    w, I = s.get(var, wunit=s.get_waveunit(), Iunit=s.units[var], copy=False)
     I += cst
     # @dev: updates the Spectrum directly because of copy=False
 
@@ -631,7 +627,7 @@ def add_array(s, a, unit=None, var=None, inplace=False):
         s = s.copy(quantity=var)
 
     # Add inplace       ( @dev: we have copied already if needed )
-    w, I = s.get(var, wunit=s.get_waveunit(), copy=False)
+    w, I = s.get(var, wunit=s.get_waveunit(), Iunit=s.units[var], copy=False)
     I += a
     # @dev: updates the Spectrum directly because of copy=False
 
@@ -715,10 +711,8 @@ def sub_baseline(s, left, right, unit=None, var=None, inplace=False):
     if not inplace:
         s = s.copy(quantity=var)
 
-    # @EP:
-
     # Substract inplace       ( @dev: we have copied already if needed )
-    w, I = s.get(var, wunit=s.get_waveunit(), copy=False)
+    w, I = s.get(var, wunit=s.get_waveunit(), Iunit=s.units[var], copy=False)
     I -= np.linspace(left, right, num=np.size(I))
     # @dev: updates the Spectrum directly because of copy=False
 
diff --git a/radis/spectrum/rescale.py b/radis/spectrum/rescale.py
index 17ce2b0bd..6d5408033 100644
--- a/radis/spectrum/rescale.py
+++ b/radis/spectrum/rescale.py
@@ -14,15 +14,15 @@
 
 from warnings import warn
 
+import astropy.units as u
 import numpy as np
-from astropy import units as unit
 from numpy import exp
 from numpy import log as ln
 
 from radis.misc.basics import all_in, any_in, compare_lists
 from radis.misc.debug import printdbg
 from radis.phys.units_astropy import convert_and_strip_units
-from radis.spectrum.equations import calc_radiance
+from radis.spectrum.equations import abscoeff2xsection, calc_radiance
 from radis.spectrum.utils import CONVOLUTED_QUANTITIES, NON_CONVOLUTED_QUANTITIES
 
 # List of all spectral variables sorted by priority during recomputation
@@ -99,13 +99,11 @@ def _build_update_graph(
 
     Returns
     -------
-
     derivation: dict
         {spectral_quantity: [list of combinations of spectral quantities needed to calculate it]}
 
     Examples
     --------
-
     to recompute a Spectrum under nonequilibrium, non optically thin case
     (note that all paths are not there yet)::
 
@@ -143,6 +141,70 @@ def _build_update_graph(
                       ['radiance_noslit'],
                       ['transmittance_noslit']],
          etc. }
+
+    Try it ::
+
+        import radis
+        s = radis.test_spectrum()
+        _build_update_graph(s)
+
+    Outputs ::
+        {'transmittance_noslit': [['absorbance'],
+          ['abscoeff'],
+          ['absorbance'],
+          ['emisscoeff'],
+          ['emissivity_noslit'],
+          ['transmittance'],
+          ['radiance'],
+          ['radiance_noslit'],
+          ['xsection']],
+         'absorbance': [['transmittance_noslit'],
+          ['abscoeff'],
+          ['abscoeff'],
+          ['emisscoeff'],
+          ['emissivity_noslit'],
+          ['transmittance'],
+          ['radiance'],
+          ['radiance_noslit'],
+          ['transmittance_noslit'],
+          ['xsection']],
+         'abscoeff': [['absorbance'],
+          ['xsection'],
+          ['absorbance'],
+          ['emisscoeff'],
+          ['emissivity_noslit'],
+          ['transmittance'],
+          ['radiance'],
+          ['radiance_noslit'],
+          ['transmittance_noslit'],
+          ['xsection']],
+         'radiance_noslit': [['emisscoeff', 'abscoeff'],
+          ['abscoeff'],
+          ['absorbance'],
+          ['emisscoeff'],
+          ['emissivity_noslit'],
+          ['transmittance'],
+          ['radiance'],
+          ['transmittance_noslit'],
+          ['xsection']],
+         'emisscoeff': [['radiance_noslit', 'abscoeff'],
+          ['abscoeff'],
+          ['absorbance'],
+          ['emissivity_noslit'],
+          ['transmittance'],
+          ['radiance'],
+          ['radiance_noslit'],
+          ['transmittance_noslit'],
+          ['xsection']],
+         'xsection': [['abscoeff']],
+         'emissivity_noslit': [['abscoeff'],
+          ['absorbance'],
+          ['emisscoeff'],
+          ['transmittance'],
+          ['radiance'],
+          ['radiance_noslit'],
+          ['transmittance_noslit'],
+          ['xsection']]}
     """
     # Get defaults
     if path_length is None:
@@ -258,9 +320,9 @@ def derives_from(what, *from_keys):
 
     # ------------------------------------------------------------
 
-    if __debug__:
-        printdbg("_build_update_graph: dependence/equivalence tree:")
-        printdbg(derivation)
+    # if __debug__:
+    #     printdbg("_build_update_graph: dependence/equivalence tree:")
+    #     printdbg(derivation)
 
     return derivation
 
@@ -528,10 +590,10 @@ def update(
 
     # Get path length
     if "path_length" in list(spec.conditions.keys()):
-        path_length = spec.conditions["path_length"]
+        path_length_cm = spec.conditions["path_length"]
         true_path_length = True
     else:
-        path_length = 1  # some stuff can still be updated.
+        path_length_cm = 1  # we dont know the exact path length, but some stuff can still be updated.
         true_path_length = False
 
     # Update thermal equilibrium
@@ -564,8 +626,8 @@ def update(
         spec,
         quantity,
         # same path length (absolute value can matter)
-        path_length,
-        path_length,
+        path_length_cm,
+        path_length_cm,
         # same mole fractions (only the ratio matters)
         1,
         1,
@@ -603,11 +665,12 @@ def rescale_abscoeff(
     initial,
     old_mole_fraction,
     new_mole_fraction,
-    old_path_length,
+    old_path_length_cm,
     wunit,
     units,
     extra,
     true_path_length,
+    assume_equilibrium,
 ):
     r"""
 
@@ -621,15 +684,12 @@ def rescale_abscoeff(
 
     References
     ----------
-
     Scale absorption coefficient :
-
     .. math::
 
         k_2=k_1 \frac{x_2}{x_1}
 
     Or recompute from the absorbance:
-
     .. math::
 
         k_2=\frac{A_1}{L_1} \cdot \frac{x_2}{x_1}
@@ -640,7 +700,6 @@ def rescale_abscoeff(
         and non-resonant broadening coefficients are different
 
     Or from the transmittance :
-
     .. math::
 
         k_2=\frac{-\ln(T_1)}{L_1} \cdot \frac{x_2}{x_1}
@@ -652,6 +711,10 @@ def rescale_abscoeff(
         \k_2 =\sigma_1 * \frac{x p}{k_b T}  \cdot \frac{x_2}{x_1}
 
     """
+    if __debug__:
+        printdbg(
+            f"recomputing `abscoeff` from initial {initial} and already rescaled {list(rescaled.keys())}, knowing `old_mole_fraction={old_mole_fraction}`, `new_mole_fraction={new_mole_fraction}`, `old_path_length_cm={old_path_length_cm}`, `true_path_length={true_path_length}`, and `assume_equilibrium={assume_equilibrium}`"
+        )
 
     unit = None
 
@@ -665,18 +728,30 @@ def rescale_abscoeff(
     # First get initial abscoeff
     # ---------------------
     if "abscoeff" in initial:
-        _, abscoeff_init = spec.get("abscoeff", wunit=wunit)
+        _, abscoeff_init = spec.get(
+            "abscoeff", wunit=wunit, Iunit=spec.units["abscoeff"], copy=False
+        )
     elif "absorbance" in initial and true_path_length:  # aka: true path_lengths given
         if __debug__:
             printdbg("... rescale: abscoeff k_1 = A_1/L_1")
-        _, A = spec.get("absorbance", wunit=wunit)
-        abscoeff_init = A / old_path_length  # recalculate initial
+        _, A = spec.get(
+            "absorbance", wunit=wunit, Iunit=spec.units["absorbance"], copy=False
+        )
+        assert (
+            units["absorbance"] == ""
+        )  # TODO: allow different units, setting Iunit above?
+        abscoeff_init = A / old_path_length_cm  # recalculate initial
         unit = "cm-1"
     elif "transmittance_noslit" in initial and true_path_length:
         if __debug__:
             printdbg("... rescale: abscoeff k_1 = -ln(T_1)/L_1")
         # Get abscoeff from transmittance
-        _, T1 = spec.get("transmittance_noslit", wunit=wunit)
+        _, T1 = spec.get(
+            "transmittance_noslit",
+            wunit=wunit,
+            Iunit=spec.units["transmittance_noslit"],
+            copy=False,
+        )
 
         # We'll have a problem if the spectrum is optically thick
         b = T1 == 0  # no transmittance: optically thick mask
@@ -690,7 +765,7 @@ def rescale_abscoeff(
                 raise ValueError(msg)
 
         # Else, let's calculate it
-        abscoeff_init = -ln(T1) / old_path_length  # recalculate initial
+        abscoeff_init = -ln(T1) / old_path_length_cm  # recalculate initial
         unit = "cm-1"
     elif (
         "xsection" in initial
@@ -713,6 +788,18 @@ def rescale_abscoeff(
 
     elif "abscoeff" in extra:  # cant calculate this one but let it go
         abscoeff_init = None
+    elif assume_equilibrium and "emissivity_noslit" in initial:
+        raise NotImplementedError(
+            "Recompute `abscoeff` from `emissivity_noslit` at equilibrium is possible but not implemented"
+        )
+    elif assume_equilibrium and "radiance_noslit" in initial:
+        raise NotImplementedError(
+            "Recompute `abscoeff` from `radiance_noslit` at equilibrium is possible but not implemented"
+        )
+    elif assume_equilibrium and "transmittance_noslit" in initial:
+        raise NotImplementedError(
+            "Recompute `abscoeff` from `transmittance_noslit` at equilibrium is possible but not implemented"
+        )
     else:
         raise ValueError(
             "Can't rescale abscoeff if not all of the following are given : transmittance_noslit ({0}) ".format(
@@ -739,21 +826,27 @@ def rescale_abscoeff(
 
 
 def _recompute_all_at_equilibrium(
-    spec, rescaled, wavenumber, Tgas, new_path_length, true_path_length, units
+    spec, rescaled, wavenumber, Tgas, new_path_length_cm, true_path_length, units
 ):
     """
 
     Parameters
     ----------
-
     rescaled: dict
         abscoeff must be rescaled already
     """
 
     def get_unit_radiance():
-        return spec.units.get("radiance_noslit", "mW/cm2/sr/nm")
+        """In which unit to store radiance : use the existing one by default,
+        or use "mW/cm2/sr/nm if spectrum stored in wavelength; 'mW/cm2/sr/cm-1'
+        if Spectrum stored in wavenumber"""
+        if spec.get_waveunit() == "cm-1":
+            return spec.units.get("radiance_noslit", "mW/cm2/sr/cm-1")
+        else:
+            return spec.units.get("radiance_noslit", "mW/cm2/sr/nm")
 
     def get_unit_emisscoeff(unit_radiance):
+        """Basically, unit of ``unit_radiance`` divided by ``cm``"""
         if "/cm2" in unit_radiance:
             return unit_radiance.replace("/cm2", "/cm3")
         else:
@@ -762,19 +855,20 @@ def get_unit_emisscoeff(unit_radiance):
     assert true_path_length
 
     abscoeff = rescaled["abscoeff"]
-    path_length = new_path_length
+    path_length_cm = new_path_length_cm
 
-    absorbance = abscoeff * path_length
+    absorbance = abscoeff * path_length_cm
 
     # Generate output quantities
+    Iunit_radiance = get_unit_radiance()
     transmittance_noslit = exp(-absorbance)
     emissivity_noslit = 1 - transmittance_noslit
     radiance_noslit = calc_radiance(
-        wavenumber, emissivity_noslit, Tgas, unit=get_unit_radiance()
+        wavenumber, emissivity_noslit, Tgas, unit=Iunit_radiance
     )
     b = transmittance_noslit == 1  # optically thin mask
     emisscoeff = np.empty_like(abscoeff)
-    emisscoeff[b] = radiance_noslit[b] / path_length  # recalculate (opt thin)
+    emisscoeff[b] = radiance_noslit[b] / path_length_cm  # recalculate (opt thin)
     emisscoeff[~b] = (
         radiance_noslit[~b] / (1 - transmittance_noslit[~b]) * abscoeff[~b]
     )  # recalculate (non opt thin)
@@ -791,7 +885,7 @@ def get_unit_emisscoeff(unit_radiance):
     units["absorbance"] = ""
     units["transmittance_noslit"] = ""
     units["emissivity_noslit"] = ""
-    units["radiance_noslit"] = get_unit_radiance()
+    units["radiance_noslit"] = Iunit_radiance
     units["emisscoeff"] = get_unit_emisscoeff(units["radiance_noslit"])
 
     return rescaled, units
@@ -806,7 +900,7 @@ def rescale_emisscoeff(
     initial,
     old_mole_fraction,
     new_mole_fraction,
-    old_path_length,
+    old_path_length_cm,
     optically_thin,
     wunit,
     units,
@@ -817,13 +911,11 @@ def rescale_emisscoeff(
 
     Parameters
     ----------
-
     spec: Spectrum
 
 
     References
     ----------
-
     If optically thin, compute from the radiance :
 
     .. math::
@@ -854,11 +946,14 @@ def rescale_emisscoeff(
 
     unit = None
 
-    def get_unit(unit_radiance):
+    def get_emisscoeff_unit(unit_radiance):
+        """Basically, units of ``unit_radiance`` divided by ``cm``"""
         if "/cm2" in unit_radiance:
             return unit_radiance.replace("/cm2", "/cm3")
         else:
-            return unit_radiance + "/cm"  # will be simplified by Pint afterwards
+            return (
+                unit_radiance + "/cm"
+            )  # will be simplified by Pint/Astropy.units afterwards
 
     # case where we recomputed it already (somehow... ex: no_change signaled)
     if "emisscoeff" in rescaled:
@@ -874,37 +969,39 @@ def get_unit(unit_radiance):
         if __debug__:
             printdbg("... rescale: emisscoeff j1 = j1")
         _, emisscoeff_init = spec.get(
-            "emisscoeff", wunit=wunit, Iunit=units["emisscoeff"]
+            "emisscoeff", wunit=wunit, Iunit=units["emisscoeff"], copy=False
         )
 
     elif "radiance_noslit" in initial and true_path_length and optically_thin:
         if __debug__:
             printdbg("... rescale: emisscoeff j_1 = I_1/L_1")
-        _, I = spec.get("radiance_noslit", wunit=wunit, Iunit=units["radiance_noslit"])
-        emisscoeff_init = I / old_path_length  # recalculate initial
-        unit = get_unit(units["radiance_noslit"])
+        _, I = spec.get(
+            "radiance_noslit", wunit=wunit, Iunit=units["radiance_noslit"], copy=False
+        )
+        emisscoeff_init = I / old_path_length_cm  # recalculate initial
+        unit = get_emisscoeff_unit(units["radiance_noslit"])
 
     elif "radiance_noslit" in initial and true_path_length and "abscoeff" in initial:
         if __debug__:
             printdbg("... rescale: emisscoeff j_1 = k_1*I_1/(1-exp(-k_1*L_1))")
         # get emisscoeff from (initial) abscoeff and (initial) radiance
         _, I = spec.get("radiance_noslit", wunit=wunit, Iunit=units["radiance_noslit"])
-        _, k = spec.get("abscoeff", wunit=wunit, Iunit=units["abscoeff"])
+        _, k = spec.get("abscoeff", wunit=wunit, Iunit=units["abscoeff"], copy=False)
 
         # Recalculate in the optically thin range (T=1) and elsewhere
         b = k == 0  # optically thin mask
         emisscoeff_init = np.empty_like(k)
         # ... optically thin case
         # recalculate (opt thin)
-        emisscoeff_init[b] = I[b] / old_path_length
+        emisscoeff_init[b] = I[b] / old_path_length_cm
 
         # ... non optically thin case:
         # ... recalculate transmittance from abscoeff
-        T_b = exp(-k[~b] * old_path_length)  # recalculate initial
+        T_b = exp(-k[~b] * old_path_length_cm)  # recalculate initial
         # ... and solve the RTE on an homogeneous slab
         # recalculate (non opt thin)
         emisscoeff_init[~b] = k[~b] * I[~b] / (1 - T_b)
-        unit = get_unit(units["radiance_noslit"])
+        unit = get_emisscoeff_unit(unit_radiance=units["radiance_noslit"])
 
     elif (
         "radiance_noslit" in initial
@@ -924,16 +1021,16 @@ def get_unit(unit_radiance):
         emisscoeff_init = np.empty_like(T)
         # ... optically thin case
         # recalculate (opt thin)
-        emisscoeff_init[b] = I[b] / old_path_length
+        emisscoeff_init[b] = I[b] / old_path_length_cm
 
         # ... non optically thin case:
         # ... recalculate abscoeff from transmittance
         T_b = T[~b]
-        k_b = -ln(T_b) / old_path_length
+        k_b = -ln(T_b) / old_path_length_cm  # cm-1
         # ... and solve the RTE on an homogeneous slab
         # recalculate (non opt thin)
         emisscoeff_init[~b] = k_b * I[~b] / (1 - T_b)
-        unit = get_unit(units["radiance_noslit"])
+        unit = get_emisscoeff_unit(unit_radiance=units["radiance_noslit"])
 
     else:
         if optically_thin:
@@ -990,10 +1087,10 @@ def rescale_absorbance(
     initial,
     old_mole_fraction,
     new_mole_fraction,
-    old_path_length,
-    new_path_length,
+    old_path_length_cm,
+    new_path_length_cm,
     waveunit,
-    units,
+    rescaled_units,
     extra,
     true_path_length,
 ):
@@ -1001,12 +1098,10 @@ def rescale_absorbance(
 
     Parameters
     ----------
-
     spec: Spectrum
 
     References
     ----------
-
     Rescale the absorbance:
 
     .. math::
@@ -1044,8 +1139,8 @@ def rescale_absorbance(
     if "absorbance" in rescaled:
         if __debug__:
             printdbg("... rescale: absorbance was scaled already")
-        assert "absorbance" in units
-        return rescaled, units
+        assert "absorbance" in rescaled_units
+        return rescaled, rescaled_units
 
     # Get scaled absorbance directly A1
     # ---------------------------------
@@ -1054,24 +1149,29 @@ def rescale_absorbance(
         if __debug__:
             printdbg("... rescale: absorbance A_2 = A_1*(x_2/x_1)*(L_2/L_1)")
         _, absorbance = spec.get(
-            "absorbance", wunit=waveunit, Iunit=units["absorbance"]
+            "absorbance", wunit=waveunit, Iunit=spec.units["absorbance"], copy=True
         )
+        assert spec.units["absorbance"] == ""
         absorbance *= new_mole_fraction / old_mole_fraction  # rescale x
-        absorbance *= new_path_length / old_path_length  # rescale L
+        absorbance *= new_path_length_cm / old_path_length_cm  # rescale L
         unit = ""
     elif "abscoeff" in rescaled and true_path_length:  # in cm
         if __debug__:
             printdbg("... rescale: absorbance A_2 = k_2*L_2")
         abscoeff = rescaled["abscoeff"]  # x already scaled
-        absorbance = abscoeff * new_path_length  # calculate L
-        if "abscoeff" in spec.units:
-            assert spec.units["abscoeff"] == "cm-1"
+        absorbance = abscoeff * new_path_length_cm  # calculate L
+        assert rescaled_units["abscoeff"] == "cm-1"
         unit = ""
     elif "transmittance_noslit" in initial and true_path_length:
         if __debug__:
             printdbg("... rescale: absorbance A_2 = -ln(T1)*(x_2/x_1)*(L_2/L_1)")
         # Get absorbance from transmittance
-        _, T1 = spec.get("transmittance_noslit", wunit=waveunit)
+        _, T1 = spec.get(
+            "transmittance_noslit",
+            wunit=waveunit,
+            Iunit=spec.units["transmittance_noslit"],
+            copy=False,
+        )
 
         # We'll have a problem if the spectrum is optically thick
         b = T1 == 0  # no transmittance: optically thick mask
@@ -1087,7 +1187,7 @@ def rescale_absorbance(
         # Else, let's calculate it
         absorbance = -ln(T1)
         absorbance *= new_mole_fraction / old_mole_fraction  # rescale x
-        absorbance *= new_path_length / old_path_length  # rescale L
+        absorbance *= new_path_length_cm / old_path_length_cm  # rescale L
         unit = ""
     else:
         msg = (
@@ -1109,10 +1209,10 @@ def rescale_absorbance(
     # Export rescaled value
     if absorbance is not None:
         rescaled["absorbance"] = absorbance
-    if unit is not None:
-        units["absorbance"] = unit
+    if rescaled_units is not None:
+        rescaled_units["absorbance"] = unit
 
-    return rescaled, units
+    return rescaled, rescaled_units
 
 
 # ... transmittance
@@ -1122,8 +1222,8 @@ def rescale_transmittance_noslit(
     initial,
     old_mole_fraction,
     new_mole_fraction,
-    old_path_length,
-    new_path_length,
+    old_path_length_cm,
+    new_path_length_cm,
     waveunit,
     units,
     extra,
@@ -1133,13 +1233,11 @@ def rescale_transmittance_noslit(
 
     Parameters
     ----------
-
     spec: Spectrum
 
 
     References
     ----------
-
     Compute from the absorbance :
 
     .. math::
@@ -1173,7 +1271,7 @@ def rescale_transmittance_noslit(
 
     unit = None
 
-    def get_unit():
+    def get_transmittance_unit():
         return ""
 
     # case where we recomputed it already (somehow... ex: no_change signaled)
@@ -1192,14 +1290,15 @@ def get_unit():
             printdbg("... rescale: transmittance_noslit T_2 = exp(-A_2)")
         absorbance = rescaled["absorbance"]  # x and L already scaled
         transmittance_noslit = exp(-absorbance)  # recalculate
-        unit = get_unit()
+        unit = get_transmittance_unit()
     elif "abscoeff" in rescaled and true_path_length:
         if __debug__:
             printdbg("... rescale: transmittance_noslit T_2 = exp(-k_2*L_2)")
         abscoeff = rescaled["abscoeff"]  # x already scaled
-        absorbance = abscoeff * new_path_length  # calculate
+        assert units["abscoeff"] == "cm-1"
+        absorbance = abscoeff * new_path_length_cm  # calculate
         transmittance_noslit = exp(-absorbance)  # recalculate
-        unit = get_unit()
+        unit = get_transmittance_unit()
     elif "transmittance_noslit" in initial:
         if __debug__:
             printdbg(
@@ -1208,13 +1307,17 @@ def get_unit():
             )
         # get transmittance from initial transmittance
         _, T1 = spec.get(
-            "transmittance_noslit", wunit=waveunit, Iunit=units["transmittance_noslit"]
+            "transmittance_noslit",
+            wunit=waveunit,
+            Iunit=units["transmittance_noslit"],
+            copy=False,
         )
 
         # We'll have a problem if the spectrum is optically thick
         b = T1 == 0  # optically thick mask
         if b.sum() > 0 and (
-            new_mole_fraction < old_mole_fraction or new_path_length < old_path_length
+            new_mole_fraction < old_mole_fraction
+            or new_path_length_cm < old_path_length_cm
         ):
             # decreasing mole fractions/ path length could increase the transmittance
             # but this information was lost in the saturation
@@ -1228,8 +1331,9 @@ def get_unit():
         # Else, just get absorbance
         absorbance = -ln(T1)
         absorbance *= new_mole_fraction / old_mole_fraction  # rescale x
-        absorbance *= new_path_length / old_path_length  # rescale L
+        absorbance *= new_path_length_cm / old_path_length_cm  # rescale L
         transmittance_noslit = exp(-absorbance)
+        unit = get_transmittance_unit()
     else:
         msg = "Missing data to rescale transmittance. Expected scaled absorbance ({0})".format(
             "absorbance" in rescaled
@@ -1258,8 +1362,8 @@ def rescale_transmittance(
     initial,
     old_mole_fraction,
     new_mole_fraction,
-    old_path_length,
-    new_path_length,
+    old_path_length_cm,
+    new_path_length_cm,
     waveunit,
     units,
     extra,
@@ -1268,7 +1372,6 @@ def rescale_transmittance(
 
     Parameters
     ----------
-
     spec: Spectrum
     """
 
@@ -1301,11 +1404,11 @@ def rescale_radiance_noslit(
     initial,
     old_mole_fraction,
     new_mole_fraction,
-    old_path_length,
-    new_path_length,
+    old_path_length_cm,
+    new_path_length_cm,
     optically_thin,
     waveunit,
-    units,
+    rescaled_units,
     extra,
     true_path_length,
 ):
@@ -1313,13 +1416,11 @@ def rescale_radiance_noslit(
 
     Parameters
     ----------
-
     spec: Spectrum
 
 
     References
     ----------
-
     If optically thin, calculate from the emission coefficient :
 
     .. math::
@@ -1359,11 +1460,16 @@ def rescale_radiance_noslit(
     :py:attr:`~radis.spectrum.rescale.rescale_radiance_noslit`,
 
     """
+    if __debug__:
+        printdbg(
+            f"recomputing `radiance_noslit` from initial {initial} and already rescaled {list(rescaled.keys())}, knowing `old_mole_fraction={old_mole_fraction}`, `new_mole_fraction={new_mole_fraction}`, `old_path_length_cm={old_path_length_cm}`, `new_path_length_cm={new_path_length_cm}`, `true_path_length={true_path_length}`"
+        )
 
     unit = None
 
     def get_radiance_unit(unit_emisscoeff):
-        """get radiance_noslit unit from emisscoeff unit."""
+        """get radiance_noslit unit from emisscoeff unit.
+        Basically just multiply by ``cm``"""
         if "/cm3" in unit_emisscoeff:
             return unit_emisscoeff.replace("/cm3", "/cm2")
         else:
@@ -1373,8 +1479,8 @@ def get_radiance_unit(unit_emisscoeff):
     if "radiance_noslit" in rescaled:
         if __debug__:
             printdbg("... rescale: radiance_noslit was scaled already")
-        assert "radiance_noslit" in units
-        return rescaled, units
+        assert "radiance_noslit" in rescaled_units
+        return rescaled, rescaled_units
 
     # Rescale!
     if "emisscoeff" in rescaled and true_path_length and optically_thin:
@@ -1383,8 +1489,8 @@ def get_radiance_unit(unit_emisscoeff):
                 "... rescale: radiance_noslit I_2 = j_2 * L_2 " + "(optically thin)"
             )
         emisscoeff = rescaled["emisscoeff"]  # x already scaled
-        radiance_noslit = emisscoeff * new_path_length  # recalculate L
-        unit = get_radiance_unit(units["emisscoeff"])
+        radiance_noslit = emisscoeff * new_path_length_cm  # recalculate L
+        unit = get_radiance_unit(rescaled_units["emisscoeff"])
 
     elif (
         "emisscoeff" in rescaled
@@ -1397,6 +1503,7 @@ def get_radiance_unit(unit_emisscoeff):
             printdbg("... rescale: radiance_noslit I_2 = j_2*(1-T_2)/k_2")
         emisscoeff = rescaled["emisscoeff"]  # x already scaled
         abscoeff = rescaled["abscoeff"]  # x already scaled
+        assert rescaled_units["abscoeff"] == "cm-1"
         # mole_fraction, path_length already scaled
         transmittance_noslit = rescaled["transmittance_noslit"]
         b = transmittance_noslit == 1  # optically thin mask
@@ -1404,8 +1511,8 @@ def get_radiance_unit(unit_emisscoeff):
         radiance_noslit[~b] = (
             emisscoeff[~b] / abscoeff[~b] * (1 - transmittance_noslit[~b])
         )
-        radiance_noslit[b] = emisscoeff[b] * new_path_length  # optically thin limit
-        unit = get_radiance_unit(units["emisscoeff"])
+        radiance_noslit[b] = emisscoeff[b] * new_path_length_cm  # optically thin limit
+        unit = get_radiance_unit(unit_emisscoeff=rescaled_units["emisscoeff"])
 
     elif (
         "emisscoeff" in rescaled
@@ -1417,13 +1524,16 @@ def get_radiance_unit(unit_emisscoeff):
             printdbg("... rescale: radiance_noslit I_2 = j_2*(1-exp(-k_2*L_2))/k_2")
         emisscoeff = rescaled["emisscoeff"]  # x already scaled
         abscoeff = rescaled["abscoeff"]  # x already scaled
+        assert rescaled_units["abscoeff"] == "cm-1"
         b = abscoeff == 0  # optically thin mask
         radiance_noslit = np.empty_like(emisscoeff)  # calculate
         radiance_noslit[~b] = (
-            emisscoeff[~b] / abscoeff[~b] * (1 - exp(-abscoeff[~b] * new_path_length))
+            emisscoeff[~b]
+            / abscoeff[~b]
+            * (1 - exp(-abscoeff[~b] * new_path_length_cm))
         )
-        radiance_noslit[b] = emisscoeff[b] * new_path_length  # optically thin limit
-        unit = get_radiance_unit(units["emisscoeff"])
+        radiance_noslit[b] = emisscoeff[b] * new_path_length_cm  # optically thin limit
+        unit = get_radiance_unit(unit_emisscoeff=rescaled_units["emisscoeff"])
 
     elif "radiance_noslit" in initial and optically_thin:
         if __debug__:
@@ -1432,10 +1542,14 @@ def get_radiance_unit(unit_emisscoeff):
                 + "(optically thin)"
             )
         _, radiance_noslit = spec.get(
-            "radiance_noslit", wunit=waveunit, Iunit=units["radiance_noslit"]
+            "radiance_noslit",
+            wunit=waveunit,
+            Iunit=rescaled_units["radiance_noslit"],
+            copy=True,
         )
         radiance_noslit *= new_mole_fraction / old_mole_fraction  # rescale
-        radiance_noslit *= new_path_length / old_path_length  # rescale
+        radiance_noslit *= new_path_length_cm / old_path_length_cm  # rescale
+        unit = rescaled_units["radiance_noslit"]
 
     else:
         if optically_thin:
@@ -1455,7 +1569,7 @@ def get_radiance_unit(unit_emisscoeff):
                 raise ValueError(msg)
         else:
             msg = (
-                "Missing data to recalculate radiance_noslit for a non-optically thin column with thermal_equilibrium={0}. You need at least "
+                "Missing data to recalculate radiance_noslit for a non-optically thin column. You need at least "
                 + "scaled emisscoeff ({0}), scaled transmittance_noslit ({1}), ".format(
                     "emisscoeff" in rescaled, "transmittance_noslit" in rescaled
                 )
@@ -1474,9 +1588,9 @@ def get_radiance_unit(unit_emisscoeff):
     if radiance_noslit is not None:
         rescaled["radiance_noslit"] = radiance_noslit
     if unit is not None:
-        units["radiance_noslit"] = unit
+        rescaled_units["radiance_noslit"] = unit
 
-    return rescaled, units
+    return rescaled, rescaled_units
 
 
 # ... radiance
@@ -1488,11 +1602,11 @@ def rescale_radiance(
     initial,
     old_mole_fraction,
     new_mole_fraction,
-    old_path_length,
-    new_path_length,
+    old_path_length_cm,
+    new_path_length_cm,
     optically_thin,
     waveunit,
-    units,
+    rescaled_units,
     extra,
     true_path_length,
 ):
@@ -1516,8 +1630,8 @@ def rescale_radiance(
     if "radiance" in rescaled:
         if __debug__:
             printdbg("... rescale: radiance was scaled already")
-        assert "radiance" in units
-        return rescaled, units, apply_slit
+        assert "radiance" in rescaled_units
+        return rescaled, rescaled_units, apply_slit
 
     # Rescale!
     if "radiance_noslit" in rescaled:
@@ -1525,13 +1639,13 @@ def rescale_radiance(
     else:
         raise NotImplementedError("rescale radiance not implemented yet")
 
-    return rescaled, units, apply_slit
+    return rescaled, rescaled_units, apply_slit
 
 
 # ... emissivity_noslit
 
 
-def rescale_emissivity_noslit(spec, rescaled, units, extra, true_path_length):
+def rescale_emissivity_noslit(spec, rescaled, rescaled_units, extra, true_path_length):
     r"""
 
     Parameters
@@ -1554,13 +1668,15 @@ def rescale_emissivity_noslit(spec, rescaled, units, extra, true_path_length):
     """
 
     # case where we recomputed it already (somehow... ex: no_change signaled)
+    # -------------------------
     if "emissivity_noslit" in rescaled:
         if __debug__:
             printdbg("... rescale: emissivity_noslit was scaled already")
-        assert "emissivity_noslit" in units
-        return rescaled, units
+        assert "emissivity_noslit" in rescaled_units
+        return rescaled, rescaled_units
 
-    # Rescale!
+    # Or: rescale!
+    # -------------
     if "transmittance_noslit" in rescaled:
         if __debug__:
             printdbg("... rescale: emissivity_noslit e_2 = 1 - T_2")
@@ -1579,9 +1695,9 @@ def rescale_emissivity_noslit(spec, rescaled, units, extra, true_path_length):
     # Export rescaled value
     if emissivity_noslit is not None:
         rescaled["emissivity_noslit"] = emissivity_noslit
-        units["emissivity_noslit"] = ""
+        rescaled_units["emissivity_noslit"] = ""
 
-    return rescaled, units
+    return rescaled, rescaled_units
 
 
 # ... cross sections
@@ -1593,10 +1709,10 @@ def rescale_xsection(
     initial,
     old_mole_fraction,
     new_mole_fraction,
-    old_path_length,
-    new_path_length,
+    old_path_length_cm,
+    new_path_length_cm,
     waveunit,
-    units,
+    rescaled_units,
     extra,
     true_path_length,
 ):
@@ -1638,8 +1754,8 @@ def rescale_xsection(
     if "xsection" in rescaled:
         if __debug__:
             printdbg("... rescale: xsection was scaled already")
-        assert "xsection" in units
-        return rescaled, units
+        assert "xsection" in rescaled_units
+        return rescaled, rescaled_units
 
     # Get scaled xsection directly from XS1
     # -------------------------------------
@@ -1648,8 +1764,10 @@ def rescale_xsection(
         # cross-sections are unchanged
         if __debug__:
             printdbg("... rescale: xsection XS_2 = XS_1")
-        _, xsection = spec.get("xsection", wunit=waveunit, Iunit=units["xsection"])
-        unit = units["xsection"]  # units unchanged
+        _, xsection = spec.get(
+            "xsection", wunit=waveunit, Iunit=rescaled_units["xsection"]
+        )
+        unit = rescaled_units["xsection"]  # units unchanged
     elif (
         "abscoeff" in rescaled
         and "Tgas" in spec.conditions
@@ -1659,14 +1777,14 @@ def rescale_xsection(
         if __debug__:  # cm-1
             printdbg("... rescale: xsection XS_2 = k_2 * (k_b * T / x / p)")
         abscoeff = rescaled["abscoeff"]  # x already scaled
-        pressure_Pa = spec.conditions["pressure_mbar"] * 1e2
-        x = spec.conditions["mole_fraction"]
-        Tgas = spec.conditions["Tgas"]  # K
-        from radis.phys.constants import k_b
+        assert spec.units["abscoeff"] == "cm-1"
 
-        xsection = abscoeff * (k_b * Tgas / x / pressure_Pa) * 1e6  # cm2
-        if "abscoeff" in spec.units:
-            assert spec.units["abscoeff"] == "cm-1"
+        xsection = abscoeff2xsection(
+            abscoeff_cm1=abscoeff,
+            Tgas_K=spec.conditions["Tgas"],
+            mole_fraction=spec.conditions["mole_fraction"],
+            pressure_Pa=spec.conditions["pressure_mbar"] * 1e2,
+        )
         unit = "cm2"
     else:
         msg = (
@@ -1691,16 +1809,16 @@ def rescale_xsection(
     if xsection is not None:
         rescaled["xsection"] = xsection
     if unit is not None:
-        units["xsection"] = unit
+        rescaled_units["xsection"] = unit
 
-    return rescaled, units
+    return rescaled, rescaled_units
 
 
 def _recalculate(
     spec,
     quantity,
-    new_path_length,
-    old_path_length,
+    new_path_length_cm,
+    old_path_length_cm,
     new_mole_fraction,
     old_mole_fraction,
     true_path_length=True,
@@ -1805,7 +1923,8 @@ def _recalculate(
     # There are two cases: either we are actually rescaling to another length /
     # mole fraction, or we are just updating() without changing length / mole fraction
     no_change = (
-        new_mole_fraction == old_mole_fraction and new_path_length == old_path_length
+        new_mole_fraction == old_mole_fraction
+        and new_path_length_cm == old_path_length_cm
     )
 
     # Quickly stop if no change
@@ -1843,7 +1962,7 @@ def _recalculate(
     recompute = set(recompute)  # remove duplicates
 
     # Get units
-    units = spec.units.copy()
+    rescaled_units = spec.units.copy()
 
     # Recompute!
     # ----------
@@ -1853,22 +1972,25 @@ def _recalculate(
     # If no_change, just set everyone as rescaled already
     if no_change:
         for k in initial:
-            rescaled[k] = spec.get(k)[1]  # note: creates a copy
+            rescaled[k] = spec.get(k, wunit=spec.get_waveunit(), Iunit=spec.units[k])[
+                1
+            ]  # note: creates a copy
 
     # Start with abscoeff
 
     if "abscoeff" in recompute:
-        rescaled, units = rescale_abscoeff(
+        rescaled, rescaled_units = rescale_abscoeff(
             spec,
             rescaled,
             initial,  # Todo: remove rescaled = ... Dict is mutable no?
             old_mole_fraction,
             new_mole_fraction,
-            old_path_length,
+            old_path_length_cm,
             waveunit,
-            units,
+            rescaled_units,
             extra,
             true_path_length,
+            assume_equilibrium,
         )
 
     if (
@@ -1887,122 +2009,128 @@ def _recalculate(
             )
         wavenumber = spec.get_wavenumber()
         Tgas = spec.conditions["Tgas"]
-        rescaled, units = _recompute_all_at_equilibrium(
-            spec, rescaled, wavenumber, Tgas, new_path_length, true_path_length, units
+        rescaled, rescaled_units = _recompute_all_at_equilibrium(
+            spec,
+            rescaled,
+            wavenumber,
+            Tgas,
+            new_path_length_cm,
+            true_path_length,
+            rescaled_units,
         )
         apply_slit = "radiance" in recompute or "transmittance" in recompute
 
     else:
 
         if "emisscoeff" in recompute:
-            rescaled, units = rescale_emisscoeff(
+            rescaled, rescaled_units = rescale_emisscoeff(
                 spec,
                 rescaled,
                 initial,
                 old_mole_fraction,
                 new_mole_fraction,
-                old_path_length,
+                old_path_length_cm,
                 optically_thin,
                 waveunit,
-                units,
+                rescaled_units,
                 extra,
                 true_path_length,
             )
 
         if "absorbance" in recompute:
-            rescaled, units = rescale_absorbance(
+            rescaled, rescaled_units = rescale_absorbance(
                 spec,
                 rescaled,
                 initial,
                 old_mole_fraction,
                 new_mole_fraction,
-                old_path_length,
-                new_path_length,
+                old_path_length_cm,
+                new_path_length_cm,
                 waveunit,
-                units,
+                rescaled_units,
                 extra,
                 true_path_length,
             )
 
         if "transmittance_noslit" in recompute:
-            rescaled, units = rescale_transmittance_noslit(
+            rescaled, rescaled_units = rescale_transmittance_noslit(
                 spec,
                 rescaled,
                 initial,
                 old_mole_fraction,
                 new_mole_fraction,
-                old_path_length,
-                new_path_length,
+                old_path_length_cm,
+                new_path_length_cm,
                 waveunit,
-                units,
+                rescaled_units,
                 extra,
                 true_path_length,
             )
 
         if "radiance_noslit" in recompute:
-            rescaled, units = rescale_radiance_noslit(
+            rescaled, rescaled_units = rescale_radiance_noslit(
                 spec,
                 rescaled,
                 initial,
                 old_mole_fraction,
                 new_mole_fraction,
-                old_path_length,
-                new_path_length,
+                old_path_length_cm,
+                new_path_length_cm,
                 optically_thin,
                 waveunit,
-                units,
+                rescaled_units,
                 extra,
                 true_path_length,
             )
 
         if "emissivity_noslit" in recompute:
-            rescaled, units = rescale_emissivity_noslit(
-                spec, rescaled, units, extra, true_path_length
+            rescaled, rescaled_units = rescale_emissivity_noslit(
+                spec, rescaled, rescaled_units, extra, true_path_length
             )
 
         if "radiance" in recompute:
-            rescaled, units, slit_needed = rescale_radiance(
+            rescaled, rescaled_units, slit_needed = rescale_radiance(
                 spec,
                 rescaled,
                 initial,
                 old_mole_fraction,
                 new_mole_fraction,
-                old_path_length,
-                new_path_length,
+                old_path_length_cm,
+                new_path_length_cm,
                 optically_thin,
                 waveunit,
-                units,
+                rescaled_units,
                 extra,
                 true_path_length,
             )
             apply_slit = apply_slit or slit_needed
 
         if "transmittance" in recompute:
-            rescaled, units, slit_needed = rescale_transmittance(
+            rescaled, rescaled_units, slit_needed = rescale_transmittance(
                 spec,
                 rescaled,
                 initial,
                 old_mole_fraction,
                 new_mole_fraction,
-                old_path_length,
-                new_path_length,
+                old_path_length_cm,
+                new_path_length_cm,
                 waveunit,
-                units,
+                rescaled_units,
                 extra,
             )
             apply_slit = apply_slit or slit_needed
 
     if "xsection" in recompute:
-        rescaled, units = rescale_xsection(
+        rescaled, rescaled_units = rescale_xsection(
             spec,
             rescaled,
             initial,
             old_mole_fraction,
             new_mole_fraction,
-            old_path_length,
+            old_path_length_cm,
             optically_thin,
             waveunit,
-            units,
+            rescaled_units,
             extra,
             true_path_length,
         )
@@ -2021,8 +2149,7 @@ def _recalculate(
             spec._q[q] = rescaled[q]
 
     # Update units
-    for k, u in units.items():
-        spec.units[k] = u
+    spec.units.update(rescaled_units)
 
     # Reapply slit if needed
     # TODO: replace with directly convolving with slit stored in conditions
@@ -2095,6 +2222,9 @@ def _recalculate(
                 + "rescaled. This can lead to error. Rescaled spectrum "
                 + "contains: {0}".format(rescaled_list)
             )
+    # ... all values have units
+    for k in spec.get_vars():
+        assert k in spec.units
 
 
 def rescale_path_length(
@@ -2106,30 +2236,24 @@ def rescale_path_length(
 
     Parameters
     ----------
-
     spec: Spectrum
-
-    new_path_length: float
+    new_path_length: float (cm) or `~astropy.units.quantity.Quantity`
         new path length
-
-    old_path_length: float, or None
+    old_path_length: float (cm) or `~astropy.units.quantity.Quantity`, or None
         if None, current path length (conditions['path_length']) is used
 
 
     Other Parameters
     ----------------
-
     inplace: boolean
         if ``True``, modifies the Spectrum object directly. Else, returns
         a copy. Default ``False``.
-
     force: boolean
         if False, won't allow rescaling to 0 (not to loose information).
         Default ``False``
 
     Returns
     -------
-
     s_rescaled: Spectrum
         a rescaled Spectrum.
         if ``inplace=True``, then ``s`` has been rescaled already and
@@ -2137,7 +2261,6 @@ def rescale_path_length(
 
     Notes
     -----
-
     Implementation:
 
         To deal with all the input cases, we first make a list of what has to
@@ -2178,8 +2301,8 @@ def rescale_path_length(
             "Rescaling to 0 will loose information. Choose force " "= True"
         )
     # Convert units
-    new_path_length = convert_and_strip_units(new_path_length, unit.cm)
-    old_path_length = convert_and_strip_units(old_path_length, unit.cm)
+    new_path_length = convert_and_strip_units(new_path_length, u.cm)
+    old_path_length = convert_and_strip_units(old_path_length, u.cm)
 
     for q in ["transmittance", "radiance"]:
         qns = q + "_noslit"
@@ -2195,11 +2318,11 @@ def rescale_path_length(
     assume_equilibrium = spec.conditions.get("thermal_equilibrium", False)
     _recalculate(
         spec,
-        "same",
-        new_path_length,
-        old_path_length,
-        1,
-        1,
+        quantity="same",
+        new_path_length_cm=new_path_length,
+        old_path_length_cm=old_path_length,
+        new_mole_fraction=1,
+        old_mole_fraction=1,
         assume_equilibrium=assume_equilibrium,
     )
 
@@ -2223,30 +2346,24 @@ def rescale_mole_fraction(
 
     Parameters
     ----------
-
     spec: Spectrum
-
     new_mole_fraction: float
         new mole fraction
-
     old_mole_fraction: float, or None
         if None, current mole fraction (conditions['mole_fraction']) is used
 
 
     Other Parameters
     ----------------
-
     inplace: boolean
         if ``True``, modifies the Spectrum object directly. Else, returns
         a copy. Default ``False``.
-
     force: boolean
         if False, won't allow rescaling to 0 (not to loose information).
         Default ``False``
 
     Returns
     -------
-
     s_rescaled: Spectrum
         a rescaled Spectrum.
         if ``inplace=True``, then ``s`` has been rescaled already and
@@ -2254,7 +2371,6 @@ def rescale_mole_fraction(
 
     Notes
     -----
-
     Implementation:
 
         similar to rescale_path_length() but we have to scale abscoeff & emisscoeff
@@ -2329,18 +2445,20 @@ def rescale_mole_fraction(
         path_length = spec.conditions["path_length"]
         true_path_length = True
     else:
-        path_length = 1
+        path_length = (
+            1  # path length not given, but we can still recompute a few spectral arrays
+        )
         true_path_length = False
 
     # Rescale
     assume_equilibrium = spec.conditions.get("thermal_equilibrium", False)
     _recalculate(
         spec,
-        "same",
-        path_length,
-        path_length,
-        new_mole_fraction,
-        old_mole_fraction,
+        quantity="same",
+        new_path_length_cm=path_length,
+        old_path_length_cm=path_length,
+        new_mole_fraction=new_mole_fraction,
+        old_mole_fraction=old_mole_fraction,
         true_path_length=true_path_length,
         assume_equilibrium=assume_equilibrium,
         verbose=verbose,
diff --git a/radis/spectrum/spectrum.py b/radis/spectrum/spectrum.py
index d7e36f102..2c688190a 100644
--- a/radis/spectrum/spectrum.py
+++ b/radis/spectrum/spectrum.py
@@ -69,7 +69,7 @@
 )
 from radis.misc.debug import printdbg
 from radis.misc.plot import split_and_plot_by_parts
-from radis.misc.signal import resample
+from radis.misc.signal import resample, resample_even
 from radis.phys.air import air2vacuum, vacuum2air
 from radis.phys.convert import cm2nm, conv2, nm2cm
 from radis.phys.units import Unit, convert_universal
@@ -919,7 +919,15 @@ def from_specutils(self, spectrum, var="radiance"):
     # ----------------
     # XXX =====================================================================
 
-    def get(self, var, wunit="default", Iunit="default", copy=True, trim_nan=False):
+    def get(
+        self,
+        var,
+        wunit="default",
+        Iunit="default",
+        copy=True,
+        trim_nan=False,
+        return_units=False,
+    ):
         """Retrieve a spectral quantity from a Spectrum object. You can select
         wavespace unit, intensity unit, or propagation medium.
 
@@ -928,7 +936,7 @@ def get(self, var, wunit="default", Iunit="default", copy=True, trim_nan=False):
         var: variable (``'absorbance'``, ``'transmittance'``, ``'xsection'`` etc.)
             Should be a defined quantity among :data:`~radis.spectrum.utils.CONVOLUTED_QUANTITIES`
             or :data:`~radis.spectrum.utils.NON_CONVOLUTED_QUANTITIES`.
-            To get the full list of quantities defined in this Spectrum object use
+            To get the full list of spectral arrays defined in this Spectrum object use
             the :meth:`~radis.spectrum.spectrum.Spectrum.get_vars` method.
         wunit: ``'nm'``, ``'cm'``, ``'nm_vac'``.
             wavespace unit: wavelength in air (``'nm'``), wavenumber
@@ -941,6 +949,14 @@ def get(self, var, wunit="default", Iunit="default", copy=True, trim_nan=False):
             For ``var="radiance"``, one can use per wavelength (~ 'W/m2/sr/nm')
             or per wavenumber (~ 'W/m2/sr/cm-1') units
 
+            .. note::
+
+                if using ``default``, the returned unit may be different from
+                Spectrum.units[var].  e.g, for a Spectrum with radiance stored
+                in 'mW/cm2/sr/nm', getting it with `wunit='cm-1', Iunit='default'`
+                will return Iunit in 'mW/cm2/sr/cm-1' for consistency.
+                Use ``return_units`` to be sure about which units are used.
+
         Other Parameters
         ----------------
         copy: bool
@@ -949,6 +965,11 @@ def get(self, var, wunit="default", Iunit="default", copy=True, trim_nan=False):
         trim_nan: bool
             if ``True``, removes ``nan`` on the sides of the spectral array
             (and corresponding wavespace). Default ``False``.
+        return_units: bool
+            if ``True``, return dimensioned Astropy arrays.
+            Is using ``return_units = 'as_str'``, return `wunit` and `Iunit`
+            as extra variables, i.e, ``w, I, wunit, Iunit = s.get(..., return_units='as_str')``
+            Default ``False``
 
         Returns
         -------
@@ -968,14 +989,16 @@ def get(self, var, wunit="default", Iunit="default", copy=True, trim_nan=False):
 
             _, R = s.get('radiance_noslit', wunit='nm', Iunit='W/cm2/sr/nm')
 
+        Use with `return_units` to get dimensioned Astropy Quantities ::
+
+            w, R  = s.get('radiance_noslit', return_units=True)
+
         See Also
         --------
         :meth:`~radis.spectrum.spectrum.Spectrum.get_radiance`,
         :meth:`~radis.spectrum.spectrum.Spectrum.get_radiance_noslit`,
         :ref:`the Spectrum page <label_spectrum>`
         """
-        # TODO: allow to get radiance in (W/sr/cm2/nm) or (W/sr/cm2) multiplying
-        # by FWHM.
 
         # check input
         if not var in self.get_vars():
@@ -1002,7 +1025,10 @@ def get(self, var, wunit="default", Iunit="default", copy=True, trim_nan=False):
 
         # Get wavespace (in correct unit, and correct medium)
         if wunit == "default":
-            wunit = self.get_waveunit()
+            if "default_output_unit" in self.conditions:
+                wunit = self.c["default_output_unit"]
+            else:
+                wunit = self.get_waveunit()
         wunit = cast_waveunit(wunit)
         if wunit == "cm-1":
             w = self.get_wavenumber(copy=copy)
@@ -1025,28 +1051,55 @@ def get(self, var, wunit="default", Iunit="default", copy=True, trim_nan=False):
 
         # Convert y unit if necessary
         Iunit0 = self.units[var]
+        if Iunit == "default":
+            # Add unit consistent with the namespace
+            # https://github.com/radis/radis/issues/456
+            if wunit in ["nm", "nm_air"]:
+                if var in [
+                    "radiance",
+                    "radiance_noslit",
+                    "emisscoeff",
+                    "emisscoeff_continuum",
+                ]:
+                    Iunit = self.units[var].replace("/cm-1", "/nm")
+            elif wunit in ["cm-1"]:
+                if var in [
+                    "radiance",
+                    "radiance_noslit",
+                    "emisscoeff",
+                    "emisscoeff_continuum",
+                ]:
+                    Iunit = self.units[var].replace("/nm", "/cm-1")
+
+        # Retrieve data (with correct unit)
         if Iunit != "default" and Iunit != Iunit0:
             if var in ["radiance", "radiance_noslit"]:
                 # deal with the case where we want to get a radiance in per
                 # wavelength unit (~ W/sr/cm2/nm) in wavenumber units (~ W/sr/cm2/cm-1),
                 # or the other way round
+                w_cm = self.get_wavenumber(copy=False)
+                if trim_nan:
+                    w_cm = w_cm[m : M + 1]
                 I = convert_universal(
                     I,
                     Iunit0,
                     Iunit,
-                    self,
-                    per_nm_is_like="mW/sr/cm2/nm",
-                    per_cm_is_like="mW/sr/cm2/cm-1",
+                    w_cm,
+                    per_nm_is_like="mW/cm2/sr/nm",
+                    per_cm_is_like="mW/cm2/sr/cm-1",
                 )
             elif var in ["emisscoeff"]:
-                # idem for emisscoeff in (~ W/sr/cm3/nm) or (~ /sr/cm3/cm-1)
+                # idem for emisscoeff in (~ W/sr/cm3/nm) or (~ W/sr/cm3/cm-1)
+                w_cm = self.get_wavenumber(copy=False)
+                if trim_nan:
+                    w_cm = w_cm[m : M + 1]
                 I = convert_universal(
                     I,
                     Iunit0,
                     Iunit,
-                    self,
-                    per_nm_is_like="mW/sr/cm3/nm",
-                    per_cm_is_like="mW/sr/cm3/cm-1",
+                    w_cm,
+                    per_nm_is_like="mW/cm3/sr/nm",
+                    per_cm_is_like="mW/cm3/sr/cm-1",
                 )
             elif var in ["absorbance"]:  # no unit
                 assert Iunit in ["", "1"]
@@ -1059,7 +1112,12 @@ def get(self, var, wunit="default", Iunit="default", copy=True, trim_nan=False):
         else:
             pass
 
-        return w, I
+        if return_units == "as_str":  # only used internally
+            return w, I, wunit, Iunit
+        elif return_units:
+            return w * Unit(wunit), I * Unit(Iunit)
+        else:
+            return w, I
 
     def _get_wavespace(self, copy=True):
         """Return wavespace (if the same for all quantities)
@@ -1249,9 +1307,12 @@ def get_name(self):
         elif self.file is not None:
             name = "{0}".format(basename(self.file))
         else:
-            name = "{0}spectrum{1}".format(
-                self.conditions.get("molecule", ""), id(self)
-            )
+            name_params = []
+            for (key, unit) in [("molecule", ""), ("dbformat", ""), ("Tgas", "K")]:
+                if key in self.conditions:
+                    name_params.append(f"{self.conditions[key]}{unit}")
+
+            name = "-".join(name_params) + "-#" + str(id(self))[-4:]
 
         return name
 
@@ -1348,7 +1409,9 @@ def update(self, quantity="all", optically_thin="default", verbose=True):
                                     [2.71414065e-06, 2.88341489e-06, 3.06942277e-06, 3.27445689e-06,
                                      3.50121831e-06, 3.75290756e-06, 4.03334037e-06, 4.34709612e-06,
                                      4.69971017e-06, 5.09792551e-06],
-                                    'abscoeff', wunit='cm-1', Iunit='cm-1',
+                                    'abscoeff',
+                                    wunit='cm-1',
+                                    Iunit='cm-1',
                                     conditions={'path_length':1, # cm
                                                 'thermal_equilibrium':True,
                                                 'Tgas':700,  # K
@@ -1721,6 +1784,7 @@ def get_power(self, unit="mW/cm2/sr"):
         Examples
         --------
         ::
+
             s.get_power('W/cm2/sr')
 
         .. minigallery:: radis.spectrum.spectrum.Spectrum.offset
@@ -1938,20 +2002,13 @@ def plot(
                 if var.replace("_noslit", "") in params:  # favour convolved quantities
                     var = var.replace("_noslit", "")
 
-        if wunit == "default":
-            wunit = self.get_waveunit()
-        wunit = cast_waveunit(wunit)
         # Get variable
-        x, y = self.get(var, wunit=wunit, Iunit=Iunit)
+        x, y, wunit, Iunit = self.get(
+            var, wunit=wunit, Iunit=Iunit, return_units="as_str"
+        )
 
         # Get labels
         xlabel = format_xlabel(wunit, show_medium)
-        if Iunit == "default":
-            try:
-                Iunit0 = self.units[var]
-            except KeyError:  # unit not defined in dictionary
-                Iunit0 = "a.u"
-            Iunit = Iunit0
 
         # cosmetic changes
         ylabel = "{0} ({1})".format(make_up(var), make_up_unit(Iunit, var))
@@ -1987,20 +2044,20 @@ def clean_error_msg(string):
         if not force and (fig.gca().get_xlabel().lower() not in ["", xlabel.lower()]):
             raise ValueError(
                 "Error while plotting {0}. Cannot plot ".format(var)
-                + "on a same figure with different xlabel: {0}, {1}".format(
+                + "on a same figure with different xlabel: Current: {0}, New: {1}".format(
                     clean_error_msg(fig.gca().get_xlabel()), clean_error_msg(xlabel)
                 )
-                + "Use force=True if you really want to plot"
+                + ". If it's a unit problem, change unit with `s.plot(..., wunit=...)`. Use force=True if you really want to plot as is."
             )
         label1 = clean_error_msg(fig.gca().get_ylabel().lower())
         label2 = clean_error_msg(ylabel.lower())
         if not force and (label1 not in ["", label2]):
             raise ValueError(
                 "Error while plotting {0}. Cannot plot ".format(var)
-                + "on a same figure with different ylabel: \n{0}\n{1}".format(
+                + "on a same figure with different ylabel: \nCurrent: {0}\nNew: {1}".format(
                     clean_error_msg(fig.gca().get_ylabel()), clean_error_msg(ylabel)
                 )
-                + "\nIf it's a unit problem, you can change the unit with `Iunit=`. Use `force=True` if you really want to plot with these units."
+                + "\nIf it's a unit problem, you can change the unit with `s.plot(..., Iunit=)`. Use `force=True` if you really want to plot with these units."
             )
 
         # Add extra plotting parameters
@@ -2376,6 +2433,7 @@ def apply_slit(
         slit_dispersion=None,
         slit_dispersion_threshold=0.01,
         auto_recenter_crop=True,
+        assert_evenly_spaced=True,
         verbose=True,
         inplace=True,
         *args,
@@ -2436,6 +2494,11 @@ def apply_slit(
 
         Other Parameters
         ----------------
+        assert_evenly_spaced: boolean, or ``'resample'``
+            for the convolution to be accurate, ``w`` should be evenly spaced. If
+            ``assert_evenly_spaced=True``, then we check this is the case, and raise
+            an error if arrays is not evenly spaced. If ``'resample'``, then we resample
+            ``w`` and ``I`` if needed. Recommended, but it takes some time.
         auto_recenter_crop: bool
             if ``True``, recenter slit and crop zeros on the side when importing
             an experimental slit. Default ``True``.
@@ -2704,7 +2767,7 @@ def apply_slit(
                     mode=mode,
                     wunit=waveunit,
                     verbose=verbose,
-                    assert_evenly_spaced=False,
+                    assert_evenly_spaced=assert_evenly_spaced,
                     # assumes Spectrum is correct by construction
                     **kwargsconvolve,
                 )
@@ -3059,9 +3122,8 @@ def get_overlay(overlay):
             transmittance)"""
 
             if isinstance(overlay, str):  # either get it from the Spectrum
-                w, I = self.get(overlay, wunit=wunit)
+                w, I, _, units = self.get(overlay, wunit=wunit, return_units="as_str")
                 name = overlay
-                units = self.units[overlay]
                 return (w, I, name, units)
 
             else:  # Or use a given tuple or arrays
@@ -3418,7 +3480,6 @@ def resample(
 
                 s1.resample(s2.get_wavenumber())
                 s1.resample(s2)            # also valid
-
         unit: ``'same'``, ``'nm'``, ``'cm-1'``, ``'nm_vac'``
             unit of new wavespace. It ``'same'`` it is assumed to be the current
             waveunit. Default ``'same'``. The spectrum waveunit is changed to this
@@ -3456,11 +3517,15 @@ def resample(
         Examples
         --------
 
+        Convert a Spectrum in wavenumber to wavelengths::
+
+            s_nm = s.resample(s.get_wavelength(), "nm", inplace=False)
+
         .. minigallery:: radis.spectrum.spectrum.Spectrum.resample
 
         See Also
         --------
-        :func:`radis.misc.signal.resample`
+        :func:`radis.misc.signal.resample`, :py:meth:`radis.spectrum.spectrum.Spectrum.resample_even`
         """
         # Check inputs (check for deprecated)
         if if_conflict_drop is not None:
@@ -3553,6 +3618,63 @@ def get_filling(variable):
 
         return s
 
+    def resample_even(
+        self, energy_threshold=5e-3, print_conservation=False, inplace=True
+    ):
+        """Resample spectrum over the same waverange, but evenly spaced.
+
+        .. warning::
+            This may result in information loss. Resampling is done with
+            oversampling and spline interpolation. These parameters can be adjusted,
+            and energy conservation ensured with the appropriate parameters.
+
+        Parameters
+        ----------
+        energy_threshold: float or ``None``
+            if energy conservation (integrals on the intersecting range) is above
+            this threshold, raise an error. If ``None``, dont check for energy conservation
+            Default 5e-3 (0.5%)
+        print_conservation: boolean
+            if ``True``, prints energy conservation. Default ``False``.
+        inplace: boolean
+            if ``True``, modifies the Spectrum object directly. Else, returns
+            a copy. Default ``True``.
+        **kwargs: **dict
+            all other arguments are sent to :func:`radis.misc.signal.resample`
+
+        Returns
+        -------
+        Spectrum : resampled Spectrum object. If using ``inplace=True``, the Spectrum
+            object has been modified anyway.
+
+        Examples
+        --------
+
+        .. minigallery:: radis.spectrum.spectrum.Spectrum.resample
+
+        See Also
+        --------
+        :func:`radis.misc.signal.resample_even`, :py:meth:`radis.spectrum.spectrum.Spectrum.resample`
+        """
+        if inplace:
+            s = self
+        else:
+            s = self.copy()
+
+        w = s._q["wavespace"]  # wavelength or wavenumber range
+        for (k, I) in s._q.items():
+            if k == "wavespace":
+                continue
+            w_new, Inew = resample_even(
+                w, I, resfactor=2, print_conservation=print_conservation
+            )
+
+            s._q[k] = Inew
+        # update wavespace
+        s._q["wavespace"] = w_new
+
+        return s
+
     # %% ======================================================================
     # Semi public functions
     # ----------------
@@ -3720,10 +3842,27 @@ def __copy__(self, copy_lines=True, quantity="all"):
         #                                        # Spectrum() recreates a copy anyway
         if quantity == "all":
             quantities = dict(self._get_items())
+        elif quantity in self.get_vars():
+            quantities = {
+                quantity: self.get(
+                    quantity,
+                    wunit=self.get_waveunit(),
+                    Iunit=self.units[quantity],
+                    copy=True,
+                )
+            }
         else:
+            self.update(quantity)  # needed to get self.units[quantity]
             quantities = {
-                quantity: self.get(quantity, wunit=self.get_waveunit())
-            }  # dict(self._get_items())
+                quantity: self.get(
+                    quantity,
+                    wunit=self.get_waveunit(),
+                    Iunit=self.units[quantity],
+                    copy=True,
+                )  # copy=True still needed as the waverange is the one of the original array
+            }
+            del self._q[quantity]  # no need to keep it in this Spectrum
+            # del self.units[quantity]   # we can keep the units (+ they'll be copied below)
 
         try:
             units = deepcopy(self.units)
@@ -3765,6 +3904,7 @@ def __copy__(self, copy_lines=True, quantity="all"):
             wunit=waveunit,
             references=references,
             name=name,
+            check_wavespace=False,  # no need, as Spectrum was normally already properly defined
             warnings=False,  # saves about 3.5 ms on the Performance test object
         )
 
@@ -4095,12 +4235,12 @@ def _add_quantity(self, name, w, I, check_wavespace=True):
             else:
                 # Check Wavelength/wavenumber is evently spaced
                 if check_wavespace:
-                    if not evenly_distributed(w, tolerance=1e-5):
+                    if not evenly_distributed(w, atolerance=1e-5):
                         warn(
                             "Wavespace is not evenly spaced ({0:.3f}%) for {1}.".format(
                                 np.abs(np.diff(w)).max() / w.mean() * 100, name
                             )
-                            + " This may create problems when convolving with slit function"
+                            + " This may create problems if later convolving with slit function (`s.apply_slit()`). You can use `s.resample_even()`"
                         )
                 self._q["wavespace"] = np.array(w)  # copy
 
@@ -4159,6 +4299,9 @@ def __str__(self):
             print(
                 " " * 2,
                 k,
+                f"\t[{self.units[k]}]"
+                if (k in self.units and self.units[k] != "")
+                else "",
                 "\t({0:,d} points{1})".format(
                     len(v),
                     ", {0} nans".format(count_nans(v)) if anynan(v) else "",
@@ -4268,7 +4411,9 @@ def max(self, value_only=False):
         """
 
         var = self._get_unique_var(operation_name="max")
-        w, I = self.get(var, wunit=self.get_waveunit(), copy=False)
+        w, I = self.get(
+            var, wunit=self.get_waveunit(), Iunit=self.units[var], copy=False
+        )
         if value_only:
             return I[~np.isnan(I)].max()
         else:
@@ -4291,7 +4436,9 @@ def min(self, value_only=True):
         """
 
         var = self._get_unique_var(operation_name="min")
-        w, I = self.get(var, wunit=self.get_waveunit(), copy=False)
+        w, I = self.get(
+            var, wunit=self.get_waveunit(), Iunit=self.units[var], copy=False
+        )
         if value_only:
             return I[~np.isnan(I)].min()
         else:
@@ -4361,7 +4508,9 @@ def normalize(
 
         if wrange is not None and len(wrange) > 0:
             wmin, wmax = wrange
-            w, I = s.get(var, wunit=wunit, copy=False)  # (faster not to copy)
+            w, I = s.get(
+                var, wunit=wunit, Iunit=s.units[var], copy=False
+            )  # (faster not to copy)
             b = (w > wmin) & (w < wmax)
             if normalize_how == "max":
                 norm = np.nanmax(I[b])
@@ -4381,13 +4530,17 @@ def normalize(
 
         else:
             if normalize_how == "max":
-                norm = np.nanmax(s.get(var, copy=False)[1])
+                norm = np.nanmax(
+                    s.get(var, wunit=wunit, Iunit=s.units[var], copy=False)[1]
+                )
                 norm_unit = s.units[var]
             elif normalize_how == "mean":
-                norm = np.nanmean(s.get(var, copy=False)[1])
+                norm = np.nanmean(
+                    s.get(var, wunit=wunit, Iunit=s.units[var], copy=False)[1]
+                )
                 norm_unit = s.units[var]
             elif normalize_how == "area":
-                w, I = s.get(var, wunit=wunit, copy=False)
+                w, I = s.get(var, wunit=wunit, Iunit=s.units[var], copy=False)
                 norm = nantrapz(I, w)
                 norm_unit = Unit(s.units[var]) * Unit(wunit)
 
diff --git a/radis/spectrum/utils.py b/radis/spectrum/utils.py
index 3f3b10d47..4d8065fd6 100644
--- a/radis/spectrum/utils.py
+++ b/radis/spectrum/utils.py
@@ -16,8 +16,8 @@
 
 # Waverange, wavespaces
 # ... aliases:
-WAVENUM_UNITS = ["cm", "cm-1", "wavenumber"]
-WAVELEN_UNITS = ["nm", "wavelength"]
+WAVENUM_UNITS = ["cm", "cm-1"]
+WAVELEN_UNITS = ["nm", "nm_air"]
 WAVELENVAC_UNITS = ["nm_vac", "nm_vacuum"]
 # ... internal names
 WAVESPACE = ["nm", "nm_vac", "cm-1"]
diff --git a/radis/tools/slit.py b/radis/tools/slit.py
index b7bff4fea..f0347704a 100644
--- a/radis/tools/slit.py
+++ b/radis/tools/slit.py
@@ -50,7 +50,7 @@
 from scipy.interpolate import splev, splrep
 from scipy.signal import oaconvolve
 
-from radis.misc.arrays import anynan, evenly_distributed
+from radis.misc.arrays import anynan, evenly_distributed, evenly_distributed_fast
 from radis.misc.basics import is_float
 from radis.misc.debug import printdbg
 from radis.misc.signal import resample_even
@@ -546,18 +546,15 @@ def convolve_with_slit(
     ----------------
     k: int
         order of spline interpolation. 3: cubic, 1: linear. Default 1.
-
     bplot: boolean
         if ``True``, plots the convolve slit function (for debugging)
-
     verbose: boolean
         more blabla
-
-    assert_evenly_spaced: boolean
+    assert_evenly_spaced: boolean, or ``'resample'``
         for the convolution to be accurate, ``w`` should be evenly spaced. If
-        ``assert_evenly_spaced=True``, then we check this is the case, and resample
+        ``assert_evenly_spaced=True``, then we check this is the case, and raise
+        an error if arrays is not evenly spaced. If ``'resample'``, then we resample
         ``w`` and ``I`` if needed. Recommended, but it takes some time.
-
     wunit: ``'nm'``, ``'cm-1'``
         just for printing messages. However, ``w`` and ``w_slit`` should be in the
         same wavespace.
@@ -677,16 +674,24 @@ def convolve_with_slit(
     #    evenly spaced
     # --------------
 
-    # ... Resample if not evenly spaced
-    # TODO: add a criteria based based on FWHM rather than absolute?
-    wstep = abs(np.diff(w)).min()  # spectrum wavelength spacing
-    if assert_evenly_spaced:
-        if not evenly_distributed(w, tolerance=wstep * 1e-3):
-            # TODO: automatically find a resampling factor?
-            warn("Spectrum not evenly spaced. Resampling")
+    # check first & last spacing, always :
+    is_evenly_spaced = evenly_distributed_fast(w, rtolerance=1e-3)
+    wstep = w[1] - w[0]
+
+    # ... check with details
+    if is_evenly_spaced and assert_evenly_spaced:
+        wstep = abs(np.diff(w)).min()  # spectrum wavelength spacing
+        is_evenly_spaced = evenly_distributed(w, atolerance=wstep * 1e-3)
+
+    if not is_evenly_spaced:
+        if assert_evenly_spaced == "resample":
             w, I = resample_even(w, I, resfactor=2, print_conservation=True)
             wstep = abs(np.diff(w)).min()  # new spectrum wavelength spacing
 
+        raise ValueError(
+            "Spectrum is not evenly spaced. Cannot apply slit function. Please resample with `convolve_with_slit(..., convolve_with_slit='resample')` ; or if using a Spectrum object `s`, with `s.resample_even()`"
+        )
+
     # ... Check that the slit is not reversed (interpolation requires objects are sorted)
     reverse = w_slit[-1] < w_slit[0]
     if reverse:


Test patch:
diff --git a/radis/test/lbl/test_bands.py b/radis/test/lbl/test_bands.py
index bdc832f07..c426e8c13 100644
--- a/radis/test/lbl/test_bands.py
+++ b/radis/test/lbl/test_bands.py
@@ -72,23 +72,27 @@ def test_plot_all_CO2_bandheads(verbose=True, plot=False, *args, **kwargs):
 
     assert len(s_bands) == 3
 
-    # Ensure that recombining gives the same
     s_merged = MergeSlabs(*list(s_bands.values()))
-    assert s_tot.compare_with(s_merged, "radiance_noslit", plot=False)
-
-    if verbose:
-        printm("Recombining bands give the same Spectrum")
+    s_merged.name = "Merged bands"
 
     # %%
     if plot:
         s_tot.apply_slit(1, "nm")
         s_tot.name = "Full spectrum"
-        s_tot.plot(wunit="nm", lw=3)
+        s_tot.plot("radiance", wunit="nm", lw=3)
         for band, s in s_bands.items():
             s.plot(wunit="nm", nfig="same")
+        s_merged.apply_slit(1, "nm")
+        s_merged.plot("radiance", wunit="nm", nfig="same", lw=2)
         plt.legend(loc="upper left")
         plt.ylim(ymax=0.25)
 
+    # Ensure that recombining gives the same
+    assert s_tot.compare_with(s_merged, "radiance_noslit", plot=False)
+
+    if verbose:
+        printm("Recombining bands give the same Spectrum")
+
     # %% Compare with equilibrium bands now
     s_bands_eq = sf.eq_bands(Tgas)
     s_merged_eq = MergeSlabs(*list(s_bands_eq.values()))
diff --git a/radis/test/lbl/test_base.py b/radis/test/lbl/test_base.py
index d4403e34e..75bab4e51 100644
--- a/radis/test/lbl/test_base.py
+++ b/radis/test/lbl/test_base.py
@@ -13,7 +13,7 @@
 import radis
 from radis import get_residual, sPlanck
 from radis.lbl import SpectrumFactory
-from radis.lbl.base import get_waverange
+from radis.lbl.base import get_wavenumber_range
 from radis.misc.printer import printm
 from radis.misc.progress_bar import ProgressBar
 from radis.misc.utils import Default
@@ -367,7 +367,7 @@ def test_optically_thick_limit_1iso(verbose=True, plot=True, *args, **kwargs):
             s.rescale_path_length(1e6)
 
             if plot:
-                s.plot(wunit="nm", nfig="same", lw=4)
+                s.plot(wunit="nm", Iunit="mW/cm2/sr/nm", nfig="same", lw=4)
 
             if verbose:
                 printm(
@@ -482,7 +482,7 @@ def test_optically_thick_limit_2iso(verbose=True, plot=True, *args, **kwargs):
 
                 nfig = "test_opt_thick_limit_2iso {0}".format(s.name)
                 plt.figure(nfig).clear()
-                s.plot(wunit="nm", nfig=nfig, lw=4)
+                s.plot(wunit="nm", Iunit="mW/cm2/sr/nm", nfig=nfig, lw=4)
                 s_plck.plot(wunit="nm", nfig=nfig, Iunit="mW/cm2/sr/nm", lw=2)
                 plt.legend()
 
@@ -505,29 +505,34 @@ def test_optically_thick_limit_2iso(verbose=True, plot=True, *args, **kwargs):
         radis.config["DEBUG_MODE"] = DEBUG_MODE
 
 
-def test_get_waverange(*args, **kwargs):
+def test_get_wavenumber_range(*args, **kwargs):
+    """Test waverange conversions for various inputs"""
 
-    assert get_waverange(1 * u.um, 10 * u.um, medium="vacuum") == (1000, 10000)
+    assert get_wavenumber_range(
+        1 * u.um, 10 * u.um, medium="vacuum", return_input_wunit=True
+    ) == (1000, 10000, "nm_vac")
 
     # 'wunit' is none
     # ...'wmin/wmax' is none
     # ...... wavenumber is passed > should return the same wavenumbers
-    assert get_waverange(wavenum_min=10, wavenum_max=20, wunit=Default("cm-1")) == (
-        10,
-        20,
-    )
+    assert get_wavenumber_range(
+        wavenum_min=10, wavenum_max=20, wunit=Default("cm-1"), return_input_wunit=True
+    ) == (10, 20, "cm-1")
 
     # ...... wavelength is passed > should convert and return the wavenumbers
     assert np.isclose(
-        get_waverange(
-            wavelength_min=1, wavelength_max=2, medium="vacuum", wunit=Default("cm-1")
+        get_wavenumber_range(
+            wavelength_min=1,
+            wavelength_max=2,
+            medium="vacuum",
+            wunit=Default("cm-1"),
         ),
         (5000000.0, 10000000.0),
     ).all()
 
     # ....... passed both wavenumber and wavelength > should throw error
     with pytest.raises(ValueError):
-        get_waverange(
+        get_wavenumber_range(
             wavenum_min=1000,
             wavenum_max=2000,
             wavelength_min=1,
@@ -538,22 +543,22 @@ def test_get_waverange(*args, **kwargs):
 
     # ...... passed neither wavenumber nor wavlength > should throw error
     with pytest.raises(ValueError):
-        get_waverange(wunit=Default("cm-1"))
+        get_wavenumber_range(wunit=Default("cm-1"))
 
     # ... 'wmin/wmax' are passed as values without accompanying units
     # ...... wavenumber is passed > should throw error due to multiple input wavespace values
     with pytest.raises(ValueError):
-        get_waverange(
+        get_wavenumber_range(
             wavenum_min=1, wavenum_max=2, wmin=10, wmax=20, wunit=Default("cm-1")
         )
 
     # ...... wavelength is passed > should throw error due to multiple input wavespace values
     with pytest.raises(ValueError):
-        get_waverange(wmin=10, wmax=20, wavelength_min=1, wavelength_max=2)
+        get_wavenumber_range(wmin=10, wmax=20, wavelength_min=1, wavelength_max=2)
 
     # ...... passed both wavenumber and wavelength > should throw error due to multiple input wavespace values
     with pytest.raises(ValueError):
-        get_waverange(
+        get_wavenumber_range(
             wmin=1,
             wmax=2,
             wavenum_min=10,
@@ -564,12 +569,14 @@ def test_get_waverange(*args, **kwargs):
         )
 
     # ...... passed neither wavenumber nor wavelength > should return wavenumber after converting wmin/wmax assuming default units
-    assert get_waverange(wmin=10, wmax=20, wunit=Default("cm-1")) == (10.0, 20.0)
+    assert get_wavenumber_range(
+        wmin=10, wmax=20, wunit=Default("cm-1"), return_input_wunit=True
+    ) == (10.0, 20.0, "cm-1")
 
     # ... 'wmin/wmax' are passed as values with accompanying units
     # ...... wavenumber is passed > should throw error due to multiple input wavespace values
     with pytest.raises(ValueError):
-        get_waverange(
+        get_wavenumber_range(
             wavenum_min=10,
             wavenum_max=20,
             wmin=100 * (1 / u.cm),
@@ -579,7 +586,7 @@ def test_get_waverange(*args, **kwargs):
 
     # ...... wavelength is passed > should throw error due to multiple input wavespace values
     with pytest.raises(ValueError):
-        get_waverange(
+        get_wavenumber_range(
             wmin=100 * (1 / u.cm),
             wmax=200 * (1 / u.cm),
             wavelength_min=10,
@@ -589,7 +596,7 @@ def test_get_waverange(*args, **kwargs):
 
     # ...... passed both wavenumber and wavelength > should throw error due to multiple input wavespace values
     with pytest.raises(ValueError):
-        get_waverange(
+        get_wavenumber_range(
             wavenum_min=10,
             wavenum_max=20,
             wmin=100 * (1 / u.cm),
@@ -600,12 +607,12 @@ def test_get_waverange(*args, **kwargs):
         )
 
     # ...... passed neither wavenumber nor wavelength > should return wavenumber after converting wmin/wmax from accompanying unit
-    assert get_waverange(
+    assert get_wavenumber_range(
         wmin=100 * (1 / u.cm), wmax=200 * (1 / u.cm), wunit=Default("cm-1")
     ) == (100.0, 200.0)
 
     assert np.isclose(
-        get_waverange(
+        get_wavenumber_range(
             wmin=1 * u.cm, wmax=2 * u.cm, medium="vacuum", wunit=Default("cm-1")
         ),
         (0.5, 1.0),
@@ -615,15 +622,15 @@ def test_get_waverange(*args, **kwargs):
     # ... 'wmin/wmax' is none
     # ...... wavenumber is passed > should throw error as wunit can only be passed with wmin/wmax
     with pytest.raises(ValueError):
-        get_waverange(wavenum_min=1, wavenum_max=2, wunit="cm")
+        get_wavenumber_range(wavenum_min=1, wavenum_max=2, wunit="cm")
 
     # ...... wavelength is passed > should throw error as wunit can only be passed with wmin/wmax
     with pytest.raises(ValueError):
-        get_waverange(wavelength_min=1, wavelength_max=2, wunit="cm")
+        get_wavenumber_range(wavelength_min=1, wavelength_max=2, wunit="cm")
 
     # ...... passed both wavenumber and wavelength > should throw error as wunit can only be passed with wmin/wmax
     with pytest.raises(ValueError):
-        get_waverange(
+        get_wavenumber_range(
             wavenum_min=1,
             wavenum_max=2,
             wavelength_min=10,
@@ -633,20 +640,22 @@ def test_get_waverange(*args, **kwargs):
 
     # ...... passed neither wavenumber nor wavlength > should throw error
     with pytest.raises(ValueError):
-        get_waverange(wunit="cm")
+        get_wavenumber_range(wunit="cm")
 
     # ... 'wmin/wmax' are passed as values without accompanying units
     # ...... wavenumber is passed > should throw error as only one set of wavespace parameters can be passed
     with pytest.raises(ValueError):
-        get_waverange(wmin=1, wmax=2, wavenum_min=10, wavenum_max=20, wunit="cm")
+        get_wavenumber_range(wmin=1, wmax=2, wavenum_min=10, wavenum_max=20, wunit="cm")
 
     # ...... wavelength is passed > should throw error as only one set of wavespace parameters can be passed
     with pytest.raises(ValueError):
-        get_waverange(wmin=1, wmax=2, wavelength_min=10, wavelength_max=20, wunit="cm")
+        get_wavenumber_range(
+            wmin=1, wmax=2, wavelength_min=10, wavelength_max=20, wunit="cm"
+        )
 
     # ...... passed both wavenumber and wavelength > should throw error as only one set of wavespace parameters can be passed
     with pytest.raises(ValueError):
-        get_waverange(
+        get_wavenumber_range(
             wmin=1,
             wmax=2,
             wavenum_min=10,
@@ -658,22 +667,22 @@ def test_get_waverange(*args, **kwargs):
 
     # ...... passed neither wavenumber nor wavlength > should return wavenumber after converting wmin/wmax with given wunit
     assert np.isclose(
-        get_waverange(wmin=1, wmax=2, wunit="cm"),
+        get_wavenumber_range(wmin=1, wmax=2, wunit="cm"),
         (0.4998637271242577, 0.9997274542485038),
     ).all()
 
-    assert get_waverange(wmin=1, wmax=2, wunit="cm-1") == (1.0, 2.0)
+    assert get_wavenumber_range(wmin=1, wmax=2, wunit="cm-1") == (1.0, 2.0)
 
     # ... 'wmin/wmax' are passed as values with accompanying units
     # ...... wavenumber is passed > should throw error as only one set of wavespace parameters can be passed
     with pytest.raises(ValueError):
-        get_waverange(
+        get_wavenumber_range(
             wmin=1 * u.cm, wmax=2 * u.cm, wavenum_min=10, wavenum_max=20, wunit="cm"
         )
 
     # ...... wavelength is passed > should throw error as only one set of wavespace parameters can be passed
     with pytest.raises(ValueError):
-        get_waverange(
+        get_wavenumber_range(
             wmin=1 * u.cm,
             wmax=2 * u.cm,
             wavelength_min=10,
@@ -683,7 +692,7 @@ def test_get_waverange(*args, **kwargs):
 
     # ...... passed both wavenumber and wavelength > should throw error as only one set of wavespace parameters can be passed
     with pytest.raises(ValueError):
-        get_waverange(
+        get_wavenumber_range(
             wmin=1 * u.cm,
             wmax=2 * u.cm,
             wavenum_min=1,
@@ -695,20 +704,22 @@ def test_get_waverange(*args, **kwargs):
 
     # ...... passed neither wavenumber nor wavlength > should return wavenumber after converting wmin/wmax with given wunit
     assert np.isclose(
-        get_waverange(wmin=1 * u.cm, wmax=2 * u.cm, wunit="cm"),
+        get_wavenumber_range(wmin=1 * u.cm, wmax=2 * u.cm, wunit="cm"),
         (0.4998637271242577, 0.9997274542485038),
     ).all()
 
-    assert get_waverange(wmin=1 * (1 / u.cm), wmax=2 * (1 / u.cm), wunit="cm-1") == (
+    assert get_wavenumber_range(
+        wmin=1 * (1 / u.cm), wmax=2 * (1 / u.cm), wunit="cm-1"
+    ) == (
         1.0,
         2.0,
     )
     assert np.isclose(
-        get_waverange(wmin=1 * u.cm, wmax=2 * u.cm, wunit="cm"),
+        get_wavenumber_range(wmin=1 * u.cm, wmax=2 * u.cm, wunit="cm"),
         (0.4998637271242577, 0.9997274542485038),
     ).all()
     assert np.isclose(
-        get_waverange(
+        get_wavenumber_range(
             wavelength_min=1 * u.cm, wavelength_max=2 * u.cm, wunit=Default("cm-1")
         ),
         (0.4998637271242577, 0.9997274542485038),
@@ -716,18 +727,91 @@ def test_get_waverange(*args, **kwargs):
 
     # ... passed wmin/wmax with units different from wunit > should throw error due to conflicting units
     with pytest.raises(ValueError):
-        get_waverange(wmin=1 * u.cm, wmax=2 * u.cm, wunit="cm-1")
+        get_wavenumber_range(wmin=1 * u.cm, wmax=2 * u.cm, wunit="cm-1")
+
+
+def test_input_wunit(plot=True, *args, **kwargs):
+    """Test spectrum default units are properly when giving inputs in wavenumber
+    or wavelength ; and test that intensity units are consistent too"""
+
+    if not plot:
+        raise ValueError("This test requires plotting")
+
+    import matplotlib.pyplot as plt
+
+    from radis import test_spectrum
+
+    s_from_wavenumber = test_spectrum()
+    s_from_wavenumber.plot()
+
+    assert s_from_wavenumber.c["waveunit"] == "cm-1"
+    assert plt.gca().get_xlabel() == "Wavenumber (cm⁻¹)"
+    assert plt.gca().get_ylabel() == "radiance (mW/cm²/sr/cm⁻¹)"
+
+    assert s_from_wavenumber.units["radiance_noslit"] == "mW/cm2/sr/cm-1"
+
+    #%%
+    # Compute the same spectrum, but with wavelength input (and compare that
+    # areas are the same at the end)
+    wmin = s_from_wavenumber.get_wavelength().min()
+    wmax = s_from_wavenumber.get_wavelength().max()
+
+    import astropy.units as u
+
+    s_fromwl = test_spectrum(wmin=wmin * u.nm, wmax=wmax * u.nm)
+    s_fromwl.plot()
+
+    import matplotlib.pyplot as plt
+
+    assert s_fromwl.c["waveunit"] == "cm-1"  # calculation unit
+    assert s_fromwl.c["default_output_unit"] == "nm"  # default get/plot unit
+    assert plt.gca().get_xlabel() == "Wavelength (nm)"
+    assert plt.gca().get_ylabel() == "radiance (mW/cm²/sr/nm)"
+
+    # Check Python technical caveats do not apply :
+    # ... Check that units have not changed :
+    assert s_from_wavenumber.units["radiance_noslit"] == "mW/cm2/sr/cm-1"
+    # ... Check we're not sharing the same object:
+    assert s_from_wavenumber.units is not s_fromwl.units  # Python problems
+
+    #%%
+    import numpy as np
+
+    print(s_from_wavenumber.get_power())
+    print(s_fromwl.get_power())
+
+    print(s_from_wavenumber.get("radiance_noslit", Iunit="mW/cm2/sr/cm-1")[1])
+    print(s_fromwl.get("radiance_noslit", Iunit="mW/cm2/sr/cm-1")[1])
+
+    # Check power is approximately the same
+    # (note : should be exactly the same actually; see https://github.com/radis/radis/issues/460)
+    # TODO: change assertion with equality once #460 is fixed
+    assert np.isclose(s_from_wavenumber.get_power(), s_fromwl.get_power(), rtol=1e-4)
+
+    # Check algebra & line-of-sight work
+    (s_from_wavenumber // s_fromwl).plot()
+    (s_from_wavenumber > s_fromwl).plot()
+    (
+        s_from_wavenumber.take("radiance_noslit") + s_fromwl.take("radiance_noslit")
+    ).plot()
+    assert np.isclose(
+        (
+            s_from_wavenumber.take("radiance_noslit") + s_fromwl.take("radiance_noslit")
+        ).max(),
+        2 * s_from_wavenumber.take("radiance_noslit").max(),
+    )
 
 
 def _run_testcases(verbose=True, plot=True):
 
+    test_input_wunit()
     test_linestrength_calculations()
     test_export_populations(plot=plot, verbose=verbose)
     test_export_rovib_fractions(plot=plot, verbose=verbose)
     test_populations_CO2_hamiltonian(plot=plot, verbose=verbose)
     test_optically_thick_limit_1iso(plot=plot, verbose=verbose)
     test_optically_thick_limit_2iso(plot=plot, verbose=verbose)
-    test_get_waverange()
+    test_get_wavenumber_range()
 
 
 if __name__ == "__main__":
diff --git a/radis/test/lbl/test_broadening.py b/radis/test/lbl/test_broadening.py
index 07b7ede2c..b03898ef6 100644
--- a/radis/test/lbl/test_broadening.py
+++ b/radis/test/lbl/test_broadening.py
@@ -493,7 +493,7 @@ def test_broadening_LDM_noneq(verbose=True, plot=False, *args, **kwargs):
     )
 
     # Compare
-    res = get_residual(s_ldm_eq, s_ldm_noneq, "radiance_noslit")
+    res = get_residual(s_ldm_eq, s_ldm_noneq, "radiance_noslit", Iunit="mW/cm2/sr/nm")
 
     if verbose:
         print("Residual:", res)
@@ -962,4 +962,5 @@ def _run_testcases(plot=False, verbose=True, *args, **kwargs):
 
 
 if __name__ == "__main__":
-    printm("test_broadening: ", _run_testcases(plot=True, verbose=True, debug=False))
+    # printm("test_broadening: ", _run_testcases(plot=True, verbose=True, debug=False))
+    printm("Testing broadening:", pytest.main(["test_broadening.py", "--pdb"]))
diff --git a/radis/test/lbl/test_calc.py b/radis/test/lbl/test_calc.py
index 8be7e396e..883f30c1d 100644
--- a/radis/test/lbl/test_calc.py
+++ b/radis/test/lbl/test_calc.py
@@ -70,13 +70,9 @@ def test_calc_spectrum(verbose=True, plot=True, warnings=True, *args, **kwargs):
     (just used here as a performance monitoring)
 
     - neq 0.9.20: 18.7s
-
     - neq 0.9.20*: 15.4s   (removed 2nd loop on 1st isotope because of groupby().apply())
-
     - neq 0.9.20**: 11.7s  (after replacing fill_Evib with map() )
-
     - neq 0.9.21: 9.4s     (improve Qrot / nrot fetching performance)
-
     - neq 0.9.22: 8.4s
 
     Starting from RADIS 1.0.1, the test is run on [HITRAN-2020]_, which
@@ -88,6 +84,7 @@ def test_calc_spectrum(verbose=True, plot=True, warnings=True, *args, **kwargs):
 
     - radis 0.9.20 : 2.49 s    on [HITRAN-2020]
                      4.05 s    on [CDSD-HITEMP-JMIN]
+    - radis 0.12.2 : 1.1s      on [HITRAN-2020]   (most time spent in the Evib, Erot look-up)
 
     """
     if verbose:
@@ -101,6 +98,8 @@ def test_calc_spectrum(verbose=True, plot=True, warnings=True, *args, **kwargs):
     s = calc_spectrum(
         wavelength_min=4165,
         wavelength_max=4200,
+        # wavenum_min=1e7/4200,
+        # wavenum_max=1e7/4165,
         #                          databank='CDSD-HITEMP-JMIN',
         databank="hitran",  # not appropriate for these temperatures, but convenient for automatic testing
         # databank="HITRAN-CO2-TEST",  # to use, but only has 1 isotope. TODO add new test file with 2 isotopes
@@ -128,9 +127,9 @@ def test_calc_spectrum(verbose=True, plot=True, warnings=True, *args, **kwargs):
     s.apply_slit((2, 2.5), "nm", shape="trapezoidal")
 
     if plot:
-        s.plot(wunit="nm")
+        s.plot(wunit="nm", Iunit="mW/cm2/sr/nm")
 
-    w, I = s.get("radiance", wunit="nm")
+    w, I = s.get("radiance", wunit="nm", Iunit="mW/cm2/sr/nm")
     # Compare against hardcoded results (neq 0.9.22, 28/06/18)
     #        I_ref = np.array([0.28694463, 0.29141711, 0.32461613, 0.32909566, 0.21939511, 0.18606445,
     #                          0.19740763, 0.16948599, 0.16780345, 0.15572173, 0.16770853, 0.14966064,
@@ -792,6 +791,34 @@ def test_calc_spectrum_multiple_molecules_wstep_auto(
     assert s.get_conditions()["wstep"] in ("N/A", 0.013)
 
 
+def test_check_wavelength_range(verbose=True, warnings=True, *args, **kwargs):
+    """Check that input wavelength is correctly taken into account.
+    See https://github.com/radis/radis/issues/214
+    """
+    if verbose:
+        printm("Testing calc_spectrum wavelength range")
+
+    wstep = 0.01
+
+    s = calc_spectrum(
+        wavelength_min=4348,  # nm
+        wavelength_max=5000,
+        molecule="CO",
+        isotope="1,2,3",
+        pressure=1.01325,  # bar
+        Tvib=1700,  # K
+        Trot=1700,  # K
+        databank="HITRAN-CO-TEST",
+        wstep=wstep,
+    )
+    w, I = s.get("radiance_noslit", wunit="nm", Iunit="mW/sr/cm2/nm")
+
+    assert np.isclose(w.min(), 4348, atol=wstep)
+    assert np.isclose(w.max(), 5000, atol=wstep)
+
+    return True
+
+
 def _run_testcases(plot=True, verbose=True, warnings=True, *args, **kwargs):
 
     # Test sPlanck and conversion functions
@@ -824,33 +851,7 @@ def _run_testcases(plot=True, verbose=True, warnings=True, *args, **kwargs):
     test_calc_spectrum_multiple_molecules_inputerror()
     test_calc_spectrum_multiple_molecules_wstep_auto()
 
-    return True
-
-
-def check_wavelength_range(verbose=True, warnings=True, *args, **kwargs):
-    """Check that input wavelength is correctly taken into account.
-    See https://github.com/radis/radis/issues/214
-    """
-    if verbose:
-        printm("Testing calc_spectrum wavelength range")
-
-    wstep = 0.01
-
-    s = calc_spectrum(
-        wavelength_min=4348,  # nm
-        wavelength_max=5000,
-        molecule="CO",
-        isotope="1,2,3",
-        pressure=1.01325,  # bar
-        Tvib=1700,  # K
-        Trot=1700,  # K
-        databank="HITRAN-CO-TEST",
-        wstep=wstep,
-    )
-    w, I = s.get("radiance_noslit", wunit="nm", Iunit="mW/sr/cm2/nm")
-
-    assert np.isclose(w.min(), 4348, atol=wstep)
-    assert np.isclose(w.max(), 5000, atol=wstep)
+    test_check_wavelength_range()
 
     return True
 
diff --git a/radis/test/levels/test_partfunc.py b/radis/test/levels/test_partfunc.py
index 795d37173..b400eef19 100644
--- a/radis/test/levels/test_partfunc.py
+++ b/radis/test/levels/test_partfunc.py
@@ -1049,4 +1049,4 @@ def _run_testcases(verbose=True, warnings=True, *args, **kwargs):
 
 if __name__ == "__main__":
     # printm("Testing parfunc: {0}".format(_run_testcases()))
-    printm("Testing factory:", pytest.main(["test_partfunc.py", "--pdb"]))
+    printm("Testing partfunc.py:", pytest.main(["test_partfunc.py", "--pdb"]))
diff --git a/radis/test/spectrum/test_compare.py b/radis/test/spectrum/test_compare.py
index 15b53ee57..087b7974d 100644
--- a/radis/test/spectrum/test_compare.py
+++ b/radis/test/spectrum/test_compare.py
@@ -115,7 +115,9 @@ def test_get_residual():
         mole_fraction=0.1,
     )
     s1.apply_slit(1, "nm")
-    w1, I1 = s1.get("radiance", copy=False, wunit=s1.get_waveunit())
+    w1, I1 = s1.get(
+        "radiance", copy=False, wunit=s1.get_waveunit(), Iunit="mW/cm2/sr/nm"
+    )
 
     # Fake experimental spectrum (identical)
     s_expe_1 = experimental_spectrum(w1, I1, Iunit="mW/cm2/sr/nm", wunit="cm-1")
@@ -126,12 +128,13 @@ def test_get_residual():
         # ignore_nan=False,
         normalize=True,
     )
-    assert residual1 == 0
+    assert residual1 < 1e-4
 
     # Fake experimental spectrum with a new unit (but identical)
     I2 = I1 * 1e-3
     s_expe_2 = experimental_spectrum(w1, I2, Iunit="W/cm2/sr/nm", wunit="cm-1")
 
+    assert np.isclose(s_expe_1.max(), s_expe_2.max())
     residual2 = get_residual(
         s1,
         s_expe_2,
@@ -149,7 +152,7 @@ def test_get_residual():
             normalize_how=bool_how,
         )
         # print('residu2 ={}'.format(residual2))
-        assert residual2_bis < 1e-13
+        assert residual2_bis < 1e-3
 
     # Fake experimental spectrum with a nan
     I3 = I1.copy()
@@ -157,7 +160,7 @@ def test_get_residual():
     s_expe_3 = experimental_spectrum(w1, I3, Iunit="W/cm2/sr/nm", wunit="cm-1")
 
     # Intrestingly, the "mean" method seams to be not adapted for spectra with nans
-    criterion = [1e-10, 2e-6, 2e-2]
+    criterion = [7e-5, 2e-6, 2e-2]
     for index, bool_how in enumerate(["max", "area", "mean"]):
         residual3 = get_residual(
             s1,
@@ -183,4 +186,5 @@ def _run_testcases(plot=True, verbose=True, warnings=True, *args, **kwargs):
 
 
 if __name__ == "__main__":
-    print("Test_compare.py: ", _run_testcases())
+    # print("Test_compare.py: ", _run_testcases())
+    test_get_residual()
diff --git a/radis/test/spectrum/test_operations.py b/radis/test/spectrum/test_operations.py
index bf0fffc9b..4bae6cddd 100644
--- a/radis/test/spectrum/test_operations.py
+++ b/radis/test/spectrum/test_operations.py
@@ -10,16 +10,7 @@
 
 from radis.los import MergeSlabs, SerialSlabs
 from radis.spectrum.compare import get_diff, plot_diff
-from radis.spectrum.operations import (
-    Radiance,
-    Radiance_noslit,
-    Transmittance_noslit,
-    add_constant,
-    crop,
-    multiply,
-    offset,
-    sub_baseline,
-)
+from radis.spectrum.operations import add_constant, crop, multiply, offset, sub_baseline
 from radis.test.utils import getTestFile
 from radis.tools.database import load_spec
 
@@ -78,35 +69,34 @@ def test_cut_recombine(verbose=True, *args, **kwargs):
 @pytest.mark.fast
 def test_invariants(*args, **kwargs):
     """Ensures adding 0 or multiplying by 1 does not change the spectra"""
-    from radis import load_spec
-    from radis.test.utils import getTestFile
 
     s = load_spec(getTestFile("CO_Tgas1500K_mole_fraction0.01.spec"))
-    s.update()
-    s = Radiance_noslit(s)
+    s = s.take("radiance_noslit")
 
     assert s.compare_with(
         add_constant(s, 0, "W/cm2/sr/nm"), plot=False, spectra_only="radiance_noslit"
     )
     assert s.compare_with(multiply(s, 1), plot=False, spectra_only="radiance_noslit")
+    assert (s + 0) == s
 
     assert 3 * s / 3 == s
     assert (1 + s) - 1 == s
 
+    assert s.max().value == s._q["radiance_noslit"].max()
+
 
 @pytest.mark.fast
 def test_operations_inplace(verbose=True, *args, **kwargs):
 
-    from radis.spectrum.operations import Radiance_noslit
+    from radis.phys.units import Unit
 
     s = load_spec(getTestFile("CO_Tgas1500K_mole_fraction0.01.spec"), binary=True)
-    s.update("radiance_noslit", verbose=False)
-    s = Radiance_noslit(s)
+    s = s.take("radiance_noslit")
 
     # Add 1, make sure it worked
-    I_max = s.get("radiance_noslit")[1].max()
-    s += 1
-    assert s.get("radiance_noslit")[1].max() == I_max + 1
+    I_max = s.get("radiance_noslit", Iunit="mW/cm2/sr/nm")[1].max()
+    s += 1 * Unit("mW/cm2/sr/nm")
+    assert s.get("radiance_noslit", Iunit="mW/cm2/sr/nm")[1].max() == I_max + 1
     if verbose:
         print("test_operations: s += 1: OK")
 
@@ -149,11 +139,11 @@ def test_multiplyAndAddition(verbose=True, plot=False, *args, **kwargs):
     s = load_spec(getTestFile("CO_Tgas1500K_mole_fraction0.01.spec"), binary=True)
     s.update("radiance_noslit", verbose=False)
     s.apply_slit(0.1)
-    s = Radiance(s)
+    s = s.take("radiance")
     assert s.units["radiance"] == "mW/cm2/sr/nm"
 
     s_bis = add_constant(s, 1, "mW/cm2/sr/nm")
-    w, Idiff = get_diff(s_bis, s, "radiance")
+    w, Idiff = get_diff(s_bis, s, "radiance", Iunit="mW/cm2/sr/nm")
     test = Idiff[1] - 1
     assert np.all(test < 1e-10)
 
@@ -194,7 +184,7 @@ def test_other_algebraic_operations(verbose=True, plot=False, *args, **kwargs):
     (s // s).plot(nfig="same")
 
     # Test substraction of Spectra
-    s_tr = Transmittance_noslit(s)
+    s_tr = s.take("transmittance_noslit")
     assert (s_tr - 1.0 * s_tr).get_integral("transmittance_noslit") == 0
 
     # TODO: add test
@@ -208,7 +198,7 @@ def test_other_algebraic_operations(verbose=True, plot=False, *args, **kwargs):
         2 * s
 
     s.apply_slit(0.1, "nm")
-    s_rad = Radiance(s)
+    s_rad = s.take("radiance")
 
     # Test multiplication with float
     s.plot(lw=2, nfig="Multiplication (by scalar): 2*s", wunit="nm")
@@ -233,7 +223,7 @@ def test_TestBaseline(plot=False, *args, **kwargs):
     s = load_spec(getTestFile("CO_Tgas1500K_mole_fraction0.01.spec"), binary=True)
     s.update("radiance_noslit", verbose=False)
     s.apply_slit(0.1)
-    s = Radiance_noslit(s)
+    s = s.take("radiance_noslit")
     assert s.units["radiance"] == "mW/cm2/sr/nm"
 
     s2 = sub_baseline(s, 2e-4, -2e-4)
@@ -246,32 +236,33 @@ def test_TestBaseline(plot=False, *args, **kwargs):
 @pytest.mark.fast
 def test_dimensioned_operations(*args, **kwargs):
 
-    import astropy.units as u
     import numpy as np
 
-    from radis import Radiance, load_spec
-    from radis.spectrum import sub_baseline
-    from radis.test.utils import getTestFile
+    import radis.phys.units as u
 
     # Generate the equivalent of an experimental spectrum
     s = load_spec(getTestFile(r"CO_Tgas1500K_mole_fraction0.01.spec"), binary=True)
     s.update()  # add radiance, etc.
     s.apply_slit(0.5)  # nm
-    s = Radiance(s)
+    s = s.take("radiance")
 
     # Test
     assert s.units["radiance"] == "mW/cm2/sr/nm"
-    Imax = s.get("radiance", trim_nan=True)[1].max()
+    Imax = s.get("radiance", Iunit="mW/cm2/sr/nm", trim_nan=True)[1].max()
 
     # add a baseline
     s += 0.1 * u.Unit("W/cm2/sr/nm")
 
-    assert np.isclose(s.get("radiance", trim_nan=True)[1].max(), Imax + 100)
+    assert np.isclose(
+        s.get("radiance", Iunit="mW/cm2/sr/nm", trim_nan=True)[1].max(), Imax + 100
+    )
 
     # remove a baseline (we could also have used s=-0.1, but we're testing another function here)
     s = sub_baseline(s, 0.1 * u.Unit("W/cm2/sr/nm"), 0.1 * u.Unit("W/cm2/sr/nm"))
 
-    assert np.isclose(s.get("radiance", trim_nan=True)[1].max(), Imax)
+    assert np.isclose(
+        s.get("radiance", Iunit="mW/cm2/sr/nm", trim_nan=True)[1].max(), Imax
+    )
 
     # Test division
     # Example : a manual normalization
@@ -294,9 +285,6 @@ def test_dimensioned_operations(*args, **kwargs):
 @pytest.mark.fast
 def test_get_baseline(plot=False, *args, **kwargs):
 
-    from radis import load_spec
-    from radis.test.utils import getTestFile
-
     s = load_spec(getTestFile(r"CO2_measured_spectrum_4-5um.spec"), binary=True)
 
     from radis.spectrum.operations import get_baseline
@@ -315,6 +303,33 @@ def test_get_baseline(plot=False, *args, **kwargs):
         plt.legend()
 
 
+@pytest.mark.fast
+def test_resample_even(*args, **kwargs):
+    """Test :py:meth:`~radis.spectrum.spectrum.Spectrum.resample_even` ,
+    and that :py:meth:`~radis.spectrum.spectrum.Spectrum.apply_slit` properly fails
+    if Spectrum is not evenly distributed"""
+
+    s = load_spec(getTestFile(r"CO_Tgas1500K_mole_fraction0.01.spec"), binary=True)
+    s.update()
+
+    s_nm = s.resample(s.get_wavelength(), "nm", inplace=False)
+    assert s_nm.c["waveunit"] == "nm"
+
+    # except failure if applying a slit function:
+    with pytest.raises(ValueError) as err:
+        s_nm.apply_slit(1.5, "nm")
+        assert "Spectrum is not evenly spaced" in str(err)
+
+    # test that it works after resampling :
+    s_nm.resample_even()
+    s_nm.apply_slit(1.5, "nm")
+
+    from radis.misc.arrays import evenly_distributed, evenly_distributed_fast
+
+    assert evenly_distributed(s_nm.get_wavelength())
+    assert evenly_distributed_fast(s_nm.get_wavelength())
+
+
 def _run_testcases(verbose=True, plot=False, *args, **kwargs):
     """Test procedures"""
 
@@ -329,6 +344,7 @@ def _run_testcases(verbose=True, plot=False, *args, **kwargs):
     test_TestBaseline(verbose=verbose, plot=plot, *args, **kwargs)
     test_dimensioned_operations(*args, **kwargs)
     test_get_baseline(plot=plot, *args, **kwargs)
+    test_resample_even()
 
     return True
 
diff --git a/radis/test/spectrum/test_rescale.py b/radis/test/spectrum/test_rescale.py
index 6dc923736..f3d52b504 100644
--- a/radis/test/spectrum/test_rescale.py
+++ b/radis/test/spectrum/test_rescale.py
@@ -124,6 +124,231 @@ def test_get_recompute(verbose=True, *args, **kwargs):
     )
 
 
+def test_rescale_vs_direct_computation(verbose=True, *args, **kwargs):
+    """Compare spectral arrays initially computed with recomputed spectral arrays
+
+    Notes
+    -----
+    Use verbose=2 to get DEBUG_MODE output
+    """
+    if verbose >= 2:
+        DEBUG_MODE = radis.config["DEBUG_MODE"]
+        radis.config["DEBUG_MODE"] = True
+
+    # %%
+    # Equilibrium Spectrum
+    # -------------------
+    from radis import calc_spectrum
+
+    s = calc_spectrum(
+        wavelength_min=4165,
+        wavelength_max=4200,
+        mole_fraction=1e-3,
+        path_length=0.3,
+        molecule="CO2",
+        isotope="1",
+        databank="HITRAN-CO2-TEST",
+        Tgas=600,
+        verbose=False,
+    )
+
+    s2 = s.copy()
+
+    import numpy as np
+
+    # 1. Compare initially computed with recomputed
+    for var in ["abscoeff", "radiance_noslit"]:
+
+        # default
+        Imax_computed = s.take(var).max()
+        # delete existing variables
+        del s2._q[var]
+        # recompute:
+        Imax_recomputed = s2.take(var).max()
+        # compare:
+        if verbose:
+            print("==>", var, Imax_computed)
+            print("==>", var, Imax_recomputed)
+        assert np.isclose(Imax_computed, Imax_recomputed)
+
+    # 2. Do this systematically, for all configurations where `s` can be recomputed from.
+    from radis.spectrum.rescale import _build_update_graph
+
+    """Typical :py:func:`radis.spectrum.rescale._build_update_graph(s)` result :
+    ::
+        {
+            "transmittance_noslit": [
+                ["absorbance"],
+                ["abscoeff"],
+                ["absorbance"],
+                ["emisscoeff"],
+                ["emissivity_noslit"],
+                ["transmittance"],
+                ["radiance"],
+                ["radiance_noslit"],
+                ["xsection"],
+            ],
+            "absorbance": [
+                ["transmittance_noslit"],
+                ["abscoeff"],
+                ["abscoeff"],
+                ["emisscoeff"],
+                ["emissivity_noslit"],
+                ["transmittance"],
+                ["radiance"],
+                ["radiance_noslit"],
+                ["transmittance_noslit"],
+                ["xsection"],
+            ],
+            "abscoeff": [
+                ["absorbance"],
+                ["xsection"],
+                ["absorbance"],
+                ["emisscoeff"],
+                ["emissivity_noslit"],
+                ["transmittance"],
+                ["radiance"],
+                ["radiance_noslit"],
+                ["transmittance_noslit"],
+                ["xsection"],
+            ],
+            "radiance_noslit": [
+                ["emisscoeff", "abscoeff"],
+                ["abscoeff"],
+                ["absorbance"],
+                ["emisscoeff"],
+                ["emissivity_noslit"],
+                ["transmittance"],
+                ["radiance"],
+                ["transmittance_noslit"],
+                ["xsection"],
+            ],
+            "emisscoeff": [
+                ["radiance_noslit", "abscoeff"],
+                ["abscoeff"],
+                ["absorbance"],
+                ["emissivity_noslit"],
+                ["transmittance"],
+                ["radiance"],
+                ["radiance_noslit"],
+                ["transmittance_noslit"],
+                ["xsection"],
+            ],
+            "xsection": [["abscoeff"]],
+            "emissivity_noslit": [
+                ["abscoeff"],
+                ["absorbance"],
+                ["emisscoeff"],
+                ["transmittance"],
+                ["radiance"],
+                ["radiance_noslit"],
+                ["transmittance_noslit"],
+                ["xsection"],
+            ],
+        }
+    """
+    from radis.misc.basics import all_in
+
+    for var, vars_needed_list in _build_update_graph(s).items():
+        for vars_needed in vars_needed_list:
+            if all_in(vars_needed, s.get_vars()):
+                if verbose:
+                    print(f"Recomputing {var} from {vars_needed} only")
+                s2 = s.copy()
+                for v in list(s2.get_vars()):
+                    if v not in vars_needed:
+                        del s2._q[v]
+                # recompute & compare:
+                try:
+                    s2.get(var)
+                except NotImplementedError:
+                    pass
+                else:
+                    assert np.isclose(s.take(var).max(), s2.take(var).max())
+                    if verbose:
+                        print(
+                            f"Checked {var} recomputed from {vars_needed} is the same"
+                        )
+
+    #%%
+    # Nonequilibrium spectrum, computed in wavenumber
+    # -----------------------------------------------
+
+    from radis import calc_spectrum
+
+    s = calc_spectrum(
+        wavenum_max=1e7 / 4165,
+        wavenum_min=1e7 / 4200,
+        mole_fraction=1e-3,
+        path_length=0.3,
+        molecule="CO2",
+        isotope="1",
+        databank="HITRAN-CO2-TEST",
+        Tvib=2000,
+        Trot=1000,
+        verbose=False,
+    )
+
+    import numpy as np
+
+    # Compare initially computed with recomputed
+    for var in ["transmittance_noslit", "emisscoeff", "radiance_noslit"]:
+
+        s2 = s.copy()
+
+        # default
+        Imax_computed = s.take(var).max()
+        # delete existing variables
+        del s2._q[var]
+        # recompute:
+        Imax_recomputed = s2.take(var).max()
+        # compare:
+        if verbose:
+            print("==>", var, Imax_computed)
+            print("==>", var, Imax_recomputed)
+        assert np.isclose(Imax_computed, Imax_recomputed)
+
+    """Typical :py:func:`radis.spectrum.rescale._build_update_graph(s)` result :
+    ::
+
+        {
+            "transmittance_noslit": [["absorbance"]],
+            "absorbance": [["transmittance_noslit"], ["abscoeff"]],
+            "abscoeff": [["absorbance"], ["xsection"]],
+            "radiance_noslit": [["emisscoeff", "abscoeff"]],
+            "emisscoeff": [["radiance_noslit", "abscoeff"]],
+            "xsection": [["abscoeff"]],
+        }
+    """
+
+    for var, vars_needed_list in _build_update_graph(s).items():
+        for vars_needed in vars_needed_list:
+            if all_in(vars_needed, s.get_vars()):
+                if verbose:
+                    print(f"Recomputing {var} from {vars_needed} only")
+                s2 = s.copy()
+                for v in list(s2.get_vars()):
+                    if v not in vars_needed:
+                        del s2._q[v]
+                # recompute & compare:
+                try:
+                    s2.get(var)
+                except NotImplementedError:
+                    pass
+                else:
+                    assert np.isclose(s.take(var).max(), s2.take(var).max())
+                    if verbose:
+                        print(
+                            f"Checked {var} recomputed from {vars_needed} is the same"
+                        )
+
+    #
+    if verbose >= 2:
+        radis.config["DEBUG_MODE"] = DEBUG_MODE
+
+    #%%
+
+
 def test_recompute_equilibrium(verbose=True, warnings=True, plot=True, *args, **kwargs):
     """Test that spectral quantities recomputed under equilibrium assumption
     yields the same output as with non equilibrium routines when Tvib = Trot"""
@@ -331,6 +556,7 @@ def _run_all_tests(verbose=True, warnings=True, *args, **kwargs):
     test_compression(verbose=verbose, warnings=warnings, *args, **kwargs)
     test_update_transmittance(verbose=verbose, warnings=warnings, *args, **kwargs)
     test_get_recompute(verbose=verbose, warnings=warnings, *args, **kwargs)
+    test_rescale_vs_direct_computation(verbose=verbose, *args, **kwargs)
     test_recompute_equilibrium(verbose=verbose, warnings=warnings, *args, **kwargs)
     test_rescale_all_quantities(verbose=verbose, *args, **kwargs)
     test_xsections(*args, **kwargs)
diff --git a/radis/test/spectrum/test_spectrum.py b/radis/test/spectrum/test_spectrum.py
index 7d0e1d0bd..631b69ea4 100644
--- a/radis/test/spectrum/test_spectrum.py
+++ b/radis/test/spectrum/test_spectrum.py
@@ -621,4 +621,23 @@ def _run_testcases(
 
 
 if __name__ == "__main__":
-    print(("Test spectrum: ", _run_testcases(debug=False, close_plots=False)))
+    # print(("Test spectrum: ", _run_testcases(debug=False, close_plots=False)))
+
+    import radis
+
+    s = radis.test_spectrum()
+    assert s.c["default_output_unit"] == "cm-1"
+
+    assert (s.get("abscoeff")[0] == s.get_wavenumber()).all()
+
+    s.c["default_output_unit"] = "nm"
+
+    assert (s.get("abscoeff")[0] == s.get_wavelength()).all()
+    s.c["default_output_unit"] = "nm_vac"
+    assert (s.get("abscoeff")[0] == s.get_wavelength(medium="vacuum")).all()
+
+    #%%
+    from radis import load_spec
+    from radis.test.utils import getTestFile
+
+    s_exp = load_spec(getTestFile("CO2_measured_spectrum_4-5um.spec"), binary=True)
diff --git a/radis/test/test_examples.py b/radis/test/test_examples.py
index 90b132c53..516eb76d3 100644
--- a/radis/test/test_examples.py
+++ b/radis/test/test_examples.py
@@ -19,6 +19,7 @@
 ]
 
 
+@pytest.fixture(scope="session")
 @pytest.mark.parametrize("script", scripts)
 def test_script_execution(script):
     """Run all examples in `radis/examples <https://github.com/radis/radis/tree/develop/examples>`__"""
diff --git a/radis/test/tools/test_slit.py b/radis/test/tools/test_slit.py
index f7a05e65b..1823cd12d 100644
--- a/radis/test/tools/test_slit.py
+++ b/radis/test/tools/test_slit.py
@@ -402,7 +402,7 @@ def test_normalisation_mode(plot=True, close_plots=True, verbose=True, *args, **
     s.update()
     # spectrum convolved with area=1
     s.apply_slit(FWHM, norm_by="area", plot_slit=plot)
-    w_area, I_area = s.get("radiance")
+    w_area, I_area = s.get("radiance", wunit="nm")
     if plot:
         fig = plt.figure(fig_prefix + "Spectrum in cm-1 + slit in nm")
         fig.clear()
@@ -715,6 +715,7 @@ def test_resampling(rtol=1e-2, verbose=True, plot=True, warnings=True, *args, **
         w_nm, T_nm, wunit="nm"
     )  # a new spectrum stored in nm
     # sCO_nm = theoretical_spectrum(w_nm, I_nm, wunit='nm', Iunit='mW/cm2/sr/nm') #  a new spectrum stored in nm
+    sCO_nm.resample_even()
 
     if plot:
         fig = plt.figure(fig_prefix + "auto-resampling")
@@ -840,3 +841,4 @@ def _run_testcases(plot=True, close_plots=False, verbose=True, *args, **kwargs):
 
 if __name__ == "__main__":
     print(("Testing slit.py: ", _run_testcases(plot=True)))
+    # printm("Testing slit.py:", pytest.main(["test_slit.py", "--pdb"]))
diff --git a/radis/test/utils.py b/radis/test/utils.py
index 228689c2d..395703e27 100644
--- a/radis/test/utils.py
+++ b/radis/test/utils.py
@@ -180,6 +180,10 @@ def test_spectrum(**kwargs):
 
     conditions.update(kwargs)
 
+    if "wmin" in kwargs and "wmax" in kwargs:
+        conditions.pop("wavenum_min")
+        conditions.pop("wavenum_max")
+
     s = calc_spectrum(**conditions)
     return s