Repository: radis/radis (language: python)

Issue:
Implement all new fitting modules and accompanied gallery examples into RADIS
This pull request is to implement all of my modules, developed separately in my repo [RADIS-Spectrum-Fitting-Benchmark](https://github.com/TranHuuNhatHuy/RADIS-Spectrum-Fitting-Benchmark), into RADIS codebase. This pull request features:

- Add `new_fitting.py`, which contains my new fitting module and models, into radis/tools/.
- Add gallery examples to `radis/examples`, including:
   - `plot_newfitting_Tgas.py` : the most basic example of how to use new fitting module, including the formats and so on.
   - `plot_newfitting_Trot-Tvib-molfrac.py` : real-life fitting case provided by Mr. @CorentinGrimaldi , featuring non-LTE CO spectrum in which we will fit `Trot`, `Tvib` and `mole_fraction`.
- Add spectrum and other auxiliary files into `radis/test/files` that are necessary for those gallery examples above.

More illustrative cases, such as comparison between the old and new fitting modules, comparison between different algorithms, or Mr. @minouHub's case, will be pushed in another PR after this PR has been merged.

[This pull request is a newer version of the old one, #519 ].


Discussion / hints:
# [Codecov](https://codecov.io/gh/radis/radis/pull/520?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radis) Report
> Merging [#520](https://codecov.io/gh/radis/radis/pull/520?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radis) (9541fb6) into [develop](https://codecov.io/gh/radis/radis/commit/3465180a2483f7529721917be0372ed8ad1b5651?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=radis) (3465180) will **decrease** coverage by `1.07%`.
> The diff coverage is `0.00%`.

```diff
@@             Coverage Diff             @@
##           develop     #520      +/-   ##
===========================================
- Coverage    74.53%   73.46%   -1.08%     
===========================================
  Files          140      141       +1     
  Lines        19141    19421     +280     
===========================================
  Hits         14267    14267              
- Misses        4874     5154     +280     
```

Looks good to me. You can merge 

Fix patch:
diff --git a/docs/conf.py b/docs/conf.py
index c8dad57b4..b224f4e06 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -177,6 +177,7 @@ def setup(app):
     "fitroom": ("https://fitroom.readthedocs.io/en/latest/", None),
     "habanero": ("https://habanero.readthedocs.io/en/latest/", None),
     "joblib": ("https://joblib.readthedocs.io/en/latest/", None),
+    "lmfit": ("https://lmfit.github.io/lmfit-py/", None),
     "numpy": ("https://numpy.org/doc/stable/", None),
     "matplotlib": ("https://matplotlib.org/stable/", None),
     "pandas": ("https://pandas.pydata.org/pandas-docs/stable/", None),
diff --git a/examples/plot_newfitting_Tgas-molfrac.py b/examples/plot_newfitting_Tgas-molfrac.py
new file mode 100644
index 000000000..78badbce8
--- /dev/null
+++ b/examples/plot_newfitting_Tgas-molfrac.py
@@ -0,0 +1,122 @@
+# -*- coding: utf-8 -*-
+"""
+
+================================================================================
+Fit an LTE spectrum with multiple fit parameters using new fitting module
+================================================================================
+
+With the new fitting module introduced in :py:func:`~radis.tools.new_fitting.fit_spectrum` function,
+and in the example of Tgas fitting using new fitting module, we can see its 1-temperature fitting
+performance for equilibrium condition.
+
+This example features how new fitting module can fit an equilibrium spectrum, with multiple fit
+parameters, including gas temperature, mole fraction, and wavelength offset.
+
+This is a real fitting case introduced by Mr. Nicolas Minesi, featuring CO spectrum with absorbance
+as spectral quantity to be fitted. As we can see, he stored the experimental result in a MATLAB file,
+and from there a Spectrum object is generated. It is worth noticing that, the result seems to differ
+slightly from ground-truth, due to the fact that currently RADIS uses air broadening parameteres for
+calculation, while this experiment was originally conducted in Argon. Future updates on other molecules'
+broadening coefficients will increase the accuracy of these cases with non-air diluents.
+
+"""
+
+import astropy.units as u
+import scipy.io
+
+from radis import Spectrum
+from radis.test.utils import getTestFile
+from radis.tools.new_fitting import fit_spectrum
+
+# ------------------------------------ Step 1. Load experimental spectrum ------------------------------------ #
+
+
+data_file = "trimmed_1857_VoigtCO_Minesi.mat"
+data = scipy.io.loadmat(getTestFile(data_file), simplify_cells=True)["CO_resu_Voigt"]
+index = 20
+s_experimental = Spectrum.from_array(
+    data["nu"], data["A_exp"][:, index], "absorbance", wunit="cm-1", unit=""
+)  # adimensioned
+
+
+# ------------------------------------ Step 2. Fill ground-truths and data ------------------------------------ #
+
+
+# Experimental conditions which will be used for spectrum modeling. Basically, these are known ground-truths.
+experimental_conditions = {
+    "molecule": "CO",  # Molecule ID
+    "isotope": "1",  # Isotope ID, can have multiple at once
+    "wmin": 2010.6,  # Starting wavelength/wavenumber to be cropped out from the original experimental spectrum.
+    "wmax": 2011.6,  # Ending wavelength/wavenumber for the cropping range.
+    "wunit": "cm-1",  # Accompanying unit of those 2 wavelengths/wavenumbers above.
+    "pressure": 1
+    * u.bar,  # Total pressure of gas, in "bar" unit by default, but you can use Astropy units too.
+    "path_length": 10
+    * u.cm,  # Experimental path length, in "cm" unit by default, but you can use Astropy units too.
+    "databank": "hitemp",  # Databank used for calculation. Must be stated.
+}
+
+# List of parameters to be fitted.
+fit_parameters = {
+    "Tgas": 5000,  # Fit parameter, accompanied by its initial value.
+    "mole_fraction": 0.05,  # Species mole fraction, from 0 to 1.
+    "offset": "0 cm-1",  # Experimental offset, must be a blank space separating offset amount and unit.
+}
+
+# List of bounding ranges applied for those fit parameters above.
+bounding_ranges = {
+    "Tgas": [
+        2000,
+        9000,
+    ],  # Bounding ranges for each fit parameter stated above. You can skip this step, but not recommended.
+    "mole_fraction": [0, 1],  # Species mole fraction, from 0 to 1.
+    "offset": [
+        -0.1,
+        0.1,
+    ],  # Experimental offset, must be a blank space separating offset amount and unit
+}
+
+# Fitting pipeline setups.
+fit_properties = {
+    "method": "lbfgsb",  # Preferred fitting method. By default, "leastsq".
+    "fit_var": "absorbance",  # Spectral quantity to be extracted for fitting process, such as "radiance", "absorbance", etc.
+    "normalize": False,  # Either applying normalization on both spectra or not.
+    "max_loop": 300,  # Max number of loops allowed. By default, 100.
+    "tol": 1e-20,  # Fitting tolerance, only applicable for "lbfgsb" method.
+}
+
+"""
+
+For the fitting method, you can try one among 17 different fitting methods and algorithms of LMFIT,
+introduced in `LMFIT method list <https://lmfit.github.io/lmfit-py/fitting.html#choosing-different-fitting-methods>`.
+
+You can see the benchmark result of these algorithms here:
+`RADIS Newfitting Algorithm Benchmark <https://github.com/radis/radis-benchmark/blob/master/manual_benchmarks/plot_newfitting_comparison_algorithm.py>`.
+
+"""
+
+
+# ------------------------------------ Step 3. Run the fitting and retrieve results ------------------------------------ #
+
+
+# Conduct the fitting process!
+s_best, result, log = fit_spectrum(
+    s_exp=s_experimental,
+    fit_params=fit_parameters,
+    bounds=bounding_ranges,
+    model=experimental_conditions,
+    pipeline=fit_properties,
+)
+
+
+# Now investigate the log
+
+print("\nResidual history: \n")
+print(log["residual"])
+
+print("\nFitted values history: \n")
+for fit_val in log["fit_vals"]:
+    print(fit_val)
+
+print("\nTotal fitting time: ")
+print(log["time_fitting"], end=" s\n")
diff --git a/examples/plot_newfitting_Tgas.py b/examples/plot_newfitting_Tgas.py
index e930d448f..ebcedc43a 100644
--- a/examples/plot_newfitting_Tgas.py
+++ b/examples/plot_newfitting_Tgas.py
@@ -2,50 +2,57 @@
 """
 
 ================================================================================
-Fit an equilibrium spectrum with Tgas as fit parameter using new fitting modules
+New fitting module introduction - simple 1-temperature LTE case
 ================================================================================
 
-RADIS has its own fitting feature, as shown in 1 temperature fit example, where you have
-to manually create the spectrum model, input the experimental spectrum and other ground-truths
-into numerous RADIS native functions, as well as adjust the fitting pipeline yourself.
+RADIS has its own fitting feature, as shown in
+`1 temperature fit example <https://radis.readthedocs.io/en/latest/auto_examples/plot_1T_fit.html>`,
+where you have to manually create the spectrum model, input the experimental spectrum and other
+ground-truths into numerous RADIS native functions, and adjust the fitting pipeline yourself.
 
-Now with the new fitting module released, all you have to do is to prepare a .spec file
-containing your experimental spectrum, fill some JSON forms describing the ground-truth conditions
-just like how you fill your medical checkup paper, call the function fit_spectrum() and let it do
-all the work! If you are not satisfied with the result, you can simply adjust the parameters in your
-JSON, such as slit and path_length , then recall the function again, until the results are satisfied.
-
-This is way easier and more convenient than dwelling into RADIS documentation to find out and
-learn how to use the current example, especially for new RADIS users. This new fitting module
-aims to provide an end-to-end fitting experience, with minimum amount of RADIS knowledge needed.
+Now with the new fitting module released, all you have to do is to prepare a .spec file containing
+your experimental spectrum, fill some JSON forms describing the ground-truth conditions just like how
+you fill a medical checkup paper, call the function :py:func:`~radis.tools.new_fitting.fit_spectrum`
+and let it do all the work! If you are not satisfied with the result, you can simply adjust the
+parameters in your JSON such as slit and path_length, recall the function until the results are good.
 
 Instruction:
 
 - Step 1: prepare a .spec file. Create a .spec file containing your experimental spectrum. You can
-  do it with RADIS by saving a Spectrum object with s.store(). If your current data is not a Spectrum
-  object, you can convert it to a Spectrum object from Python arrays or from text files, and then save
-  it as .spec file using s.store() as mentioned above.
+  do it with RADIS by saving a Spectrum object with :py:meth:`~radis.spectrum.spectrum.Spectrum.store`.
+  If your current data is not a Spectrum object, you can convert it to a Spectrum object from Python
+  arrays or from text files, and then save it as .spec file as above.
 
 - Step 2: fill the JSON forms. There are 4 JSON forms you need to fill: `experimental_conditions` with
   ground-truth data about your experimental environment, `fit_parameters` with the parameters you need
   to fit (such as Tgas, mole fraction, etc.), `bounding_ranges` with fitting ranges for parameters you
   listed in `fit_parameters`, and `fit_properties` for some fitting pipeline references.
 
-- Step 3: call `fit_spectrum()` with the experimental spectrum and 4 JSON forms and see the result.
+- Step 3: call :py:func:`~radis.tools.new_fitting.fit_spectrum` with the experimental spectrum and 4
+  JSON forms, then see the result.
 
 This example features fitting an experimental spectrum with Tgas, using new fitting modules.
 
 
 """
 
+import astropy.units as u
+
 from radis import load_spec
 from radis.test.utils import getTestFile
 from radis.tools.new_fitting import fit_spectrum
 
+# ------------------------------------ Step 1. Load experimental spectrum ------------------------------------ #
+
+
 # Load an experimental spectrum. You can prepare yours, or fetch one of them in the radis/test/files directory.
 my_spec = getTestFile("synth-NH3-1-500-2000cm-P10-mf0.01-p1.spec")
 s_experimental = load_spec(my_spec)
 
+
+# ------------------------------------ Step 2. Fill ground-truths and data ------------------------------------ #
+
+
 # Experimental conditions which will be used for spectrum modeling. Basically, these are known ground-truths.
 experimental_conditions = {
     "molecule": "NH3",  # Molecule ID
@@ -54,11 +61,12 @@
     "wmax": 1050,  # Ending wavelength/wavenumber for the cropping range.
     "wunit": "cm-1",  # Accompanying unit of those 2 wavelengths/wavenumbers above.
     "mole_fraction": 0.01,  # Species mole fraction, from 0 to 1.
-    "pressure": 10,  # Partial pressure of gas, in "bar" unit.
-    "path_length": 1,  # Experimental path length, in "cm" unit.
+    "pressure": 1e6
+    * u.Pa,  # Total pressure of gas, in "bar" unit by default, but you can also use Astropy units.
+    "path_length": 10
+    * u.mm,  # Experimental path length, in "cm" unit by default, but you can also use Astropy units.
     "slit": "1 nm",  # Experimental slit, must be a blank space separating slit amount and unit.
     "offset": "-0.2 nm",  # Experimental offset, must be a blank space separating offset amount and unit.
-    "wstep": 0.003,  # Resolution of wavenumber grid, in cm-1.
     "databank": "hitran",  # Databank used for the spectrum calculation. Must be stated.
 }
 
@@ -79,13 +87,26 @@
 
 # Fitting pipeline setups.
 fit_properties = {
-    "method": "leastsq",  # Preferred fitting method from the 17 confirmed methods of LMFIT stated in week 4 blog. By default, "leastsq".
+    "method": "least_squares",  # Preferred fitting method from the 17 confirmed methods of LMFIT stated in week 4 blog. By default, "leastsq".
     "fit_var": "radiance",  # Spectral quantity to be extracted for fitting process, such as "radiance", "absorbance", etc.
     "normalize": False,  # Either applying normalization on both spectra or not.
     "max_loop": 150,  # Max number of loops allowed. By default, 200.
     "tol": 1e-15,  # Fitting tolerance, only applicable for "lbfgsb" method.
 }
 
+"""
+
+For the fitting method, you can try one among 17 different fitting methods and algorithms of LMFIT,
+introduced in `LMFIT method list <https://lmfit.github.io/lmfit-py/fitting.html#choosing-different-fitting-methods>`.
+
+You can see the benchmark result of these algorithms here:
+`RADIS Newfitting Algorithm Benchmark <https://github.com/radis/radis-benchmark/blob/master/manual_benchmarks/plot_newfitting_comparison_algorithm.py>`.
+
+"""
+
+
+# ------------------------------------ Step 3. Run the fitting and retrieve results ------------------------------------ #
+
 
 # Conduct the fitting process!
 s_best, result, log = fit_spectrum(
@@ -94,7 +115,7 @@
     bounds=bounding_ranges,  # Bounding ranges for those fit parameters.
     model=experimental_conditions,  # Experimental ground-truths conditions.
     pipeline=fit_properties,  # Fitting pipeline references.
-    verbose=False,  # If you want a clean result, stay False. If you want to see more about each loop, True.
+    fit_kws={"gtol": 1e-12},
 )
 
 
diff --git a/examples/plot_newfitting_Trot-Tvib-molfrac.py b/examples/plot_newfitting_Trot-Tvib-molfrac.py
index 10512b001..964e86e0d 100644
--- a/examples/plot_newfitting_Trot-Tvib-molfrac.py
+++ b/examples/plot_newfitting_Trot-Tvib-molfrac.py
@@ -2,7 +2,7 @@
 """
 
 ================================================================================
-Fit a non-LTE spectrum with multiple fit parameters
+Fit a non-LTE spectrum with multiple fit parameters using new fitting module
 ================================================================================
 
 With the new fitting module introduced in :py:func:`~radis.tools.new_fitting.fit_spectrum` function,
@@ -12,13 +12,15 @@
 This example features how new fitting module can fit non-equilibrium spectrum, with multiple fit
 parameters, such as vibrational/rotational temperatures, or mole fraction, etc.
 
-This is a real fitting case introduced by Mr. Corentin Grimaldi, featuring CO molecule emitting on
-a wide range of spectrum. This case also features a user-defined trapezoid slit function, which is
-accounted by the distorion (dispersion) of the slit, as a result of the influences from experimental
-parameters and spectrometer dispersion parameters during the experiment.
+This is a real fitting case introduced by Mr. Corentin Grimaldi, a scientist and member of RADIS
+community. This case features CO molecule emitting on a wide range of spectrum. This case also includes
+a user-defined trapezoid slit function, which is accounted by the distorion (dispersion) of the slit,
+as a result of the influences from experimental parameters and spectrometer dispersion parameters
+during the experiment.
 
 """
 
+import astropy.units as u
 import numpy as np
 
 from radis import load_spec
@@ -58,21 +60,30 @@ def slit_dispersion(w):
 # ------------------------ End of Mr. Grimaldi's script ------------------------ #
 
 
-# Load experimental spectrum.
+# ------------------------------------ Step 1. Load experimental spectrum ------------------------------------ #
+
+
 specName = (
     "Corentin_0_100cm_DownSampled_20cm_10pctCO2_1-wc-gw450-gr300-sl1500-acc5000-.spec"
 )
 s_experimental = load_spec(getTestFile(specName)).offset(-0.6, "nm")
 
+
+# ------------------------------------ Step 2. Fill ground-truths and data ------------------------------------ #
+
+
 # Experimental conditions which will be used for spectrum modeling. Basically, these are known ground-truths.
 experimental_conditions = {
     "molecule": "CO",  # Molecule ID
     "isotope": "1,2,3",  # Isotope ID, can have multiple at once
-    "wmin": 2270,  # Starting wavelength/wavenumber to be cropped out from the original experimental spectrum.
-    "wmax": 2400,  # Ending wavelength/wavenumber for the cropping range.
-    "wunit": "nm",  # Accompanying unit of those 2 wavelengths/wavenumbers above.
-    "pressure": 1.01325,  # Partial pressure of gas, in "bar" unit.
-    "path_length": 1 / 195,  # Experimental path length, in "cm" unit.
+    "wmin": 2270
+    * u.nm,  # Starting wavelength/wavenumber to be cropped out from the original experimental spectrum.
+    "wmax": 2400 * u.nm,  # Ending wavelength/wavenumber for the cropping range.
+    "pressure": 1.01325
+    * u.bar,  # Total pressure of gas, in "bar" unit by default, but you can also use Astropy units.
+    "path_length": 1
+    / 195
+    * u.cm,  # Experimental path length, in "cm" unit by default, but you can also use Astropy units.
     "slit": {  # Experimental slit. In simple form: "[value] [unit]", i.e. "-0.2 nm". In complex form: a dict with parameters of apply_slit()
         "slit_function": (top_slit_nm, base_slit_nm),
         "unit": "nm",
@@ -81,7 +92,6 @@ def slit_dispersion(w):
         "slit_dispersion": slit_dispersion,
         "inplace": False,
     },
-    "wstep": 0.001,  # Resolution of wavenumber grid, in cm-1.
     "cutoff": 0,  # (RADIS native) Discard linestrengths that are lower that this to reduce calculation time, in cm-1.
     "databank": "hitemp",  # Databank used for calculation. Must be stated.
 }
@@ -104,13 +114,26 @@ def slit_dispersion(w):
 
 # Fitting pipeline setups.
 fit_properties = {
-    "method": "lbfgsb",  # Preferred fitting method from the 17 confirmed methods of LMFIT stated in week 4 blog. By default, "leastsq".
+    "method": "lbfgsb",  # Preferred fitting method. By default, "leastsq".
     "fit_var": "radiance",  # Spectral quantity to be extracted for fitting process, such as "radiance", "absorbance", etc.
     "normalize": False,  # Either applying normalization on both spectra or not.
     "max_loop": 300,  # Max number of loops allowed. By default, 200.
     "tol": 1e-20,  # Fitting tolerance, only applicable for "lbfgsb" method.
 }
 
+"""
+
+For the fitting method, you can try one among 17 different fitting methods and algorithms of LMFIT,
+introduced in `LMFIT method list <https://lmfit.github.io/lmfit-py/fitting.html#choosing-different-fitting-methods>`.
+
+You can see the benchmark result of these algorithms here:
+`RADIS Newfitting Algorithm Benchmark <https://github.com/radis/radis-benchmark/blob/master/manual_benchmarks/plot_newfitting_comparison_algorithm.py>`.
+
+"""
+
+
+# ------------------------------------ Step 3. Run the fitting and retrieve results ------------------------------------ #
+
 
 # Conduct the fitting process!
 s_best, result, log = fit_spectrum(
diff --git a/examples/plot_newfitting_comparison_oldnew.py b/examples/plot_newfitting_comparison_oldnew.py
new file mode 100644
index 000000000..e55da2b14
--- /dev/null
+++ b/examples/plot_newfitting_comparison_oldnew.py
@@ -0,0 +1,229 @@
+# -*- coding: utf-8 -*-
+"""
+
+================================================================================
+Compare performance between old 1T fitting example and new fitting module
+================================================================================
+
+This example will try to fit an experimental spectrum using two fitting pipelines:
+
+- Old fitting module of :py:func:`~radis.tools.fitting.fit_spectrum`, which is the
+  current RADIS fitting module. You can find the gallery example featuring it at
+  `1 temperature fit <https://radis.readthedocs.io/en/latest/auto_examples/plot_1T_fit.html>`.
+- New fitting module of :py:func:`~radis.tools.new_fitting.fit_spectrum`, a new
+  fitting interface for a more practical and interactive fitting experience. See
+  `1 temperature fit using new module <https://radis.readthedocs.io/en/latest/auto_examples/plot_newfitting_Tgas.html>`
+
+With this, you can compare their overall performance, including number of loops,
+fitting time, and final residuals between experimental and best-fit spectra, as
+an indicator of fit accuracy.
+
+"""
+
+import time
+
+import astropy.units as u
+
+from radis import SpectrumFactory, load_spec, plot_diff
+from radis.test.utils import getTestFile, setup_test_line_databases
+from radis.tools.fitting import LTEModel
+from radis.tools.new_fitting import fit_spectrum
+
+# ------------------------------------ OLD 1-TEMPERATURE FITTING EXAMPLE ------------------------------------ #
+
+setup_test_line_databases()
+
+# fit range
+wlmin = 4167
+wlmax = 4180
+
+s_exp = (
+    load_spec(getTestFile("CO2_measured_spectrum_4-5um.spec"))
+    .crop(wlmin, wlmax, "nm")
+    .normalize()
+    .sort()
+    .offset(-0.2, "nm")
+)
+
+
+def LTEModel_withslitnorm(factory, fit_parameters, fixed_parameters):
+    s = LTEModel(factory, fit_parameters, fixed_parameters)
+    # we could also have added a fittable parameter, such as an offset,
+    # or made the slit width a fittable parameter.
+    # ... any paramter in model_input will be fitted.
+    # s.offset(model_input["offset"], 'nm')
+    s.apply_slit(1.4, "nm")
+    return s.take("radiance").normalize()
+
+
+sf = SpectrumFactory(
+    wlmin * u.nm,
+    wlmax * u.nm,
+    wstep=0.001,  # cm-1
+    pressure=1 * 1e-3,  # bar
+    cutoff=1e-25,
+    isotope="1,2",
+    path_length=10,  # cm-1
+    mole_fraction=1,
+    truncation=1,  # cm-1
+)
+sf.warnings["MissingSelfBroadeningWarning"] = "ignore"
+sf.warnings["HighTemperatureWarning"] = "ignore"
+sf.load_databank("HITRAN-CO2-TEST")
+
+begin_time_mark = time.time()
+
+s_best, best = sf.fit_spectrum(
+    s_exp.take("radiance"),
+    model=LTEModel_withslitnorm,
+    fit_parameters={
+        "Tgas": 300,
+        # "offset": 0
+    },
+    bounds={
+        "Tgas": [300, 2000],
+        # "offset": [-1, 1],
+    },
+    plot=True,
+    solver_options={
+        "maxiter": 15,  # 👈 increase to let the fit converge
+        "ftol": 1e-15,
+    },
+    verbose=2,
+)
+
+end_time_mark = time.time()
+
+plot_diff(s_exp, s_best)
+
+# Information to report later in comparison result
+oldfitting_residual = best.fun
+oldfitting_loops = best.nfev
+oldfitting_time = end_time_mark - begin_time_mark
+
+# ------------------------------------------------------------------------------------------------------------ #
+
+
+# -------------------------------------------- NEW FITTING MODULE -------------------------------------------- #
+
+
+# ------------------------------------ Step 1. Load experimental spectrum ------------------------------------ #
+
+
+# Load an experimental spectrum. You can prepare yours, or fetch one of them in the radis/test/files directory.
+my_spec = getTestFile("CO2_measured_spectrum_4-5um.spec")
+s_experimental = load_spec(my_spec).offset(-0.2, "nm")
+
+
+# ------------------------------------ Step 2. Fill ground-truths and data ------------------------------------ #
+
+
+# Experimental conditions which will be used for spectrum modeling. Basically, these are known ground-truths.
+experimental_conditions = {
+    "molecule": "CO2",  # Molecule ID
+    "isotope": "1,2",  # Isotope ID, can have multiple at once
+    "wmin": 4167
+    * u.nm,  # Starting wavelength/wavenumber to be cropped out from the original experimental spectrum.
+    "wmax": 4180 * u.nm,  # Ending wavelength/wavenumber for the cropping range.
+    "mole_fraction": 1,  # Species mole fraction, from 0 to 1.
+    "pressure": 1
+    * 1e-3
+    * u.bar,  # Total pressure of gas, in "bar" unit by default, but you can also use Astropy units.
+    "path_length": 10
+    * u.cm,  # Experimental path length, in "cm" unit by default, but you can also use Astropy units.
+    "slit": "1.4 nm",  # Experimental slit, must be a blank space separating slit amount and unit.
+    "wstep": 0.001,  # Resolution of wavenumber grid, in cm-1.
+    "databank": "hitran",  # Databank used for the spectrum calculation. Must be stated.
+}
+
+# List of parameters to be fitted, accompanied by their initial values.
+fit_parameters = {
+    "Tgas": 1150,  # Gas temperature, in K.
+}
+
+# List of bounding ranges applied for those fit parameters above.
+# You can skip this step and let it use default bounding ranges, but this is not recommended.
+# Bounding range must be at format [<lower bound>, <upper bound>].
+bounding_ranges = {
+    "Tgas": [
+        300,
+        2000,
+    ],
+}
+
+# Fitting pipeline setups.
+fit_properties = {
+    "method": "lbfgsb",  # Preferred fitting method. By default, "leastsq".
+    "fit_var": "radiance",  # Spectral quantity to be extracted for fitting process, such as "radiance", "absorbance", etc.
+    "normalize": True,  # Either applying normalization on both spectra or not.
+    "max_loop": 150,  # Max number of loops allowed. By default, 200.
+    "tol": 1e-15,  # Fitting tolerance, only applicable for "lbfgsb" method.
+}
+
+"""
+
+For the fitting method, you can try one among 17 different fitting methods and algorithms of LMFIT,
+introduced in `LMFIT method list <https://lmfit.github.io/lmfit-py/fitting.html#choosing-different-fitting-methods>`.
+
+You can see the benchmark result of these algorithms here:
+`RADIS Newfitting Algorithm Benchmark <https://github.com/radis/radis-benchmark/blob/master/manual_benchmarks/plot_newfitting_comparison_algorithm.py>`.
+
+"""
+
+
+# ------------------------------------ Step 3. Run the fitting and retrieve results ------------------------------------ #
+
+
+# Conduct the fitting process!
+s_bestfit, result, log = fit_spectrum(
+    s_exp=s_experimental,  # Experimental spectrum.
+    fit_params=fit_parameters,  # Fit parameters.
+    bounds=bounding_ranges,  # Bounding ranges for those fit parameters.
+    model=experimental_conditions,  # Experimental ground-truths conditions.
+    pipeline=fit_properties,  # Fitting pipeline references.
+    verbose=False,  # If you want a clean result, stay False. If you want to see more about each loop, True.
+)
+
+# Information to report later in comparison result
+newfitting_residual = log["residual"][-1]
+newfitting_loops = result.nfev
+newfitting_time = log["time_fitting"]
+
+
+# ---------------------------------------------------------------------------------------------------------------------- #
+
+
+# ---------------------------------- PERFORMANCE COMPARISON BETWEEN 2 FITTING METHODS ---------------------------------- #
+
+print(
+    "\n\n\n====================  PERFORMANCE COMPARISON BETWEEN 2 FITTING METHODS  ===================="
+)
+
+print("\n1. LAST RESIDUAL\n")
+print(f"- Old 1T fitting example:   \t{oldfitting_residual}")
+print(f"- New fitting module:       \t{newfitting_residual}")
+
+print("\n2. NUMBER OF FITTING LOOPS\n")
+print(f"- Old 1T fitting example:   \t{oldfitting_loops} loops")
+print(f"- New fitting module:       \t{newfitting_loops} loops")
+
+print("\n3. TOTAL TIME TAKEN (s)\n")
+print(f"- Old 1T fitting example:   \t{oldfitting_time} s")
+print(f"- New fitting module:       \t{newfitting_time} s")
+
+print(
+    "\n==========================================================================================\n"
+)
+
+
+"""
+
+From the comparative result above, which includes last residual, total time fitting and total number of
+fitting loops of each fitting method, we can see that under exactly the same ground-truth conditions,
+fitting method (L-BFGS-B) and fitting pipeline, we can see that the new fitting module provides a better
+accuracy, while requiring less fitting loops and time for execution.
+
+You are free to change the experimental spectrum and its accompanied ground-truth conditions, but please
+make sure to keep the same inputs between the two for a transparent comparative result.
+
+"""
diff --git a/radis/tools/fitting.py b/radis/tools/fitting.py
index 59c3a0646..25d65b269 100644
--- a/radis/tools/fitting.py
+++ b/radis/tools/fitting.py
@@ -49,6 +49,7 @@ def LTEModel(factory, fittable_parameters, fixed_parameters={}) -> Spectrum:
     See Also
     --------
     :py:func:`~radis.tools.fitting.Tvib12Tvib3Trot_NonLTEModel`
+    :py:func:`~radis.tools.new_fitting.residual_LTE`
 
     """
 
@@ -103,6 +104,7 @@ def Tvib12Tvib3Trot_NonLTEModel(
     See Also
     --------
     :py:func:`~radis.tools.fitting.LTEModel`
+    :py:func:`~radis.tools.new_fitting.residual_NonLTE`
     """
 
     fittable_parameters = fittable_parameters.copy()
@@ -207,6 +209,7 @@ def fit_spectrum(
     See Also
     --------
     :py:meth:`~radis.lbl.factory.SpectrumFactory.fit_spectrum`,
+    :py:func:`~radis.tools.new_fitting.fit_spectrum`,
     :py:mod:`fitroom`
     """
 
diff --git a/radis/tools/new_fitting.py b/radis/tools/new_fitting.py
index 6bcf1aeb9..14cca0b60 100644
--- a/radis/tools/new_fitting.py
+++ b/radis/tools/new_fitting.py
@@ -68,6 +68,16 @@ def residual_LTE(params, conditions, s_data, sf, log, verbose):
     residual: float
         residuals of the two spectra, using RADIS's get_residual().
 
+    Examples
+    --------
+
+    .. minigallery:: radis.tools.new_fitting.residual_LTE
+
+    See Also
+    --------
+    :py:func:`~radis.tools.new_fitting.residual_NonLTE`
+    :py:func:`~radis.tools.fitting.LTEModel`
+
     """
 
     # GENERATE LTE SPECTRUM BASED ON THOSE PARAMETERS
@@ -147,12 +157,6 @@ def residual_LTE(params, conditions, s_data, sf, log, verbose):
         current_fitvals.append(float(params[param]))
     log["fit_vals"].append(current_fitvals)
 
-    # Print information of fitting process
-    if verbose:
-        for param in params:
-            print(f"{param} = {float(params[param])}")
-        print(f"\nResidual = {residual}\n")
-
     return residual
 
 
@@ -189,6 +193,16 @@ def residual_NonLTE(params, conditions, s_data, sf, log, verbose):
     residual: float
         residuals of the two spectra, using RADIS's get_residual().
 
+    Examples
+    --------
+
+    .. minigallery:: radis.tools.new_fitting.residual_NonLTE
+
+    See Also
+    --------
+    :py:func:`~radis.tools.new_fitting.residual_LTE`
+    :py:func:`~radis.tools.fitting.radis.tools.fitting.Tvib12Tvib3Trot_NonLTEModel`
+
     """
 
     # GENERATE NON-LTE SPECTRUM BASED ON THOSE PARAMETERS
@@ -277,12 +291,6 @@ def residual_NonLTE(params, conditions, s_data, sf, log, verbose):
         current_fitvals.append(float(params[param]))
     log["fit_vals"].append(current_fitvals)
 
-    # Print information of fitting process
-    if verbose:
-        for param in params:
-            print(f"{param} = {float(params[param])}")
-        print(f"\nResidual = {residual}\n")
-
     return residual
 
 
@@ -313,12 +321,10 @@ def get_conditions(
     verbose: bool
         by default, True, print details about the fitting progress.
 
-    Returns
-    -------
-    conditions: dict
-        Dictionary object containing all ground-truth information.
-    params: Parameters
-        LMFIT.Parameters object that contains fit parameters.
+    Examples
+    --------
+
+    .. minigallery:: radis.tools.new_fitting.get_conditions
 
     """
 
@@ -457,7 +463,7 @@ def get_conditions(
 # -------------------------------------------------------------------------------------------- #
 
 
-def spectrum_refinement(s_data, conditions, verbose=True) -> Union[Spectrum, dict]:
+def spectrum_refinement(s_data, conditions, verbose) -> Union[Spectrum, dict]:
     """Receive an experimental spectrum and further refine it according to the provided pipeline
     and ground-truths. Refinement process includes extracting the desired spectrum quantity,
     removing NaN values, and other additional convolutions. Finally, a refined Spectrum object is
@@ -484,6 +490,11 @@ def spectrum_refinement(s_data, conditions, verbose=True) -> Union[Spectrum, dic
     conditions : dict
         the input conditions that might have been added fit_var (in case user didn't).
 
+    Examples
+    --------
+
+    .. minigallery:: radis.tools.new_fitting.spectrum_refinement
+
     """
 
     # Extract spectrum
@@ -517,7 +528,7 @@ def spectrum_refinement(s_data, conditions, verbose=True) -> Union[Spectrum, dic
     s_refined = (
         Spectrum(
             {fit_var: (s_data_mtr[0], s_data_mtr[1])},
-            wunit=conditions["model"]["wunit"],
+            wunit=s_data.get_waveunit(),
             units={fit_var: s_data.units[fit_var]},
         )
         .take(fit_var)
@@ -548,7 +559,9 @@ def fit_spectrum(
     pipeline=None,
     bounds=None,
     input_file=None,
-    verbose=False,
+    verbose=True,
+    show_plot=True,
+    fit_kws={},
 ) -> Union[Spectrum, MinimizerResult, dict]:
     """Fit an experimental spectrum (from here referred as "data spectrum") with a modeled one,
     then derive the fit results. Data spectrum is loaded from the path stated in the JSON file,
@@ -577,6 +590,12 @@ def fit_spectrum(
     ----------
     verbose : bool
         print details about each loop of the fitting process. By default, False.
+    show_plot : bool
+        show experimental spectrum before and after being fitted. By default, True. Showing the plot
+        will pause the runtime, so if you are running benchmark which involves multiple fitting cases,
+        set it to False.
+    fit_kws: dict
+        Options to pass to the minimizer being used. See :py:func:`~lmfit.minimizer.minimize`
 
     Returns
     -------
@@ -589,10 +608,85 @@ def fit_spectrum(
         not quite covered by the Minimizer, including: residual and fit
         values after each fitting loop, and total time elapsed.
 
+    Examples
+    --------
+
+    ::
+        from radis import load_spec
+        from radis.test.utils import getTestFile
+        from radis.tools.new_fitting import fit_spectrum
+
+        # Load an experimental spectrum. You can prepare yours, or fetch one of them in the radis/test/files directory.
+        my_spec = getTestFile("synth-NH3-1-500-2000cm-P10-mf0.01-p1.spec")
+        s_experimental = load_spec(my_spec)
+
+        # Experimental conditions which will be used for spectrum modeling. Basically, these are known ground-truths.
+        experimental_conditions = {
+            "molecule": "NH3",  # Molecule ID
+            "isotope": "1",  # Isotope ID, can have multiple at once
+            "wmin": 1000,  # Starting wavelength/wavenumber to be cropped out from the original experimental spectrum.
+            "wmax": 1050,  # Ending wavelength/wavenumber for the cropping range.
+            "wunit": "cm-1",  # Accompanying unit of those 2 wavelengths/wavenumbers above.
+            "mole_fraction": 0.01,  # Species mole fraction, from 0 to 1.
+            "pressure": 10,  # Total pressure of gas, in "bar" unit by default, but you can also use Astropy units.
+            "path_length": 1,  # Experimental path length, in "cm" unit by default, but you can also use Astropy units.
+            "slit": "1 nm",  # Experimental slit, must be a blank space separating slit amount and unit.
+            "offset": "-0.2 nm",  # Experimental offset, must be a blank space separating offset amount and unit.
+            "wstep": 0.003,  # Resolution of wavenumber grid, in cm-1.
+            "databank": "hitran",  # Databank used for the spectrum calculation. Must be stated.
+        }
+
+        # List of parameters to be fitted, accompanied by their initial values.
+        fit_parameters = {
+            "Tgas": 700,  # Gas temperature, in K.
+        }
+
+        # List of bounding ranges applied for those fit parameters above.
+        bounding_ranges = {
+            "Tgas": [
+                500,
+                2000,
+            ],
+        }
+
+        # Fitting pipeline setups.
+        fit_properties = {
+            "method": "leastsq",  # Preferred fitting method. By default, "leastsq".
+            "fit_var": "radiance",  # Spectral quantity to be extracted for fitting process, such as "radiance", "absorbance", etc.
+            "normalize": False,  # Either applying normalization on both spectra or not.
+            "max_loop": 150,  # Max number of loops allowed. By default, 200.
+            "tol": 1e-15,  # Fitting tolerance, only applicable for "lbfgsb" method.
+        }
+
+        # Conduct the fitting process!
+        s_best, result, log = fit_spectrum(
+            s_exp=s_experimental,  # Experimental spectrum.
+            fit_params=fit_parameters,  # Fit parameters.
+            bounds=bounding_ranges,  # Bounding ranges for those fit parameters.
+            model=experimental_conditions,  # Experimental ground-truths conditions.
+            pipeline=fit_properties,  # Fitting pipeline references.
+            verbose=False,  # If you want a clean result, stay False. If you want to see more about each loop, True.
+        )
+    .. minigallery:: radis.tools.new_fitting.fit_spectrum
+    Notes
+    -----
+    For the fitting method, you can try one among 17 different fitting methods and algorithms of LMFIT, introduced
+    in `LMFIT method list <https://lmfit.github.io/lmfit-py/fitting.html#choosing-different-fitting-methods>`.
+    See Also
+    --------
+    :py:func:`~radis.tools.fitting.fit_spectrum()`,
+    :py:meth:`~radis.lbl.factory.SpectrumFactory.fit_spectrum`,
+    :py:mod:`fitroom`
+
     """
 
     begin = time.time()  # Start the fitting time counter
 
+    if verbose:
+        print(
+            "\n======================= COMMENCE FITTING PROCESS ======================="
+        )
+
     # ACQUIRE AND REFINE EXPERIMENTAL SPECTRUM s_data
 
     if input_file != None:  # The user uses JSON file
@@ -621,37 +715,39 @@ def fit_spectrum(
 
         # Merge all separated dict inputs into a single dict first
 
-        conditions = {"model": {}, "fit": {}, "bounds": {}, "pipeline": {}}
-
-        for key in model:
-            conditions["model"][key] = model[key]
-        for key in fit_params:
-            conditions["fit"][key] = fit_params[key]
-        if bounds != None:
-            for key in bounds:
-                conditions["bounds"][key] = bounds[key]
-        for key in pipeline:
-            conditions["pipeline"][key] = pipeline[key]
+        conditions = {
+            "model": model,
+            "fit": fit_params,
+            "bounds": bounds if bounds is not None else {},
+            "pipeline": pipeline,
+        }
 
         # Load data, create Parameters object, from the above dict
         conditions, params = get_conditions(conditions)
 
         # Load and crop the experimental spectrum
-        s_data = s_exp.crop(
-            conditions["model"]["wmin"],
-            conditions["model"]["wmax"],
-            conditions["model"]["wunit"],
-        )
+        if "wunit" in conditions["model"]:  # Users state wunit directly
+            s_data = s_exp.crop(
+                conditions["model"]["wmin"],
+                conditions["model"]["wmax"],
+                conditions["model"]["wunit"],
+            )
+        else:  # Users don't state wunit directly, probably astropy units
+            s_data = s_exp.crop(
+                conditions["model"]["wmin"],
+                conditions["model"]["wmax"],
+            )
 
     # Log the time mark when experimental data is retrieved
     end_exp_load = time.time()
     time_exp_load = end_exp_load - begin
 
-    print(f"\nSuccessfully retrieved the experimental data in {time_exp_load}s.\n")
+    if verbose:
+        print(f"\nSuccessfully retrieved the experimental data in {time_exp_load}s.")
 
     # Further refine the data spectrum before calculating diff
 
-    s_data, conditions = spectrum_refinement(s_data, conditions)
+    s_data, conditions = spectrum_refinement(s_data, conditions, verbose)
 
     pipeline = conditions["pipeline"]
     fit_var = pipeline["fit_var"]
@@ -660,8 +756,9 @@ def fit_spectrum(
     end_exp_refine = time.time()
     time_exp_refine = end_exp_refine - end_exp_load
 
-    print(f"Successfully refined the experimental data in {time_exp_refine}s.")
     if verbose:
+        print(f"Successfully refined the experimental data in {time_exp_refine}s.")
+    if show_plot:
         s_data.plot(show=True)
 
     # PRE-MINIMIZATION SETUP
@@ -706,11 +803,13 @@ def fit_spectrum(
 
     # Determine additional fitting conditions for the minimizer
 
-    fit_kws = {
-        "max_nfev": pipeline["max_loop"]
-        if "max_loop" in pipeline
-        else 200,  # Default max number of loops: 200
-    }
+    fit_kws.update(
+        {
+            "max_nfev": pipeline["max_loop"]
+            if "max_loop" in pipeline
+            else 200,  # Default max number of loops: 200
+        }
+    )
 
     if "tol" in pipeline:
         if pipeline["method"] == "lbfgsb":
@@ -733,7 +832,9 @@ def fit_spectrum(
 
         sf.fetch_databank(databank, load_columns="equilibrium")
 
-        print("\nCommence fitting process for LTE spectrum!")
+        if verbose:
+            print("\nCommence fitting process for LTE spectrum!\n")
+
         result = minimize(
             residual_LTE,
             params,
@@ -747,7 +848,9 @@ def fit_spectrum(
 
         sf.fetch_databank(databank, load_columns="noneq")
 
-        print("\nCommence fitting process for non-LTE spectrum!")
+        if verbose:
+            print("\nCommence fitting process for non-LTE spectrum!")
+
         result = minimize(
             residual_NonLTE,
             params,
@@ -758,7 +861,7 @@ def fit_spectrum(
 
     # POST-MINIMIZATION REPORT
 
-    print(fit_report(result))  # Report the fitting result
+    print(fit_report(result), end="\n")  # Report the fitting result
 
     # REGENERATE THE BEST-FIT SPECTRUM, AS A RESULT FROM THE MINIMIZER
 
@@ -820,62 +923,26 @@ def fit_spectrum(
     end_fitting = time.time()
     time_fitting = float(end_fitting - begin_fitting)
 
-    print(f"\nSuccesfully finished the fitting process in {time_fitting}s.")
+    if verbose:
+        print(f"\nSuccessfully finished the fitting process in {time_fitting}s.")
+
     log["time_fitting"] = time_fitting
 
     # PLOT THE DIFFERENCE BETWEEN THE TWO
 
-    # We must show this plot regardless of verbose
-    plot_diff(
-        s_data,
-        s_result,
-        fit_var,
-        method=["diff", "ratio"],
-        show=True,
-    )
-
-    return s_result, result, log
-
-
-# ----------------------------------------------------------------------- #
+    # If show_plot = True by default, show the diff plot
+    if show_plot:
+        plot_diff(
+            s_data,
+            s_result,
+            fit_var,
+            method=["diff", "ratio"],
+            show=True,
+        )
 
+    if verbose:
+        print(
+            "\n======================== END OF FITTING PROCESS ========================\n"
+        )
 
-# if __name__ == "__main__":
-
-#     spec_list = [
-#         "CO2_measured_spectrum_4-5um",                                      # 0
-#         "synth-CO-1-1800-2300-cm-1-P3-t1500-v-r-mf0.1-p1-sl1nm",            # 1
-#         "synth-CO2-1-500-1100-cm-1-P2-t900-v-r-mf0.5-p1-sl1nm",             # 2
-#         "synth-CO2-1-500-3000-cm-1-P93-t740-v-r-mf0.96-p1-sl1nm",           # 3
-#         "synth-CO2-1-3300-3700-cm-1-P0.005-t3000-v-r-mf0.01-p1-sl1.4nm",    # 4
-#         "synth-H2O-1-1000-2500-cm-1-P0.5-t1500-v-r-mf0.5-p1-sl1nm",         # 5
-#         "synth-NH3-1-500-2000-cm-1-P10-t1000-v-r-mf0.01-p1-sl1nm",          # 6
-#         "synth-O2-1-7500-8000-cm-1-P1.01325-t298.15-v-r-mf0.21-p1-sl1nm",   # 7
-#     ]
-
-#     for i in range(len(spec_list)):
-#         input_path = f"../data/LTE/ground-truth/{spec_list[i]}.json"
-#         _, result, log, pipeline = fit_spectrum(input_path)
-#         json_data = {
-#             "fileName": f"{spec_list[i]}.spec",
-#             "pipeline": {
-#                 "method": pipeline["method"],
-# 	            "fit_var": "radiance",
-# 	            "normalize": pipeline["normalize"],
-# 	            "max_loop": 100
-#             },
-#             "result": {
-#                 "last_residual": log["residual"][-1],
-#                 "loops": result.nfev,
-#                 "time": log["time_fitting"]
-#             }
-#         }
-
-#         with open(f"../data/LTE/result/{spec_list[i]}/best_fit/pipeline.json", 'w') as f:
-#             json.dump(json_data, f, indent = 2)
-#             print("JSON file successfully created.")
-
-#         with open(f"../data/LTE/result/{spec_list[i]}/best_fit/log_residuals.txt", 'w') as f:
-#             for resi in log["residual"]:
-#                 f.write(f"{resi}\n")
-#             f.close()
+    return s_result, result, log


Test patch:
diff --git a/radis/test/files/trimmed_1857_VoigtCO_Minesi.mat b/radis/test/files/trimmed_1857_VoigtCO_Minesi.mat
new file mode 100644
index 000000000..67eac2100
Binary files /dev/null and b/radis/test/files/trimmed_1857_VoigtCO_Minesi.mat differ