Repository: radis/radis (language: python)

Issue:
"emissivity_noslit" convolved by the slit "on the go"
### 🎯 Describe what you'd like
<!--- A clear and concise description of what you want to happen. -->
The request is to have the possibility to get (and plot) the variable "emissivity_noslit" convoluted by the slit "on the go".

According to @erwanp , in the future Radis will only store the non convoluted variables and the slit function, and will calculate the convoluted variables "on the go" when asked by the user (through "get" or "plot" methods for example).
This could be the first step toward the future handling of the slit by Radis. 

### Why "emissivity" convoluted by the slit ?
After the abel inversion of a calibrated experimental spectrum (in W/(cm2.sr.nm) representing the power coming from the source received on the surface of the camera by nanometer, we obtain an abel inverted spectrum in W/(cm3.sr.nm).
This is to be compared to the emissivity of the source.

In radis the emissivity is not convoluted by the slit, but the abel inverted spectrum has a slit applied by (re)construction. 

At this point in Radis, to compare the abel inverted spectrum to a simulation, we can ask for a spectrum convoluted by the slit (in W/(cm2.sr.nm) for a path length of 1 cm, thus virtually having the result in W/(cm3.sr.nm).
However we can not plot the 2 spectra on the same figure without using the 'force' mode...
And the method to get what we want is a bit tricky...

### 💡 Possible implementations
<!--- Not obligatory, but suggest an idea for implementing addition or change -->

---
<sub>*Tip : you can vote for [New Features](https://github.com/radis/radis/issues/79) on https://feathub.com/radis/radis*</sub>



Discussion / hints:
Note that as early as this morning, some users using RADIS for heat transfer also wanted to be able to use a rectangular slit to create a moving-average of a high-resolution absorption coefficient spectrum (`abscoeff`), and generate a narrow-band-model spectrum. This would be very useful for them too !!

I'm labelling it as a Good First Issue, as most methods/functions are already implemented, and should just be properly re-ordered.

Fix patch:
diff --git a/docs/spectrum/spectrum.rst b/docs/spectrum/spectrum.rst
index 10072b675..959cdfbd7 100644
--- a/docs/spectrum/spectrum.rst
+++ b/docs/spectrum/spectrum.rst
@@ -53,7 +53,6 @@ community chats on `Gitter <https://gitter.im/radis-radiation/community>`__ or
 How to generate a Spectrum?
 ===========================
 
-
 Calculate a Spectrum
 --------------------
 
@@ -80,25 +79,30 @@ or calculates it and stores it if you didn't. Very useful for spectra that
 require long computational times!
 
 
-Initialize from text file or arrays
------------------------------------
-
-Spectrum objects can also be generated from NumPy arrays or files.
+Initialize from Python arrays
+-----------------------------
 
-From numpy arrays, use :py:meth:`~radis.spectrum.spectrum.Spectrum.from_array` ::
+The standard way to build a Radis Spectrum is from a dictionary of
+:py:mod:`numpy` arrays::
 
-    # w, T are two numpy arrays
+    # w, k, I are numpy arrays for wavenumbers, absorption coefficient, and radiance.
     from radis import Spectrum
-    s = Spectrum.from_array(w, T, 'transmittance_noslit',
-                               waveunit='nm', unit='') # adimensioned
+    s = Spectrum({"wavenumber":w, "abscoeff":k, "radiance_noslit":I},
+                 units={"radiance_noslit":"mW/cm2/sr/nm", "abscoeff":"cm-1"})
+Or::
 
+    s = Spectrum({"abscoeff":(w,k), "radiance_noslit":(w,I)},
+                 wunit="cm-1"
+                 units={"radiance_noslit":"mW/cm2/sr/nm", "abscoeff":"cm-1"})
 
-From a file, use :py:meth:`~radis.spectrum.spectrum.Spectrum.from_txt` ::
+You can also use the :py:meth:`~radis.spectrum.spectrum.Spectrum.from_array`
+convenience function::
 
-    # 'exp_spectrum.txt' contains a spectrum
+    # w, T are two numpy arrays
     from radis import Spectrum
-    s = Spectrum.from_txt('exp_spectrum.txt', 'radiance',
-                               waveunit='nm', unit='mW/cm2/sr/nm')
+    s = Spectrum.from_array(w, T, 'transmittance_noslit',
+                               wunit='nm', unit='') # adimensioned
+
 
 Convenience functions have been added to handle the usual cases:
 :func:`~radis.spectrum.models.calculated_spectrum`,
@@ -112,25 +116,20 @@ Convenience functions have been added to handle the usual cases:
     s3 = experimental_spectrum(w, I, wunit='nm', Iunit='W/cm2/sr/nm')   # creates 'radiance'
 
 
-Generate a Blackbody (Planck) function object
----------------------------------------------
-
-In RADIS you can either use the :func:`~radis.phys.blackbody.planck` and
-:func:`~radis.phys.blackbody.planck_wn` functions that generate Planck
-radiation arrays for wavelength and wavenumber, respectively.
+Initialize from a text file
+---------------------------
 
-Or, you can use the :func:`~radis.phys.blackbody.sPlanck` function that
-returns a :class:`~radis.spectrum.spectrum.Spectrum` object, with all
-the associated methods (add in a line-of-sight, compare, etc.)
+Spectrum objects can also be generated directly from a text file.
 
-Example::
+From a file, use :py:meth:`~radis.spectrum.spectrum.Spectrum.from_txt` ::
 
-    s = sPlanck(wavelength_min=3000, wavelength_max=50000,
-                T=288, eps=1)
-    s.plot()
+    # 'exp_spectrum.txt' contains a spectrum
+    from radis import Spectrum
+    s = Spectrum.from_txt('exp_spectrum.txt', 'radiance',
+                               wunit='nm', unit='mW/cm2/sr/nm')
 
-
-.. minigallery:: radis.sPlanck
+It is, however, recommended to use the RADIS ``.spec`` json format to store
+and load arrays :
 
 
 Load from a .spec file
@@ -157,6 +156,43 @@ Spectrum conditions and get Spectrum that suits specific parameters
 .. minigallery:: radis.load_spec
 
 
+Want a quick spectrum
+---------------------
+
+You need a spectrum immediatly, to run some tests ? Use :py:func:`~radis.test.utils.test_spectrum` ::
+
+    s = radis.test_spectrum()
+    s.apply_slit(0.5, 'nm')
+    s.plot('radiance')
+
+This returns the CO spectrum from the :ref:`first documentation example <label_first_example>`
+
+.. figure:: examples/co_spectrum_700K.png
+:scale: 60 %
+
+
+Generate a Blackbody (Planck) function object
+---------------------------------------------
+
+In RADIS you can either use the :func:`~radis.phys.blackbody.planck` and
+:func:`~radis.phys.blackbody.planck_wn` functions that generate Planck
+radiation arrays for wavelength and wavenumber, respectively.
+
+Or, you can use the :func:`~radis.phys.blackbody.sPlanck` function that
+returns a :class:`~radis.spectrum.spectrum.Spectrum` object, with all
+the associated methods (add in a line-of-sight, compare, etc.)
+
+Example::
+
+    s = sPlanck(wavelength_min=3000, wavelength_max=50000,
+                T=288, eps=1)
+    s.plot()
+
+.. minigallery:: radis.sPlanck
+
+
+
+
 .. _label_spectral_arrays:
 
 Spectral Arrays
@@ -211,7 +247,7 @@ Custom spectral arrays with arbitrary units can be defined when creating a Spect
 
     # w, I are two numpy arrays
     s = Spectrum.from_array(w, I, 'irradiance',
-                               waveunit='nm', unit='w/cm2/nm')
+                               wunit='nm', unit='w/cm2/nm')
 
 Although not recommended, it is also possible to directly edit the dictionary containing the objects.
 For instance, this is done in
@@ -362,11 +398,11 @@ Example of output::
 Know if a spectrum has nan
 --------------------------
 
-:py:func:`~radis.spectrum.utils.has_nan` looks over all spectral
+:py:meth:`~radis.spectrum.spectrum.Spectrum.has_nan` looks over all spectral
 quantities. ``print(s)`` will also show the number of nans per quantity ::
 
-    from radis.spectrum.utils import has_nan
-    has_nan(s)
+    s = radis.test_spectrum()
+    s.has_nan()
 
 
 How to export ?
diff --git a/examples/plot_multi_temperature_fit.py b/examples/plot_multi_temperature_fit.py
index dcdba2b91..76e828e2d 100644
--- a/examples/plot_multi_temperature_fit.py
+++ b/examples/plot_multi_temperature_fit.py
@@ -39,7 +39,7 @@
         join("test_CO2_3Tvib_vs_klarenaar_data", "klarenaar_2017_digitized_data.csv")
     ),
     "transmittance_noslit",
-    waveunit="cm-1",
+    wunit="cm-1",
     unit="",
     delimiter=",",
     name="Klarenaar 2017",
diff --git a/examples/plot_slit_function.py b/examples/plot_slit_function.py
index 231478782..b55846d0d 100644
--- a/examples/plot_slit_function.py
+++ b/examples/plot_slit_function.py
@@ -13,4 +13,4 @@
 from radis.test.utils import getTestFile
 
 my_slit = getTestFile("slitfunction.txt")  # for the example here
-plot_slit(my_slit, waveunit="nm")
+plot_slit(my_slit, wunit="nm")
diff --git a/radis/lbl/bands.py b/radis/lbl/bands.py
index f583d3fb0..66c7826dc 100644
--- a/radis/lbl/bands.py
+++ b/radis/lbl/bands.py
@@ -368,10 +368,10 @@ def add_attr(attr):
                 lines=lines,
                 units=self.units,
                 cond_units=self.cond_units,
-                waveunit=self.params.waveunit,  # cm-1
+                wunit=self.params.waveunit,  # cm-1
                 name=band,
                 # dont check input (much faster, and Spectrum
-                warnings=False,
+                check_wavespace=False,
                 # is freshly baken so probably in a good format
             )
 
@@ -725,10 +725,10 @@ def add_attr(attr):
                 lines=lines,
                 units=self.units,
                 cond_units=self.cond_units,
-                waveunit=self.params.waveunit,  # cm-1
+                wunit=self.params.waveunit,  # cm-1
                 name=band,
                 # dont check input (much faster, and Spectrum
-                warnings=False,
+                check_wavespace=False,
                 # is freshly baken so probably in a good format
             )
 
diff --git a/radis/lbl/base.py b/radis/lbl/base.py
index 0c12922a1..2c0c8b9cc 100644
--- a/radis/lbl/base.py
+++ b/radis/lbl/base.py
@@ -79,7 +79,7 @@
     from .loader import KNOWN_LVLFORMAT, DatabankLoader, df_metadata
 except ImportError:  # if ran from here
     from radis.lbl.loader import KNOWN_LVLFORMAT, DatabankLoader, df_metadata
-
+from radis.misc.arrays import anynan
 from radis.misc.basics import all_in, transfer_metadata
 from radis.misc.debug import printdbg
 from radis.misc.log import printwarn
@@ -291,7 +291,7 @@ def assert_no_nan(self, df, column):
         from radis.misc.printer import get_print_full
 
         try:
-            assert np.isnan(df[column]).sum() == 0
+            assert not anynan(df[column])
         except AssertionError as err:
             index = np.isnan(df[column]).idxmax()
             if self.input.molecule == "CO2":
diff --git a/radis/lbl/broadening.py b/radis/lbl/broadening.py
index 56ca2b290..ed3fb38f7 100644
--- a/radis/lbl/broadening.py
+++ b/radis/lbl/broadening.py
@@ -587,7 +587,7 @@ def voigt_lineshape(w_centered, hwhm_lorentz, hwhm_voigt, jit=True):
     # Normalize
     lineshape /= integral
 
-    assert not np.isnan(lineshape).any()
+    # assert not anynan(lineshape).any()
 
     return lineshape
 
