Repository: radis/radis (language: python)

Issue:
Automatic wstep attribution in calc_spectrum()
### What is currently done 
Currently the default "wstep" value is 0.01, so when executing the code below:
```
spec_h2o = calc_spectrum(wavelength_min=4195,
                         wavelength_max=4215,
                         Tgas=300,
                         pressure=1.01325,
                         mole_fraction=0.001,
                         molecule='H2O',
                         isotope='1,2,3,4',
                         path_length=1,
                         databank='hitran',
                         )
```
An AccuracyWarning is raised: 
"Some lines are too narrow (FWHM ~ 0.027 cm⁻¹) for the current spectral grid (wstep=0.01). Please reduce wstep to below 0.009 cm⁻¹. You can also ignore by setting `warnings={'AccuracyWarning':'ignore'}` (if you know what you're doing!)"

### 🎯 Describe what you'd like
<!--- A clear and concise description of what you want to happen. -->
The option to choose "auto" for wstep in the calc_spectrum() function. (Maybe update the default behavior to it)

Since the warning already give the value that is pertinent for the current case, the wstep="auto" option should just take it (0.009 in the example).




---
<sub>*Tip : you can vote for [New Features](https://github.com/radis/radis/issues/79) on https://feathub.com/radis/radis*</sub>



Discussion / hints:
Yes !! It's already in #184 ,   but your description is better.

I'll close **the other one** as a duplicate ;) 

The only thing to choose is : **which number of points per line to enforce** ? 

To answer this question, eventually, we'll want to ensure a certain numerical error (note : eventually RADIS will have `performance/auto/precision` modes. The accepted error would be different for each. see the Discussion in https://github.com/radis/radis/discussions/197 )

1. For a perfect Voigt (numerical convolution), use different `wstep` and measure the error 
2. At runtime, multiply this by the linestrength to have an idea of the total error 
3. Determine the proper wstep at runtime ! 

Hey @erwanp, the  `number of points per line to enforce` will directly relate to the adequate `wstep` for the line having minimum FWHM such that for a particular `number of wstep`  the `minimum_length` is sufficient enough to get good accuracy. For the precisions mode, we are talking about having an error threshold for each mode.
Suppose we are talking about `auto` mode. What we will want to do is find the numerical convolution of it using the broadening method = 'convolve' instead of the `olivero_1977` Voigt FWHM approximation, and get the error? 
Not sure how we will get the error from the Gaussian and Lorentzian HWHM.
After that, we multiply it with linestrength `S` and get the total error for different `wstep`. Depending on the value having the least error we will select that `wstep`.

All the necessary checks can be added in the [_check_accuracy()](https://github.com/radis/radis/blob/180f63e4129c1307cd151f132f4342813c8ffbcd/radis/lbl/broadening.py#L876) function itself to recompute `wstep`. 
A simple approach will be manually defining the number of grid points per line for different modes like for 
`performance`  - `2 wstep /line`, for `auto` - `5 wstep/line`, etc. 
The method you have proposed will be more accurate per condition but we will have to make sure that all the checks aren't that computationally expensive.
Now about your idea, I have some doubts regarding the process. Depending on the `broadening method` selected we will have either Voigt HWHM or both Lorentzian and Gaussian HWHM. In default cases, we choose the `Voigt` method so sticking to that how we will calculate the error based on `wstep`?

Fix patch:
diff --git a/docs/lbl/lbl.rst b/docs/lbl/lbl.rst
index ed8e238f6..9b62cf04c 100644
--- a/docs/lbl/lbl.rst
+++ b/docs/lbl/lbl.rst
@@ -44,6 +44,7 @@ from the latest HITRAN database, and plot the transmittance: ::
         		molecule='CO',
         		mole_fraction=0.5,
         		isotope=1,
+            wstep=0.01,
             databank='hitran'   # or 'hitemp'
     		  	)
 	s.plot('transmittance_noslit')
@@ -62,7 +63,8 @@ transmittance: ::
         	Tgas=700,
         	path_length=0.1,
         	mole_fraction={'CO2':0.5, 'CO':0.5},
-        	isotope=1,
+        	wstep=0.01,
+          isotope=1,
     		)
 	s.plot('transmittance_noslit')
 
@@ -81,6 +83,7 @@ following commands give the same result: ::
             molecule=["CO2", "CO"],
             mole_fraction=1,
             isotope={"CO2": "1,2", "CO": "1,2,3"},
+            wstep=0.01,
             verbose=verbose,
       )
 
@@ -133,7 +136,6 @@ Be careful to be consistent and not to give partial or contradictory inputs. ::
         )
 
 
-
 Flow Chart
 ----------
 
@@ -207,10 +209,10 @@ before the spectrum calculation begins.
 # TODO: perform timing test to see how much time calculating log_2gs separately takes
 
 Currently, GPU-powered spectra calculations are supported only at thermal equilibrium
