Feature request : Allow multiple diluents calculation for multiple molecules in calc_spectrum
### 🐛 Describe the bug

Currently, we can computer spectrum for a single molecule in non-air diluents (implemented in [PR: 495](https://github.com/radis/radis/pull/495). For example: 
```
calc_spectrum(
        wavenum_min=2245,
        wavenum_max=2250,
        Tgas=700,
        path_length=0.1,
        molecule="CO",
        mole_fraction=0.2,
        isotope=1,
        databank="hitran",
        diluent={"CO2": 0.2, "H2O": 0.1, "air": 0.5}
)
```

If we want to compute spectrum for multiple molecules for multiple diluents, that is not currently implemented and will raise a `NotImplementError`:

```
calc_spectrum(
        wavenum_min=2245,
        wavenum_max=2250,
        Tgas=700,
        path_length=0.1,
        mole_fraction={"CO":0.4, "CO2":0.4},
        isotope=1,
        databank="hitran",
        export_lines=True,
)
```




### 💡 Possible solutions

Modify the [calc_spectrum](https://github.com/radis/radis/blob/develop/radis/lbl/calc.py) class, and check if molefraction is a dictionary, if yes then for every molecule in molefraction, use the diluent to create spectrum using SpectrumFactory(which is already done), and merge all the spectrums.

Note: Make sure to evaluate all the molefractions of diluents and molecules properly.

### 🎲 Radis version

0.13.1

### 💻 Operating system

Linux
diff --git a/docs/dev/_test.rst b/docs/dev/_test.rst
index 6d24882c..93d40460 100644
--- a/docs/dev/_test.rst
+++ b/docs/dev/_test.rst
@@ -108,6 +108,9 @@ adding the following lines within your test function::
         import matplotlib.pyplot as plt
         plt.ion()   # dont get stuck with Matplotlib if executing through pytest
 
+You can refer to Pull Request: https://github.com/radis/radis/pull/495 to see
+how test cases are written when a new feature is added.
+
 See: https://github.com/statsmodels/statsmodels/issues/3697
 
 
diff --git a/radis/lbl/bands.py b/radis/lbl/bands.py
index 0df41e63..ded83841 100644
--- a/radis/lbl/bands.py
+++ b/radis/lbl/bands.py
@@ -105,6 +105,7 @@ class BandFactory(BroadenFactory):
         self,
         Tgas,
         mole_fraction=None,
+        diluent=None,
         path_length=None,
         pressure=None,
         levels="all",
@@ -233,6 +234,9 @@ class BandFactory(BroadenFactory):
 
         # Line broadening
 
+        # ... generates molefraction for diluents
+        self._generate_diluent_molefraction(mole_fraction, diluent)
+
         # ... calculate broadening  HWHM
         self._calc_broadening_HWHM()
 
@@ -395,6 +399,7 @@ class BandFactory(BroadenFactory):
         Trot,
         Ttrans=None,
         mole_fraction=None,
+        diluent=None,
         path_length=None,
         pressure=None,
         vib_distribution="boltzmann",
@@ -556,6 +561,9 @@ class BandFactory(BroadenFactory):
 
         # Line broadening
 
+        # ... generates molefraction for diluents
+        self._generate_diluent_molefraction(mole_fraction, diluent)
+
         # ... calculate broadening  HWHM
         self._calc_broadening_HWHM()
 
diff --git a/radis/lbl/broadening.py b/radis/lbl/broadening.py
index 4e73d020..ad162c14 100644
--- a/radis/lbl/broadening.py
+++ b/radis/lbl/broadening.py
@@ -248,7 +248,16 @@ def gaussian_FT(w_centered, hwhm):
 
 
 def pressure_broadening_HWHM(
-    airbrd, selbrd, Tdpair, Tdpsel, pressure_atm, mole_fraction, Tgas, Tref
+    airbrd,
+    selbrd,
+    Tdpair,
+    Tdpsel,
+    pressure_atm,
+    mole_fraction,
+    Tgas,
+    Tref,
+    diluent,
+    diluent_broadening_coeff,
 ):
     """Calculates collisional broadening HWHM over all lines by scaling
     tabulated HWHM for new pressure and mole fractions conditions [1]_
@@ -278,6 +287,11 @@ def pressure_broadening_HWHM(
     Tref: float [K]
         reference temperature at which tabulated HWHM pressure
         broadening coefficients were tabulated
+    diluent: dictionary
+        contains diluent and their mole fraction
+    diluent_broadening_coeff: dictionary
+        contains all non air diluents broadening coefficients
+
 
     Returns
     -------
@@ -314,15 +328,37 @@ def pressure_broadening_HWHM(
     # | dev note: in that case we simplify the expression by calculation the
     # | power function once only.
 
-    if Tdpsel is None:
-        gamma_lb = ((Tref / Tgas) ** Tdpair) * (
-            (airbrd * pressure_atm * (1 - mole_fraction))
-            + (selbrd * pressure_atm * mole_fraction)
-        )
-    else:
-        gamma_lb = ((Tref / Tgas) ** Tdpair) * (
-            airbrd * pressure_atm * (1 - mole_fraction)
-        ) + ((Tref / Tgas) ** Tdpsel) * (selbrd * pressure_atm * mole_fraction)
+    diluent_molecules = diluent.keys()
+
+    # check if gamma_diluent and n_diluent exists or not
+    try:
+        gamma_lb = 0
+        for diluent_molecule, diluent_mole_fraction in diluent.items():
+            if diluent_molecule != "air":
+                gamma_lb += (
+                    (Tref / Tgas)
+                    ** diluent_broadening_coeff["n_" + diluent_molecule.lower()]
+                ) * (
+                    diluent_broadening_coeff["gamma_" + diluent_molecule.lower()]
+                    * pressure_atm
+                    * diluent_mole_fraction
+                )
+        # Adding air coefficient
+        if "air" in diluent_molecules:
+            gamma_lb += ((Tref / Tgas) ** Tdpair) * (
+                (airbrd * pressure_atm * diluent["air"])
+            )
+        # Adding self coefficient
+        if Tdpsel is None:  # use Tdpair instead
+            gamma_lb += ((Tref / Tgas) ** Tdpair) * (
+                (selbrd * pressure_atm * mole_fraction)
+            )
+        else:
+            gamma_lb += ((Tref / Tgas) ** Tdpsel) * (
+                selbrd * pressure_atm * mole_fraction
+            )
+    except KeyError as err:
+        raise KeyError("Column not found {0}".format(err))
 
     return gamma_lb
 
@@ -408,6 +444,8 @@ def voigt_broadening_HWHM(
     mole_fraction,
     Tgas,
     Tref,
+    diluent,
+    diluent_broadening_coeff,
 ):
     """Calculate Voigt profile half-width at half-maximum (HWHM) from the
     Gaussian and Collisional broadening with the empirical formula of [Olivero-1977]_
@@ -474,7 +512,16 @@ def voigt_broadening_HWHM(
 
     # Collisional broadening HWHM
     gamma_lb = pressure_broadening_HWHM(
-        airbrd, selbrd, Tdpair, Tdpsel, pressure_atm, mole_fraction, Tgas, Tref
+        airbrd,
+        selbrd,
+        Tdpair,
+        Tdpsel,
+        pressure_atm,
+        mole_fraction,
+        Tgas,
+        Tref,
+        diluent,
+        diluent_broadening_coeff,
     )
 
     # Doppler Broadening HWHM:
@@ -815,8 +862,10 @@ class BroadenFactory(BaseFactory):
 
         Run this method before using `_calc_lineshape`
         """
+
         # Init variables
         df = self.df1
+        diluent = self._diluent
 
         if len(df) == 0:
             return  # no lines
@@ -838,14 +887,37 @@ class BroadenFactory(BaseFactory):
             self.params.broadening_method
         )  # Lineshape broadening algorithm
 
+        # diluent and their broadening coeff dictionary
+        diluent_broadening_coeff = {}
+        for key in diluent:
+            if key != "air":
+                diluent_broadening_coeff["gamma_" + key.lower()] = df[
+                    "gamma_" + key.lower()
+                ]
+                diluent_broadening_coeff["n_" + key.lower()] = df["n_" + key.lower()]
+
         # Get broadenings
         if broadening_method == "voigt":
             # Adds hwhm_voigt, hwhm_gauss, hwhm_lorentz:
-            self._add_voigt_broadening_HWHM(df, pressure_atm, mole_fraction, Tgas, Tref)
+            self._add_voigt_broadening_HWHM(
+                df,
+                pressure_atm,
+                mole_fraction,
+                Tgas,
+                Tref,
+                diluent,
+                diluent_broadening_coeff,
+            )
         elif broadening_method in ["convolve", "fft"]:
             # Adds hwhm_lorentz:
             self._add_collisional_broadening_HWHM(
-                df, pressure_atm, mole_fraction, Tgas, Tref
+                df,
+                pressure_atm,
+                mole_fraction,
+                Tgas,
+                Tref,
+                diluent,
+                diluent_broadening_coeff,
             )
             # Add hwhm_gauss:
             self._add_doppler_broadening_HWHM(df, Tgas)
@@ -946,7 +1018,16 @@ class BroadenFactory(BaseFactory):
 
         return
 
-    def _add_voigt_broadening_HWHM(self, df, pressure_atm, mole_fraction, Tgas, Tref):
+    def _add_voigt_broadening_HWHM(
+        self,
+        df,
+        pressure_atm,
+        mole_fraction,
+        Tgas,
+        Tref,
+        diluent,
+        diluent_broadening_coeff,
+    ):
         """Update dataframe with Voigt HWHM.
 
         Returns
@@ -997,6 +1078,8 @@ class BroadenFactory(BaseFactory):
             mole_fraction,
             Tgas,
             Tref,
+            diluent,
+            diluent_broadening_coeff,
         )
 
         # Update dataframe
@@ -1007,7 +1090,14 @@ class BroadenFactory(BaseFactory):
         return
 
     def _add_collisional_broadening_HWHM(
-        self, df, pressure_atm, mole_fraction, Tgas, Tref
+        self,
+        df,
+        pressure_atm,
+        mole_fraction,
+        Tgas,
+        Tref,
+        diluent,
+        diluent_broadening_coeff,
     ):
         """Update dataframe with collisional HWHM [1]_
 
@@ -1063,11 +1153,15 @@ class BroadenFactory(BaseFactory):
             mole_fraction,
             Tgas,
             Tref,
+            diluent,
+            diluent_broadening_coeff,
         )
 
         # Update dataframe
         df["hwhm_lorentz"] = wl
 
+        #####
+
         return
 
     def _add_doppler_broadening_HWHM(self, df, Tgas):
diff --git a/radis/lbl/calc.py b/radis/lbl/calc.py
index 0786794b..b64458d8 100644
--- a/radis/lbl/calc.py
+++ b/radis/lbl/calc.py
@@ -46,6 +46,7 @@ def calc_spectrum(
     molecule=None,
     isotope="all",
     mole_fraction=1,
+    diluent="air",
     path_length=1,
     databank="hitran",
     medium="air",
@@ -129,6 +130,17 @@ def calc_spectrum(
 
             mole_fraction={'CO2': 0.8, 'CO':0.2}​
 
+    diluent: str or dictionary
+            can be a string of a single diluent or a dictionary containing diluent
+            name as key and its mole_fraction as value
+            For single diluent ::
+
+                diluent = 'CO2'
+
+            For multiple diluents ::
+
+                diluent = { 'CO2': 0.6, 'H2O':0.2}
+
     path_length: float [:math:`cm`] or `~astropy.units.quantity.Quantity`
         slab size. Default ``1`` cm​. Use arbitrary units::
 
@@ -399,6 +411,12 @@ def calc_spectrum(
     # (dealt with separately because we cannot use them to guess what are the input molecules)
     DICT_INPUT_DICT_ARGUMENTS = {"overpopulation": overpopulation}
 
+    # Check if mole_fraction is dictionary of multiple molecules and have non air diluent, if yes then raise NotImplementedError
+    if isinstance(mole_fraction, dict) and isinstance(diluent, dict):
+        raise NotImplementedError(
+            "Spectrum calculation for multiple molecules in multiple diluents is not yet implemented."
+        )
+
     def _check_molecules_are_consistent(
         molecule_reference_set, reference_name, new_argument, new_argument_name
     ):
@@ -544,6 +562,7 @@ def calc_spectrum(
             mode=mode,
             export_lines=export_lines,
             return_factory=return_factory,
+            diluent=diluent,
             **kwargs_molecule,
         )
         if return_factory:
@@ -601,6 +620,7 @@ def _calc_spectrum_one_molecule(
     mode,
     export_lines,
     return_factory=False,
+    diluent="air",
     **kwargs,
 ) -> Spectrum:
     """See :py:func:`~radis.lbl.calc.calc_spectrum`
@@ -668,6 +688,7 @@ def _calc_spectrum_one_molecule(
         pressure=pressure,
         wstep=wstep,
         truncation=truncation,
+        mole_fraction=mole_fraction,
         neighbour_lines=neighbour_lines,
         cutoff=cutoff,
         parsum_mode=parsum_mode,
@@ -676,11 +697,23 @@ def _calc_spectrum_one_molecule(
         chunksize=chunksize,
         broadening_method=broadening_method,
         export_lines=export_lines,
+        diluent=diluent,
         **kwargs,
     )
     # Have consistent output units
     sf.input_wunit = input_wunit
 
+    # Checking diluent other than air present
+    if isinstance(diluent, str):
+        if diluent == "air":
+            diluent_other_than_air = False
+        else:
+            diluent_other_than_air = True
+    else:
+        diluent_other_than_air = len(diluent) > 1 or (
+            len(diluent) == 1 and "air" not in diluent
+        )
+
     # Load databank
     # -------------
 
@@ -746,18 +779,25 @@ def _calc_spectrum_one_molecule(
             # constants (not all molecules are supported!)
             conditions["levelsfmt"] = "radis"
             conditions["lvl_use_cached"] = use_cached
-
         # Columns to load
         if export_lines:
             conditions["load_columns"] = "all"
+            conditions["extra_params"] = "all"
+        elif not _equilibrium and diluent_other_than_air:
+            conditions["load_columns"] = ["noneq", "diluent"]
+            conditions["extra_params"] = "all"
         elif not _equilibrium:
             conditions["load_columns"] = "noneq"
+        elif diluent_other_than_air:
+            conditions["load_columns"] = ["equilibrium", "diluent"]
+            conditions["extra_params"] = "all"
         else:
             conditions["load_columns"] = "equilibrium"
 
         conditions["load_energies"] = not _equilibrium
         # Details to identify lines
         conditions["parse_local_global_quanta"] = (not _equilibrium) or export_lines
+
         # Finally, LOAD :
         sf.fetch_databank(**conditions)
     elif exists(databank):
@@ -802,8 +842,12 @@ def _calc_spectrum_one_molecule(
         # Columns to load
         if export_lines:
             conditions["load_columns"] = "all"
+        elif not _equilibrium and diluent_other_than_air:
+            conditions["load_columns"] = ["noneq", "diluent"]
         elif not _equilibrium:
             conditions["load_columns"] = "noneq"
+        elif diluent_other_than_air:
+            conditions["load_columns"] = ["equilibrium", "diluent"]
         else:
             conditions["load_columns"] = "equilibrium"
 
@@ -816,8 +860,12 @@ def _calc_spectrum_one_molecule(
         # Columns to load
         if export_lines:
             load_columns = "all"
+        elif not _equilibrium and diluent_other_than_air:
+            load_columns = ["noneq", "diluent"]
         elif not _equilibrium:
             load_columns = "noneq"
+        elif diluent_other_than_air:
+            load_columns = ["equilibrium", "diluent"]
         else:
             load_columns = "equilibrium"
 
diff --git a/radis/lbl/factory.py b/radis/lbl/factory.py
index 4459d6e8..3bd21b36 100644
--- a/radis/lbl/factory.py
+++ b/radis/lbl/factory.py
@@ -155,6 +155,9 @@ class SpectrumFactory(BandFactory):
     medium: ``'air'``, ``'vacuum'``
         propagating medium when giving inputs with ``'wavenum_min'``, ``'wavenum_max'``.
         Does not change anything when giving inputs in wavenumber. Default ``'air'``
+    diluent: ``str`` or ``dictionary``
+            can be a string of a single diluent or a dictionary containing diluent
+            name as key and its mole_fraction as value. Default ``air``.
 
     Other Parameters
     ----------------
@@ -398,6 +401,7 @@ class SpectrumFactory(BandFactory):
         export_populations=None,
         export_lines=False,
         emulate_gpu=False,
+        diluent="air",
         **kwargs,
     ):
 
@@ -514,6 +518,16 @@ class SpectrumFactory(BandFactory):
         if not isinstance(isotope, str):
             isotope = ",".join([str(k) for k in list_if_float(isotope)])
 
+        # If molecule present in diluent, raise error
+        if (isinstance(diluent, str) and diluent == molecule) or (
+            isinstance(diluent, dict) and molecule in diluent.keys()
+        ):
+            raise KeyError(
+                "{0} is being called as molecule and diluent, please remove it from diluent.".format(
+                    molecule
+                )
+            )
+
         # Initialize input conditions
         self.input.wavenum_min = wavenum_min
         self.input.wavenum_max = wavenum_max
@@ -535,6 +549,8 @@ class SpectrumFactory(BandFactory):
         # Initialize computation variables
         self.params.wstep = wstep
         self.params.pseudo_continuum_threshold = pseudo_continuum_threshold
+        self.params.diluent = diluent
+        self._diluent = None
 
         if cutoff is None:
             # If None, use no cutoff : https://github.com/radis/radis/pull/259
@@ -648,7 +664,13 @@ class SpectrumFactory(BandFactory):
     # XXX =====================================================================
 
     def eq_spectrum(
-        self, Tgas, mole_fraction=None, path_length=None, pressure=None, name=None
+        self,
+        Tgas,
+        mole_fraction=None,
+        path_length=None,
+        diluent=None,
+        pressure=None,
+        name=None,
     ) -> Spectrum:
         """Generate a spectrum at equilibrium.
 
@@ -658,6 +680,9 @@ class SpectrumFactory(BandFactory):
             Gas temperature (K)
         mole_fraction: float
             database species mole fraction. If None, Factory mole fraction is used.
+        diluent: str or dictionary
+            can be a string of a single diluent or a dictionary containing diluent
+            name as key and its mole_fraction as value
         path_length: float or `~astropy.units.quantity.Quantity`
             slab size (cm). If ``None``, the default Factory
             :py:attr:`~radis.lbl.factory.SpectrumFactor.input.path_length` is used.
@@ -705,7 +730,6 @@ class SpectrumFactory(BandFactory):
         --------
         :meth:`~radis.lbl.factory.SpectrumFactory.non_eq_spectrum`
         """
-
         # %% Preprocessing
         # --------------------------------------------------------------------
 
@@ -786,6 +810,9 @@ class SpectrumFactory(BandFactory):
         # ----------------------------------------------------------------------
         # Line broadening
 
+        # ... generates molefraction for diluents
+        self._generate_diluent_molefraction(mole_fraction, diluent)
+
         # ... calculate broadening  HWHM
         self._calc_broadening_HWHM()
 
@@ -850,6 +877,7 @@ class SpectrumFactory(BandFactory):
                 "lines_cutoff": self._Nlines_cutoff,
                 "lines_in_continuum": self._Nlines_in_continuum,
                 "thermal_equilibrium": True,
+                "diluents": self._diluent,
                 "radis_version": version,
                 "spectral_points": (
                     self.params.wavenum_max_calc - self.params.wavenum_min_calc
@@ -925,6 +953,7 @@ class SpectrumFactory(BandFactory):
         self,
         Tgas,
         mole_fraction=None,
+        diluent=None,
         path_length=None,
         pressure=None,
         name=None,
@@ -1193,6 +1222,7 @@ class SpectrumFactory(BandFactory):
                 ],
                 "lines_calculated": _Nlines_calculated,
                 "thermal_equilibrium": True,
+                "diluents": self._diluent,
                 "radis_version": version,
                 "emulate_gpu": emulate,
                 "spectral_points": (
@@ -1411,6 +1441,7 @@ class SpectrumFactory(BandFactory):
         Ttrans=None,
         mole_fraction=None,
         path_length=None,
+        diluent=None,
         pressure=None,
         vib_distribution="boltzmann",
         rot_distribution="boltzmann",
@@ -1435,6 +1466,9 @@ class SpectrumFactory(BandFactory):
             which is the RT characteristic time)
         mole_fraction: float
             database species mole fraction. If None, Factory mole fraction is used.
+        diluent: str or dictionary
+            can be a string of a single diluent or a dictionary containing diluent
+            name as key and its mole_fraction as value
         path_length: float or `~astropy.units.quantity.Quantity`
             slab size (cm). If ``None``, the default Factory
             :py:attr:`~radis.lbl.factory.SpectrumFactor.input.path_length` is used.
@@ -1518,7 +1552,6 @@ class SpectrumFactory(BandFactory):
 
         # %% Preprocessing
         # --------------------------------------------------------------------
-
         # Convert units
         Tvib = convert_and_strip_units(Tvib, u.K)
         Trot = convert_and_strip_units(Trot, u.K)
@@ -1549,7 +1582,6 @@ class SpectrumFactory(BandFactory):
         self.input.overpopulation = overpopulation
         self.input.rot_distribution = rot_distribution
         self.input.vib_distribution = vib_distribution
-
         # Get translational temperature
         Tgas = Ttrans
         if Tgas is None:
@@ -1566,7 +1598,6 @@ class SpectrumFactory(BandFactory):
 
         # New Profiler object
         self._reset_profiler(verbose)
-
         # Check variables
         self._check_inputs(mole_fraction, max(flatten(Tgas, Tvib, Trot)))
 
@@ -1628,6 +1659,9 @@ class SpectrumFactory(BandFactory):
 
         # Line broadening
 
+        # ... generates molefraction for diluents
+        self._generate_diluent_molefraction(mole_fraction, diluent)
+
         # ... calculate broadening  HWHM
         self._calc_broadening_HWHM()
 
@@ -1725,6 +1759,7 @@ class SpectrumFactory(BandFactory):
                 "lines_cutoff": self._Nlines_cutoff,
                 "lines_in_continuum": self._Nlines_in_continuum,
                 "thermal_equilibrium": False,  # dont even try to guess if it's at equilibrium
+                "diluents": self._diluent,
                 "radis_version": version,
                 "spectral_points": (
                     self.params.wavenum_max_calc - self.params.wavenum_min_calc
@@ -1884,6 +1919,50 @@ class SpectrumFactory(BandFactory):
 
         return
 
+    def _generate_diluent_molefraction(self, mole_fraction, diluent):
+        from radis.misc.warning import MoleFractionError
+
+        molecule = self.input.molecule
+        # Check if diluent is same as molecule or not
+        if (isinstance(diluent, str) and diluent == molecule) or (
+            isinstance(diluent, dict) and molecule in diluent.keys()
+        ):
+            raise KeyError(
+                "{0} is being called as molecule and diluent, please remove it from diluent.".format(
+                    molecule
+                )
+            )
+
+        # Using Spectrum Factory values of diluents
+        if diluent is None:
+            # If diluent is air, then add remaining molefraction as 'air'
+            if isinstance(self.params.diluent, str):
+                diluents = {self.params.diluent: 1 - mole_fraction}
+            else:
+                diluents = self.params.diluent.copy()
+        else:
+            if isinstance(diluent, str):
+                diluents = {diluent: 1 - mole_fraction}
+            else:
+                diluents = diluent.copy()
+
+        # Checking mole_fraction of molecule and diluent
+        total_mole_fraction = mole_fraction + sum(list(diluents.values()))
+        if total_mole_fraction > 1:
+            raise MoleFractionError(
+                "Total molefraction = {0} of molecule and diluents greater than 1. Please set appropriate molefraction value of molecule and diluents.".format(
+                    total_mole_fraction
+                )
+            )
+        elif total_mole_fraction < 1:
+            raise MoleFractionError(
+                "Total molefraction = {0} of molecule and diluents less than 1. Please set appropriate molefraction value of molecule and diluents".format(
+                    total_mole_fraction
+                )
+            )
+
+        self._diluent = diluents
+
     def predict_time(self):
         """predict_time(self) uses the input parameters like Spectral Range, Number of lines, wstep,
         truncation to predict the estimated calculation time for the Spectrum
diff --git a/radis/lbl/loader.py b/radis/lbl/loader.py
index d9a66b7d..beda1ff9 100644
--- a/radis/lbl/loader.py
+++ b/radis/lbl/loader.py
@@ -245,6 +245,18 @@ required_non_eq = [
 ]
 """list: column names required for non-equilibrium calculations.
 See load_column= key of fetch_databank() and load_databank() """
+
+broadening_coeff = [
+    "gamma_co2",
+    "n_co2",
+    "gamma_h2o",
+    "n_h2o",
+    "gamma_he",
+    "n_he",
+    "gamma_h2",
+    "n_h2",
+]
+"""list: column names required for non-air diluent calculations."""
 # TODO refactor : directly go & parse the identifications names (globu, locu, etc.)
 # in radis.io.hitran ?
 # For the moment we just try to be exhaustive
@@ -474,6 +486,7 @@ class Parameters(ConditionDict):
         "wavenum_min_calc",
         "waveunit",
         "wstep",
+        "diluent",
     ]
 
     def __init__(self):
@@ -501,6 +514,7 @@ class Parameters(ConditionDict):
         self.wavenum_min_calc = None  #: float: minimum calculated wavenumber (cm-1) initialized by SpectrumFactory
         self.waveunit = "cm-1"  #: waverange unit: should be cm-1.
         self.wstep = None  #: float: spectral resolution (cm-1)
+        self.diluent = {}  # dict: molecule : mole fraction
         self.dxL = _lorentzian_step(
             0.01
         )  #: float : Lorentzian step for LDM lineshape database. Default _lorentzian_step(0.01)
@@ -891,6 +905,20 @@ class DatabankLoader(object):
         self.df0 = None  # type : pd.DataFrame
         self._reset_references()  # bibliographic references
 
+    def columns_list_to_load(self, load_columns_type):
+        # Which columns to load
+        if load_columns_type == "equilibrium":
+            columns = list(drop_all_but_these)
+        elif load_columns_type == "noneq":
+            columns = list(set(drop_all_but_these) | set(required_non_eq))
+        elif load_columns_type == "diluent":
+            columns = list(broadening_coeff)
+        else:
+            raise ValueError(
+                f"Expected a list or 'all' for `load_columns`, got `load_columns={load_columns_type}"
+            )
+        return columns
+
     def fetch_databank(
         self,
         source="hitran",
@@ -987,7 +1015,7 @@ class DatabankLoader(object):
         parallel: bool
             if ``True``, uses joblib.parallel to load database with multiple processes
             (works only for HITEMP files)
-        load_columns: list, ``'all'``, ``'equilibrium'``, ``'noneq'``
+        load_columns: list, ``'all'``, ``'equilibrium'``, ``'noneq'``, ``diluent``,
             columns names to load.
             If ``'equilibrium'``, only load the columns required for equilibrium
             calculations. If ``'noneq'``, also load the columns required for
@@ -995,6 +1023,8 @@ class DatabankLoader(object):
             If ``'all'``, load everything. Note that for performances, it is
             better to load only certain columsn rather than loading them all
             and dropping them with ``drop_columns``.
+            If ``diluent`` then all additional columns required for calculating spectrum
+            in that diluent is loaded.
             Default ``'equilibrium'``.
 
             .. warning::
@@ -1092,20 +1122,32 @@ class DatabankLoader(object):
         self.params.db_use_cached = db_use_cached
         self.params.lvl_use_cached = lvl_use_cached
 
-        # %% Which columns to load
-        if load_columns == "equilibrium":
-            columns = list(drop_all_but_these)
-        elif load_columns == "noneq":
-            columns = list(set(drop_all_but_these) | set(required_non_eq))
-        elif load_columns == "all":
+        # Which columns to load
+        columns = []
+        if "all" in load_columns:
             columns = None  # see fetch_hitemp, fetch_hitran, etc.
-        elif isinstance(load_columns, list):
-            columns = list(set(drop_all_but_these) | set(load_columns))
-        else:
+        elif isinstance(load_columns, str) and load_columns in ["equilibrium", "noneq"]:
+            columns = self.columns_list_to_load(load_columns)
+        elif load_columns == "diluent":
             raise ValueError(
-                f"Expected a list or 'all' for `load_columns`, got `load_columns={load_columns}"
+                "Please use diluent along with 'equilibrium' or 'noneq' in a list like ['diluent','noneq']"
             )
 
+        elif isinstance(load_columns, list) and "all" not in load_columns:
+            for load_columns_type in load_columns:
+                if load_columns_type in ["equilibrium", "noneq", "diluent"]:
+                    for col in self.columns_list_to_load(load_columns_type):
+                        columns.append(col)
+                elif load_columns_type in list(
+                    set(drop_all_but_these)
+                    | set(required_non_eq)
+                    | set(broadening_coeff)
+                ):
+                    columns.append(load_columns_type)
+                else:
+                    raise ValueError("invalid column name provided")
+            columns = list(set(columns))
+
         # %% Init Line database
         # ---------------------
         self._reset_references()  # bibliographic references
@@ -2074,6 +2116,7 @@ class DatabankLoader(object):
         load_columns,
         include_neighbouring_lines=True,
     ) -> pd.DataFrame:
+
         """Loads all available database files and keep the relevant one.
         Returns a Pandas dataframe.
 
@@ -2120,7 +2163,6 @@ class DatabankLoader(object):
             because of lineshape broadening. The ``neighbour_lines``
             parameter is used to determine the limit. Default ``True``.
         """
-
         # Check inputs
         assert db_use_cached in [True, False, "regen", "force"]
 
@@ -2144,20 +2186,32 @@ class DatabankLoader(object):
                 + drop_auto_columns_for_levelsfmt[levelsfmt]
             )
 
-        # Which columns to load
-        if load_columns == "equilibrium":
-            columns = list(drop_all_but_these)
-        elif load_columns == "noneq":
-            columns = list(set(drop_all_but_these) | set(required_non_eq))
-        elif load_columns == "all":
+        # which columns to load
+        columns = []
+        if "all" in load_columns:
             columns = None  # see fetch_hitemp, fetch_hitran, etc.
-        elif isinstance(load_columns, list):
-            columns = list(set(drop_all_but_these) | set(load_columns))
-        else:
+        elif isinstance(load_columns, str) and load_columns in ["equilibrium", "noneq"]:
+            columns = self.columns_list_to_load(load_columns)
+        elif load_columns == "diluent":
             raise ValueError(
-                f"Expected a list or 'all' for `load_columns`, got `load_columns={load_columns}"
+                "Please use diluent along with 'equilibrium' or 'noneq' in a list like ['diluent','noneq']"
             )
 
+        elif isinstance(load_columns, list) and "all" not in load_columns:
+            for load_columns_type in load_columns:
+                if load_columns_type in ["equilibrium", "noneq", "diluent"]:
+                    for col in self.columns_list_to_load(load_columns_type):
+                        columns.append(col)
+                elif load_columns_type in list(
+                    set(drop_all_but_these)
+                    | set(required_non_eq)
+                    | set(broadening_coeff)
+                ):
+                    columns.append(load_columns_type)
+                else:
+                    raise ValueError("invalid column name provided")
+            columns = list(set(columns))
+
         # subroutine load_and_concat
         # --------------------------------------
         def load_and_concat(files):
diff --git a/radis/misc/warning.py b/radis/misc/warning.py
index 3b4e6fe5..7c84baad 100644
--- a/radis/misc/warning.py
+++ b/radis/misc/warning.py
@@ -79,6 +79,10 @@ class InconsistentDatabaseError(ValueError):
     pass
 
 
+class MoleFractionError(ValueError):
+    pass
+
+
 # Warnings
 
 
diff --git a/radis/test/lbl/test_calc.py b/radis/test/lbl/test_calc.py
index 883f30c1..64262a9f 100644
--- a/radis/test/lbl/test_calc.py
+++ b/radis/test/lbl/test_calc.py
@@ -819,6 +819,84 @@ def test_check_wavelength_range(verbose=True, warnings=True, *args, **kwargs):
     return True
 
 
+def test_non_air_diluent_calc(verbose=True, plot=False, warnings=True, *args, **kwargs):
+    from radis import plot_diff
+
+    if verbose:
+        printm("Regenerating database with all params, and calculating non_air diluent")
+
+    # Regenerating database of CO and calculating non_air_diluent spectrum
+    s1 = calc_spectrum(
+        wavenum_min=2245,
+        wavenum_max=2250,
+        Tgas=700,
+        path_length=0.1,
+        molecule="CO",
+        mole_fraction=0.2,
+        isotope=1,
+        databank="hitran",
+        use_cached="regen",
+        diluent={"CO2": 0.2, "H2O": 0.1, "air": 0.5},
+        export_lines=True,
+    )
+
+    # Calculating spectrum with air as diluent
+    s2 = calc_spectrum(
+        wavenum_min=2245,
+        wavenum_max=2250,
+        Tgas=700,
+        path_length=0.1,
+        molecule="CO",
+        mole_fraction=0.2,
+        isotope=1,
+        databank="hitran",
+        export_lines=True,
+    )
+
+    assert s1.get_conditions()["diluents"] == {"CO2": 0.2, "H2O": 0.1, "air": 0.5}
+    assert s2.get_conditions()["diluents"] == {"air": 0.8}
+
+    hwhm_voigt_s1 = s1.lines.hwhm_voigt
+    hwhm_voigt_s2 = s2.lines.hwhm_voigt
+
+    assert (hwhm_voigt_s2 < hwhm_voigt_s1).all()
+    assert np.isclose(hwhm_voigt_s2[0], 0.022718389546218788)
+    assert np.isclose(hwhm_voigt_s1[0], 0.02530070148135749)
+
+    if plot:
+        plot_diff(
+            s1,
+            s2,
+            method=["diff"],
+            label1="Diluent CO2, H2O, air",
+            label2="Diluent air",
+            title="Diluent CO2:0.2, H2O:0.1, air:0.5",
+        )
+
+    return True
+
+
+def test_diluent_invalid(verbose=True, plot=False, *args, **kwargs):
+
+    with pytest.raises(KeyError) as err:
+        calc_spectrum(
+            wavenum_min=2245,
+            wavenum_max=2250,
+            Tgas=700,
+            path_length=0.1,
+            molecule="CO",
+            mole_fraction=0.2,
+            isotope=1,
+            databank="hitran",
+            diluent="CO",
+        )
+
+    assert (
+        "is being called as molecule and diluent, please remove it from diluent."
+        in str(err.value)
+    )
+
+
 def _run_testcases(plot=True, verbose=True, warnings=True, *args, **kwargs):
 
     # Test sPlanck and conversion functions
@@ -852,6 +930,7 @@ def _run_testcases(plot=True, verbose=True, warnings=True, *args, **kwargs):
     test_calc_spectrum_multiple_molecules_wstep_auto()
 
     test_check_wavelength_range()
+    test_non_air_diluent_calc()
 
     return True
 
diff --git a/radis/test/lbl/test_factory.py b/radis/test/lbl/test_factory.py
index be5dccf6..90bdce7a 100644
--- a/radis/test/lbl/test_factory.py
+++ b/radis/test/lbl/test_factory.py
@@ -644,6 +644,87 @@ def test_all_spectrum_using_wstep_auto(verbose=True, plot=False, *args, **kwargs
     assert sf._wstep == "auto"
 
 
+@pytest.mark.fast
+def test_non_air_diluent(verbose=True, plot=False, *args, **kwargs):
+
+    sf = SpectrumFactory(
+        wavelength_min=4200,
+        wavelength_max=4500,
+        cutoff=1e-23,
+        molecule="CO",
+        isotope="1,2",
+        truncation=5,
+        neighbour_lines=10,
+        path_length=0.1,
+        mole_fraction=0.1,
+        medium="vacuum",
+        optimization=None,
+        verbose=verbose,
+    )
+    sf.warnings.update(
+        {
+            "MissingSelfBroadeningWarning": "ignore",
+            "NegativeEnergiesWarning": "ignore",
+            "LinestrengthCutoffWarning": "ignore",
+            "HighTemperatureWarning": "ignore",
+            "AccuracyWarning": "ignore",
+            "PerformanceWarning": "ignore",
+        }
+    )
+
+    sf.load_databank("HITRAN-CO", load_columns=["diluent", "equilibrium"])
+
+    # Calculating spectrum for different diluents
+    sf.eq_spectrum(Tgas=2000)
+    wl1 = sf.df1["hwhm_lorentz"]
+    assert sf._diluent == {"air": 0.9}
+
+    sf.eq_spectrum(Tgas=2000, diluent={"CO2": 0.4, "air": 0.5})
+    wl2 = sf.df1["hwhm_lorentz"]
+    assert sf._diluent == {"CO2": 0.4, "air": 0.5}
+
+    sf.eq_spectrum(Tgas=2000, diluent="CO2")
+    wl3 = sf.df1["hwhm_lorentz"]
+    assert sf._diluent == {"CO2": 0.9}
+
+    assert (wl1 < wl2).all() and (wl2 < wl3).all()
+
+
+@pytest.mark.fast
+def test_diluents_molefraction(verbose=True, plot=False, *args, **kwargs):
+    from radis.misc.warning import MoleFractionError
+
+    sf = SpectrumFactory(
+        wavelength_min=4300,
+        wavelength_max=4500,
+        wstep=0.01,
+        cutoff=1e-30,
+        pressure=1,
+        isotope=[1],
+        verbose=verbose,
+        diluent={"CO2": 0.4, "air": 0.2},
+    )
+    sf.load_databank("HITRAN-CO", load_columns=["diluent", "equilibrium"])
+    # Molefraction (molecule + diluent) < 1
+    with pytest.raises(MoleFractionError) as err:
+        sf.eq_spectrum(Tgas=300, mole_fraction=0.3)
+    assert (
+        "of molecule and diluents less than 1. Please set appropriate molefraction value of molecule and diluents"
+        in str(err.value)
+    )
+
+    # Molefraction (molecule + diluent) > 1
+    with pytest.raises(MoleFractionError) as err:
+        sf.eq_spectrum(Tgas=300, mole_fraction=0.6)
+    assert (
+        "of molecule and diluents greater than 1. Please set appropriate molefraction value of molecule and diluents."
+        in str(err.value)
+    )
+
+    # Molefraction (molecule + diluent) == 1
+    sf.eq_spectrum(Tgas=300, mole_fraction=0.4)
+
+
 # --------------------------
 if __name__ == "__main__":
 
['docs/dev/_test.rst', 'radis/lbl/loader.py', 'radis/test/lbl/test_factory.py', 'radis/lbl/calc.py', 'radis/lbl/bands.py', 'radis/lbl/factory.py', 'radis/lbl/broadening.py', 'radis/misc/warning.py', 'radis/test/lbl/test_calc.py']