diff --git a/radis/lbl/factory.py b/radis/lbl/factory.py
index bac774990..7257e781a 100644
--- a/radis/lbl/factory.py
+++ b/radis/lbl/factory.py
@@ -827,9 +827,9 @@ def eq_spectrum(
             lines=lines,
             units=self.units,
             cond_units=self.cond_units,
-            waveunit=self.params.waveunit,  # cm-1
+            wunit=self.params.waveunit,  # cm-1
             # dont check input (much faster, and Spectrum
-            warnings=False,
+            check_wavespace=False,
             # is freshly baken so probably in a good format
             name=name,
             references=dict(self.reftracker),
@@ -1127,9 +1127,9 @@ def eq_spectrum_gpu(
             conditions=conditions,
             lines=lines,
             cond_units=self.cond_units,
-            waveunit=self.params.waveunit,  # cm-1
+            wunit=self.params.waveunit,  # cm-1
             # dont check input (much faster, and Spectrum
-            warnings=False,
+            check_wavespace=False,
             # is freshly baken so probably in a good format
             name=name,
             references=dict(self.reftracker),
@@ -1482,9 +1482,9 @@ def non_eq_spectrum(
             lines=lines,
             units=self.units,
             cond_units=self.cond_units,
-            waveunit=self.params.waveunit,  # cm-1
+            wunit=self.params.waveunit,  # cm-1
             # dont check input (much faster, and Spectrum
-            warnings=False,
+            check_wavespace=False,
             # is freshly baken so probably in a good format
             name=name,
             references=dict(self.reftracker),
diff --git a/radis/lbl/overp.py b/radis/lbl/overp.py
index e4777a37c..b3a5a451f 100644
--- a/radis/lbl/overp.py
+++ b/radis/lbl/overp.py
@@ -38,8 +38,7 @@
     rescale_radiance_noslit,
     rescale_transmittance_noslit,
 )
-from radis.spectrum.spectrum import is_spectrum
-from radis.spectrum.utils import NON_CONVOLUTED_QUANTITIES
+from radis.spectrum.spectrum import Spectrum
 
 # keys that may not be equals for different bands
 _IGNORE_KEYS = ["band", "band_htrn", "viblvl_u", "viblvl_l"]
@@ -93,7 +92,7 @@ def __init__(
 
         # Check inputs
         for br, s in bands.items():
-            if not is_spectrum(s):
+            if not isinstance(s, Spectrum):
                 raise ValueError(
                     "`bands` must be a list of Spectrum objects. "
                     + "Got {0}".format(type(s))
@@ -845,10 +844,7 @@ def rescale_updown_levels(
     # Save (only) the ones that were in the spectrum initially
     for q in rescaled:
         if q in initial:
-            if q in NON_CONVOLUTED_QUANTITIES:
-                spec._q[q] = rescaled[q]
-            else:
-                spec._q_conv[q] = rescaled[q]
+            spec._q[q] = rescaled[q]
 
     # Update units
     for k, u in units.items():
@@ -856,8 +852,8 @@ def rescale_updown_levels(
 
     # Drop convoluted values
     for q in ["transmittance", "radiance"]:
-        if q in list(spec._q_conv.keys()):
-            del spec._q_conv[q]
+        if q in list(spec._q.keys()):
+            del spec._q[q]
 
     # Reapply slit if possible
     if (
diff --git a/radis/los/slabs.py b/radis/los/slabs.py
index 0c940ac01..69f1f901c 100644
--- a/radis/los/slabs.py
+++ b/radis/los/slabs.py
@@ -24,9 +24,10 @@
 import numpy as np
 from numpy import abs, allclose, arange, diff
 
+from radis.misc.arrays import anynan
 from radis.misc.basics import in_all, merge_lists
 from radis.misc.debug import printdbg
-from radis.spectrum.spectrum import Spectrum, is_spectrum
+from radis.spectrum.spectrum import Spectrum
 
 # %% Slabs / Multi-layers / Radiative Transfer Equation (RTE)
 # ----------------------------------------------------------------------
@@ -159,7 +160,7 @@ def SerialSlabs(*slabs, **kwargs) -> Spectrum:
         raise ValueError("Empty list of slabs")
 
     elif len(slabs) == 1:
-        if not is_spectrum(slabs[0]):
+        if not isinstance(slabs[0], Spectrum):
             raise TypeError(
                 "SerialSlabs takes an unfolded list of Spectrum as "
                 + "argument: *list (got {0})".format(type(slabs[0]))
@@ -331,16 +332,14 @@ def _check_valid(s):
     - quantities used for solving the LOS have nan
     """
 
-    if not is_spectrum(s):
+    if not isinstance(s, Spectrum):
         raise TypeError(
             "All inputs must be Spectrum objects (got: {0})".format(type(s))
         )
-    sdict = s._get_items()
-    for k in sdict.keys():
+    for k, v in s._q.items():
         if (
             k in ["transmittance_noslit", "radiance_noslit", "abscoeff", "emisscoeff"]
-            and np.isnan(sdict.get(k)[1]).any()
-        ):
+        ) and anynan(v):
             warn(
                 "Nans detected in Spectrum object for multi-slab operation. "
                 + "Results may be wrong!"
@@ -454,7 +453,6 @@ def same_wavespace(wl):
                     s.resample(
                         wnew,
                         unit=waveunit,
-                        if_conflict_drop="convoluted",
                         out_of_bounds=out_of_bounds,
                         inplace=True,  # we already copied 'if not modify_inputs'
                     )
@@ -473,7 +471,6 @@ def same_wavespace(wl):
                     s.resample(
                         wnew,
                         unit=waveunit,
-                        if_conflict_drop="convoluted",
                         out_of_bounds=out_of_bounds,
                         inplace=True,  # we already copied 'if not modify_inputs'
                     )
@@ -624,7 +621,7 @@ def MergeSlabs(*slabs, **kwargs) -> Spectrum:
         raise ValueError("Empty list of slabs")
 
     elif len(slabs) == 1:
-        if not is_spectrum(slabs[0]):
+        if not isinstance(slabs[0], Spectrum):
             raise TypeError(
                 "MergeSlabs takes an unfolded list of Spectrum as "
                 + "argument: (got {0})".format(type(slabs[0]))
diff --git a/radis/misc/arrays.py b/radis/misc/arrays.py
index 5adb3d1c2..d5af9dd38 100644
--- a/radis/misc/arrays.py
+++ b/radis/misc/arrays.py
@@ -106,7 +106,7 @@ def array_allclose(a, b, rtol=1e-5, atol=1e-8, equal_nan=True):
 
 
 def nantrapz(I, w, dx=1.0, axis=-1):
-    """Returns np.trapz(I, w) discarding nan."""
+    """Returns :py:func:`~numpy.trapz` (I, w) discarding nan."""
     b = ~np.isnan(I)
     return np.trapz(I[b], w[b], dx=dx, axis=axis)
 
@@ -275,7 +275,50 @@ def evenly_distributed(w, tolerance=1e-5):
     return (np.abs((np.diff(w) - mean_step)) > tolerance).sum() == 0
 
 
-@numba.jit
+def anynan(a):
+    """Returns whether ``a`` has at least one :py:attr:`~numpy.nan`
+
+    Fastest implementation for arrays with >10^4 elements
+    https://stackoverflow.com/a/45011547/5622825
+    """
+    return np.isnan(np.dot(a, a))
+
+
+@numba.njit
+def first_nonnan_index(a):
+    """Returns index of first non-nan value in ``a``
+
+    Returns None is all values are :py:attr:`~numpy.nan`
+
+    See Also
+    --------
+
+    :func:`~radis.misc.arrays.last_nonnan_index`
+    """
+    for i in range(a.size):
+        if not np.isnan(a[i]):
+            return i
+    return None
+
+
+@numba.njit
+def last_nonnan_index(a):
+    """Returns index of first non-nan value in ``a``
+
+    Returns None is all values are :py:attr:`~numpy.nan`
+
+    See Also
+    --------
+
+    :func:`~radis.misc.arrays.first_nonnan_index`
+    """
+    for i in range(a.size - 1, 0, -1):
+        if not np.isnan(a[i]):
+            return i
+    return None
+
+
+@numba.njit
 def is_sorted(a):
     """Returns whether ``a`` is sorted in ascending order.
 
@@ -292,7 +335,7 @@ def is_sorted(a):
     return True
 
 
-@numba.jit
+@numba.njit
 def is_sorted_backward(a):
     """Returns whether ``a`` is sorted in descending order.
 
diff --git a/radis/misc/curve.py b/radis/misc/curve.py
index e86d364b4..59297173e 100644
--- a/radis/misc/curve.py
+++ b/radis/misc/curve.py
@@ -58,13 +58,10 @@ def curve_distance(w1, I1, w2, I2, discard_out_of_bounds=True):
 
     Parameters
     ----------
-
     w1, I1: array
         range and values for first curve
-
     w2, I2: array
         range and values for 2nd curve
-
     discard_out_of_bounds: boolean
         if True, distance for out of bound values is set to nan. Else, it
         will be the distance from the last point.
@@ -101,19 +98,19 @@ def curve_add(w1, I1, w2, I2, is_sorted=False, kind="linear"):
 
     Similar to OriginPro's "Simple Curve Math Substract"
 
+    .. note::
+        for higher accuracy, choose (w2, I2) has the curve with the highest
+        resolution
+
 
     Parameters
     ----------
-
     w1, I1: array
         range and values for first curve
-
     w2, I2: array
         range and values for 2nd curve
-
     is_sorted: boolean
         (optional) if True, doesnt sort input arrays
-
     kind: str
         interpolation kind. Default 'linear'. See scipy.interpolate.interp1d
 
@@ -136,19 +133,19 @@ def curve_substract(w1, I1, w2, I2, is_sorted=False, kind="linear"):
 
     Similar to OriginPro's "Simple Curve Math Substract"
 
+    .. note::
+        for higher accuracy, choose (w2, I2) has the curve with the highest
+        resolution
+
 
     Parameters
     ----------
-
     w1, I1: array
         range and values for first curve
-
     w2, I2: array
         range and values for 2nd curve
-
     is_sorted: boolean
         (optional) if True, doesnt sort input arrays
-
     kind: str
         interpolation kind. Default 'linear'. See scipy.interpolate.interp1d
 
@@ -171,26 +168,25 @@ def curve_multiply(w1, I1, w2, I2, is_sorted=False, kind="linear"):
 
     Similar to OriginPro's "Simple Curve Math Substract"
 
+    .. note::
+        for higher accuracy, choose (w2, I2) has the curve with the highest
+        resolution
+
 
     Parameters
     ----------
-
     w1, I1: array
         range and values for first curve
-
     w2, I2: array
         range and values for 2nd curve
-
     is_sorted: boolean
         (optional) if True, doesnt sort input arrays
-
     kind: str
         interpolation kind. Default 'linear'. See scipy.interpolate.interp1d
 
 
     Returns
     -------
-
     w1, Iproduct: array
         product ``I1 * I2`` interpolated on the first range ``w1``
     """
@@ -206,33 +202,31 @@ def curve_divide(w1, I1, w2, I2, is_sorted=False, kind="linear", interpolation=1
 
     Similar to OriginPro's "Simple Curve Math Substract"
 
+    .. note::
+        for higher accuracy, choose (w2, I2) has the curve with the highest
+        resolution
+
 
     Parameters
     ----------
-
     w1, I1: array
         range and values for first curve
-
     w2, I2: array
         range and values for 2nd curve
 
     Other Parameters
     ----------------
-
     is_sorted: boolean
         (optional) if ``True``, assumes that both input arrays are sorted already.
         Default ``False``.
-
     kind: str
         interpolation kind. Default 'linear'. See scipy.interpolate.interp1d
-
     interpolation: int, optional
         If 1, interpolate on w1, I1. Else, on w2, I2.
         Default 1
 
     Returns
     -------
-
     w1, Idiv: array
         Division ``I1 / I2`` interpolated on the first or second range according to reverseInterpolation
     """
@@ -257,15 +251,24 @@ def _curve_interpolate(w1, I1, w2, I2, is_sorted=False, kind="linear"):
     For the interpolation to work, this requires to first sort ``w1`` and ``w2``,
     interpolate, then revert to the initial ``w1``
 
+    .. note::
+        for higher accuracy, choose (w2, I2) has the curve with the highest
+        resolution
+
     Returns
     -------
 
     I2_inter: np.array
         ``I2`` interpolated on ``w1``.
     """
+    if (abs((w1[1] - w1[0])) < 2 * abs((w2[1] - w2[0]))) and (
+        abs((w1[-1] - w1[-2])) < 2 * abs((w2[-1] - w2[-2]))
+    ):
+        # w1 looks like the spectra with highest accuracy
+        warnings.warn("Reverse your spectra in interpolation for a better accuracy")
 
     # First sort both
-    if not is_sorted:
+    if not is_sorted:  # TODO : add first check with `radis.misc.is_sorted(a) `
         b = np.argsort(w1)
         w1, I1 = w1[b], I1[b]
         # Get reverse sorting array:
diff --git a/radis/misc/signal.py b/radis/misc/signal.py
index 04412f9ae..8f5ade720 100644
--- a/radis/misc/signal.py
+++ b/radis/misc/signal.py
@@ -13,12 +13,16 @@
 """
 
 
-from warnings import warn
-
-from numpy import abs, isnan, linspace, nan, trapz
+from numpy import abs, isnan, linspace, nan, trapz, zeros_like
 from scipy.interpolate import splev, splrep
 
-from radis.misc.arrays import is_sorted, is_sorted_backward
+from radis.misc.arrays import (
+    anynan,
+    first_nonnan_index,
+    is_sorted,
+    is_sorted_backward,
+    last_nonnan_index,
+)
 from radis.misc.debug import printdbg
 
 
@@ -28,7 +32,7 @@ def resample(
     xspace_new,
     k=1,
     ext="error",
-    energy_threshold=1e-3,
+    energy_threshold=5e-3,
     print_conservation=True,
 ):
     """Resample (xspace, vector) on a new space (xspace_new) of evenly
@@ -62,7 +66,7 @@ def resample(
     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 1e-3 (0.1%)
+        Default 5e-3 (0.5%)
     print_conservation: boolean
         if True, prints energy conservation
 
@@ -73,7 +77,6 @@ def resample(
 
     Notes
     -----
-
     Note that depending upon the from_space > to_space operation, sorting may
     be reversed.
 
@@ -119,9 +122,22 @@ def resample(
     else:
         raise ValueError("Unexpected value for `ext`: {0}".format(ext))
 
-    if isnan(vector).sum() > 0:
+    # Handle nans in vector :
+    # FITPACK will fail if there are nans in the vector being interpolated.
+    # here we trim the sides; and raise an error if there are nans in the middle
+    m = first_nonnan_index(vector)
+    M = last_nonnan_index(vector)
+    if m is None:
+        raise ValueError("Vector being resample has only nans. Check your data")
+    nanmask = False
+    if m > 0 or M < len(vector) - 1:  # else, no change has to be made
+        nanmask = zeros_like(vector, dtype=bool)
+        nanmask[m : M + 1] = True
+        xspace = xspace[nanmask]
+        vector = vector[nanmask]
+    if anynan(vector):
         raise ValueError(
-            "Resampled vector has {0} nans. Interpolation will fail".format(
+            "Vector being resample has {0} nans. Interpolation will fail".format(
                 isnan(vector).sum()
             )
         )
@@ -295,65 +311,11 @@ def resample_even(
     return xspace_new, vector_new
 
 
-# %% Test routines
-
-
-def _test(verbose=True, debug=False, plot=True, warnings=True, *args, **kwargs):
-    """Test procedures.
-
-    Parameters
-    ----------
-
-    debug: boolean
-        swamps the console namespace with local variables. Default ``False``
-    """
-
-    import matplotlib.pyplot as plt
-    from numpy import linspace, loadtxt
-
-    from radis.phys.convert import cm2nm, nm2cm
-    from radis.test.utils import getTestFile
-
-    # Test even resampling
-
-    w_nm, I_nm = loadtxt(getTestFile("spectrum.txt")).T
-    w_cm, I_cm = resample_even(
-        nm2cm(w_nm),
-        I_nm,
-        resfactor=2,
-        energy_threshold=1e-3,
-        print_conservation=verbose,
+if __name__ == "__main__":
+    from radis.test.spectrum.test_spectrum import (
+        test_resampling_function,
+        test_resampling_nan_function,
     )
 
-    if plot:
-        plt.figure()
-        plt.xlabel("Wavelength (nm)")
-        plt.ylabel("Intensity")
-        plt.plot(w_nm, I_nm, "-ok", label="original")
-        plt.plot(cm2nm(w_cm), I_cm, "-or", label="resampled")
-        plt.legend()
-
-    # Test resampling
-
-    w_crop = linspace(376, 381, 100)
-    I_crop = resample(w_nm, I_nm, w_crop, energy_threshold=0.01)
-
-    if plot:
-        plt.figure()
-        plt.xlabel("Wavelength (nm)")
-        plt.ylabel("Intensity")
-        plt.plot(w_nm, I_nm, "-ok", label="original")
-        plt.plot(w_crop, I_crop, "-or", label="resampled")
-        plt.legend()
-
-    if debug:
-        globals().update(locals())
-
-    if warnings:
-        warn("Testing resampling: no quantitative tests defined yet")
-
-    return True  # no standard tests yet
-
-
-if __name__ == "__main__":
-    _test(debug=False)
+    test_resampling_function()
+    test_resampling_nan_function()
diff --git a/radis/phys/blackbody.py b/radis/phys/blackbody.py
index 69e847d8b..7fb025f6f 100644
--- a/radis/phys/blackbody.py
+++ b/radis/phys/blackbody.py
@@ -241,7 +241,7 @@ def sPlanck(
         conditions=conditions,
         units={"radiance_noslit": Iunit, "transmittance_noslit": "", "absorbance": ""},
         cond_units={"wstep": waveunit},
-        waveunit=waveunit,
+        wunit=waveunit,
         name="Planck {0}K, eps={1:.2g}".format(T, eps),
     )
 
diff --git a/radis/spectrum/__init__.py b/radis/spectrum/__init__.py
index 11aaf6098..0e63fbc23 100644
--- a/radis/spectrum/__init__.py
+++ b/radis/spectrum/__init__.py
@@ -1,6 +1,7 @@
 # -*- coding: utf-8 -*-
 """Spectrum-class module for post-processing."""
 
+from . import rescale, spectrum, utils
 from .compare import (
     get_diff,
     get_distance,
@@ -19,7 +20,7 @@
     get_baseline,
     sub_baseline,
 )
-from .spectrum import Spectrum, is_spectrum
+from .spectrum import Spectrum
 
 __all__ = [
     "get_diff",
diff --git a/radis/spectrum/compare.py b/radis/spectrum/compare.py
index 5b82c6a46..a8cffdcb2 100644
--- a/radis/spectrum/compare.py
+++ b/radis/spectrum/compare.py
@@ -31,7 +31,7 @@
 from radis.misc.arrays import array_allclose
 from radis.misc.basics import compare_dict, compare_lists
 from radis.misc.curve import curve_distance, curve_divide, curve_substract
-from radis.spectrum.spectrum import Spectrum, is_spectrum
+from radis.spectrum.spectrum import Spectrum
 from radis.spectrum.utils import cast_waveunit, format_xlabel, make_up, make_up_unit
 
 # %% ======================================================================
@@ -41,7 +41,13 @@
 
 
 def get_diff(
-    s1, s2, var, wunit="default", Iunit="default", resample=True, diff_window=0
+    s1: Spectrum,
+    s2: Spectrum,
+    var,
+    wunit="default",
+    Iunit="default",
+    resample=True,
+    diff_window=0,
 ):
     # type: (Spectrum, Spectrum, str, str, str, str, bool, int) -> np.array, np.array
     """Get the difference between 2 spectra.
@@ -97,13 +103,14 @@ def get_diff(
     :func:`~radis.spectrum.compare.plot_diff`,
     :meth:`~radis.spectrum.spectrum.compare_with`
     """
+    # TODO : ensure objects are sorted and if so; forward to curve_substract()
 
     var, wunit, Iunit = get_default_units(s1, s2, var=var, wunit=wunit, Iunit=Iunit)
 
     # Get data
     # ----
-    w1, I1 = s1.get(var, wunit=wunit, Iunit=Iunit)
-    w2, I2 = s2.get(var, wunit=wunit, Iunit=Iunit)
+    w1, I1 = s1.get(var, wunit=wunit, Iunit=Iunit, trim_nan=True)
+    w2, I2 = s2.get(var, wunit=wunit, Iunit=Iunit, trim_nan=True)
 
     if not resample:
         if not array_allclose(w1, w2):
@@ -125,7 +132,9 @@ def get_diff(
     return w1, Idiff
 
 
-def get_ratio(s1, s2, var, wunit="default", Iunit="default", resample=True):
+def get_ratio(
+    s1: Spectrum, s2: Spectrum, var, wunit="default", Iunit="default", resample=True
+):
     # type: (Spectrum, Spectrum, str, str, str, str, bool) -> np.array, np.array
     """Get the ratio between two spectra Basically returns w1, I1 / I2 where
     (w1, I1) and (w2, I2) are the values of s1 and s2 for variable var. (w2,
@@ -177,7 +186,9 @@ def get_ratio(s1, s2, var, wunit="default", Iunit="default", resample=True):
     return curve_divide(w1, I1, w2, I2)
 
 
-def get_distance(s1, s2, var, wunit="default", Iunit="default", resample=True):
+def get_distance(
+    s1: Spectrum, s2: Spectrum, var, wunit="default", Iunit="default", resample=True
+):
     # type: (Spectrum, Spectrum, str, str, str, str, bool) -> np.array, np.array
     """Get a regularized Euclidian distance between two spectra ``s1`` and
     ``s2``
@@ -250,8 +261,8 @@ def get_distance(s1, s2, var, wunit="default", Iunit="default", resample=True):
 
 
 def get_residual(
-    s1,
-    s2,
+    s1: Spectrum,
+    s2: Spectrum,
     var,
     norm="L2",
     ignore_nan=False,
@@ -291,9 +302,9 @@ def get_residual(
     Other Parameters
     ----------------
     ignore_nan: boolean
-        if ``True``, ignore nan in the difference between s1 and s2 (ex: out of bound)
+        if ``True``, ignore ``nan`` in the difference between s1 and s2 (ex: out of bound)
         when calculating residual. Default ``False``. Note: ``get_residual`` will still
-        fail if there are nan in initial Spectrum.
+        fail if there are ``nan`` in initial Spectrum.
     normalize: bool, or tuple
         if ``True``, normalize the two spectra before calculating the residual.
         If a tuple (ex: ``(4168, 4180)``), normalize on this range only. The unit
@@ -363,24 +374,24 @@ def get_residual(
     if ignore_nan:
         b = np.isnan(dI)
         wdiff, dI = wdiff[~b], dI[~b]
-    warningText = (
-        'NaN output in residual. You should use "ignore_nan=True". Read the help.'
-    )
+
     if norm == "L2":
         output = np.sqrt((dI ** 2).sum()) / len(dI)
-        if np.isnan(output):
-            warn(warningText, UserWarning)
-        return output
     elif norm == "L1":
         output = (np.abs(dI)).sum() / len(dI)
-        if np.isnan(output):
-            warn.warning(warningText, UserWarning)
-        return output
     else:
         raise ValueError("unexpected value for norm")
 
+    if np.isnan(output):
+        warn(
+            'NaN output in residual. You should use "ignore_nan=True". Read the help of get_residual.',
+            UserWarning,
+        )
+
+    return output
+
 
-def get_residual_integral(s1, s2, var, ignore_nan=False) -> float:
+def get_residual_integral(s1: Spectrum, s2: Spectrum, var, ignore_nan=False) -> 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.
@@ -457,7 +468,9 @@ def get_residual_integral(s1, s2, var, ignore_nan=False) -> float:
     return np.abs(np.trapz(dI, wdiff)) / norm
 
 
-def get_default_units(s1, s2, var=None, wunit="default", Iunit="default"):
+def get_default_units(
+    s1: Spectrum, s2: Spectrum, var=None, wunit="default", Iunit="default"
+):
     """Get ``wunit``, Iunit for var; compatible with both spectra s1 and s2
 
     Parameters
@@ -527,8 +540,8 @@ def get_default_units(s1, s2, var=None, wunit="default", Iunit="default"):
 
 
 def _get_wdiff_Idiff(
-    s1,
-    s2,
+    s1: Spectrum,
+    s2: Spectrum,
     var=None,
     wunit="default",
     Iunit="default",
@@ -632,8 +645,8 @@ def _get_wdiff_Idiff(
 
 
 def plot_diff(
-    s1,
-    s2,
+    s1: Spectrum,
+    s2: Spectrum,
     var=None,
     wunit="default",
     Iunit="default",
@@ -1013,7 +1026,7 @@ def plot_diff(
     return fig, axes
 
 
-def averageDistance(s1, s2, var="radiance"):
+def averageDistance(s1: Spectrum, s2: Spectrum, var="radiance"):
     """Return the average distance between two spectra.
     It's important to note for fitting that if averageDistance(s1, s2)==0
     then s1 = s2
@@ -1056,8 +1069,8 @@ def averageDistance(s1, s2, var="radiance"):
 
 
 def compare_spectra(
-    first,
-    other,
+    first: Spectrum,
+    other: Spectrum,
     spectra_only=False,
     plot=True,
     wunit="default",
@@ -1068,7 +1081,7 @@ def compare_spectra(
     ignore_conditions=["calculation_time"],
     normalize=False,
     **kwargs
-):
+) -> bool:
     """Compare Spectrum with another Spectrum object.
 
     Parameters
@@ -1089,7 +1102,8 @@ def compare_spectra(
         in which wavespace to compare (and plot). If default, natural wavespace
         of first Spectrum is taken
     rtol: float
-        relative difference to use for spectral quantities comparison
+        relative difference to use for spectral quantities comparison.
+        absolute difference ``atol`` is set to 0.
     ignore_nan: boolean
         if ``True``, nans are ignored when comparing spectral quantities
     ignore_outliers: boolean, or float
@@ -1130,7 +1144,7 @@ def compare_spectra(
     # Check inputs
     if not 0 <= ignore_outliers < 1:
         raise ValueError("ignore_outliers should be < 1, or False")
-    if not is_spectrum(other):
+    if not isinstance(other, Spectrum):
         raise TypeError(
             "2nd object is not a Spectrum: got class {0}".format(other.__class__)
         )
@@ -1174,21 +1188,6 @@ def _compare_dataframes(df1, df2, name):
             for error message
         """
 
-        #            if compare_lists(df1.keys(), df2.keys(), verbose=False) != 1:
-        #                if verbose: print('... keys in {0} dont match:'.format(name))
-        #                compare_lists(list(df1.keys()), list(df2.keys()),
-        #                              verbose=True)
-        #                out = False
-        #            elif compare_lists(df1.index, df2.index, verbose=False) != 1:
-        #                if verbose: print('... index in {0} dont match:'.format(name))
-        #                compare_lists(list(df1.index), list(df2.index),
-        #                              verbose=True)
-        #                out = False
-        #            else:
-        #                out = (df1 == df2).all().all()
-        #
-        #            return out
-
         from pandas.util.testing import assert_frame_equal
 
         try:
@@ -1208,7 +1207,7 @@ def _compare_dataframes(df1, df2, name):
 
         return out
 
-    def _compare_variables(I, Ie):
+    def _compare_arrays(I, Ie):
         """Compare spectral quantities I and Ie."""
 
         if ignore_nan:
@@ -1216,12 +1215,16 @@ def _compare_variables(I, Ie):
             I = I[b]
             Ie = Ie[b]
 
+        # in all cases, we ignore if there are nans at the same
+        # on both arrays (happens with convolved values in 'valid' mode)
+        equal_nan = True
+
         if ignore_outliers:
-            out = (~np.isclose(I, Ie, rtol=rtol, atol=0)).sum() / len(
-                I
-            ) < ignore_outliers
+            out = (
+                ~np.isclose(I, Ie, rtol=rtol, atol=0, equal_nan=equal_nan)
+            ).sum() / len(I) < ignore_outliers
         else:
-            out = np.allclose(I, Ie, rtol=rtol, atol=0)
+            out = np.allclose(I, Ie, rtol=rtol, atol=0, equal_nan=equal_nan)
 
         return bool(out)
 
@@ -1249,7 +1252,7 @@ def _display_difference(q, q0):
         vars = [k for k in first.get_vars() if k in other.get_vars()]
 
     if spectra_only:
-        # Compare spectral variables
+        # Compare spectral arrays
         # -----------
         for k in vars:
             w, q = first.get(k, wunit=wunit)
@@ -1267,19 +1270,19 @@ def _display_difference(q, q0):
                 if len(q) > len(q0):
                     f = interp1d(w, q, kind="cubic")
                     new_q = f(w0)
-                    b1 = _compare_variables(new_q, q0)
+                    b1 = _compare_arrays(new_q, q0)
                     if not b1 and verbose:
                         _display_difference(new_q, q0)
                 else:
                     f = interp1d(w0, q0, kind="cubic")
                     new_q0 = f(w)
-                    b1 = _compare_variables(q, new_q0)
+                    b1 = _compare_arrays(q, new_q0)
                     if not b1 and verbose:
                         _display_difference(q, new_q0)
 
             else:  # no need to interpolate
                 b1 = np.allclose(w, w0, rtol=rtol, atol=0)
-                b1 *= _compare_variables(q, q0)
+                b1 *= _compare_arrays(q, q0)
 
                 if not b1 and verbose:
                     _display_difference(q, q0)
@@ -1314,7 +1317,7 @@ def _display_difference(q, q0):
                 b1 = False
             else:
                 b1 = np.allclose(w, w0, rtol=rtol, atol=0)
-                b1 *= _compare_variables(q, q0)
+                b1 *= _compare_arrays(q, q0)
                 if not b1 and verbose:
                     error = np.nanmax(abs(q / q0 - 1))
                     avgerr = np.nanmean(abs(q / q0 - 1))
@@ -1476,7 +1479,7 @@ def _display_difference(q, q0):
                 b6 = False
             else:
                 b6 = np.allclose(ws, ws0, rtol=rtol, atol=0)
-                b6 *= _compare_variables(Is, Is0)
+                b6 *= _compare_arrays(Is, Is0)
                 if not b6 and verbose:
                     print("Slit functions dont match")
         b *= b6
diff --git a/radis/spectrum/models.py b/radis/spectrum/models.py
index bc0dd6889..e7107c490 100644
--- a/radis/spectrum/models.py
+++ b/radis/spectrum/models.py
@@ -110,7 +110,7 @@ def calculated_spectrum(
         np.array(w),
         np.array(I),
         "radiance_noslit",
-        waveunit=wunit,
+        wunit=wunit,
         unit=Iunit,
         conditions=conditions,
         cond_units=cond_units,
@@ -181,7 +181,7 @@ def transmittance_spectrum(
         np.array(w),
         np.array(T),
         "transmittance_noslit",
-        waveunit=wunit,
+        wunit=wunit,
         unit=Tunit,
         conditions=conditions,
         cond_units=cond_units,
@@ -258,7 +258,7 @@ def experimental_spectrum(
         np.array(w),
         np.array(I),
         "radiance",
-        waveunit=wunit,
+        wunit=wunit,
         unit=Iunit,
         conditions=conditions,
         cond_units=cond_units,
diff --git a/radis/spectrum/operations.py b/radis/spectrum/operations.py
index 3f392126c..516bb7b0f 100644
--- a/radis/spectrum/operations.py
+++ b/radis/spectrum/operations.py
@@ -77,8 +77,7 @@
 # %% Filter Spectra
 
 
-def Transmittance(s):
-    # type: (Spectrum) -> Spectrum
+def Transmittance(s: Spectrum) -> Spectrum:
     """Returns a new Spectrum with only the ``transmittance`` component of ``s``
 
     Parameters
@@ -118,7 +117,7 @@ def Transmittance(s):
     return s.copy(copy_lines=True, quantity="transmittance")
 
 
-def Transmittance_noslit(s):
+def Transmittance_noslit(s: Spectrum) -> Spectrum:
     """Returns a new Spectrum with only the ``transmittance_noslit`` component of ``s``
 
     Parameters
@@ -157,7 +156,7 @@ def Transmittance_noslit(s):
     return s.copy(copy_lines=True, quantity="transmittance_noslit")
 
 
-def Radiance(s):
+def Radiance(s: Spectrum) -> Spectrum:
     """Returns a new Spectrum with only the ``radiance`` component of ``s``
 
     Parameters
@@ -196,7 +195,7 @@ def Radiance(s):
     return s.copy(copy_lines=True, quantity="radiance")
 
 
-def Radiance_noslit(s):
+def Radiance_noslit(s: Spectrum) -> Spectrum:
     """Returns a new Spectrum with only the ``radiance_noslit`` component of ``s``
 
     Parameters
@@ -238,7 +237,7 @@ def Radiance_noslit(s):
 # Useful for determining line-of-sight contributions:
 
 
-def PerfectAbsorber(s):
+def PerfectAbsorber(s: Spectrum) -> Spectrum:
     """Makes a new Spectrum with the same transmittance/absorbance as Spectrum
     ``s``, but with radiance set to 0.
     Useful to get contribution of different slabs in line-of-sight
@@ -289,8 +288,8 @@ def PerfectAbsorber(s):
             s_tr._q[k] *= 0
 
     for k in ["radiance", "emissivity"]:
-        if k in s_tr._q_conv:
-            s_tr._q_conv[k] *= 0
+        if k in s_tr._q:
+            s_tr._q[k] *= 0
 
     # Deactivate equilibrium conditions (so Kirchoff's law cannot be used anymore)
     s_tr.conditions["thermal_equilibrium"] = False
@@ -303,7 +302,7 @@ def PerfectAbsorber(s):
 # %% Change wavelength
 
 
-def crop(s, wmin=None, wmax=None, wunit=None, inplace=False):
+def crop(s: Spectrum, wmin=None, wmax=None, wunit=None, inplace=False) -> Spectrum:
     # type: (Spectrum, float, float, str, str, bool) -> Spectrum
     """Crop spectrum to ``wmin-wmax`` range in ``wunit``
 
@@ -358,17 +357,6 @@ def crop(s, wmin=None, wmax=None, wunit=None, inplace=False):
             "wmin should be < wmax (Got: {0:.2f}, {1:.2f})".format(wmin, wmax)
         )
 
-    if len(s._q) > 0 and len(s._q_conv) > 0:
-        raise NotImplementedError(
-            "Cant crop this Spectrum as there are both convoluted "
-            + "and not convoluted quantities stored"
-        )
-        # Could bring unexpected errors... For instance, if cropping both
-        # with slit and without slit  quantities to the same waverange,
-        # reapplying the slit in 'valid' mode would reduce the wavelength range
-        # of the convoluted quantities
-        # Implementation: better ask User to drop some of the quantities themselves
-
     if not inplace:
         s = s.copy()
 
@@ -428,7 +416,7 @@ def crop(s, wmin=None, wmax=None, wunit=None, inplace=False):
     else:
         raise ValueError(stored_waveunit)
 
-    # Crop non convoluted
+    # Crop
     if len(s._q) > 0:
         b = ones_like(s._q["wavespace"], dtype=bool)
         if wmin:
@@ -438,16 +426,6 @@ def crop(s, wmin=None, wmax=None, wunit=None, inplace=False):
         for k, v in s._q.items():
             s._q[k] = v[b]
 
-    # Crop convoluted
-    if len(s._q_conv) > 0:
-        b = ones_like(s._q_conv["wavespace"], dtype=bool)
-        if wmin:
-            b *= wmin <= s._q_conv["wavespace"]
-        if wmax:
-            b *= s._q_conv["wavespace"] <= wmax
-        for k, v in s._q_conv.items():
-            s._q_conv[k] = v[b]
-
     return s
 
 
@@ -819,7 +797,7 @@ def add_spectra(s1, s2, var=None, force=False):
 
     name = s1.get_name() + "+" + s2.get_name()
 
-    sub = Spectrum.from_array(w1, I1 + I2, var, waveunit=wunit1, unit=Iunit1, name=name)
+    sub = Spectrum.from_array(w1, I1 + I2, var, wunit=wunit1, unit=Iunit1, name=name)
     #    warn("Conditions of the left spectrum were copied in the substraction.", Warning)
     return sub
 
@@ -883,7 +861,7 @@ def substract_spectra(s1, s2, var=None):
 
     name = s1.get_name() + "-" + s2.get_name()
 
-    sub = Spectrum.from_array(w1, I1 - I2, var, waveunit=wunit1, unit=Iunit1, name=name)
+    sub = Spectrum.from_array(w1, I1 - I2, var, wunit=wunit1, unit=Iunit1, name=name)
     #    warn("Conditions of the left spectrum were copied in the substraction.", Warning)
     return sub
 
@@ -978,7 +956,7 @@ def concat_spectra(s1, s2, var=None):
     name = s1.get_name() + "&" + s2.get_name()  # use "&" instead of "+"
 
     concat = Spectrum.from_array(
-        w_tot, I_tot, var, waveunit=wunit1, unit=Iunit1, name=name
+        w_tot, I_tot, var, wunit=wunit1, unit=Iunit1, name=name
     )
 
     return concat
@@ -996,6 +974,9 @@ def offset(s, offset, unit, name=None, inplace=False):
         Constant to add to all quantities in the Spectrum.
     unit: 'nm' or 'cm-1'
         unit for ``offset``.
+
+    Other Parameters
+    ----------------
     name: str
         name of output spectrum
     inplace: bool
@@ -1010,13 +991,9 @@ def offset(s, offset, unit, name=None, inplace=False):
 
     See Also
     --------
-
     call as a Spectrum method directly: :py:meth:`~radis.spectrum.spectrum.Spectrum.offset`
     """
 
-    has_var = len(s._q) > 0
-    has_conv_var = len(s._q_conv) > 0
-
     stored_waveunit = s.get_waveunit()
 
     # Convert offset to correct unit:
@@ -1024,76 +1001,33 @@ def offset(s, offset, unit, name=None, inplace=False):
         if unit == "nm":
             # Note @EP: here we're offsetting by a constant value in 'nm', which is
             # not a constant value in 'cm-1'. The offset is an array
-            if has_var:
-                offset_q = -dnm_air2dcm(
-                    offset, s.get_wavelength(which="non_convoluted")
-                )  # this is an array
-            if has_conv_var:
-                offset_qconv = -dnm_air2dcm(
-                    offset, s.get_wavelength(which="convoluted")
-                )  # this is an array
+            offset_q = -dnm_air2dcm(offset, s.get_wavelength())  # this is an array
         elif unit == "nm_vac":
-            if has_var:
-                offset_q = -dnm2dcm(
-                    offset, s.get_wavelength(which="non_convoluted")
-                )  # this is an array
-            if has_conv_var:
-                offset_qconv = -dnm2dcm(
-                    offset, s.get_wavelength(which="convoluted")
-                )  # this is an array
+            offset_q = -dnm2dcm(offset, s.get_wavelength())  # this is an array
         elif unit == "cm-1":
-            if has_var:
-                offset_q = offset
-            if has_conv_var:
-                offset_qconv = offset
+            offset_q = offset
         else:
             raise ValueError(unit)
     elif stored_waveunit == "nm":  # wavelength air
         if unit == "nm":
-            if has_var:
-                offset_q = offset
-            if has_conv_var:
-                offset_qconv = offset
+            offset_q = offset
         elif unit == "nm_vac":
             # Note @EP: strictly speaking, the offset should change a little bit
             # as nm_vac > nm depends on the wavelength. Neglected here. # TODO ?
-            if has_var:
-                offset_q = offset
-            if has_conv_var:
-                offset_qconv = offset
+            offset_q = offset
         elif unit == "cm-1":
-            if has_var:
-                offset_q = -dcm2dnm_air(
-                    offset, s.get_wavenumber(which="non_convoluted")
-                )  # this is an array
-            if has_conv_var:
-                offset_qconv = -dcm2dnm_air(
-                    offset, s.get_wavenumber(which="convoluted")
-                )  # this is an array
+            offset_q = -dcm2dnm_air(offset, s.get_wavenumber())  # this is an array
         else:
             raise ValueError(unit)
     elif stored_waveunit == "nm_vac":  # wavelength vacuum
         if unit == "nm":
             # Note @EP: strictly speaking, the offset should change a little bit
             # as nm > nm_vac depends on the wavelength. Neglected here. # TODO ?
-            if has_var:
-                offset_q = offset
-            if has_conv_var:
-                offset_qconv = offset
+            offset_q = offset
         elif unit == "nm_vac":
-            if has_var:
-                offset_q = offset
-            if has_conv_var:
-                offset_qconv = offset
+            offset_q = offset
         elif unit == "cm-1":
-            if has_var:
-                offset_q = -dcm2dnm(
-                    offset, s.get_wavenumber(which="non_convoluted")
-                )  # this is an array
-            if has_conv_var:
-                offset_qconv = -dcm2dnm(
-                    offset, s.get_wavenumber(which="convoluted")
-                )  # this is an array
+            offset_q = -dcm2dnm(offset, s.get_wavenumber())  # this is an array
         else:
             raise ValueError(unit)
     else:
@@ -1103,11 +1037,8 @@ def offset(s, offset, unit, name=None, inplace=False):
         s = s.copy()
 
     # Update all variables
-    if has_var:
-        s._q["wavespace"] += offset_q
-        # @dev: updates the Spectrum directly because of copy=False
-    if has_conv_var:
-        s._q_conv["wavespace"] += offset_qconv
+    s._q["wavespace"] += offset_q
+    # @dev: updates the Spectrum directly because of copy=False
 
     if name:
         s.name = name
diff --git a/radis/spectrum/rescale.py b/radis/spectrum/rescale.py
index 1deaaf924..35f5b34a3 100644
--- a/radis/spectrum/rescale.py
+++ b/radis/spectrum/rescale.py
@@ -594,7 +594,7 @@ def rescale_abscoeff(
     old_mole_fraction,
     new_mole_fraction,
     old_path_length,
-    waveunit,
+    wunit,
     units,
     extra,
     true_path_length,
@@ -623,18 +623,18 @@ def rescale_abscoeff(
     # First get initial abscoeff
     # ---------------------
     if "abscoeff" in initial:
-        _, abscoeff_init = spec.get("abscoeff", wunit=waveunit)
+        _, abscoeff_init = spec.get("abscoeff", wunit=wunit)
     elif "absorbance" in initial and true_path_length:  # aka: true path_lengths given
         if __debug__:
             printdbg("... rescale: abscoeff k1 = A/L1")
-        _, A = spec.get("absorbance", wunit=waveunit)
+        _, A = spec.get("absorbance", wunit=wunit)
         abscoeff_init = A / old_path_length  # recalculate initial
         unit = "cm-1"
     elif "transmittance_noslit" in initial and true_path_length:
         if __debug__:
             printdbg("... rescale: abscoeff k1 = -ln(T1)/L1")
         # Get abscoeff from transmittance
-        _, T1 = spec.get("transmittance_noslit", wunit=waveunit)
+        _, T1 = spec.get("transmittance_noslit", wunit=wunit)
 
         # We'll have a problem if the spectrum is optically thick
         b = T1 == 0  # no transmittance: optically thick mask
@@ -748,7 +748,7 @@ def rescale_emisscoeff(
     new_mole_fraction,
     old_path_length,
     optically_thin,
-    waveunit,
+    wunit,
     units,
     extra,
     true_path_length,
@@ -783,15 +783,13 @@ def get_unit(unit_radiance):
         if __debug__:
             printdbg("... rescale: emisscoeff j1 = j1")
         _, emisscoeff_init = spec.get(
-            "emisscoeff", wunit=waveunit, Iunit=units["emisscoeff"]
+            "emisscoeff", wunit=wunit, Iunit=units["emisscoeff"]
         )
 
     elif "radiance_noslit" in initial and true_path_length and optically_thin:
         if __debug__:
             printdbg("... rescale: emisscoeff j1 = I1/L1")
-        _, I = spec.get(
-            "radiance_noslit", wunit=waveunit, Iunit=units["radiance_noslit"]
-        )
+        _, 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"])
 
@@ -799,10 +797,8 @@ def get_unit(unit_radiance):
         if __debug__:
             printdbg("... rescale: emisscoeff j1 = k1*I1/(1-exp(-k1*L1))")
         # get emisscoeff from (initial) abscoeff and (initial) radiance
-        _, I = spec.get(
-            "radiance_noslit", wunit=waveunit, Iunit=units["radiance_noslit"]
-        )
-        _, k = spec.get("abscoeff", wunit=waveunit, Iunit=units["abscoeff"])
+        _, I = spec.get("radiance_noslit", wunit=wunit, Iunit=units["radiance_noslit"])
+        _, k = spec.get("abscoeff", wunit=wunit, Iunit=units["abscoeff"])
 
         # Recalculate in the optically thin range (T=1) and elsewhere
         b = k == 0  # optically thin mask
@@ -827,11 +823,9 @@ def get_unit(unit_radiance):
         if __debug__:
             printdbg("... rescale: emisscoeff j1 = k1*I1/(1-T1)")
         # get emisscoeff from (initial) transmittance and (initial) radiance
-        _, I = spec.get(
-            "radiance_noslit", wunit=waveunit, Iunit=units["radiance_noslit"]
-        )
+        _, I = spec.get("radiance_noslit", wunit=wunit, Iunit=units["radiance_noslit"])
         _, T = spec.get(
-            "transmittance_noslit", wunit=waveunit, Iunit=units["transmittance_noslit"]
+            "transmittance_noslit", wunit=wunit, Iunit=units["transmittance_noslit"]
         )
 
         # Recalculate in the optically thin range (T=1) and elsewhere
@@ -1559,7 +1553,7 @@ def _recalculate(
                     + "not look at equilibrium. Check spectrum conditions"
                 )
             )
-        wavenumber = spec.get_wavenumber("non_convoluted")
+        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
@@ -1669,18 +1663,15 @@ def _recalculate(
     # delete former convoluted value if apply_slit will be used (just to be sure
     # we arent keeping a non rescaled value if something goes wrong)
     if apply_slit:
-        for k in list(spec._q_conv.keys()):
-            if k != "wavespace":
-                del spec._q_conv[k]
+        for k in list(spec._q.keys()):
+            if k != "wavespace" and k in CONVOLUTED_QUANTITIES:
+                del spec._q[k]
                 del spec.units[k]
 
     # Save (only) the ones that we want, unless we want everything ('greedy')
     for q in rescaled:
         if q in wanted:  # or greedy:
-            if q in NON_CONVOLUTED_QUANTITIES:
-                spec._q[q] = rescaled[q]
-            else:
-                spec._q_conv[q] = rescaled[q]
+            spec._q[q] = rescaled[q]
 
     # Update units
     for k, u in units.items():
@@ -1732,7 +1723,9 @@ def _recalculate(
     # ... "everyone is here": check we didnt miss anyone
     rescaled_list = list(rescaled)
     # add the new quantities added by apply_slit
-    rescaled_list = rescaled_list + spec.get_vars("convoluted")
+    rescaled_list = rescaled_list + [
+        k for k in spec.get_vars() if k in CONVOLUTED_QUANTITIES
+    ]
     for q in wanted:
         if not q in rescaled_list:
             raise AssertionError(
diff --git a/radis/spectrum/spectrum.py b/radis/spectrum/spectrum.py
index 68612f873..2eedee164 100644
--- a/radis/spectrum/spectrum.py
+++ b/radis/spectrum/spectrum.py
@@ -59,7 +59,14 @@
 from radis.db.references import doi
 
 # from radis.lbl.base import print_conditions
-from radis.misc.arrays import count_nans, evenly_distributed, nantrapz
+from radis.misc.arrays import (
+    anynan,
+    count_nans,
+    evenly_distributed,
+    first_nonnan_index,
+    last_nonnan_index,
+    nantrapz,
+)
 from radis.misc.debug import printdbg
 from radis.misc.plot import split_and_plot_by_parts
 from radis.misc.signal import resample
@@ -69,7 +76,9 @@
 from radis.spectrum.rescale import rescale_mole_fraction, rescale_path_length, update
 from radis.spectrum.utils import (
     CONVOLUTED_QUANTITIES,
-    NON_CONVOLUTED_QUANTITIES,
+    WAVELEN_UNITS,
+    WAVELENVAC_UNITS,
+    WAVENUM_UNITS,
     cast_waveunit,
     format_xlabel,
     make_up,
@@ -94,13 +103,23 @@ class Spectrum(object):
 
     Parameters
     ----------
-    quantities: dict of tuples   {'quantity':(wavenum, quantity)}
+    quantities: dict of tuples   ``{'quantity':(w, a)}``  or dict  ``{'wavelength/wavenumber': w, quantity': a}``
         where quantities are spectral quantities (absorbance, radiance, etc.)
-        and wavenum is in :math:`cm^{-1}`
-        example::
+        and wavenum is in :math:`cm^{-1}` or :math:`nm` (see ``waveunit``)
+        Example::
+
+            # w, k, I are numpy arrays for wavenumbers, absorption coefficient, and radiance.
+            from radis import Spectrum
+            s = Spectrum({"wavenumber":w, "abscoeff":k, "radiance_noslit":I}, wunit='cm-1',
+                         units={"radiance_noslit":"mW/cm2/sr/nm", "abscoeff":"cm-1"})
+        Or::
+
+            s = Spectrum({"abscoeff":(w,k), "radiance_noslit":(w,I)},
+                         wunit="cm-1"
+                         units={"radiance_noslit":"mW/cm2/sr/nm", "abscoeff":"cm-1"})
 
-        {'radiance_noslit':(wavenum, radiance_noslit),
-        'absorbance':(wavenum, absorbance)}
+        See also: :py:meth:`~radis.spectrum.spectrum.Spectrum.from_array`
+        and :py:meth:`~radis.spectrum.spectrum.Spectrum.from_txt`
 
     units: dict
         units for quantities
@@ -109,16 +128,13 @@ class Spectrum(object):
     ----------------
     conditions: dict
         physical conditions and calculation parameters
-    cond_units: dict
-        units for conditions
-    wavespace: ``'nm'``, ``'cm-1'``, ``'nm_vac'`` or ``None``
+    wunit: ``'nm'``, ``'cm-1'``, ``'nm_vac'`` or ``None``
         wavelength in air (``'nm'``), wavenumber (``'cm-1'``), or wavelength in vacuum (``'nm_vac'``).
+        If ``None``, ``'wavespace'`` must be defined in ``conditions``.
         Quantities should be evenly distributed along this space for fast
         convolution with the slit function
-        If ``None``, ``'wavespace'`` must be defined in ``conditions``.
-        (non-uniform slit function is not implemented anyway... )
-        Defaults None (but raises an error if wavespace is not defined in
-        conditions neither)
+    cond_units: dict
+        units for conditions
 
 
     Other Parameters
@@ -240,10 +256,10 @@ class Spectrum(object):
         from radis import Spectrum, calculated_spectrum
         s1 = calculated_spectrum(w, I, wunit='nm', Iunit='mW/cm2/sr/nm')
         s2 = Spectrum.from_array(w, I, 'radiance_noslit',
-                               waveunit='nm', unit='mW/cm2/sr/nm')
+                               wunit='nm', unit='mW/cm2/sr/nm')
         s3 = Spectrum({'radiance_noslit': (w, I)},
                       units={'radiance_noslit':'mW/cm2/sr/nm'},
-                      waveunit='nm')
+                      wunit='nm')
 
     See more examples in the [Spectrum]_ page.
 
@@ -262,7 +278,7 @@ class Spectrum(object):
     -----
     Implementation:
 
-        quantities are stored in ``self._q`` and ``self._q_conv`` dictionaries.
+        quantities are stored in the ``self._q`` dictionary.
         They are better accessed with the :meth:`~radis.spectrum.spectrum.Spectrum.get`
         method that deals with units and wavespace
 
@@ -278,6 +294,12 @@ class Spectrum(object):
     ----------
     conditions : dict
         Stores computation / measurement conditions
+    c: dict
+        convenience wrapper to ``conditions``::
+
+            s.c["calculation_time"] is s.conditions["calculation_time"]
+            >> True
+
     populations: dict
         Stores molecules, isotopes, electronic states and vibrational or
         rovibrational populations
@@ -302,8 +324,8 @@ class Spectrum(object):
     # hardcode attribute names, but can save a lot of memory if hundreds of spectra
     __slots__ = [
         "_q",
-        "_q_conv",
         "units",
+        "c",
         "conditions",
         "cond_units",
         "populations",
@@ -322,17 +344,16 @@ def __init__(
         cond_units=None,
         populations=None,
         lines=None,
-        waveunit=None,
+        wunit=None,
         name=None,
         references={},
-        warnings=True,
+        check_wavespace=True,
+        **kwargs,
     ):
-        # TODO: make it possible to give quantities={'wavespace':w, 'radiance':I,
-        # 'transmittance':T, etc.} directly too (instead of {'radiance':(w,I), etc.})
-
         # TODO: add help on creating a Spectrum from a dictionary
 
         # Check inputs
+        # ---------------
         # ... Replace None attributes with dictionaries
         if conditions is None:
             conditions = {}
@@ -346,26 +367,56 @@ def __init__(
             references = {}
         self._init_annotations()  # typing hints for get() and plot()
 
+        if "warnings" in kwargs:
+            check_wavespace = kwargs.pop("warnings")
+            warn(
+                "Spectrum(warning=) renamed Spectrum(check_wavespace=) in 0.9.30",
+                DeprecationWarning,
+            )
+        # ... waveunit in Spectrum renamed wunit
+        if "waveunit" in kwargs:
+            assert wunit is None  # not given twice
+            warn(
+                "waveunit parameter in Spectrum(waveunit=) is now named `wunit`",
+                DeprecationWarning,
+            )
+            wunit = kwargs.pop("waveunit")
+        waveunit = wunit  # (but not renamed within the code or conditions yet)
+        if len(kwargs) > 0:
+            raise ValueError("Unexpected input: {0}".format(list(kwargs.keys())))
+
         # Deal with deprecated inputs
         # ... wavespace renamed waveunit
         if "wavespace" in conditions:
             warn(
-                DeprecationWarning("wavespace key in conditions is now named waveunit")
+                "wavespace key in conditions is now named waveunit", DeprecationWarning
             )
             conditions["waveunit"] = conditions["wavespace"]
             del conditions["wavespace"]
 
-        # Waveunit
-        # ... Cast in standard waveunit format
+        # ... make sure waveunit is given. Either implicitely, or explicitely:
         if waveunit is None and not "waveunit" in conditions:
-            raise AssertionError(
-                "waveunit ('nm', 'cm-1'?) has to be defined in `conditions`"
-                + "or with waveunit="
-            )
+            if "wavelength" in quantities:
+                raise AssertionError(
+                    "waveunit ('nm', 'nm_vac' ?) has to be defined in `conditions`"
+                    + "or with `wunit='nm'` or `wunit='nm_vac'`"
+                )
+            elif "wavenumber" in quantities:
+                raise AssertionError(
+                    "waveunit ('cm-1'?) has to be defined in `conditions`"
+                    + "or with `wunit='cm-1'`"
+                )
+            else:
+                raise AssertionError(
+                    "waveunit ('nm', 'cm-1'?) has to be defined in `conditions`"
+                    + "or with `wunit=`"
+                )
+        # ... Cast in standard waveunit format
         if waveunit is not None:
             waveunit = cast_waveunit(waveunit)
         if "waveunit" in conditions:
             conditions["waveunit"] = cast_waveunit(conditions["waveunit"])
+
         # ... Make sure unit match
         if "waveunit" in conditions:
             if waveunit is not None and conditions["waveunit"] != waveunit:
@@ -374,38 +425,81 @@ def __init__(
                         conditions["waveunit"], waveunit
                     )
                 )
-        else:  # ... or define them in dictionary
+        elif waveunit is not None:  # ... or define them in dictionary
             conditions["waveunit"] = waveunit
+        if "wavelength" in quantities:
+            assert "wavenumber" not in quantities and "wavespace" not in quantities
+            assert conditions["waveunit"] in WAVELEN_UNITS + WAVELENVAC_UNITS
+        if "wavenumber" in quantities:
+            assert "wavelength" not in quantities and "wavespace" not in quantities
+            assert conditions["waveunit"] in WAVENUM_UNITS
+        if "wavespace" in quantities:
+            assert "wavenumber" not in quantities and "wavelength" not in quantities
 
         # Check quantities format
         if len(quantities) == 0:
             raise AssertionError(
                 "Spectrum is created with no quantities. Add "
                 + "`radiance`, `transmittance` etc... with dict "
-                + "format. e.g: {`radiance`: (w,I)}"
+                + "format. e.g: {`wavenumber`:w, `radiance`: I}"
             )
 
-        self._q = {}  #: dict: stores non convoluted quantities
-        self._q_conv = {}  #: dict: stores convoluted (slit) quantities
+        # Create the arrays
+        # -----------------
+
+        self._q = {}  #: dict: stores spectral arrays
 
         self._slit = {}  #: dict: hold slit function
 
+        # infer format:
+        tuple_format = (
+            "wavelength" not in quantities
+            and "wavenumber" not in quantities
+            and "wavespace" not in quantities
+        )
+        k0 = list(quantities.keys())[0]
+        v0 = list(quantities.values())[0]
         for k, v in quantities.items():
-            try:
-                assert len(v) == 2
-            except AssertionError:
+            if len(v) != len(v0):
+                raise AssertionError(
+                    f"Input arrays should have the same length. Got : {k} : len {len(v)}, {k0} : len {len(v0)}"
+                )
+            if tuple_format and len(v) != 2:
                 raise AssertionError(
-                    "Attributes should have format (wavenumber, "
-                    + "quantity) . Error with `{0}`".format(k)
+                    "Input arrays should have format `{'quantity':(wavespace, array)} or {'wavelength/wavenumber':wavespace, 'quantity':array}` but not both. Got :"
+                    + "{0}".format(quantities),
                 )
 
-            w, I = v
+        if tuple_format:
+            for k, (w, I) in quantities.items():
 
-            # creates a copy of w,I
-            self._add_quantity(k, w, I, warnings=warnings)
+                # creates and stores a copy of w,I
+                self._add_quantity(k, w, I, check_wavespace=check_wavespace)
+        else:
+            for k, I in quantities.items():
+                if k in ["wavelength", "wavenumber", "wavespace"]:
+                    continue
+                # add:
+                if not "wavespace" in self._q:
+                    # adding wavenumber space
+                    if "wavespace" in quantities.keys():
+                        w = quantities["wavespace"]
+                    elif "wavelength" in quantities.keys():
+                        w = quantities["wavelength"]
+                    elif "wavenumber" in quantities.keys():
+                        w = quantities["wavenumber"]
+
+                    # creates and stores a copy of w,I :
+                    self._add_quantity(k, w, I, check_wavespace=True)
+                else:
+                    # this also checks that all arrays have same length
+                    self._add_quantity(k, w, I, check_wavespace=False)
 
         # Finally, add our attributes
-        self.conditions = conditions
+        self.conditions = self.c = conditions
+        """ dict: computation conditions, or experimetnal parameters, or
+        any metadata you need to store with the Spectrum object.
+        """
         self.populations = populations
         self.lines = lines
         self.units = units
@@ -425,17 +519,16 @@ def __init__(
     # %% Constructors
 
     @classmethod
-    def from_array(self, w, I, quantity, waveunit, unit, *args, **kwargs):
+    def from_array(self, w, I, quantity, wunit, unit, waveunit=None, *args, **kwargs):
         """Construct Spectrum from 2 arrays.
 
         Parameters
         ----------
-
         w, I: array
             waverange and vector
         quantity: str
             spectral quantity name
-        waveunit: ``'nm'``, ``'cm-1'``, ``'nm_vac'``
+        wunit: ``'nm'``, ``'cm-1'``, ``'nm_vac'``
             unit of waverange:         wavelength in air (``'nm'``), wavenumber
             (``'cm-1'``), or wavelength in vacuum (``'nm_vac'``).
         unit: str
@@ -445,7 +538,6 @@ def from_array(self, w, I, quantity, waveunit, unit, *args, **kwargs):
 
         Other Parameters
         ----------------
-
         conditions: dict
             physical conditions and calculation parameters
         cond_units: dict
@@ -480,7 +572,7 @@ def from_array(self, w, I, quantity, waveunit, unit, *args, **kwargs):
 
             from radis import Spectrum
             s = Spectrum.from_array(w, I, 'radiance_noslit',
-                                   waveunit='nm', unit='mW/cm2/sr/nm')
+                                   wunit='nm', unit='mW/cm2/sr/nm')
 
         To create a spectrum with absorption and emission components
         (e.g: ``radiance_noslit`` and ``transmittance_noslit``, or ``emisscoeff``
@@ -490,7 +582,7 @@ class directly. Ex::
             from radis import Spectrum
             s = Spectrum({'abscoeff': (w, A), 'emisscoeff': (w, E)},
                          units={'abscoeff': 'cm-1', 'emisscoeff':'W/cm2/sr/nm'},
-                         waveunit='nm')
+                         wunit='nm')
 
         .. minigallery:: radis.spectrum.spectrum.Spectrum.from_array
 
@@ -505,6 +597,13 @@ class directly. Ex::
         :func:`~radis.tools.database.load_spec`,
         :ref:`the Spectrum page <label_spectrum>`
         """
+        # Deprecated inputs
+        if waveunit is not None:
+            warn(
+                "`waveunit=` parameter in from_array is now named `wunit=`",
+                DeprecationWarning,
+            )
+            wunit = waveunit
 
         quantities = {quantity: (w, I)}
         units = {quantity: unit}
@@ -516,11 +615,11 @@ class directly. Ex::
                 conditions[k] = kwargs.pop(k)
 
         return self(
-            quantities, units, waveunit=waveunit, conditions=conditions, *args, **kwargs
+            quantities, units, wunit=wunit, conditions=conditions, *args, **kwargs
         )
 
     @classmethod
-    def from_txt(self, file, quantity, waveunit, unit, *args, **kwargs):
+    def from_txt(self, file, quantity, wunit, unit, waveunit=None, *args, **kwargs):
         """Construct Spectrum from txt file.
 
         Parameters
@@ -529,7 +628,7 @@ def from_txt(self, file, quantity, waveunit, unit, *args, **kwargs):
             file name
         quantity: str
             spectral quantity name
-        waveunit: ``'nm'``, ``'cm-1'``, ``'nm_vac'``
+        wunit: ``'nm'``, ``'cm-1'``, ``'nm_vac'``
             unit of waverange: wavelength in air (``'nm'``), wavenumber
             (``'cm-1'``), or wavelength in vacuum (``'nm_vac'``).
         unit: str
@@ -585,7 +684,7 @@ def from_txt(self, file, quantity, waveunit, unit, *args, **kwargs):
         ``delimiter`` key is forwarded to :py:func:`~numpy.loadtxt`::
 
             from radis import Spectrum
-            s = Spectrum.from_txt('spectrum.csv', 'radiance', waveunit='nm',
+            s = Spectrum.from_txt('spectrum.csv', 'radiance', wunit='nm',
                                       unit='W/cm2/sr/nm', delimiter=',')
 
 
@@ -597,7 +696,7 @@ class directly. Ex::
             from radis import Spectrum
             s = Spectrum({'abscoeff': (w, A), 'emisscoeff': (w, E)},
                          units={'abscoeff': 'cm-1', 'emisscoeff':'W/cm2/sr/nm'},
-                         waveunit='nm')
+                         wunit='nm')
 
         .. minigallery:: radis.spectrum.spectrum.Spectrum.from_txt
 
@@ -618,6 +717,13 @@ class directly. Ex::
         :func:`~radis.tools.database.load_spec`,
         :ref:`the Spectrum page <label_spectrum>`
         """
+        # Deprecated inputs
+        if waveunit is not None:
+            warn(
+                "`waveunit=` parameter in from_array is now named `wunit=`",
+                DeprecationWarning,
+            )
+            wunit = waveunit
 
         # Get input for loadtxt
         kwloadtxt = {}
@@ -639,9 +745,7 @@ class directly. Ex::
             if k in kwargs:
                 conditions[k] = kwargs.pop(k)
 
-        s = self(
-            quantities, units, waveunit=waveunit, conditions=conditions, *args, **kwargs
-        )
+        s = self(quantities, units, wunit=wunit, conditions=conditions, *args, **kwargs)
 
         # Store filename
         s.file = file
@@ -652,7 +756,7 @@ class directly. Ex::
     # ----------------
     # XXX =====================================================================
 
-    def get(self, var, wunit="default", Iunit="default", copy=True):
+    def get(self, var, wunit="default", Iunit="default", copy=True, trim_nan=False):
         """Retrieve a spectral quantity from a Spectrum object. You can select
         wavespace unit, intensity unit, or propagation medium.
 
@@ -676,9 +780,12 @@ def get(self, var, wunit="default", Iunit="default", copy=True):
 
         Other Parameters
         ----------------
-        copy: boolean
+        copy: bool
             if ``True``, returns a copy of the stored quantity (modifying it wont
             change the Spectrum object). Default ``True``.
+        trim_nan: bool
+            if ``True``, removes ``nan`` on the sides of the spectral array
+            (and corresponding wavespace). Default ``False``.
 
         Returns
         -------
@@ -720,17 +827,7 @@ def get(self, var, wunit="default", Iunit="default", copy=True):
                 )
 
         # Get quantity
-        if var in CONVOLUTED_QUANTITIES:
-            vartype = "convoluted"
-            I = self._q_conv[var]
-        elif var in NON_CONVOLUTED_QUANTITIES:  # non convoluted
-            vartype = "non_convoluted"
-            I = self._q[var]
-        elif var in self._q:  # not declared, but exists. Assumes it's non_convoluted?
-            vartype = "non_convoluted"
-            I = self._q[var]
-        else:
-            raise ValueError("Unexpected quantity: {0}".format(var))
+        I = self._q[var]
 
         # Copy
         if copy:
@@ -742,14 +839,24 @@ def get(self, var, wunit="default", Iunit="default", copy=True):
         else:
             wunit = cast_waveunit(wunit)
         if wunit == "cm-1":
-            w = self.get_wavenumber(vartype, copy=copy)
+            w = self.get_wavenumber(copy=copy)
         elif wunit == "nm":
-            w = self.get_wavelength(medium="air", which=vartype, copy=copy)
+            w = self.get_wavelength(medium="air", copy=copy)
         elif wunit == "nm_vac":
-            w = self.get_wavelength(medium="vacuum", which=vartype, copy=copy)
+            w = self.get_wavelength(medium="vacuum", copy=copy)
         else:
             raise ValueError(wunit)
 
+        # Trim nan if needed
+        if trim_nan:
+            m = first_nonnan_index(I)
+            M = last_nonnan_index(I)
+            if m is None:
+                raise ValueError("All values are nan. Check your data?")
+            if m > 0 or M < len(I) - 1:  # else, no change to be made
+                w = w[m : M + 1]
+                I = I[m : M + 1]
+
         # Convert y unit if necessary
         Iunit0 = self.units[var]
         if Iunit != "default" and Iunit != Iunit0:
@@ -788,16 +895,9 @@ def get(self, var, wunit="default", Iunit="default", copy=True):
 
         return w, I
 
-    def _get_wavespace(self, which="any", copy=True):
+    def _get_wavespace(self, copy=True):
         """Return wavespace (if the same for all quantities)
 
-        Parameters
-        ----------
-        which: 'convoluted', 'non_convoluted', ``'any'``
-            return wavelength for convoluted quantities, non convoluted quantities,
-            or any. If ``any`` and both are defined, they have to be the same else
-            an error is raised. Default ``any``.
-
         Other Parameters
         ----------------
         copy: boolean
@@ -807,60 +907,23 @@ def _get_wavespace(self, which="any", copy=True):
 
         Returns
         -------
-
         w: array_like
             (a copy of) spectrum wavespace for convoluted or non convoluted
             quantities
         """
 
-        if which == "any":
-            q_defined = "wavespace" in self._q
-            q_conv_defined = "wavespace" in self._q_conv
-            if q_defined and q_conv_defined:
-                if not len(self._q["wavespace"]) == len(self._q_conv["wavespace"]):
-                    raise ValueError(
-                        "All wavespace not equal for calculated "
-                        + "quantities. Can't use get_wavespace(). "
-                        + "Specify which=`convoluted` or `non_convoluted`"
-                    )
-                if not np.allclose(self._q["wavespace"], self._q_conv["wavespace"]):
-                    raise ValueError(
-                        "All wavespace not equal for calculated "
-                        + "quantities. Can't use get_wavespace(). "
-                        + "Specify which=`convoluted` or `non_convoluted`"
-                    )
-                w = self._q["wavespace"]
-            elif q_defined:
-                w = self._q["wavespace"]
-            else:
-                w = self._q_conv["wavespace"]
-
-        elif which == "convoluted":
-            w = self._q_conv["wavespace"]
-
-        elif which == "non_convoluted":
-            w = self._q["wavespace"]
-
-        else:
-            raise ValueError(
-                "which has to be one of 'convoluted', 'non_convoluted', "
-                + "'any'. Got {0}".format(which)
-            )
+        w = self._q["wavespace"]
 
         if copy:
             w = w.copy()
 
         return w
 
-    def get_wavelength(self, medium="air", which="any", copy=True):
+    def get_wavelength(self, medium="air", which=None, copy=True):
         """Return wavelength in defined medium.
 
         Parameters
         ----------
-        which: 'convoluted', 'non_convoluted', 'any'
-            return wavelength for convoluted quantities, non convoluted quantities,
-            or any. If any and both are defined, they have to be the same else
-            an error is raised. Default any.
         medium: ``'air'``, ``'vacuum'``
             returns wavelength as seen in air, or vacuum. Default ``'air'``.
             See :func:`~radis.phys.air.vacuum2air`, :func:`~radis.phys.air.air2vacuum`
@@ -881,13 +944,18 @@ def get_wavelength(self, medium="air", which="any", copy=True):
         --------
         :ref:`the Spectrum page <label_spectrum>`
         """
+        if which is not None:
+            raise DeprecationWarning(
+                "`which` parameter was deleted in Radis 0.9.30. Just use Spectrum.get_wavelength()"
+            )
+            # TODO: remove after 0.9.31
 
         # Check input
         if not medium in ["air", "vacuum"]:
             raise NotImplementedError("Unknown propagating medium: {0}".format(medium))
 
         # Now convert stored wavespace to the output unit
-        w = self._get_wavespace(which=which, copy=copy)
+        w = self._get_wavespace(copy=copy)
         if self.get_waveunit() == "cm-1":
             w = cm2nm(w)  # vacuum wavelength
 
@@ -914,17 +982,9 @@ def get_wavelength(self, medium="air", which="any", copy=True):
 
         return w
 
-    def get_wavenumber(self, which="any", copy=True):
+    def get_wavenumber(self, which=None, copy=True):
         """Return wavenumber (if the same for all quantities)
 
-        Parameters
-        ----------
-        which: 'convoluted', 'non_convoluted', 'any'
-            return wavenumber for convoluted quantities, non convoluted quantities,
-            or any. If any and both are defined, they have to be the same else
-            an error is raised. Default any.
-
-
         Other Parameters
         ----------------
         copy: boolean
@@ -937,7 +997,13 @@ def get_wavenumber(self, which="any", copy=True):
             (a copy of) spectrum wavenumber for convoluted or non convoluted
             quantities
         """
-        w = self._get_wavespace(which=which, copy=copy)
+        if which is not None:
+            raise DeprecationWarning(
+                "`which` parameter was deleted in Radis 0.9.30. Just use Spectrum.get_wavenumber()"
+            )
+            # TODO: remove after 0.9.31
+
+        w = self._get_wavespace(copy=copy)
 
         if self.get_waveunit() == "cm-1":  #
             pass
@@ -1296,6 +1362,41 @@ def crop(self, wmin=None, wmax=None, wunit="default", inplace=True):
 
         return crop(self, wmin=wmin, wmax=wmax, wunit=wunit, inplace=inplace)
 
+    def trim(self, inplace=True):
+        """Remove :py:attr:`~numpy.nan` common to all arrays on each side of the Spectrum.
+
+        Returns a smaller Spectrum (inplace or not).
+
+        Returns
+        -------
+        s: Spectrum : trimmed Spectrum. If ``inplace=True``, Spectrum has been updated
+              directly anyway. Allows :ref:`chaining <label_spectrum_chaining>`.
+        """
+        if inplace:
+            s = self
+        else:
+            s = self.copy()
+
+        trim_left = []
+        trim_right = []
+        for k in s.get_vars():
+            m = first_nonnan_index(s._q[k])
+            M = last_nonnan_index(s._q[k])
+            if m is None:
+                raise ValueError(
+                    f"All values are nan in {k}: Spectrum would be trimmed entirely. Are you sure?"
+                )
+            trim_left.append(m)
+            trim_right.append(M)
+        m = np.max(trim_left)
+        M = np.min(trim_right)
+
+        if m > 0 or M < len(s) - 1:  # else, no change to be made
+            for k, v in s._q.items():
+                s._q[k] = v[m : M + 1]
+
+        return s
+
     def sort(self, inplace=True):
         """Sort the Spectrum by wavelength / wavenumber.
 
@@ -1317,15 +1418,9 @@ def sort(self, inplace=True):
         else:
             s = self.copy()
 
-        if len(self._q) > 0:
-            b = np.argsort(s._q["wavespace"])
-            for k, v in s._q.items():
-                s._q[k] = v[b]
-
-        if len(s._q_conv) > 0:
-            b = np.argsort(s._q_conv["wavespace"])
-            for k, v in s._q_conv.items():
-                s._q_conv[k] = v[b]
+        b = np.argsort(s._q["wavespace"])
+        for k, v in s._q.items():
+            s._q[k] = v[b]
 
         return s
 
@@ -1405,7 +1500,7 @@ def get_integral(self, var, wunit="default", Iunit="default", **kwargs):
         """
 
         w, I = self.get(var, wunit=wunit, Iunit=Iunit, **kwargs)
-        return abs(np.trapz(I, x=w))
+        return abs(nantrapz(I, w))
 
     def get_power(self, unit="mW/cm2/sr"):
         """Returns integrated radiance (no slit) power density.
@@ -1437,42 +1532,65 @@ def get_power(self, unit="mW/cm2/sr"):
         # P is in mW/cm2/sr/nm * nm
         return conv2(P, "mW/cm2/sr", unit)
 
+    def has_nan(self, ignore_wavespace=True) -> bool:
+        """
+
+        Parameters
+        ----------
+        s : Spectrum
+            radis Spectrum.
+
+        Returns
+        -------
+        b : bool
+            returns whether Spectrum has ``nan``
+
+        Note
+        ----
+
+        ``print(s)`` will also show which spectral quantities have ````nan.
+
+        """
+
+        for k, v in self._q.items():
+            if k == "wavespace" and ignore_wavespace:
+                continue
+            if anynan(v):
+                return True
+        return False
+
     # %% Plotting routines
 
-    def get_vars(self, which="any"):
+    def get_vars(self, which=None):
         """Returns all spectral quantities stored in this object (convoluted or
         non convoluted)
 
-        Parameters
-        ----------
-        which: 'any', 'convoluted', 'non convoluted'
         """
-        if which == "any":
-            varlist = list(self._q.keys()) + list(self._q_conv.keys())
-        elif which == "convoluted":
-            varlist = list(self._q_conv.keys())
-        elif which == "non convoluted":
-            varlist = list(self._q.keys())
-        else:
-            raise ValueError("Unexpected value: {0}".format(which))
+        if which is not None:
+            raise DeprecationWarning(
+                "`which` parameter was deleted in Radis 0.9.30. Just use Spectrum.get_vars()"
+            )
+            # TODO: remove after 0.9.31
 
         # remove wavespace
-        varlist = [k for k in varlist if k != "wavespace"]
+        varlist = [k for k in self._q.keys() if k != "wavespace"]
         return varlist
 
-    def get_quantities(self, which="any"):
+    def get_quantities(self, which=None):
         """Returns all spectral quantities stored in this object (convoluted or
         non convoluted). Wrapper to
         :py:meth:`~radis.spectrum.spectrum.get_vars`
 
-        Parameters
-        ----------
-        which: 'any', 'convoluted', 'non convoluted'
         """
+        if which is not None:
+            raise DeprecationWarning(
+                "`which` parameter was deleted in Radis 0.9.30. Just use Spectrum.get_quantities()"
+            )
+            # TODO: remove after 0.9.31
 
-        return self.get_vars(which=which)
+        return self.get_vars()
 
-    def _get_items(self):
+    def _get_items(self) -> dict:
         """Return a dictionary of tuples, e.g::
 
             {'radiance':(w,I), 'transmittance_noslit':(w_ns,T)}
@@ -1486,13 +1604,6 @@ def _get_items(self):
         items = {
             k: (self._q["wavespace"], v) for k, v in self._q.items() if k != "wavespace"
         }
-        items.update(
-            {
-                k: (self._q_conv["wavespace"], v)
-                for k, v in self._q_conv.items()
-                if k != "wavespace"
-            }
-        )
         return items
 
     def plot(
@@ -1508,7 +1619,7 @@ def plot(
         force=False,
         plot_by_parts=False,
         show_ruler=False,
-        **kwargs
+        **kwargs,
     ):
         """Plot a :py:class:`~radis.spectrum.spectrum.Spectrum` object.
 
@@ -2060,8 +2171,9 @@ def apply_slit(
         slit_dispersion_threshold=0.01,
         auto_recenter_crop=True,
         verbose=True,
+        inplace=True,
         *args,
-        **kwargs
+        **kwargs,
     ):
         """Apply an instrumental slit function to all quantities in Spectrum.
         Slit function can be generated with usual shapes (see ``shape=``) or
@@ -2167,6 +2279,14 @@ def apply_slit(
         slit_dispersion_warning_threshold: float
             if not ``None``, check that slit dispersion is about constant (< ``threshold`` change)
             on the calculated range. Default 0.01 (1%). See :func:`~radis.tools.slit.offset_dilate_slit_function`
+        inplace: bool
+            if ``True``, adds convolved arrays directly in the Spectrum. If
+            ``False``, returns a new Spectrum with only the convolved arrays.
+            Note: if you want a new Spectrum with both the convolved and
+            non convolved quantities, use ::
+
+                s.copy().apply_slit()
+
         *args, **kwargs
             are forwarded to slit generation or import function
         verbose: bool
@@ -2179,7 +2299,9 @@ def apply_slit(
         Returns
         -------
         Spectrum : same Spectrum, with new spectral arrays.
-            Allows :ref:`chaining <label_spectrum_chaining>`
+            Allows :ref:`chaining <label_spectrum_chaining>`.
+            If ``inplace=False``, return a new Spectrum with the
+            new spectral arrays only.
 
 
         Notes
@@ -2257,8 +2379,9 @@ def apply_slit(
 
         if len(varlist) == 0:
             raise AssertionError(
-                "No variables to apply slit on. Variable names "
-                + "to be convolved should end with _noslit"
+                "No spectral arrays to apply slit on. Spectral arrays "
+                + "to be convolved should end with _noslit. "
+                + f"Available arrays in Spectrum : {self.get_vars()}"
             )
 
         # Forward relevant inputs to convolution instead of slit function generation
@@ -2305,7 +2428,7 @@ def apply_slit(
             verbose=verbose,
             plot=plot_slit,
             *args,
-            **kwargs
+            **kwargs,
         )
 
         # Check if dispersion is specified
@@ -2373,61 +2496,94 @@ def apply_slit(
                     wslit,
                     Islit,
                     mode=mode,
-                    waveunit=waveunit,
+                    wunit=waveunit,
                     verbose=verbose,
                     assert_evenly_spaced=False,
                     # assumes Spectrum is correct by construction
-                    **kwargsconvolve
+                    **kwargsconvolve,
                 )
 
                 if i == 0:
                     w_conv_slices.append(w_conv_window)
                 I_conv_slices[q].append(I_conv_window)
 
-        # Merge and store all variables
-        # ---------
+        # Get units
+        new_units = {}
         for q in I_conv_slices.keys():
             qns = q + "_noslit"
 
-            # Merge all slices
-            w_conv = np.hstack(w_conv_slices)
-            I_conv = np.hstack(I_conv_slices[q])
-
-            # Store
-            self._q_conv["wavespace"] = w_conv
-            self._q_conv[q] = I_conv
-
             # Get units
             if norm_by == "area":
-                self.units[q] = self.units[qns]
+                new_units[q] = self.units[qns]
             elif norm_by == "max":
                 new_unit = (Unit(unit) * Unit(self.units[qns])).to_string()
                 # because it's like if we multiplied by slit FWHM in the wavespace
                 # it was generated
-                self.units[q] = new_unit
+                new_units[q] = new_unit
             # Note: there was another mode called 'max2' where, unlike 'max',
             # unit was multiplied by [unit] not [return_unit]
             # Removed for simplification. You should stay with norm_by='area' anyway
             else:
                 raise ValueError("Unknown normalization type: {0}".format(norm_by))
 
+        # Merge and store all variables
+        # ---------
+        w_conv = np.hstack(w_conv_slices)
+        if inplace:
+            if len(self._q["wavespace"]) != len(w_conv) or not np.allclose(
+                self._q["wavespace"], w_conv
+            ):
+                raise AssertionError(
+                    "Wavespace of convolved arrays is different, cannot store it in the same Spectrum. You can use Spectrum.apply_slit(inplace=False) to return a new spectrum with only the convolved arrays"
+                )
+            for q in I_conv_slices.keys():
+                # Merge all slices
+                I_conv = np.hstack(I_conv_slices[q])
+
+                # Store
+                self._q[q] = I_conv
+
+                # Get units
+                self.units[q] = new_units[q]
+            s_out = self
+        else:
+            # make new spectrum
+            quantities = {
+                "wavespace": w_conv
+            }  # a copy will be created in Spectrum creation
+
+            for q in I_conv_slices.keys():
+
+                # Merge all slices
+                quantities[q] = np.hstack(I_conv_slices[q])
+
+            s_out = Spectrum(
+                quantities,
+                units=new_units,
+                conditions=self.conditions.copy(),
+                cond_units=self.cond_units.copy(),
+                wunit=waveunit,
+                name=self.name,
+                check_wavespace=False,
+            )
+
         # Store slit in Spectrum, in the Spectrum unit
         if store:
-            self._slit["wavespace"] = wslit0  # in 'waveunit'
-            self._slit["intensity"] = Islit0
+            s_out._slit["wavespace"] = wslit0  # in 'waveunit'
+            s_out._slit["intensity"] = Islit0
 
         # Update conditions
-        self.conditions["slit_function"] = slit_function
-        self.conditions["slit_unit"] = unit  # input slit unit
-        self.conditions["slit_dispersion"] = slit_dispersion
-        self.conditions["slit_dispersion_threshold"] = slit_dispersion_threshold
-        self.conditions["slit_shape"] = shape
+        s_out.conditions["slit_function"] = slit_function
+        s_out.conditions["slit_unit"] = unit  # input slit unit
+        s_out.conditions["slit_dispersion"] = slit_dispersion
+        s_out.conditions["slit_dispersion_threshold"] = slit_dispersion_threshold
+        s_out.conditions["slit_shape"] = shape
         # TODO: probably removed after Spectrum is stored.
-        self.conditions["norm_by"] = norm_by
+        s_out.conditions["norm_by"] = norm_by
 
-        return self  # to be able to chain: s.apply_slit().plot()
+        return s_out  # to be able to chain: s.apply_slit().plot()
 
-    def get_slit(self, unit="same"):
+    def get_slit(self, wunit="same"):
         """Get slit function that was applied to the Spectrum.
 
         Returns
@@ -2438,7 +2594,7 @@ def get_slit(self, unit="same"):
             :meth:`~radis.spectrum.spectrum.Spectrum.get_waveunit`
         """
 
-        if not unit in ["same", self.get_waveunit()]:
+        if not wunit in ["same", self.get_waveunit()]:
             raise NotImplementedError(
                 "Unit must be Spectrum waveunit: {0}".format(self.get_waveunit())
             )
@@ -2456,7 +2612,7 @@ def get_slit(self, unit="same"):
 
         return wslit, Islit
 
-    def plot_slit(self, wunit=None):
+    def plot_slit(self, wunit=None, waveunit=None):
         """Plot slit function that was applied to the Spectrum.
 
         If dispersion was used (see :meth:`~radis.spectrum.spectrum.Spectrum.apply_slit`)
@@ -2481,6 +2637,13 @@ def plot_slit(self, wunit=None):
 
         :ref:`the Spectrum page <label_spectrum>`
         """
+        # Deprecated inputs
+        if waveunit is not None:
+            warn(
+                "`waveunit=` parameter in from_array is now named `wunit=`",
+                DeprecationWarning,
+            )
+            wunit = waveunit
 
         from radis.tools.slit import (
             normalize_slit,
@@ -2514,7 +2677,7 @@ def plot_slit(self, wunit=None):
 
         # Plot in correct unit  (plot_slit deals with the conversion if needed)
         fig, ax = plot_slit(
-            wslit0, Islit0, waveunit=waveunit, plot_unit=wunit, Iunit=Iunit
+            wslit0, Islit0, wunit=waveunit, plot_unit=wunit, Iunit=Iunit
         )
 
         # Plot other slit functions if dispersion was applied:
@@ -2533,7 +2696,7 @@ def plot_slit(self, wunit=None):
                     wslit0_nm = vacuum2air(wslit0)
                 else:
                     wslit0_nm = cm2nm(wslit0)
-                w_nm = self.get_wavelength(medium="air", which="non_convoluted")
+                w_nm = self.get_wavelength(medium="air")
                 slice_windows = _cut_slices(
                     w_nm, wslit0, slit_dispersion, slit_dispersion_threshold
                 )
@@ -2565,7 +2728,7 @@ def plot_slit(self, wunit=None):
                     plot_slit(
                         wslit,
                         Islit,
-                        waveunit=waveunit,
+                        wunit=waveunit,
                         plot_unit=wunit,
                         Iunit=Iunit,
                         ls="--",
@@ -2583,7 +2746,7 @@ def line_survey(
         writefile=None,
         cutoff=None,
         *args,
-        **kwargs
+        **kwargs,
     ):
         """Plot Line Survey (all linestrengths used for calculation) Output in
         Plotly (html)
@@ -2715,7 +2878,7 @@ def get_overlay(overlay):
                 writefile=writefile,
                 cutoff=cutoff,
                 *args,
-                **kwargs
+                **kwargs,
             )
 
         else:
@@ -2889,14 +3052,13 @@ def resample(
         w_new,
         unit="same",
         out_of_bounds="nan",
-        if_conflict_drop="error",
-        energy_threshold=1e-3,
+        energy_threshold=5e-3,
         print_conservation=False,
         inplace=True,
-        **kwargs
+        if_conflict_drop=None,
+        **kwargs,
     ):
-        """Resample spectrum over a new wavelength. Fills with transparent
-        medium when out of bound (transmittance 1, radiance 0)
+        """Resample spectrum over a new wavelength/wavenumber range.
 
         .. warning::
             This may result in information loss. Resampling is done with
@@ -2908,6 +3070,9 @@ def resample(
 
             s_highres.resample(s_lowres)
 
+        Fills with ``'nan'`` or transparent medium (transmittance 1, radiance 0)
+        when out of bound  (see ``out_of_bounds``)
+
 
         Parameters
         ----------
@@ -2926,15 +3091,9 @@ def resample(
             but resampled in `nm` will be stored in `nm` from now on).
             If ``'nm'``, wavelength in air. If ``'nm_vac'``, wavelength in vacuum.
         out_of_bounds: ``'transparent'``, ``'nan'``, ``'error'``
-            what to do if resampling is out of bounds. 'transparent': fills with
-            transparent medium. 'nan': fill with nan. 'error': raises an error.
+            what to do if resampling is out of bounds. ``'transparent'``: fills with
+            transparent medium. 'nan': fill with nan. ``'error'``: raises an error.
             Default ``'nan'``
-        if_conflict_drop: ``'error'``, ``'convoluted'``, ``'non_convoluted'``
-            There is a problem if both convoluted and non convoluted (*no_slit)
-            quantities coexists, as they aren't scaled on the same wavespace
-            grid. If 'error' an error is raised. If 'convoluted', convoluted
-            quantities will be dropped. If 'non_convoluted' non convoluted quantities
-            are dropped. Default ``'error'``
         medium: ``'air'``, ``'vacuum'``, or ``'default'``
             in which medium is the new waverange is calculated if it is given
             in 'nm'. Ignored if unit='cm-1'
@@ -2942,11 +3101,10 @@ def resample(
 
         Other Parameters
         ----------------
-        *Inputs forwarded to :func:`radis.misc.signal.resample`*
-
-        energy_threshold: float
-            if energy conservation (integrals) is above this threshold, raise an
-            error.
+        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
@@ -2969,42 +3127,18 @@ def resample(
         --------
         :func:`radis.misc.signal.resample`
         """
-        # TODO (but dangerous): reapply_slit at the end of the process if slit
-        # is in conditions?
+        # Check inputs (check for deprecated)
+        if if_conflict_drop is not None:
+            raise DeprecationWarning(
+                "`if_conflict_drop` parameter was deleted in Radis 0.9.30"
+            )
+            # TODO: remove after 0.9.31
 
         if inplace:
             s = self
         else:
             s = self.copy()
 
-        # Check inputs (check for deprecated)
-
-        # ... see if convoluted / non convoluted values co-exist
-        if "wavespace" in s._q and "wavespace" in s._q_conv:
-            try:
-                assert len(s._q["wavespace"]) == len(s._q_conv["wavespace"])
-                assert np.allclose(s._q["wavespace"], s._q_conv["wavespace"])
-            except AssertionError:  # wavespaces are not the same.
-                if if_conflict_drop == "convoluted":
-                    for q in list(s._q_conv.keys()):
-                        del s._q_conv[q]
-                elif if_conflict_drop == "non_convoluted":
-                    for q in list(s._q.keys()):
-                        del s._q[q]
-                elif if_conflict_drop == "error":
-                    raise ValueError(
-                        "Cant resample as there are convoluted and non "
-                        + "convoluted quantities in the Spectrum object (and "
-                        + "wavespace are not the same). Use "
-                        + "`if_conflict_drop='convoluted' or 'non_convoluted'`"
-                    )
-                else:
-                    raise ValueError(
-                        "Unknown value for if_conflict_drop: {0}".format(
-                            if_conflict_drop
-                        )
-                    )
-
         # Get wavespace units
         stored_waveunit = s.get_waveunit()  # spectrum unit
         if unit == "same":  # resampled unit
@@ -3039,10 +3173,6 @@ def resample(
         if unit != stored_waveunit:
             s.conditions["waveunit"] = unit
 
-        # Get wavespace
-        update_q = "wavespace" in s._q
-        update_q_conv = "wavespace" in s._q_conv
-
         # Now let's resample
         def get_filling(variable):
             """Get out of bounds values for spectral quantity `variable`"""
@@ -3067,41 +3197,24 @@ def get_filling(variable):
         # ... using the (safer) .get() function because it's much faster (the
         # ... air2vacuum conversion in particular is quite slow, but has been
         # ... done once for all with get_wavelength() above )
-        if update_q:
-            for (k, I) in s._q.items():
-                if k == "wavespace":
-                    continue
-                fill_with = get_filling(k)
-                Inew = resample(
-                    w,
-                    I,
-                    w_new,
-                    ext=fill_with,
-                    energy_threshold=energy_threshold,
-                    print_conservation=False,
-                    **kwargs
-                )
-                s._q[k] = Inew
-            # update wavespace
-            s._q["wavespace"] = w_new
 
-        if update_q_conv:
-            for (k, I) in s._q_conv.items():
-                if k == "wavespace":
-                    continue
-                fill_with = get_filling(k)
-                Inew = resample(
-                    w,
-                    I,
-                    w_new,
-                    ext=fill_with,
-                    energy_threshold=energy_threshold,
-                    print_conservation=False,
-                    **kwargs
-                )
-                s._q_conv[k] = Inew
-            # update wavespace
-            s._q_conv["wavespace"] = w_new
+        for (k, I) in s._q.items():
+            if k == "wavespace":
+                continue
+            fill_with = get_filling(k)
+            Inew = resample(
+                w,
+                I,
+                w_new,
+                ext=fill_with,
+                energy_threshold=energy_threshold,
+                print_conservation=False,
+                **kwargs,
+            )
+
+            s._q[k] = Inew
+        # update wavespace
+        s._q["wavespace"] = w_new
 
         return s
 
@@ -3270,10 +3383,6 @@ def __copy__(self, copy_lines=True, quantity="all"):
         if quantity == "all":
             quantities = dict(self._get_items())
         else:
-            #            assert quantity in CONVOLUTED_QUANTITIES+NON_CONVOLUTED_QUANTITIES
-            #            if not quantity in self.get_vars():
-            #                raise ValueError("Spectrum {0} has no quantity '{1}'. Got: {2}".format(
-            #                        self.get_name(), quantity, self.get_vars()))
             quantities = {
                 quantity: self.get(quantity, wunit=self.get_waveunit())
             }  # dict(self._get_items())
@@ -3315,7 +3424,7 @@ def __copy__(self, copy_lines=True, quantity="all"):
             populations=populations,
             lines=lines,
             units=units,
-            waveunit=waveunit,
+            wunit=waveunit,
             references=references,
             name=name,
             warnings=False,  # saves about 3.5 ms on the Performance test object
@@ -3348,7 +3457,7 @@ def compare_with(
         ignore_nan=False,
         ignore_outliers=False,
         normalize=False,
-        **kwargs
+        **kwargs,
     ):
         """Compare Spectrum with another Spectrum object.
 
@@ -3426,7 +3535,7 @@ def compare_with(
             ignore_nan=ignore_nan,
             ignore_outliers=ignore_outliers,
             normalize=normalize,
-            **kwargs
+            **kwargs,
         )
 
     def cite(self, format="bibentry"):
@@ -3611,76 +3720,50 @@ def _init_annotations(self):
         except AttributeError:
             pass  # old Python version
 
-    def _add_quantity(self, name, w, I, warnings=True):
+    def _add_quantity(self, name, w, I, check_wavespace=True):
         """Add quantity.
 
         Note: creates a copy of the input array
+
+        If check_wavespace, check that array matches existing wavespace, and
+        that new wavespace is evenly spaced (required for the slit)
+        Note: this check takes a lot of time!  (few ms)
         """
 
         assert len(w) == len(I)
 
-        def check_wavespace(w):
-            """If warnings, check that array is evenly spaced. Returns a copy
-            of input array.
-
-            Note: this check takes a lot of time!  (few ms)
-            Is is not performed if warnings is False
-            """
-            if warnings:
-                # Check Wavelength/wavenumber is evently spaced
-                if not evenly_distributed(w, tolerance=1e-5):
-                    warn(
-                        "Wavespace is not evenly spaced ({0:.3f}%) for {1}.".format(
-                            np.abs(np.diff(w)).max() / w.mean() * 100, name
+        # Add wavespace
+        if "wavespace" in self._q:
+            if check_wavespace:
+                # Check new wavespace match the existing one
+                if not np.allclose(w, self._q["wavespace"]):
+                    raise ValueError(
+                        "wavespace for {0} doesnt correspond to existing wavespace".format(
+                            name
                         )
-                        + " This may create problems when convolving with slit function"
+                        + " for non convoluted quantities"
                     )
-            return np.array(w)  # copy
-
-        if name in CONVOLUTED_QUANTITIES:
-            # Add wavespace
-            if "wavespace" in self._q_conv:
-                if warnings:
-                    # Check new wavespace match the existing one
-                    if not np.allclose(w, self._q_conv["wavespace"]):
-                        raise ValueError(
-                            "wavespace for {0} doesnt correspond to existing wavespace".format(
-                                name
-                            )
-                            + " for convoluted quantities"
-                        )
             else:
-                self._q_conv["wavespace"] = np.array(w)  # copy
+                pass  # nothing to be done
+        else:
+            if name in CONVOLUTED_QUANTITIES:
+                self._q["wavespace"] = np.array(w)  # copy
                 # no need to check if wavespace is evenly spaced: we won't
                 # apply the slit function again
-
-            # Add quantity itself
-            self._q_conv[name] = np.array(I)  # copy
-
-        elif name in NON_CONVOLUTED_QUANTITIES:
-            # Add wavespace
-            if "wavespace" in self._q:
-                if warnings:
-                    # Check new wavespace match the existing one
-                    if not np.allclose(w, self._q["wavespace"]):
-                        raise ValueError(
-                            "wavespace for {0} doesnt correspond to existing wavespace".format(
-                                name
+            else:
+                # Check Wavelength/wavenumber is evently spaced
+                if check_wavespace:
+                    if not evenly_distributed(w, tolerance=1e-5):
+                        warn(
+                            "Wavespace is not evenly spaced ({0:.3f}%) for {1}.".format(
+                                np.abs(np.diff(w)).max() / w.mean() * 100, name
                             )
-                            + " for non convoluted quantities"
+                            + " This may create problems when convolving with slit function"
                         )
-            else:
-                self._q["wavespace"] = check_wavespace(w)  # copy
-
-            # Add quantity itself
-            self._q[name] = np.array(I)  # copy
+                self._q["wavespace"] = np.array(w)  # copy
 
-        else:
-            raise ValueError(
-                "Unknown quantity: {0}. Expected one of: {1}".format(
-                    name, CONVOLUTED_QUANTITIES + NON_CONVOLUTED_QUANTITIES
-                )
-            )
+        # Add quantity itself
+        self._q[name] = np.array(I)  # copy
 
         # also make the quantity accessible with s.[name] like Pandas dataframes (removed eventually)
         # setattr(self, name, quantity)   # Warning this makes another copy of it (it's a tuple!)
@@ -3727,30 +3810,32 @@ def __str__(self):
         # Print spectral quantities
         print("Spectral Quantities")
         print("-" * 40)
-        for k, v in self._get_items().items():
+        for k, v in self._q.items():
+            if k == "wavespace":
+                continue
             # print number of points with a comma separator
             print(
                 " " * 2,
                 k,
                 "\t({0:,d} points{1})".format(
-                    len(v[0]),
-                    ", {0} nans".format(count_nans(v[1]))
-                    if count_nans(v[1]) > 0
-                    else "",
+                    len(v),
+                    ", {0} nans".format(count_nans(v)) if anynan(v) else "",
                 ),
             )
 
         # Print populations
-        print("Populations Stored")
-        print("-" * 40)
-        try:
-            for k, v in self.populations.items():
-                print(" " * 2, k, "\t\t", list(v.keys()))
-        except:
-            pass
+        if self.populations:
+            print("Populations Stored")
+            print("-" * 40)
+            try:
+                for k, v in self.populations.items():
+                    print(" " * 2, k, "\t\t", list(v.keys()))
+            except:
+                pass
 
         # Print conditions
-        self.print_conditions()
+        if self.conditions:
+            self.print_conditions()
 
         return ""  # self.print_conditions()
 
@@ -3819,7 +3904,7 @@ def max(self):
 
         var = self._get_unique_var(operation_name="max")
         w, I = self.get(var, wunit=self.get_waveunit(), copy=False)
-        return I.max()
+        return I[~np.isnan(I)].max()
 
     def min(self):
         """Minimum of the Spectrum, if only one spectral quantity is available
@@ -3837,7 +3922,7 @@ def min(self):
 
         var = self._get_unique_var(operation_name="min")
         w, I = self.get(var, wunit=self.get_waveunit(), copy=False)
-        return I.min()
+        return I[~np.isnan(I)].min()
 
     def normalize(
         self,
@@ -4384,15 +4469,7 @@ def __len__(self):
         """Length of a Spectrum object = length of the wavespace if unique,
         else raises an error"""
 
-        # raises ValueError if both convolved and non convolved are defined
-        try:
-            return len(self._get_wavespace("any", copy=False))
-        except ValueError:
-            raise ValueError(
-                "All quantities do not have the same length in the Spectrum : {0}".format(
-                    {k: len(self.get(k)[0]) for k in self.get_vars()}
-                )
-            )
+        return len(self._get_wavespace(copy=False))
 
 
 # %% Private functions
@@ -4479,33 +4556,6 @@ def _cut_slices(w_spec_nm, w_slit_nm, slit_dispersion, slit_dispersion_threshold
 # -------------------
 # XXX =====================================================================
 
-# Test class
-
-
-def is_spectrum(a):
-    """Returns whether a is a Spectrum object.
-
-    Parameters
-    ----------
-    a: anything
-        a Python object
-
-    Returns
-    -------
-    bool: True if a is a Spectrum object
-
-    Notes
-    -----
-    is_spectrum compares the object class name (str): in some cases the Spectrum
-    class gets imported twice (when databases are involved, mostly), and a purely
-    isinstance() comparison fails
-    """
-    return isinstance(a, Spectrum)
-    # removed: was used initially in the early RADIS development phase. Spectrum
-    # object would not be recognized if the library was modified. The following
-    # line is more flexible :
-    # return isinstance(a, Spectrum) or repr(a.__class__) == repr(Spectrum)
-
 
 # %% Test functions
 if __name__ == "__main__":
diff --git a/radis/spectrum/utils.py b/radis/spectrum/utils.py
index fe47c24f5..fe99f4fea 100644
--- a/radis/spectrum/utils.py
+++ b/radis/spectrum/utils.py
@@ -10,8 +10,6 @@
 
 import subprocess
 
-import numpy as np
-
 from radis.misc.basics import partition
 
 # %% Definitions
@@ -153,9 +151,16 @@ def make_up(label):
     label: str
     """
 
+    # Hardfix for missing characters
+    SUPERSCRIPT_MINUS = "⁻"
+    from matplotlib.font_manager import FontProperties, findfont
+
+    if "arial" in findfont(FontProperties(family=["sans-serif"])):
+        SUPERSCRIPT_MINUS = "$^{-}$"
+
     # Improve units
-    label = label.replace(r"cm-1", r"cm⁻¹")
-    label = label.replace(r"m^-1", r"m⁻¹")
+    label = label.replace(r"cm-1", r"cm" + f"{SUPERSCRIPT_MINUS}" + r"¹")
+    label = label.replace(r"m^-1", r"m" + f"{SUPERSCRIPT_MINUS}" + r"¹")
     label = label.replace(r"m2", r"m²")
     label = label.replace(r"m3", r"m³")
     label = label.replace(r"I/I0", r"I/I₀")  # transmittance unit
@@ -345,32 +350,6 @@ def print_param(k):
     return None
 
 
-def has_nan(s):
-    """
-
-    Parameters
-    ----------
-    s : Spectrum
-        radis Spectrum.
-
-    Returns
-    -------
-    b : bool
-        returns whether Spectrum has ``nan``
-
-    Note
-    ----
-
-    ``print(s)`` will also show which spectral quantities have ````nan.
-
-    """
-
-    for k, v in s._get_items().items():
-        if np.isnan(v[1]).any():
-            return True
-    return False
-
-
 # %% Plot helper
 
 
diff --git a/radis/tools/database.py b/radis/tools/database.py
index 979b9d4b7..3e4dd17a9 100644
--- a/radis/tools/database.py
+++ b/radis/tools/database.py
@@ -73,7 +73,7 @@ class to manage them all.
 from radis.misc.basics import all_in, is_float, list_if_float
 from radis.misc.debug import printdbg
 from radis.misc.printer import printr
-from radis.spectrum.spectrum import Spectrum, is_spectrum
+from radis.spectrum.spectrum import Spectrum
 
 _scalable_inputs = ["mole_fraction", "path_length"]
 
@@ -133,7 +133,7 @@ def is_jsonable(x):
 
 
 def save(
-    s,
+    s: Spectrum,
     path,
     discard=[],
     compress=True,
@@ -232,7 +232,7 @@ def save(
     return fout  # return final name
 
 
-def _format_to_jsondict(s, discard, compress, verbose=True):
+def _format_to_jsondict(s: Spectrum, discard, compress, verbose=True):
     """Format to JSON writable dictionary.
 
     Notes
@@ -275,7 +275,6 @@ def _format_to_jsondict(s, discard, compress, verbose=True):
     # from the rest)
     if compress >= 2:
         sjson["_q"] = sjson["_q"].copy()
-        sjson["_q_conv"] = sjson["_q_conv"].copy()
         sjson = _compress(s, sjson)
 
     return sjson
@@ -375,7 +374,7 @@ def _get_fout_name(path, if_exists_then, add_date, add_info, sjson, verbose):
     return fout
 
 
-def _compress(s, sjson):
+def _compress(s: Spectrum, sjson):
     """removes all quantities that can be regenerated with s.update(), e.g,
     transmittance if abscoeff and path length are given, radiance if emisscoeff
     and abscoeff are given in non-optically thin case, etc.
@@ -398,12 +397,6 @@ def _compress(s, sjson):
         if redundant[key]:
             del sjson["_q"][key]
             discarded.append(key)
-    for key in list(sjson["_q_conv"].keys()):
-        if key == "wavespace":
-            continue
-        if redundant[key]:
-            del sjson["_q_conv"][key]
-            discarded.append(key)
 
     if len(discarded) > 0:
         print(
@@ -494,10 +487,12 @@ def _load(binary):
     return s
 
 
-def _json_to_spec(sload, file=""):
+def _json_to_spec(sload, file="") -> Spectrum:
     """Builds a Spectrum object from a JSON dictionary. Called by
     :func:`~radis.tools.database.load_spec`.
 
+    Json has been fixed from deprecating changes by _fix_format
+
     Parameters
     ----------
     sload: dict
@@ -529,13 +524,6 @@ def _json_to_spec(sload, file=""):
             for k, v in sload["_q"].items()
             if k != "wavespace"
         }
-        quantities.update(
-            {
-                k: (sload["_q_conv"]["wavespace"], v)
-                for k, v in sload["_q_conv"].items()
-                if k != "wavespace"
-            }
-        )
 
     # Generate spectrum:
     waveunit = sload["conditions"]["waveunit"]
@@ -590,7 +578,7 @@ def _json_to_spec(sload, file=""):
     s = Spectrum(
         quantities=quantities,
         conditions=conditions,
-        waveunit=waveunit,
+        wunit=waveunit,
         references=references,
         **kwargs,
     )
@@ -609,6 +597,10 @@ def _fix_format(file, sload):
     load old format precomputed spectra, and fix their attribute names.
 
     Save them again to fix the warnigns definitly.
+
+    Returns
+    -------
+    json, fixed:  fixed is True if a change was made.
     """
 
     fixed = False
@@ -629,6 +621,30 @@ def _fix_format(file, sload):
         sload["_q_conv"] = sload.pop("q_conv")
         fixed = True
 
+    # Fix _q_conv removal (0.9.30)
+    if "_q_conv" in sload:
+        if len(sload["_q_conv"]) == 0:
+            del sload["_q_conv"]
+        else:
+            printr(
+                "File {0}".format(basename(file))
+                + " has a deprecrated structure (key "
+                + "_q_conv removed in 0.9.30). Fixed this time, but regenerate "
+                + "database for faster loading."
+            )
+            if not "_q" in sload or len(sload["_q"]) == 0:
+                # only convolved quantities; just replace the dict
+                sload["_q"] = sload.pop("_q_conv")
+            else:
+                # here wavespaces arent necessarily the same
+                if sload["_q"]["wavespace"] == sload["_q_conv"]["wavespace"]:
+                    sload["_q"].update(sload.pop("_q_conv"))
+                else:
+                    # we need to interpolate
+                    raise NotImplementedError
+
+        fixed = True
+
     try:
         sload["conditions"]["waveunit"]
     except KeyError as err:
@@ -2087,7 +2103,9 @@ def find_duplicates(self, columns=None):
         onlyDuplicatedFiles = dg[dg == True]
         return onlyDuplicatedFiles
 
-    def add(self, spectrum, store_name=None, if_exists_then="increment", **kwargs):
+    def add(
+        self, spectrum: Spectrum, store_name=None, if_exists_then="increment", **kwargs
+    ):
         """Add Spectrum to database, whether it's a
         :py:class:`~radis.spectrum.spectrum.Spectrum` object or a file that
         stores one. Check it's not in database already.
@@ -2179,7 +2197,7 @@ def add(self, spectrum, store_name=None, if_exists_then="increment", **kwargs):
 
         # First, store the spectrum in a file
         # ... input is a Spectrum. Store it in database and load it from there
-        if is_spectrum(spectrum):
+        if isinstance(spectrum, Spectrum):
             # add defaults
             if store_name == None:
                 if not "add_info" in kwargs:
diff --git a/radis/tools/fitting.py b/radis/tools/fitting.py
index b547b24ab..36c2330cf 100644
--- a/radis/tools/fitting.py
+++ b/radis/tools/fitting.py
@@ -283,7 +283,7 @@ def print_fit_values(fit_values):
             ]
         )
 
-    def generate_spectrum(fit_values):
+    def generate_spectrum(fit_values) -> Spectrum:
 
         # Generate dictionary
         inputs = fit_parameters.copy()
@@ -549,7 +549,7 @@ def cost_and_plot_function(fit_values):
             )
         ),
         "transmittance_noslit",
-        waveunit="cm-1",
+        wunit="cm-1",
         unit="",
         delimiter=",",
         name="Klarenaar 2017",
diff --git a/radis/tools/slit.py b/radis/tools/slit.py
index 5cda16044..dfd6ae8c9 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 evenly_distributed
+from radis.misc.arrays import anynan, evenly_distributed
 from radis.misc.basics import is_float
 from radis.misc.debug import printdbg
 from radis.misc.signal import resample_even
@@ -275,7 +275,7 @@ def check_input_gen():
                 center=center_wavespace,
                 bplot=plot,
                 norm_by=norm_by,
-                waveunit=return_unit,
+                wunit=return_unit,
                 scale=scale_slit,
                 *args,
                 **kwargs
@@ -291,7 +291,7 @@ def check_input_gen():
                 center=center_wavespace,
                 bplot=plot,
                 norm_by=norm_by,
-                waveunit=return_unit,
+                wunit=return_unit,
                 scale=scale_slit,
                 *args,
                 **kwargs
@@ -358,7 +358,7 @@ def check_input_gen():
             center=center_wavespace,
             bplot=plot,
             norm_by=norm_by,
-            waveunit=return_unit,
+            wunit=return_unit,
             scale=scale_slit,
             *args,
             **kwargs
@@ -379,7 +379,7 @@ def check_input_gen():
         wslit, Islit = import_experimental_slit(
             slit_function,
             norm_by=norm_by,  # norm is done later anyway
-            waveunit=unit,
+            wunit=unit,
             verbose=verbose,
             auto_crop=auto_recenter_crop,
             auto_recenter=auto_recenter_crop,
@@ -461,7 +461,7 @@ def check_input_gen():
 
         if plot:  # (plot after resampling / renormalizing)
             # Plot slit
-            plot_slit(wslit, Islit, waveunit=return_unit, Iunit=Iunit)
+            plot_slit(wslit, Islit, wunit=return_unit, Iunit=Iunit)
 
     else:
         raise TypeError(
@@ -496,7 +496,8 @@ def convolve_with_slit(
     bplot=False,
     verbose=True,
     assert_evenly_spaced=True,
-    waveunit="",
+    wunit="",
+    waveunit=None,
 ):
     """Convolves spectrum (w,I) with instrumental slit function (w_slit, I_slit)
     Returns a convolved spectrum on a valid range.
@@ -547,7 +548,7 @@ def convolve_with_slit(
         ``assert_evenly_spaced=True``, then we check this is the case, and resample
         ``w`` and ``I`` if needed. Recommended, but it takes some time.
 
-    waveunit: ``'nm'``, ``'cm-1'``
+    wunit: ``'nm'``, ``'cm-1'``
         just for printing messages. However, ``w`` and ``w_slit`` should be in the
         same wavespace.
 
@@ -581,6 +582,14 @@ def convolve_with_slit(
     # 1. Check input
     # --------------
 
+    # Deprecated input:
+    if waveunit is not None:
+        warn(
+            "`waveunit=` parameter in convolve_with_slit is now named `wunit=`",
+            DeprecationWarning,
+        )
+        wunit = waveunit
+
     # Assert slit function is thin enough
     try:
         assert abs(w[-1] - w[0]) > abs(w_slit[-1] - w_slit[0])
@@ -625,7 +634,7 @@ def convolve_with_slit(
             print(w_slit)
             print(I_slit)
             print("Check figure")
-            plot_slit(w_slit, I_slit, waveunit=waveunit)
+            plot_slit(w_slit, I_slit, wunit=wunit)
             raise
 
         # can be a bug here if wstep has the wrong sign.
@@ -640,17 +649,17 @@ def convolve_with_slit(
     # -------------
 
     # Check no nan
-    if np.isnan(I_slit).sum() > 0:
+    if anynan(I_slit):
         raise ValueError("Slit has nan value")
 
     # check slit is positive
     if not (I_slit >= 0).all():
-        plot_slit(w_slit, I_slit, waveunit="")
+        plot_slit(w_slit, I_slit, wunit="")
         raise ValueError("Slit is partially negative. Check Figure")
 
     # Plot slit if asked for
     if bplot:
-        plot_slit(w_slit, I_slit, waveunit="")
+        plot_slit(w_slit, I_slit, wunit="")
 
     # 4. Convolve!
     # --------------
@@ -949,6 +958,7 @@ def remove_boundary(
 
     crop_left, crop_right: int
         number of points to discard on each side if ``mode='crop'``
+        Values are replaced with ``nan``
 
     Returns
     -------
@@ -956,6 +966,13 @@ def remove_boundary(
     w_conv, I_conv: numpy arrays
         cropped waverange and quantity
 
+    Notes
+    -----
+
+    .. note::
+        new in 0.9.30 : remove_boundary replaces off-range values with ``nan``
+        but keeps the same array size.
+
     See Also
     --------
 
@@ -969,25 +986,30 @@ def remove_boundary(
 
     """
 
-    # Remove boundary effects with the x-axis changed accordingly
+    # Remove boundary effects by adding nans; keep the same x-axis
     if mode == "valid":
         la = min(len_I, len_I_slit_interp)
         a = int((la - 1) / 2)
         b = int((la) / 2)
-        I_conv = I_conv[a:-b]
-        w_conv = w[a:-b]
+        # I_conv = I[a:-b]    # former version : we would change the array size
+        # w_conv = w[a:-b]
+        I_conv[:a] = np.nan
+        I_conv[-b:] = np.nan
+        w_conv = w
     elif mode == "same":
         I_conv = I_conv
         w_conv = w
     elif mode == "crop":
-        l = len(I_conv)
         if crop_right == 0:
             _crop_right = None
         else:
             _crop_right = -crop_right
-        I_conv = I_conv[crop_left:_crop_right]
-        w_conv = w[crop_left:_crop_right]
-        assert len(I_conv) == l - crop_left - crop_right
+        # l = len(I_conv)
+        # I_conv = I_conv[crop_left:_crop_right]
+        # w_conv = w[crop_left:_crop_right]
+        # assert len(I_conv) == l - crop_left - crop_right
+        I_conv[:crop_left] = np.nan
+        I_conv[_crop_right:] = np.nan
     else:
         raise ValueError("Unexpected mode: {0}".format(mode))
 
@@ -997,12 +1019,13 @@ def remove_boundary(
 def plot_slit(
     w,
     I=None,
-    waveunit="",
+    wunit="",
     plot_unit="same",
     Iunit=None,
     warnings=True,
     ls="-",
     title=None,
+    waveunit=None,
 ):
     """Plot slit, calculate and display FWHM, and calculate effective FWHM.
     FWHM is calculated from the limits of the range above the half width,
@@ -1046,6 +1069,14 @@ def plot_slit(
 
     """
 
+    # Deprecated input:
+    if waveunit is not None:
+        warn(
+            "`waveunit=` parameter in convolve_with_slit is now named `wunit=`",
+            DeprecationWarning,
+        )
+        wunit = waveunit
+
     import matplotlib.pyplot as plt
 
     from radis.misc.plot import set_style
@@ -1053,7 +1084,7 @@ def plot_slit(
     set_style()
 
     try:
-        from radis.plot.toolbar import add_tools  # TODO: move in publib
+        from radis.plot.toolbar import add_tools
 
         add_tools()  # includes a Ruler to measure slit
     except:
@@ -1063,26 +1094,26 @@ def plot_slit(
     if isinstance(w, str) and I is None:
         w, I = np.loadtxt(w).T
     assert len(w) == len(I)
-    if np.isnan(I).sum() > 0:
+    if anynan(I):
         warn("Slit function has nans")
         w = w[~np.isnan(I)]
         I = I[~np.isnan(I)]
     assert len(I) > 0
 
     # cast units
-    waveunit = cast_waveunit(waveunit, force_match=False)
+    wunit = cast_waveunit(wunit, force_match=False)
     plot_unit = cast_waveunit(plot_unit, force_match=False)
     if plot_unit == "same":
-        plot_unit = waveunit
+        plot_unit = wunit
 
     # Convert wavespace unit if needed
-    elif waveunit == "cm-1" and plot_unit == "nm":  # wavelength > wavenumber
+    elif wunit == "cm-1" and plot_unit == "nm":  # wavelength > wavenumber
         w = cm2nm(w)
-        waveunit = "nm"
-    elif waveunit == "nm" and plot_unit == "cm-1":  # wavenumber > wavelength
+        wunit = "nm"
+    elif wunit == "nm" and plot_unit == "cm-1":  # wavenumber > wavelength
         w = nm2cm(w)
-        waveunit = "cm-1"
-    elif waveunit == plot_unit:  # same units
+        wunit = "cm-1"
+    elif wunit == plot_unit:  # same units
         pass
     elif plot_unit == "":  # ???
         pass
@@ -1144,7 +1175,7 @@ def plot_slit(
             warn(
                 "Slit function doesnt seem centered: center measured with FWHM"
                 + " is not the array center (shift: {0:.3f}{1}): This can induce offsets!".format(
-                    abs(w[(xmin + xmax) // 2] - w[len(w) // 2]), waveunit
+                    abs(w[(xmin + xmax) // 2] - w[len(w) // 2]), wunit
                 )
             )
 
@@ -1243,12 +1274,13 @@ def import_experimental_slit(
     slit,
     norm_by="area",
     bplot=False,
-    waveunit="nm",
+    wunit="nm",
     auto_recenter=True,
     auto_crop=True,
     center=None,
     scale=1,
     verbose=True,
+    waveunit=None,
 ):
     """Import instrumental slit function and normalize it
 
@@ -1271,7 +1303,7 @@ def import_experimental_slit(
     bplot: boolean
         plot normalized slit function (for debugging). Default ``False``
 
-    waveunit: ``'nm'``, ``'cm-1'``
+    wunit: ``'nm'``, ``'cm-1'``
         used for plot only. Slit function is generated assuming you use the
         correct wavespace. No conversions are made here. Default ``'nm'``
 
@@ -1314,6 +1346,14 @@ def import_experimental_slit(
 
     """
 
+    # Deprecated input:
+    if waveunit is not None:
+        warn(
+            "`waveunit=` parameter in convolve_with_slit is now named `wunit=`",
+            DeprecationWarning,
+        )
+        wunit = waveunit
+
     # import
     if isinstance(slit, str):
         try:
@@ -1372,7 +1412,7 @@ def import_experimental_slit(
 
     # Plot slit
     if bplot:
-        plot_slit(w_slit, I_slit, waveunit=waveunit, Iunit=Iunit)
+        plot_slit(w_slit, I_slit, wunit=wunit, Iunit=Iunit)
 
     return w_slit, I_slit
 
@@ -1383,9 +1423,10 @@ def triangular_slit(
     center=0,
     norm_by="area",
     bplot=False,
-    waveunit="",
+    wunit="",
     scale=1,
     footerspacing=0,
+    waveunit=None,  # Deprecated
 ):
     r""" Generate normalized slit function
 
@@ -1409,7 +1450,7 @@ def triangular_slit(
     bplot: boolean
         plot normalized slit function (for debugging). Default ``False``
 
-    wavespace: '', 'nm', 'cm-1'
+    wunit: '', 'nm', 'cm-1'
         used for plot only. Slit function is generated assuming you use the
         correct wavespace. No conversions are made here.
 
@@ -1452,6 +1493,14 @@ def triangular_slit(
 
     """
 
+    # Deprecated input:
+    if waveunit is not None:
+        warn(
+            "`waveunit=` parameter in convolve_with_slit is now named `wunit=`",
+            DeprecationWarning,
+        )
+        wunit = waveunit
+
     # Build first half
     slope = 1 / (FWHM - wstep)
     a = int(FWHM / wstep)
@@ -1471,7 +1520,7 @@ def triangular_slit(
     # Normalize
     if norm_by == "area":  # normalize by the area
         I /= np.trapz(I, x=w)
-        Iunit = "1/{0}".format(waveunit)
+        Iunit = "1/{0}".format(wunit)
     elif norm_by == "max":  # set maximum to 1
         I /= np.max(I)
         Iunit = ""
@@ -1487,7 +1536,7 @@ def triangular_slit(
 
     # Plot slit
     if bplot:
-        plot_slit(w, I, waveunit=waveunit, Iunit=Iunit)
+        plot_slit(w, I, wunit=wunit, Iunit=Iunit)
 
     return w, I
 
@@ -1502,42 +1551,37 @@ def trapezoidal_slit(
     center=0,
     norm_by="area",
     bplot=False,
-    waveunit="",
+    wunit="",
     scale=1,
     footerspacing=0,
+    waveunit=None,
 ):
     r""" Build a trapezoidal slit. Remember that FWHM = (top + base) / 2
 
 
     Parameters
     ----------
-
     top: (nm)
         top of the trapeze
-
     base: (nm)
         base of the trapeze
-
     wstep: (nm)
         wavelength step
 
+    Other Parameters
+    ----------------
     center: (nm)
         center wavelength for the wavelength axs of the slit function
-
     norm_by: ``'area'``, ``'max'``
         normalisation type. ``'area'`` conserves energy. ``'max'`` is what is
         done in Specair and changes units. Default ``'area'``
-
     bplot: boolean
         plot normalized slit function (for debugging). Default ``False``
-
-    wavespace: '', 'nm', 'cm-1'
+    waveunit: '', 'nm', 'cm-1'
         used for plot only. Slit function is generated assuming you use the
         correct wavespace. No conversions are made here.
-
     scale: float
         multiply slit by an arbitrary factor. Default 1.
-
     footerspacing: int
         spacing (footer) on left and right. Default 10.
 
@@ -1612,7 +1656,7 @@ def trapezoidal_slit(
     # Normalize
     if norm_by == "area":  # normalize by the area
         I /= np.trapz(I, x=w)
-        Iunit = "1/{0}".format(waveunit)
+        Iunit = "1/{0}".format(wunit)
     elif norm_by == "max":  # set maximum to 1
         I /= np.max(I)
         Iunit = ""
@@ -1628,7 +1672,7 @@ def trapezoidal_slit(
 
     # Plot slit
     if bplot:
-        plot_slit(w, I, waveunit=waveunit, Iunit=Iunit)
+        plot_slit(w, I, wunit=wunit, Iunit=Iunit)
 
     return w, I
 
@@ -1639,50 +1683,44 @@ def gaussian_slit(
     center=0,
     norm_by="area",
     bplot=False,
-    waveunit="",
+    wunit="",
     calc_range=4,
     scale=1,
     footerspacing=0,
+    waveunit=None,  # Deprecated
 ):
     r""" Generate normalized slit function
 
 
     Parameters
     ----------
-
     FWHM: (nm)
         full-width at half maximum
-
     wstep: (nm)
         wavelength step
 
+    Other Parameters
+    ----------------
     center: (nm)
         center wavelength for the wavelength axs of the slit function
-
     norm_by: ``'area'``, ``'max'``
         normalisation type. ``'area'`` conserves energy. ``'max'`` is what is
         done in Specair and changes units. Default ``'area'``
-
     bplot: boolean
         plot normalized slit function (for debugging). Default ``False``
-
-    wavespace: '', 'nm', 'cm-1'
+    wunit: '', 'nm', 'cm-1'
         used for plot only. Slit function is generated assuming you use the
         correct wavespace. No conversions are made here.
-
     calc_range: float   (number of sigma)
         function broadening range expressed in number of standard deviation
-
     scale: float
         multiply slit by an arbitrary factor. Default 1.
-
     footerspacing: int
         spacing (footer) on left and right. Default 10
 
 
     Returns
     -------
-
     w, I: numpy arrays
         Slit function of shape::
 
@@ -1732,7 +1770,7 @@ def gaussian_slit(
     # Normalize
     if norm_by == "area":  # normalize by the area
         I /= np.trapz(I, x=w)
-        Iunit = "1/{0}".format(waveunit)
+        Iunit = "1/{0}".format(wunit)
     elif norm_by == "max":  # set maximum to 1
         I /= np.max(I)
         Iunit = ""
@@ -1748,7 +1786,7 @@ def gaussian_slit(
 
     # Plot slit
     if bplot:
-        plot_slit(w, I, waveunit=waveunit, Iunit=Iunit)
+        plot_slit(w, I, wunit=wunit, Iunit=Iunit)
 
     return w, I
 


Test patch:
diff --git a/radis/test/benchmark/radis_vs_hapi_CH4_full_spectrum.py b/radis/test/benchmark/radis_vs_hapi_CH4_full_spectrum.py
index 83b15b564..0ceb7fcc9 100644
--- a/radis/test/benchmark/radis_vs_hapi_CH4_full_spectrum.py
+++ b/radis/test/benchmark/radis_vs_hapi_CH4_full_spectrum.py
@@ -104,7 +104,7 @@ def calc_hapi():
     t0 = time() - t0
     print(("Calculated with HAPI in {0:.2f}s".format(t0)))
 
-    s_hapi = Spectrum.from_array(nu, coef, "abscoeff", waveunit="cm-1", unit="cm-1")
+    s_hapi = Spectrum.from_array(nu, coef, "abscoeff", wunit="cm-1", unit="cm-1")
     s_hapi.name = "HAPI ({0:.1f}s)".format(t0)
 
     plot_diff(s_hapi, s, "abscoeff")
diff --git a/radis/test/lbl/test_calc.py b/radis/test/lbl/test_calc.py
index 023e68aa0..3bd916d44 100644
--- a/radis/test/lbl/test_calc.py
+++ b/radis/test/lbl/test_calc.py
@@ -216,8 +216,10 @@ def test_calc_spectrum(verbose=True, plot=True, warnings=True, *args, **kwargs):
         plt.plot(w_ref, I_ref, "or", label="ref")
         plt.legend()
 
-    assert np.allclose(w[::100], w_ref, atol=1e-6)
-    assert np.allclose(I[::100], I_ref, atol=1e-6)
+    # [71:-71] because of the shift introduced in 0.9.30 where convolved
+    # arrays are not cropped; but np.nan values are added
+    assert np.allclose(w[71:-71][::100], w_ref, atol=1e-6)
+    assert np.allclose(I[71:-71][::100], I_ref, atol=1e-6)
 
     return True
 
@@ -284,7 +286,9 @@ def test_calc_spectrum_overpopulations(
         s.plot(wunit="nm")
 
     w, I = s.get("radiance", wunit="nm")
-    w_ref = w[::100]
+    # [71:-71] because of the shift introduced in 0.9.30 where convolved
+    # arrays are not cropped; but np.nan values are added
+    w_ref = w[71:-71][::100]
     # Compare against hardcoded results (neq 0.9.22, 28/06/18)
     #        I_ref = np.array([0.61826008, 0.65598262, 0.79760003, 0.7958013 , 0.5792486 ,
     #                          0.56727691, 0.60361258, 0.51549598, 0.51012651, 0.47133131,
@@ -355,7 +359,9 @@ def test_calc_spectrum_overpopulations(
         plt.legend()
         s.plot_populations()
 
-    assert np.allclose(I[::100], I_ref, atol=1e-6)
+    # [71:-71] because of the shift introduced in 0.9.30 where convolved
+    # arrays are not cropped; but np.nan values are added
+    assert np.allclose(I[71:-71][::100], I_ref, atol=1e-6)
 
     if verbose:
         printm("Test overpopulations: OK")
diff --git a/radis/test/misc/test_arrays.py b/radis/test/misc/test_arrays.py
index eb3958fa0..edb2ca204 100644
--- a/radis/test/misc/test_arrays.py
+++ b/radis/test/misc/test_arrays.py
@@ -19,8 +19,10 @@
     centered_diff,
     find_first,
     find_nearest,
+    first_nonnan_index,
     is_sorted,
     is_sorted_backward,
+    last_nonnan_index,
     logspace,
 )
 from radis.test.utils import setup_test_line_databases
@@ -37,6 +39,25 @@ def test_is_sorted(*args, **kwargs):
     assert not is_sorted(a[::-1])
 
 
+def test_nonnan_index(*args, **kwargs):
+
+    a = np.arange(1000) * 0.2
+
+    assert first_nonnan_index(a) == 0  # is None
+    assert last_nonnan_index(a) == 999  # is None
+
+    a[:10] = np.nan
+    a[-60:] = np.nan
+
+    assert first_nonnan_index(a) == 10
+    assert last_nonnan_index(a) == 939
+
+    a[:] = np.nan
+
+    assert first_nonnan_index(a) == None
+    assert last_nonnan_index(a) == None  # len(a)-1
+
+
 def test_find_first(*args, **kwargs):
 
     a = np.arange(10)
@@ -243,6 +264,7 @@ def test_cython_add_at_spectra(*args, **kwargs):
 
 
 if __name__ == "__main__":
+    test_nonnan_index()
     import pytest
 
     pytest.main(["test_arrays.py", "-s"])  # -s for showing console output
diff --git a/radis/test/spectrum/test_operations.py b/radis/test/spectrum/test_operations.py
index 6211157d7..7616f18ba 100644
--- a/radis/test/spectrum/test_operations.py
+++ b/radis/test/spectrum/test_operations.py
@@ -174,19 +174,11 @@ def test_offset(verbose=True, plot=False, *args, **kwargs):
     s2 = offset(s, 10, "nm", name="offset_10nm")
     if plot:
         plot_diff(s, s2)
-    assert np.allclose(
-        s2.get_wavelength(which="convoluted"), s.get_wavelength(which="convoluted") + 10
-    )
-    assert np.allclose(
-        s2.get_wavelength(which="non_convoluted"),
-        s.get_wavelength(which="non_convoluted") + 10,
-    )
+    assert np.allclose(s2.get_wavelength(), s.get_wavelength() + 10)
 
     # Test inplace version
     s.offset(10, "nm")
-    assert np.allclose(
-        s2.get_wavelength(which="convoluted"), s.get_wavelength(which="convoluted")
-    )
+    assert np.allclose(s2.get_wavelength(), s.get_wavelength())
 
 
 @pytest.mark.fast
@@ -269,17 +261,17 @@ def test_dimensioned_operations(*args, **kwargs):
 
     # Test
     assert s.units["radiance"] == "mW/cm2/sr/nm"
-    Imax = s.get("radiance")[1].max()
+    Imax = s.get("radiance", trim_nan=True)[1].max()
 
     # add a baseline
     s += 0.1 * u.Unit("W/cm2/sr/nm")
 
-    assert np.isclose(s.get("radiance")[1].max(), Imax + 100)
+    assert np.isclose(s.get("radiance", 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")[1].max(), Imax)
+    assert np.isclose(s.get("radiance", trim_nan=True)[1].max(), Imax)
 
     # Test division
     # Example : a manual normalization
diff --git a/radis/test/spectrum/test_spectrum.py b/radis/test/spectrum/test_spectrum.py
index 7ccf2faf3..58696f1e6 100644
--- a/radis/test/spectrum/test_spectrum.py
+++ b/radis/test/spectrum/test_spectrum.py
@@ -32,6 +32,61 @@
 # %% Test routines
 
 
+def test_spectrum_creation_method(*args, **kwargs):
+    import pytest
+
+    w = np.linspace(300, 700, 1000)
+    k = (np.random.rand(1000) + 0.3) ** 10
+    T = np.exp(-k * 1)
+
+    # Good inputs:
+    # ... format of quantities :
+    Spectrum({"wavelength": w, "abscoeff": k, "transmittance_noslit": T}, wunit="nm")
+    Spectrum({"abscoeff": (w, k), "transmittance_noslit": (w, T)}, wunit="nm")
+
+    Spectrum({"wavelength": w, "abscoeff": k, "transmittance_noslit": T}, wunit="nm")
+
+    # ... units:
+    Spectrum({"wavespace": w, "abscoeff": k}, wunit="nm")
+    Spectrum({"wavelength": w, "abscoeff": k}, wunit="nm")
+
+    # Bad inputs:
+    with pytest.raises(AssertionError) as err:
+        Spectrum(
+            {"wavelength": w, "abscoeff": k, "transmittance_noslit": T}
+        )  # wunit not defined
+    with pytest.raises(AssertionError) as err:
+        Spectrum(
+            {"wavenumber": w, "abscoeff": k, "transmittance_noslit": T}
+        )  # wunit not defined
+
+    # ... wavespace defined multiple times
+    with pytest.raises(AssertionError) as err:
+        Spectrum({"wavelength": w, "wavenumber": 1e7 / w, "abscoeff": k}, wunit="nm")
+    with pytest.raises(AssertionError) as err:
+        Spectrum({"wavelength": w, "wavespace": 1e7 / w, "abscoeff": k}, wunit="nm")
+
+    # ... format of quantities :
+    with pytest.raises(AssertionError) as err:
+        Spectrum({"wavenumber": w, "abscoeff": np.hstack((k, k))})
+        assert "Input arrays should have the same length" in str(err)
+
+    # ... units badly defeined :
+    with pytest.raises(AssertionError) as err:
+        Spectrum({"wavespace": w, "abscoeff": k})
+        assert "waveunit ('nm', 'cm-1'?) has to be defined" in str(err)
+
+    with pytest.raises(AssertionError) as err:
+        Spectrum({"abscoeff": (w, k)})
+        assert "waveunit ('nm', 'cm-1'?) has to be defined" in str(err)
+
+    with pytest.raises(AssertionError):
+        Spectrum({"wavenumber": w, "abscoeff": k}, wunit="nm")
+
+    with pytest.raises(AssertionError):
+        Spectrum({"wavelength": w, "abscoeff": k}, wunit="cm-1")
+
+
 def test_spectrum_get_methods(
     verbose=True, plot=True, close_plots=True, *args, **kwargs
 ):
@@ -96,6 +151,36 @@ def test_spectrum_get_methods(
         print(".. get_slit()")
 
 
+@pytest.mark.fast
+def test_trimming(verbose=True, *args, **kwargs):
+    """Test :py:meth:`radis.spectrum.spectrum.Spectrum.trim`"""
+
+    from radis.misc.arrays import count_nans
+
+    # trim both sides
+    w = np.linspace(300, 600)
+    T = np.ones_like(w)
+    T[:15] = np.nan
+    T[-10:] = np.nan
+    s = Spectrum.from_array(w, T, "transmittance_noslit", wunit="nm", unit="")
+
+    Nnans = count_nans(T)
+    assert len(s) == 50  # before trimming
+    s.trim()
+    assert len(s) == 50 - Nnans
+
+    # trim one side
+    w = np.linspace(300, 600)
+    T = np.ones_like(w)
+    T[:15] = np.nan
+    s = Spectrum.from_array(w, T, "transmittance_noslit", wunit="nm", unit="")
+
+    Nnans = count_nans(T)
+    assert len(s) == 50  # before trimming
+    s.trim()
+    assert len(s) == 50 - Nnans
+
+
 @pytest.mark.fast
 def test_copy(verbose=True, *args, **kwargs):
     """Test that a Spectrum is correctly copied
@@ -123,11 +208,8 @@ def test_copy(verbose=True, *args, **kwargs):
 
     # Test all quantities in detail
     for var in s._q.keys():
-        assert np.allclose(s._q[var], s2._q[var])
+        assert np.allclose(s._q[var], s2._q[var], equal_nan=True)
         assert not (s._q[var] is s2._q[var])
-    for var in s._q_conv.keys():
-        assert np.allclose(s._q_conv[var], s2._q_conv[var])
-        assert not (s._q_conv[var] is s2._q_conv[var])
 
     if verbose:
         print("Tested that s2 == s (but s2 is not s) after Spectrum copy")
@@ -258,7 +340,7 @@ def test_rescaling_function(verbose=True, *args, **kwargs):
     s = Spectrum.from_txt(
         getTestFile("calc_N2C_spectrum_Trot1200_Tvib3000.txt"),
         quantity="radiance_noslit",
-        waveunit="nm",
+        wunit="nm",
         unit="mW/cm2/sr/µm",  # Specair units: mW/cm2/sr/µm
         conditions={
             "Tvib": 3000,
@@ -340,6 +422,42 @@ def test_resampling_function(
     assert get_residual_integral(s, s3, "abscoeff", ignore_nan=True) < 1e-5
 
 
+def test_resampling_nan_function(verbose=True, *args, **kwargs):
+    """Test resampling functions, for Spectra with nans"""
+    from radis import get_residual
+    from radis.test.utils import getTestFile
+    from radis.tools.database import load_spec
+
+    plot = True
+
+    s = load_spec(getTestFile("CO_Tgas1500K_mole_fraction0.01.spec"), binary=True).crop(
+        2170, 2180, "cm-1"
+    )
+    s.rescale_path_length(10)  # cm
+    s.name = "original"
+
+    # We want to resample on a small range:
+    w_exp = s.get_wavenumber()[700:-700][::5]
+    # ... add a shift:
+    w_exp += np.diff(s.get_wavenumber())[0] * 2 / 3
+
+    s.update("transmittance_noslit")
+    s.apply_slit(3, "nm")
+
+    assert s.has_nan()
+
+    # s2.resample(w_exp)
+    s2 = s.resample(w_exp, inplace=False)
+    s2.name = "resampled"
+
+    assert get_residual(s2, s, "transmittance") < 1e-6
+
+    if plot:
+        from radis import plot_diff
+
+        plot_diff(s, s2, show_points=True)
+
+
 @pytest.mark.fast
 def test_noplot_different_quantities(*args, **kwargs):
     """Prevents User Errors: Ensures an error is raised if plotting different
@@ -451,6 +569,7 @@ def _run_testcases(
 
     # Test all Spectrum methods
     # -------------------------
+    test_spectrum_creation_method(*args, **kwargs)
     test_spectrum_get_methods(
         debug=debug,
         verbose=verbose,
@@ -460,6 +579,7 @@ def _run_testcases(
         **kwargs
     )
     test_copy(verbose=verbose, *args, **kwargs)
+    test_trimming(*args, **kwargs)
     test_populations(
         verbose=verbose, plot=plot, close_plots=close_plots, *args, **kwargs
     )
@@ -480,6 +600,7 @@ def _run_testcases(
     test_resampling_function(
         debug=debug, plot=plot, close_plots=close_plots, *args, **kwargs
     )
+    test_resampling_nan_function(verbose=verbose, *args, **kwargs)
 
     test_normalization(*args, **kwargs)
 
diff --git a/radis/test/tools/test_slit.py b/radis/test/tools/test_slit.py
index 47e702012..1f44fabb7 100644
--- a/radis/test/tools/test_slit.py
+++ b/radis/test/tools/test_slit.py
@@ -25,9 +25,10 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import pytest
-from numpy import abs, cos, linspace, pi, sqrt, tan, trapz
+from numpy import abs, cos, linspace, pi, sqrt, tan
 
 from radis.lbl.factory import SpectrumFactory
+from radis.misc.arrays import nantrapz
 from radis.misc.plot import fix_style, set_style
 from radis.misc.printer import printm
 from radis.phys.convert import dcm2dnm, dnm2dcm
@@ -75,7 +76,7 @@ def test_all_slit_shapes(
     s = Spectrum.from_txt(
         getTestFile("calc_N2C_spectrum_Trot1200_Tvib3000.txt"),
         quantity="radiance_noslit",
-        waveunit="nm",
+        wunit="nm",
         unit="mW/cm2/sr/µm",
     )
     wstep = np.diff(s.get_wavelength())[0]
@@ -154,7 +155,7 @@ def test_slit_unit_conversions_spectrum_in_cm(
 
         # Apply slit in nm this time
         s_nm = s_cm.copy()
-        w_cm = s_nm.get_wavenumber(which="non_convoluted")
+        w_cm = s_nm.get_wavenumber()
         slit_nm = dcm2dnm(slit_cm, w_cm[len(w_cm) // 2])
         s_nm.name = "Spec in cm-1, slit {0:.2f} nm".format(slit_nm)
         s_nm.apply_slit(slit_nm, unit="nm", shape=shape, mode="same")
@@ -196,7 +197,7 @@ def test_slit_unit_conversions_spectrum_in_nm(
     s_nm = Spectrum.from_txt(
         getTestFile("calc_N2C_spectrum_Trot1200_Tvib3000.txt"),
         quantity="radiance_noslit",
-        waveunit="nm",
+        wunit="nm",
         unit="mW/cm2/sr/µm",
         conditions={"self_absorption": False},
     )
@@ -225,7 +226,7 @@ def test_slit_unit_conversions_spectrum_in_nm(
 
         # Apply slit in nm this time
         s_cm = s_nm.copy()
-        w_nm = s_nm.get_wavelength(which="non_convoluted")
+        w_nm = s_nm.get_wavelength()
         slit_cm = dnm2dcm(slit_nm, w_nm[len(w_nm) // 2])
         s_cm.name = "Spec in nm, slit {0:.2f} cm-1".format(slit_cm)
         s_cm.apply_slit(slit_cm, unit="cm-1", shape=shape, mode="same")
@@ -325,8 +326,8 @@ def test_against_specair_convolution(
     # ... Test output is the same
     wu, Iu = s.get("radiance", Iunit="mW/cm2/sr")
 
-    As = np.trapz(Is, ws)  # Todo one day: replace with get_power() function
-    Au = np.trapz(Iu, wu)
+    As = nantrapz(Is, ws)  # Todo one day: replace with get_power() function
+    Au = nantrapz(Iu, wu)
     if verbose:
         print(
             (
@@ -338,7 +339,7 @@ def test_against_specair_convolution(
     assert np.isclose(As, Au, rtol=1e-2)
 
     # Test resampling
-    s.resample(linspace(376, 380.6, 3000), unit="nm", if_conflict_drop="convoluted")
+    s.resample(linspace(376, 380.6, 3000), unit="nm")
     s.apply_slit(slit_nm, norm_by="max")
     if plot:
         s.plot(
@@ -389,7 +390,7 @@ def test_normalisation_mode(plot=True, close_plots=True, verbose=True, *args, **
     if plot:
         ax.plot(w_max, I_max / FWHM, "r", label="(norm_by:max)/FWHM")
         ax.legend(loc="best")
-    assert np.allclose(I_area, I_max / FWHM)
+    assert np.allclose(I_area, I_max / FWHM, equal_nan=True)
     if verbose:
         print("equivalence of normalisation mode for spectrum in 'nm': OK")
 
@@ -413,7 +414,7 @@ def test_normalisation_mode(plot=True, close_plots=True, verbose=True, *args, **
     if plot:
         ax.plot(w_max, I_max / FWHM, "r", label="(norm_by:max)/FWHM")
         ax.legend(loc="best")
-    assert np.allclose(I_area, I_max / FWHM)
+    assert np.allclose(I_area, I_max / FWHM, equal_nan=True)
     if verbose:
         print("equivalence of normalisation mode for spectrum in 'cm-1': {0}: OK")
     assert is_homogeneous(s.units["radiance"], "mW/cm2/sr")
@@ -453,7 +454,7 @@ def test_slit_energy_conservation(
     P = s.get_power(unit="mW/cm2/sr")
     s.apply_slit(0.5, norm_by="area")
     w, I = s.get("radiance", wunit="nm", Iunit="mW/cm2/sr/nm")
-    Pc = abs(np.trapz(I, x=w))  # mW/cm2/sr
+    Pc = abs(np.trapz(I[~np.isnan(I)], x=w[~np.isnan(I)]))  # mW/cm2/sr
 
     b = np.isclose(P, Pc, 3e-2)
 
@@ -665,7 +666,9 @@ def slit_dispersion(w):
         # Compare 2 spectra
         s.plot(nfig="Linear dispersion effect", color="r", label="not corrected")
 
-    s2.apply_slit(slit_measured_632nm, slit_dispersion=slit_dispersion)
+    s2 = s2.apply_slit(
+        slit_measured_632nm, slit_dispersion=slit_dispersion, inplace=False
+    )
 
     if plot:
         s2.plot(nfig="same", color="k", label="corrected")
@@ -773,8 +776,8 @@ def test_resampling(rtol=1e-2, verbose=True, plot=True, warnings=True, *args, **
     w_nm_conv, T_nm_conv = sCO_nm.get("transmittance", wunit="cm-1")
 
     error = abs(
-        (trapz(1 - T_conv, w_conv) - trapz(1 - T_nm_conv, w_nm_conv))
-        / trapz(1 - T_nm_conv, w_nm_conv)
+        (nantrapz(1 - T_conv, w_conv) - nantrapz(1 - T_nm_conv, w_nm_conv))
+        / nantrapz(1 - T_nm_conv, w_nm_conv)
     )
 
     if verbose:
diff --git a/radis/test/utils.py b/radis/test/utils.py
index c94e8d12c..c2e42a9bd 100644
--- a/radis/test/utils.py
+++ b/radis/test/utils.py
@@ -117,7 +117,7 @@ def getValidationCase(file, force=False):
                 )
             ),
             "transmittance_noslit",
-            waveunit="cm-1",
+            wunit="cm-1",
             unit="",
             delimiter=",",
             name="Klarenaar 2017",
diff --git a/radis/test/validation/test_CO2_3Tvib_vs_klarenaar.py b/radis/test/validation/test_CO2_3Tvib_vs_klarenaar.py
index 1d55ceaf4..20b8302d1 100644
--- a/radis/test/validation/test_CO2_3Tvib_vs_klarenaar.py
+++ b/radis/test/validation/test_CO2_3Tvib_vs_klarenaar.py
@@ -62,7 +62,7 @@ def test_klarenaar_validation_case(
             )
         ),
         "transmittance_noslit",
-        waveunit="cm-1",
+        wunit="cm-1",
         unit="",
         delimiter=",",
         name="Klarenaar 2017",