-and therefore, the method to calculate the spectra has been named :py:func:`~radis.lbl.calc.eq_spectrum_gpu`.
+and therefore, the method to calculate the spectra has been named :py:meth:`~radis.lbl.factory.SpectrumFactory.eq_spectrum_gpu`.
 In order to use this method to calculate the spectra, follow the same steps as in the
 case of a normal equilibrium spectra, and if using :py:func:`~radis.lbl.calc.calc_spectrum`
-function set the parameter `mode` to `gpu`, or use :py:func:`~radis.lbl.calc.eq_spectrum_gpu`
+function set the parameter `mode` to `gpu`, or use :py:meth:`~radis.lbl.factory.SpectrumFactory.eq_spectrum_gpu`
 
 Consider the following example which demonstrates the above information::
 
@@ -336,6 +338,7 @@ and :py:mod:`~astropy.units` ::
                          path_length=0.1 * u.m,
                          pressure=20 * u.mbar,
                          molecule='CO2',
+                         wstep = 0.01,
                          isotope='1,2',
                          cutoff=1e-25,              # cm/molecule
                          broadening_max_width=10,   # cm-1
@@ -688,6 +691,27 @@ for an example, which can be run with (you will need the CDSD-HITEMP database in
 
     pytest radis/test/lbl/test_broadening.py -m "test_abscoeff_continuum"
 
+Choose the right wavenumber grid
+--------------------------------
+
+``wstep`` determines the wavenumber grid's resolution. Smaller the value, higher the resolution and
+vice-versa. By default **radis** uses ``wstep=0.01``. You can manually set the ``wstep`` value
+in :py:func:`~radis.lbl.calc.calc_spectrum` and :py:class:`~radis.lbl.factory.SpectrumFactory`.
+To get more accurate result you can further reduce the value, and to increase the performance you can increase the value.
+
+Based on ``wstep``, it will determine the number of gridpoints per linewidth.
+To make sure that there are enough gridpoints, Radis will raise an Accuracy Warning :py:meth:`~radis.lbl.broadening.BroadenFactory._check_accuracy`
+if number of gridpoints are less than :py:data:`~radis.params.GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD` and raises an Accuracy Error
+if number of gridpoints are less than :py:data:`~radis.params.GRIDPOINTS_PER_LINEWIDTH_ERROR_THRESHOLD`.
+
+From ``0.9.30`` a new mode ``wstep='auto'`` has been added which directly computes the optimum value of ``wstep``
+ensuring both performance and accuracy. It is ensured that there are slightly more or less than :py:data:`~radis.params.GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD`
+points for each linewidth.
+
+.. note::
+    wstep = 'auto' is optimized for performances while ensuring accuracy,
+    but is still experimental in 0.9.30. Feedback welcome!
+
 
 Database loading
 ----------------
@@ -748,7 +772,7 @@ GPU accelerated spectrum calculation
 
 Apart from the parallelization method mentioned above, RADIS also supports CUDA-native parallel computation, specifically
 for lineshape calculation and broadening. To use these GPU-accelerated methods to compute the spectra, use either :py:func:`~radis.lbl.calc.calc_spectrum`
-function with parameter `mode` set to `gpu`, or :py:func:`~radis.lbl.calc.eq_spectrum_gpu`. In order to use these methods,
+function with parameter `mode` set to `gpu`, or :py:meth:`~radis.lbl.factory.SpectrumFactory.eq_spectrum_gpu`. In order to use these methods,
 ensure that your system has an Nvidia GPU with compute capability of atleast 3.0 and CUDA Toolkit 8.0 or above. Refer to
 :ref:`GPU Spectrum Calculation on RADIS <label_radis_gpu>` to see how to setup your system to run GPU accelerated spectrum
 calculation methods, examples and performance tests.
diff --git a/radis/__init__.py b/radis/__init__.py
index 333f7da58..a35a194d8 100644
--- a/radis/__init__.py
+++ b/radis/__init__.py
@@ -98,7 +98,14 @@
 
 
 # %% Global constants
-from .params import AUTO_UPDATE_SPEC, DEBUG_MODE, OLDEST_COMPATIBLE_VERSION, USE_CYTHON
+from .params import (
+    AUTO_UPDATE_SPEC,
+    DEBUG_MODE,
+    GRIDPOINTS_PER_LINEWIDTH_ERROR_THRESHOLD,
+    GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD,
+    OLDEST_COMPATIBLE_VERSION,
+    USE_CYTHON,
+)
 
 # %% Version
 
diff --git a/radis/lbl/bands.py b/radis/lbl/bands.py
index d1c578054..ff58be276 100644
--- a/radis/lbl/bands.py
+++ b/radis/lbl/bands.py
@@ -230,6 +230,9 @@ def eq_bands(
         # ... calculate broadening  HWHM
         self._calc_broadening_HWHM()
 
+        # ... generates all wstep related entities
+        self._generate_wavenumber_arrays()
+
         # ... find weak lines and calculate semi-continuum (optional)
         I_continuum = self.calculate_pseudo_continuum()
         if I_continuum:
@@ -551,6 +554,9 @@ def non_eq_bands(
         # ... calculate broadening  HWHM
         self._calc_broadening_HWHM()
 
+        # ... generates all wstep related entities
+        self._generate_wavenumber_arrays()
+
         # ... find weak lines and calculate semi-continuum (optional)
         I_continuum = self.calculate_pseudo_continuum()
         if I_continuum:
diff --git a/radis/lbl/base.py b/radis/lbl/base.py
index aa4d89b8f..748c3cd7d 100644
--- a/radis/lbl/base.py
+++ b/radis/lbl/base.py
@@ -2830,11 +2830,12 @@ def _cutoff_linestrength(self, cutoff=None):
         Nlines_cutoff = b.sum()
 
         # Estimate time gained
-        expected_broadening_time_gain = (
-            self._broadening_time_ruleofthumb
-            * Nlines_cutoff
-            * len(self.wbroad_centered)
-        )
+        # TODO: Add a better formula to estimate time gained during broadening process
+        """ Previous method used was:
+            expected_broadening_time_gain = (
+                self._broadening_time_ruleofthumb * Nlines_cutoff * len(self.wbroad_centered)
+            )
+        """
 
         # Estimate error being made:
         if self.warnings["LinestrengthCutoffWarning"] != "ignore":
@@ -2846,9 +2847,7 @@ def _cutoff_linestrength(self, cutoff=None):
                     "Discarded {0:.2f}% of lines (linestrength<{1}cm-1/(#.cm-2))".format(
                         Nlines_cutoff / len(df.S) * 100, cutoff
                     )
-                    + " Estimated error: {0:.2f}%, Expected time saved: {1:.1f}s".format(
-                        error, expected_broadening_time_gain
-                    )
+                    + " Estimated error: {0:.2f}%".format(error)
                 )
             if error > self.misc.warning_linestrength_cutoff:
                 self.warn(
@@ -2887,24 +2886,7 @@ def _cutoff_linestrength(self, cutoff=None):
 
         time_spent = time() - t0
         if self.verbose >= 2:
-            printg(
-                "Applied linestrength cutoff in {0:.1f}s (expected time saved ~ {1:.1f}s)".format(
-                    time_spent, expected_broadening_time_gain
-                )
-            )
-
-        # Raise a warning if we dont expect any performance improvement
-        if time_spent > 0.3 and time_spent > 3 * expected_broadening_time_gain:
-            self.warn(
-                "Your linestrength cutoff is too high. "
-                + "Time spent on applying cutoff "
-                + "({0:.1f}s) is much longer than expected gain ({1:.1f}s). ".format(
-                    time_spent, expected_broadening_time_gain
-                )
-                + "If the number of lines is low, you can consider "
-                + "using cutoff=0",
-                "PerformanceWarning",
-            )
+            printg("Applied linestrength cutoff in {0:.1f}s ".format(time_spent))
 
         return
 
@@ -3044,7 +3026,10 @@ def _check_inputs(self, mole_fraction, Tmax):
         # ... that wstep and broadening_max_width were not inadvertanly changed
         # ... (would have no effect as the waverange is calculated on SpectrumFactory
         # ... initialization)
-        assert self._wstep == self.params.wstep
+
+        # Checks there if there is change in wstep value if initial wstep != "auto"
+        if self.wstep != "auto":
+            assert self.wstep == self.params.wstep
         assert self._broadening_max_width == self.params.broadening_max_width
 
     def _get_parsum(self, molecule, iso, state):
diff --git a/radis/lbl/broadening.py b/radis/lbl/broadening.py
index 28c71bffe..6d02785cb 100644
--- a/radis/lbl/broadening.py
+++ b/radis/lbl/broadening.py
@@ -79,7 +79,11 @@
 
 # from radis.misc.warning import AccuracyError, AccuracyWarning
 from radis.misc.warning import reset_warnings
-from radis.params import USE_CYTHON
+from radis.params import (
+    GRIDPOINTS_PER_LINEWIDTH_ERROR_THRESHOLD,
+    GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD,
+    USE_CYTHON,
+)
 from radis.phys.constants import Na, c_CGS, k_b_CGS
 
 # %% Broadening functions
@@ -836,13 +840,28 @@ def _calc_broadening_HWHM(self):
                 )
             )
 
-        # AccuracyWarning. Check there are enough gridpoints per line.
-        self._check_accuracy(df, self.params.wstep)
-
         if self.verbose >= 2:
             printg("Calculated broadening HWHM in {0:.2f}s".format(time() - t0))
 
-    def _check_accuracy(self, df, wstep):
+    def _calc_min_width(self, df):
+        """Calculates the minimum FWHW of the lines
+        and stores in self.min_width
+        """
+        if "hwhm_voigt" in df:
+            min_width = 2 * df.hwhm_voigt.min()
+        else:
+            min_lorentz_fwhm = 2 * df.hwhm_lorentz.min()
+            min_gauss_fwhm = 2 * df.hwhm_gauss.min()
+            # We take the max of both. Note: could also have used
+            # Olivero1977 to get the Voigt-equivlaent width of all lines,
+            # but it's quite expensive to compute
+            min_width = max(min_lorentz_fwhm, min_gauss_fwhm)
+
+        self.min_width = min_width
+
+        return
+
+    def _check_accuracy(self, wstep):
         """Check there are enough gridpoints per line.
 
         Raises
@@ -875,33 +894,25 @@ def _check_accuracy(self, df, wstep):
         # TODO: thresholds depend whether we're computing Transmittance/optically thin emission,
         # for a homogeneous slab, or self-absorbed radiance combined with other slabs.
 
-        if "hwhm_voigt" in df:
-            min_width = 2 * df.hwhm_voigt.min()
-        else:
-            min_lorentz_fwhm = 2 * df.hwhm_lorentz.min()
-            min_gauss_fwhm = 2 * df.hwhm_gauss.min()
-            # We take the max of both. Note: could also have used
-            # Olivero1977 to get the Voigt-equivlaent width of all lines,
-            # but it's quite expensive to compute
-            min_width = max(min_lorentz_fwhm, min_gauss_fwhm)
+        min_width = self.min_width
 
-        WARN_THRESHOLD = 3
-        ERROR_TRESHOLD = 1
-        if wstep > min_width / ERROR_TRESHOLD:
+        if wstep > min_width / GRIDPOINTS_PER_LINEWIDTH_ERROR_THRESHOLD:
             self.warn(
                 f"Some lines are too narrow (FWHM ~ {min_width:.2g} cm⁻¹) for "
                 + f"the current spectral grid (wstep={wstep}). Please reduce "
-                + f"wstep to (at least) below {min_width/ERROR_TRESHOLD:.2g} cm⁻¹ "
-                + f"or (suggested) {min_width/WARN_THRESHOLD:.2g} cm⁻¹. "
+                + f"wstep to (at least) below {min_width/GRIDPOINTS_PER_LINEWIDTH_ERROR_THRESHOLD:.2g} cm⁻¹ "
+                + f"or (suggested) {min_width/GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD:.2g} cm⁻¹. "
+                + "You can use wstep='auto' to get the optimal spectral grid value. "
                 + "You can also ignore by setting `warnings={'AccuracyError':'ignore'}` "
                 + "(if you know what you're doing!)",
                 "AccuracyError",
             )
-        elif wstep > min_width / WARN_THRESHOLD:
+        elif wstep > min_width / GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD:
             self.warn(
                 f"Some lines are too narrow (FWHM ~ {min_width:.2g} cm⁻¹) for "
                 + f"the current spectral grid (wstep={wstep}). Please reduce "
-                + f"wstep to below {min_width/WARN_THRESHOLD:.2g} cm⁻¹. "
+                + f"wstep to below {min_width/GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD:.2g} cm⁻¹. "
+                + "You can use wstep='auto' to get the optimal spectral grid value. "
                 + "You can also ignore by setting `warnings={'AccuracyWarning':'ignore'}` "
                 + "(if you know what you're doing!)",
                 "AccuracyWarning",
@@ -909,6 +920,8 @@ def _check_accuracy(self, df, wstep):
         else:
             pass
 
+        return
+
     def _add_voigt_broadening_HWHM(self, df, pressure_atm, mole_fraction, Tgas, Tref):
         """Update dataframe with Voigt HWHM.
 
diff --git a/radis/lbl/calc.py b/radis/lbl/calc.py
index a45205f1c..161b45bf7 100644
--- a/radis/lbl/calc.py
+++ b/radis/lbl/calc.py
@@ -150,8 +150,15 @@ def calc_spectrum(
     medium: ``'air'``, ``'vacuum'``
         propagating medium when giving inputs with ``'wavenum_min'``, ``'wavenum_max'``.
         Does not change anything when giving inputs in wavenumber. Default ``'air'``​.
-    wstep: float (:math:`cm^{-1}`)
-        Spacing of calculated spectrum. Default ``0.01 cm-1``.
+    wstep: float (:math:`cm^{-1}`)  or `'auto'`
+        Resolution of wavenumber grid. Default ``0.01`` cm-1.
+        If `'auto'`, it is ensured that there
+        are slightly more or less than :py:data:`~radis.params.GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD`
+        points for each linewidth.
+
+        .. note::
+            wstep = 'auto' is optimized for performances while ensuring accuracy,
+            but is still experimental in 0.9.30. Feedback welcome!
     broadening_max_width: float (cm-1)
         Full width over which to compute the broadening. Large values will create
         a huge performance drop (scales as :math:`~broadening_max_width^2` without LDM)
@@ -280,7 +287,6 @@ def calc_spectrum(
 
     Cite
     ----
-
     RADIS is built on the shoulders of many state-of-the-art packages and databases. If using RADIS
     to compute spectra, make sure you cite all of them, for proper reproducibility and acknowledgement of
     the work ! See :ref:`How to cite? <label_cite>`
diff --git a/radis/lbl/factory.py b/radis/lbl/factory.py
index 7f0fe011c..ee628384b 100644
--- a/radis/lbl/factory.py
+++ b/radis/lbl/factory.py
@@ -90,9 +90,10 @@
 from radis.db.molparam import MolParams
 from radis.lbl.bands import BandFactory
 from radis.lbl.base import get_waverange
-from radis.misc.basics import flatten, is_float, list_if_float
+from radis.misc.basics import flatten, is_float, list_if_float, round_off
 from radis.misc.printer import printg
 from radis.misc.utils import Default
+from radis.params import GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD
 from radis.phys.constants import k_b
 from radis.phys.convert import conv2
 from radis.phys.units import convert_emi2nm, convert_rad2nm
@@ -166,8 +167,15 @@ class SpectrumFactory(BandFactory):
         The calculated spectral range is increased (by broadening_max_width/2
         on each side) to take into account overlaps from out-of-range lines.
         Default ``10`` cm-1.
-    wstep: float (cm-1)
-        Spacing of calculated spectrum. Default ``0.01`` cm-1
+    wstep: float (cm-1) or `'auto'`
+        Resolution of wavenumber grid. Default ``0.01`` cm-1.
+        If `'auto'`, it is ensured that there
+        are slightly more or less than :py:data:`~radis.params.GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD`
+        points for each linewidth.
+
+        .. note::
+            wstep = 'auto' is optimized for performances while ensuring accuracy,
+            but is still experimental in 0.9.30. Feedback welcome!
     cutoff: float (~ unit of Linestrength: cm-1/(#.cm-2))
         discard linestrengths that are lower that this, to reduce calculation
         times. ``1e-27`` is what is generally used to generate databases such as
@@ -411,21 +419,11 @@ def __init__(
             medium,
         )
 
-        # calculated range is broader than output waverange to take into account off-range line broadening
-        wavenumber, wavenumber_calc = _generate_wavenumber_range(
-            wavenum_min, wavenum_max, wstep, broadening_max_width
-        )
-        wbroad_centered = _generate_broadening_range(wstep, broadening_max_width)
         # Store broadening max width and wstep as hidden variable (to ensure they are not changed afterwards)
-        self._wstep = wstep
         self._broadening_max_width = broadening_max_width
 
-        # Get boolean array that extracts the reduced range `wavenumber` from `wavenumber_calc`
-        woutrange = np.in1d(wavenumber_calc, wavenumber, assume_unique=True)
-        self.wbroad_centered = wbroad_centered
-        self.wavenumber = wavenumber
-        self.wavenumber_calc = wavenumber_calc
-        self.woutrange = woutrange
+        # Storing inital value of wstep if wstep != "auto"
+        self.wstep = wstep
 
         # Init variables
         # --------------
@@ -482,8 +480,8 @@ def __init__(
         self.params.broadening_max_width = broadening_max_width  # line broadening
         self.misc.export_lines = export_lines
         self.misc.export_populations = export_populations
-        self.params.wavenum_min_calc = wavenumber_calc[0]
-        self.params.wavenum_max_calc = wavenumber_calc[-1]
+        self.params.wavenum_min_calc = wavenum_min - broadening_max_width / 2
+        self.params.wavenum_max_calc = wavenum_max + broadening_max_width / 2
 
         # if optimization is ``'simple'`` or ``'min-RMS'``, or None :
         # Adjust default values of broadening method :
@@ -569,7 +567,6 @@ def eq_spectrum(
 
         Parameters
         ----------
-
         Tgas: float or `~astropy.units.quantity.Quantity`
             Gas temperature (K)
         mole_fraction: float
@@ -585,7 +582,6 @@ def eq_spectrum(
 
         Returns
         -------
-
         s : Spectrum
             Returns a :class:`~radis.spectrum.spectrum.Spectrum` object
 
@@ -597,13 +593,11 @@ def eq_spectrum(
 
         References
         ----------
-
         .. [1] RADIS doc: `Spectrum how to? <https://radis.readthedocs.io/en/latest/spectrum/spectrum.html#label-spectrum>`__
 
 
         See Also
         --------
-
         :meth:`~radis.lbl.factory.SpectrumFactory.non_eq_spectrum`
         """
 
@@ -685,6 +679,9 @@ def eq_spectrum(
         # ... calculate broadening  HWHM
         self._calc_broadening_HWHM()
 
+        # ... generates all wstep related entities
+        self._generate_wavenumber_arrays()
+
         # ... find weak lines and calculate semi-continuum (optional)
         I_continuum = self.calculate_pseudo_continuum()
         # ... apply lineshape and get absorption coefficient
@@ -811,6 +808,11 @@ def eq_spectrum_gpu(
         """Generate a spectrum at equilibrium with calculation of lineshapes
         and broadening done on the GPU.
 
+        .. note::
+            This method requires CUDA compatible hardware to execute.
+            For more information on how to setup your system to run GPU-accelerated methods
+            using CUDA and Cython, check `GPU Spectrum Calculation on RADIS <https://radis.readthedocs.io/en/latest/lbl/gpu.html>`
+
         Parameters
         ----------
         Tgas: float or `~astropy.units.quantity.Quantity`
@@ -837,14 +839,8 @@ def eq_spectrum_gpu(
                 Or directly the :meth:`~radis.spectrum.spectrum.Spectrum.plot` method
                 to plot it. See [1]_ to get an overview of all Spectrum methods
 
-        Notes
-        -----
-
-        This method requires CUDA compatible hardware to execute. For more information on how to setup your system to run GPU-accelerated methods using CUDA and Cython, check `GPU Spectrum Calculation on RADIS <https://radis.readthedocs.io/en/latest/lbl/gpu.html>`
-
         See Also
         --------
-
         :meth:`~radis.lbl.factory.SpectrumFactory.eq_spectrum`
         """
 
@@ -937,8 +933,8 @@ def eq_spectrum_gpu(
         # generate the v_arr
         v_arr = np.arange(
             self.input.wavenum_min,
-            self.input.wavenum_max + self._wstep,
-            self._wstep,
+            self.input.wavenum_max + self.params.wstep,
+            self.params.wstep,
         )
 
         # load the data
@@ -1131,16 +1127,14 @@ def non_eq_spectrum(
         rot_distribution: ``'boltzmann'``
             rotational distribution
         overpopulation: dict, or ``None``
-            add overpopulation factors for given levels:
-
-            >>> {level:overpopulation_factor}
+            add overpopulation factors for given levels::
 
+                {level:overpopulation_factor}
         name: str
             output Spectrum name (useful in batch)
 
         Returns
         -------
-
         s : Spectrum
             Returns a :class:`~radis.spectrum.spectrum.Spectrum` object
 
@@ -1152,7 +1146,6 @@ def non_eq_spectrum(
 
         References
         ----------
-
         .. [1] RADIS doc: `Spectrum how to? <https://radis.readthedocs.io/en/latest/spectrum/spectrum.html#label-spectrum>`__
 
         See Also
@@ -1273,6 +1266,9 @@ def non_eq_spectrum(
         # ... calculate broadening  HWHM
         self._calc_broadening_HWHM()
 
+        # ... generates all wstep related entities
+        self._generate_wavenumber_arrays()
+
         # ... find weak lines and calculate semi-continuum (optional)
         k_continuum, j_continuum = self.calculate_pseudo_continuum(noneq=True)
 
@@ -1422,11 +1418,50 @@ def non_eq_spectrum(
 
         return s
 
+    def _generate_wavenumber_arrays(self):
+        """define wavenumber grid vectors
+
+        `SpectrumFactory.wavenumber` is the output spectral range and
+        ``SpectrumFactory.wavenumber_calc`` the spectral range used for calculation, that
+        includes neighbour lines within ``broadening_max_width`` distance."""
+        # calculates minimum FWHM of lines
+        self._calc_min_width(self.df1)
+
+        # Setting wstep to optimal value and rounding it to a degree 3
+        if self.wstep == "auto":
+            self.params.wstep = round_off(
+                self.min_width / GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD
+            )
+            self.warnings["AccuracyWarning"] = "ignore"
+
+        wavenumber, wavenumber_calc = _generate_wavenumber_range(
+            self.input.wavenum_min,
+            self.input.wavenum_max,
+            self.params.wstep,
+            self.params.broadening_max_width,
+        )
+        wbroad_centered = _generate_broadening_range(
+            self.params.wstep, self.params.broadening_max_width
+        )
+
+        # Get boolean array that extracts the reduced range `wavenumber` from `wavenumber_calc`
+        woutrange = np.in1d(wavenumber_calc, wavenumber, assume_unique=True)
+        self.wbroad_centered = wbroad_centered
+        self.wavenumber = wavenumber
+        self.wavenumber_calc = wavenumber_calc
+        self.woutrange = woutrange
+
+        # AccuracyWarning. Check there are enough gridpoints per line.
+        self._check_accuracy(self.params.wstep)
+
+        return
+
     def _get_log_2gs(self):
         """Returns log_2gs if it already exists in the dataframe, otherwise
         computes it using gamma_air."""
         df = self.df0
         # TODO: deal with the case of gamma_self [so we don't forget]
+        # TODO (refactor) : move into BaseFactory or BroadenFactory (parent classes)
 
         # if the column already exists, then return
         if "log_2gs" in df.columns:
@@ -1547,12 +1582,10 @@ def optically_thin_power(
 
         See Also
         --------
-
         :py:meth:`~radis.lbl.factory.SpectrumFactory.eq_spectrum`,
         :py:meth:`~radis.spectrum.spectrum.Spectrum.get_power`,
         :py:meth:`~radis.spectrum.spectrum.Spectrum.get_integral`
         """
-
         # Check inputs
 
         # ... temperatures
@@ -1668,35 +1701,29 @@ def optically_thin_power(
 
 
 def _generate_wavenumber_range(wavenum_min, wavenum_max, wstep, broadening_max_width):
-    """define waverange vectors, with ``wavenumber`` the ouput spectral range
+    """define waverange vectors, with ``wavenumber`` the output spectral range
     and ``wavenumber_calc`` the spectral range used for calculation, that
     includes neighbour lines within ``broadening_max_width`` distance.
 
     Parameters
     ----------
-
     wavenum_min, wavenum_max: float
         wavenumber range limits (cm-1)
-
     wstep: float
         wavenumber step (cm-1)
-
     broadening_max_width: float
         wavenumber full width of broadening calculation: used to define which
         neighbour lines shall be included in the calculation
 
     Returns
     -------
-
     wavenumber: numpy array
         an evenly spaced array between ``wavenum_min`` and ``wavenum_max`` with
         a spacing of ``wstep``
-
     wavenumber_calc: numpy array
         an evenly spaced array between ``wavenum_min-broadening_max_width/2`` and
         ``wavenum_max+broadening_max_width/2`` with a spacing of ``wstep``
     """
-
     assert wavenum_min < wavenum_max
 
     # Output range
@@ -1733,17 +1760,14 @@ def _generate_broadening_range(wstep, broadening_max_width):
 
     Parameters
     ----------
-
     wstep: float
         wavenumber step (cm-1)
-
     broadening_max_width: float
         wavenumber full width of broadening calculation: used to define which
         neighbour lines shall be included in the calculation
 
     Returns
     -------
-
     wbroad_centered: numpy array
         an evenly spaced array, of odd-parity length, centered on 0, and of width
         ``broadening_max_width``
diff --git a/radis/misc/basics.py b/radis/misc/basics.py
index f84088be9..9ff9570da 100644
--- a/radis/misc/basics.py
+++ b/radis/misc/basics.py
@@ -451,3 +451,13 @@ def to_str(a):
         return a.decode("utf-8")
     else:
         return a
+
+
+def round_off(n):
+    # Getting rounded off value (atleast order 3)
+    for i in range(0, 10):
+        val = round(n, 3 + i)
+        if val == 0:
+            continue
+        return val
+    return 0
diff --git a/radis/params.py b/radis/params.py
index 502b81999..74f4179ff 100644
--- a/radis/params.py
+++ b/radis/params.py
@@ -39,7 +39,7 @@
 
 See Also
 --------
-:func:`~radis.tools.database._update_to_latest_format`
+:py:func:`~radis.tools.database._update_to_latest_format`
 """
 
 OLDEST_COMPATIBLE_VERSION = "0.9.1"
@@ -57,3 +57,23 @@
 --------
 :py:func:`~radis.misc.arrays.add_at`
 """
+
+GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD = 3
+"""float: to determine the optimal value
+of wstep using minimum FWHM value of spectrum.
+Makes sure there are enough gridpoints per line.
+
+See Also
+--------
+:py:meth:`~radis.lbl.broadening.BroadenFactory._check_accuracy`
+"""
+
+GRIDPOINTS_PER_LINEWIDTH_ERROR_THRESHOLD = 1
+"""float: to determine the minimum feasible value
+of wstep using minimum FWHM value of spectrum.
+Makes sure there are enough gridpoints per line.
+
+See Also
+--------
+:py:meth:`~radis.lbl.broadening.BroadenFactory._check_accuracy`
+"""
diff --git a/radis/spectrum/spectrum.py b/radis/spectrum/spectrum.py
index 611474553..c4e07b3c3 100644
--- a/radis/spectrum/spectrum.py
+++ b/radis/spectrum/spectrum.py
@@ -214,7 +214,6 @@ class Spectrum(object):
 
     See Also
     --------
-
     :func:`~radis.spectrum.models.calculated_spectrum`,
     :func:`~radis.spectrum.models.transmittance_spectrum`,
     :func:`~radis.spectrum.models.experimental_spectrum`


Test patch:
diff --git a/radis/test/lbl/test_factory.py b/radis/test/lbl/test_factory.py
index c2e4db15c..6a0b840b0 100644
--- a/radis/test/lbl/test_factory.py
+++ b/radis/test/lbl/test_factory.py
@@ -555,6 +555,95 @@ def test_temperature_units_conversion(
     assert np.isclose(s.conditions["pressure_mbar"], 20)
 
 
+@pytest.mark.fast
+def test_wstep_auto_method_sf(verbose=True, plot=False, *args, **kwargs):
+    """Test to check that on computing several spectrum from the same Spectrum
+    Factory object we get the different wstep for each case using auto method"""
+
+    from radis.misc.basics import round_off
+    from radis.params import GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD
+
+    setup_test_line_databases()  # add HITRAN-CO-TEST in ~/radis.json if not there
+
+    sf = SpectrumFactory(
+        wavelength_min=4400,
+        wavelength_max=4800,
+        mole_fraction=0.01,
+        cutoff=1e-25,
+        wstep="auto",
+        isotope=[1],
+        db_use_cached=True,
+        self_absorption=True,
+        verbose=verbose,
+    )
+
+    sf.warnings["MissingSelfBroadeningWarning"] = "ignore"
+    sf.warnings["NegativeEnergiesWarning"] = "ignore"
+    sf.warnings["HighTemperatureWarning"] = "ignore"
+
+    sf.load_databank("HITRAN-CO-TEST")
+
+    s1 = sf.eq_spectrum(300, pressure=1)
+    assert sf.wstep == "auto"
+
+    wstep_calculated = s1.get_conditions()["wstep"]
+
+    # Checking computed wstep and expected wstep are equal
+    assert wstep_calculated == round_off(
+        sf.min_width / GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD
+    )
+
+    s2 = sf.eq_spectrum(300, pressure=0.2)
+    assert sf.wstep == "auto"
+
+    s3 = sf.eq_spectrum(300, pressure=0.001)
+    assert sf.wstep == "auto"
+
+    assert (
+        s1.get_conditions()["wstep"]
+        != s2.get_conditions()["wstep"]
+        != s3.get_conditions()["wstep"]
+    )
+
+
+@pytest.mark.fast
+def test_all_spectrum_using_wstep_auto(verbose=True, plot=False, *args, **kwargs):
+    """Checks all methods to calculate Spectrum works with "auto" mode"""
+    Tgas = 1000
+
+    sf = SpectrumFactory(
+        wavelength_min=4160,
+        wavelength_max=4220,
+        mole_fraction=1,
+        path_length=0.3,
+        cutoff=1e-23,
+        molecule="CO2",
+        isotope=1,
+        wstep="auto",
+        optimization=None,
+        verbose=verbose,
+    )
+    sf.warnings["MissingSelfBroadeningWarning"] = "ignore"
+    sf.warnings["NegativeEnergiesWarning"] = "ignore"
+    sf.warnings["HighTemperatureWarning"] = "ignore"
+    sf.fetch_databank("hitran")
+
+    sf.eq_spectrum(Tgas)
+    wstep_1 = sf.params.wstep
+
+    sf.non_eq_spectrum(Tvib=Tgas, Trot=Tgas)
+    wstep_2 = sf.params.wstep
+
+    sf.eq_bands(Tgas)
+    wstep_3 = sf.params.wstep
+
+    sf.non_eq_bands(Tvib=Tgas, Trot=Tgas)
+    wstep_4 = sf.params.wstep
+
+    assert wstep_1 == wstep_2 == wstep_3 == wstep_4
+    assert sf.wstep == "auto"
+
+
 # --------------------------
 if __name__ == "__main__":
 
diff --git a/radis/test/misc/test_basics.py b/radis/test/misc/test_basics.py
index 7f2b5ad68..b14d378fd 100644
--- a/radis/test/misc/test_basics.py
+++ b/radis/test/misc/test_basics.py
@@ -14,3 +14,12 @@ def test_flatten(*args, **kwargs):
     from radis.misc.basics import flatten
 
     assert flatten(*[600, 300, (600, 600, 700)]) == [600, 300, 600, 600, 700]
+
+
+@pytest.mark.fast
+def test_round_off(*args, **kwargs):
+    from radis.misc.basics import round_off
+
+    assert round_off(0.024231242) == 0.024
+    assert round_off(0.0019932133) == 0.002
+    assert round_off(0.0000000000002193) == 0