Merge pull request #400 from radis/develop

0.11.0

README.rst

@@ -6,9 +6,9 @@
 *****************************************
 
 RADIS is a fast line-by-line code for high resolution infrared molecular spectra (emission / absorption,
-equilibrium / nonequilibrium).
+equilibrium / non-LTE) based on HITRAN/HITEMP/ExoMol.
 
-Includes post-processing tools to compare experimental spectra and spectra calculated
+It includes post-processing tools to compare experimental spectra and spectra calculated
 with RADIS or other spectral codes.
 
 User guide, install procedure and examples are available on the `RADIS Website <http://radis.readthedocs.io/>`__:

docs/lbl/lbl.rst

@@ -364,11 +364,15 @@ and :py:mod:`~astropy.units` ::
                          cutoff=1e-25,              # cm/molecule
                          broadening_max_width=10,   # cm-1
                          )
-    sf.load_databank('HITRAN-CO2-TEST')        # this database must be defined in ~/radis.json
+    sf.load_databank('HITRAN-CO2-TEST', load_columns='noneq')        # this database must be defined in ~/radis.json
     s1 = sf.eq_spectrum(Tgas=300 * u.K)
     s2 = sf.eq_spectrum(Tgas=2000 * u.K)
     s3 = sf.non_eq_spectrum(Tvib=2000 * u.K, Trot=300 * u.K)
 
+Note that for non-LTE calculations, specific columns must be loaded. This is done by using the
+``load_columns='noneq'`` parameter. See :py:meth:`~radis.lbl.loader.DatabankLoader.load_databank`
+for more information.
+
 
 .. _label_lbl_config_file:
 
@@ -737,6 +741,21 @@ points for each linewidth.
     but is still experimental in 0.9.30. Feedback welcome!
 
 
+Sparse wavenumber grid
+----------------------
+
+To compute large band spectra with a small number of lines, RADIS includes
+a sparse wavenumber implementation of the DIT algorithm, which is
+activated based on a scarcity criterion (``Nlines/Ngrid_points > 1``).
+
+The sparse version can be forced to be activated or deactivated. This behavior
+is done by setting the `SPARSE_WAVENUMBER` key of the :py:attr:`radis.config`
+dictionary, or of the ~/radis.json user file.
+
+See the :ref:`HITRAN full-range example <example_hitran_full_range>` for an
+example.
+
+
 Database loading
 ----------------
 

docs/spectrum/spectrum.rst

@@ -68,6 +68,7 @@ be used to batch-generate multiple spectra using the same line database::
 
     from radis import SpectrumFactory
     sf = SpectrumFactory(...)
+    sf.fetch_databank("hitemp", load_columns='noneq')  # or 'hitran', 'exomol', etc.
     s = sf.eq_spectrum(...)
     s2 = sf.non_eq_spectrum(...)
 
@@ -103,8 +104,14 @@ convenience function::
     s = Spectrum.from_array(w, T, 'transmittance_noslit',
                                wunit='nm', unit='') # adimensioned
 
+Dimensionned arrays can also be used directly ::
 
-Convenience functions have been added to handle the usual cases:
+    import astropy.units as u
+    w = np.linspace(200, 300) * u.nm
+    I = np.random.rand(len(w)) * u.mW/u.cm**2/u.sr/u.nm
+    s = Spectrum.from_array(w, I, 'radiance_noslit')
+
+Other convenience functions have been added to handle the usual cases:
 :func:`~radis.spectrum.models.calculated_spectrum`,
 :func:`~radis.spectrum.models.transmittance_spectrum` and
 :func:`~radis.spectrum.models.experimental_spectrum`::
@@ -116,6 +123,16 @@ Convenience functions have been added to handle the usual cases:
     s3 = experimental_spectrum(w, I, wunit='nm', Iunit='W/cm2/sr/nm')   # creates 'radiance'
 
 
+Initialize from Specutils
+-------------------------
+
+Use :py:meth:`~radis.spectrum.spectrum.Spectrum.from_specutils` to convert
+from a ``specutils`` :py:class:`specutils.spectra.spectrum1d.Spectrum1D` object ::
+
+    from radis import Spectrum
+    Spectrum.from_specutils(spectrum)
+
+
 Initialize from a text file
 ---------------------------
 
@@ -156,9 +173,21 @@ Spectrum conditions and get Spectrum that suits specific parameters
 .. minigallery:: radis.load_spec
 
 
-Want a quick spectrum
+Load from a HDF5 file
 ---------------------
 
+This is the fastest way to read a Spectrum object from disk. It keeps metadata
+and units, and you can also load only a part of a very large spectrum.
+Use :py:meth:`~radis.spectrum.spectrum.Spectrum.from_hdf5`  ::
+
+    Spectrum.from_hdf5("rad_hdf.h5", wmin=2100, wmax=2200, columns=['abscoeff', 'emisscoeff'])
+
+
+
+
+Calculate a test spectrum
+-------------------------
+
 You need a spectrum immediatly, to run some tests ? Use :py:func:`~radis.test.utils.test_spectrum` ::
 
     s = radis.test_spectrum()
@@ -446,6 +475,14 @@ automatically retrieves a spectrum from a database if you calculated it already,
 or calculates it and stores it if you didn't. Very useful for spectra that
 requires long computational times!
 
+Export to hdf5
+--------------
+
+This is the fastest way to dump a Spectrum object on disk (and also, it keeps metadata
+and therefore units !). Use :py:meth:`~radis.spectrum.spectrum.Spectrum.to_hdf5`  ::
+
+    s.to_hdf5('spec01.h5')
+
 
 Export to txt
 -------------
@@ -460,6 +497,28 @@ will export a given spectral arrays::
 
     s.savetxt('radiance_W_cm2_sr_um.csv', 'radiance_noslit', wunit='nm', Iunit='W/cm2/sr/µm')
 
+Export to Pandas
+----------------
+
+Use :py:meth:`~radis.spectrum.spectrum.Spectrum.to_pandas` ::
+
+    s.to_pandas()
+
+This will return a DataFrame with all spectral arrays as columns.
+
+
+Export to Specutils
+-------------------
+
+Use :py:meth:`~radis.spectrum.spectrum.Spectrum.to_specutils` to convert
+to a to ``specutils`` :py:class:`specutils.spectra.spectrum1d.Spectrum1D` object ::
+
+    s.to_specutils()
+
+
+.. minigallery:: radis.spectrum.spectrum.Spectrum.to_specutils
+
+
 
 How to modify a Spectrum object?
 ================================
@@ -605,13 +664,15 @@ For instance, remove a constant baseline in a given unit::
 
     s -= 0.1 * u.Unit('W/cm2/sr/nm')
 
-Here, se normalize a spectrum manually, assuming the spectrum units is in "mW/cm2/sr/nm" ::
+The :py:meth:`~radis.spectrum.spectrum.Spectrum.max` function returns a dimensionned
+value, therefore it can be used to normalize a spectrum directly :  ::
 
-    s /= s.max() * u.Unit("mW/cm2/sr/nm")
+    s /= s.max()
 
-Or below, we calibrate a Spectrum, assuming the spectrum units is in "count" ::
+Or below, we calibrate a Spectrum, assuming the spectrum units is in "count",
+and that our calibration show we have 94 :math:`mW/cm2/sr/nm` per count.  ::
 
-    s *= 100 * u.Unit("mW/cm2/sr/nm/count")
+    s *= 94 * u.Unit("mW/cm2/sr/nm") / u.Unit("count")
 
 
 .. _label_spectrum_offset_crop:
@@ -666,7 +727,7 @@ Or::
 If you want to create a new spectrum, don't forget to set ``inplace=False``
 for the first command that allows it. i.e ::
 
-    s2 = s.crop(4120, 4220, 'nm', inplace=False).apply_slit(3, 'nm').take('radiance')
+    s2 = s.crop(4120, 4220, 'nm', inplace=False).apply_slit(3, 'nm').offset(1.5, 'nm')
 
 
 .. _label_spectrum_remove_baseline:
@@ -684,6 +745,9 @@ Or remove a linear baseline with:
 - :py:func:`~radis.spectrum.operations.get_baseline`
 - :py:func:`~radis.spectrum.operations.sub_baseline`
 
+You could also use the functions available in :py:mod:`pyspecutils`,
+see :py:meth:`~radis.spectrum.spectrum.Spectrum.to_specutils`.
+
 
 Calculate transmittance from radiance with Kirchoff's law
 ---------------------------------------------------------

environment.yml

@@ -6,35 +6,35 @@ channels:
 - cantera
 dependencies:
 - python=3.8
-- numpy
-- scipy>=1.4.0
-- matplotlib
-- seaborn   # other matplotlib themes
-- cantera   # for chemical equilibrium computations
-- cython
-- pandas>=1.0.5
-- pytables  # for pandas to HDF5 export
-- plotly    # used in line-survey
 - astropy  # Unit aware calculations
-- plotly>=2.5.1  # for line survey HTML output
-- termcolor  # terminal colors
-- configparser
 - astroquery>=0.3.9  # to fetch HITRAN databases
-- pytest  # to run test suite
+- cantera   # for chemical equilibrium computations
+- configparser
+- cython
+- h5py   # load HDF5
 - joblib  # for parallel loading of SpecDatabase
+- matplotlib
+- numpy
 - numba  # just-in-time compiler
+- pandas>=1.0.5
+- plotly>=2.5.1  # for line survey HTML output
 - psutil # for getting user RAM
-- h5py   # load HDF5
+- scipy>=1.4.0
+- seaborn   # other matplotlib themes
+- termcolor  # terminal colors
+- specutils
+- vaex>=4.4.0   # load HDF5 files  (version for custom HDF5 groups)
 - pip
 - pip:
   - beautifulsoup4 # parse ExoMol website
   - lxml           # parser used for ExoMol website
-  - pyarrow       # for the feather format (temporarily needed for ExoMol)
+  - hjson          # Json with comments (for default_radis.json)
   - publib>=0.3.2  # Plotting styles for Matplotlib
   - hitran-api     # HAPI, used to access TIPS partition functions
+  - peakutils
   - json-tricks>=3.15.0  # to deal with non jsonable formats
   - mpldatacursor
+  - progressbar2    # used in vaex
   - tuna           # to generate visual/interactive performance profiles
   - tables  # for pandas to HDF5 export
-  - vaex    # load HDF5 files
   - habanero  # CrossRef API to retrieve data from doi

examples/get_molecular_parameters.py

@@ -0,0 +1,37 @@
+# -*- coding: utf-8 -*-
+"""
+========================
+Get Molecular Parameters
+========================
+
+Retrieve isotopologue abundances or molecular mass, based on HITRAN data,
+using :py:class:`~radis.db.molparam.MolParams`
+
+See also how to use custom (non-terrestrial) abundances in
+:ref:`the Custom Abundances example <example_custom_abundances>`
+
+"""
+
+#%% Access the parameters:
+from radis.db.molparam import MolParams
+
+molpar = MolParams()
+
+print("CO2 abundances for the first 2 isotopes :")
+print(molpar.get("CO2", 1, "abundance"))  # 1 for isotopologue number
+print(molpar.get("CO2", 2, "abundance"))  # 1 for isotopologue number
+
+
+print("H2O molar mass for the first 2 isotopes")
+print(molpar.get("H2O", 1, "mol_mass"))  # 1 for isotopologue number
+print(molpar.get("CO2", 2, "mol_mass"))  # 1 for isotopologue number
+
+
+print("Tabulated CO partition function at 296 K for the first 2 isotopes")
+print(molpar.get("CO", 1, "Q_296K"))  # 1 for isotopologue number
+print(molpar.get("CO", 2, "Q_296K"))  # 1 for isotopologue number
+
+
+#%% All parameters are stored in a DataFrame:
+print("All parameters: ", list(molpar.df.columns))
+print(molpar.df)

examples/plot_SpecDatabase.py

@@ -53,7 +53,7 @@
 
 # Multiple Spectrum calculation based on different values of Tgas and wstep
 for i in wstep:
-    sf.wstep = i
+    sf._wstep = i
     sf.params.wstep = i
     for j in Tgas:
         sf.eq_spectrum(Tgas=j, path_length=1)

examples/plot_compare_CO_exomol_hitemp.py

@@ -0,0 +1,54 @@
+# -*- coding: utf-8 -*-
+"""
+==============================================
+Compare CO from the ExoMol and HITEMP database
+==============================================
+
+Auto-download and calculate CO spectrum from the HITEMP database, and the
+ExoMol database. ExoMol references multiple databases for CO. Here we do not
+use the ExoMol recommended database (see :py:func:`~radis.io.exomol.get_exomol_database_list`)
+but we use the HITEMP database hosted on ExoMol.
+
+Output should be similar, but no ! By default these two databases provide
+different broadening coefficients. However, the Einstein coefficients & linestrengths
+should be the same, therefore the integrals under the lines should be similar.
+
+We verifiy this below :
+
+For further exploratino, ExoMol and HITEMP lines can be downloaded and accessed separately using
+:py:func:`~radis.io.exomol.fetch_exomol` and :py:func:`~radis.io.hitemp.fetch_hitemp`
+
+"""
+import astropy.units as u
+
+from radis import calc_spectrum, plot_diff
+
+conditions = {
+    "wmin": 2002 / u.cm,
+    "wmax": 2300 / u.cm,
+    "molecule": "CO",
+    "isotope": "2",
+    "pressure": 1.01325,  # bar
+    "Tgas": 1000,  # K
+    "mole_fraction": 0.1,
+    "path_length": 1,  # cm
+    "broadening_method": "fft",  # @ dev: Doesn't work with 'voigt'
+    "verbose": True,
+}
+
+s_exomol = calc_spectrum(
+    **conditions, databank="exomol", name="ExoMol's HITEMP (H2 broadened)"
+)
+s_hitemp = calc_spectrum(
+    **conditions,
+    databank="hitemp",
+    name="HITEMP (Air broadened)",
+)
+plot_diff(s_exomol, s_hitemp, "abscoeff")
+
+#%% Broadening coefficients are different but areas under the lines should be the same :
+import numpy as np
+
+assert np.isclose(
+    s_exomol.get_integral("abscoeff"), s_hitemp.get_integral("abscoeff"), rtol=0.001
+)

examples/plot_custom_abundances.py

@@ -1,5 +1,7 @@
 # -*- coding: utf-8 -*-
 """
+.. _example_custom_abundances:
+
 =====================
 Use Custom Abundances
 =====================

examples/plot_exomol_spectrum.py

@@ -6,6 +6,9 @@
 
 Auto-download and compute a SiO spectrum from the ExoMol database ([ExoMol-2020]_)
 
+ExoMol lines can be downloaded and accessed separately using
+:py:func:`~radis.io.exomol.fetch_exomol`
+
 """
 
 from radis import calc_spectrum
@@ -20,7 +23,24 @@
     mole_fraction=0.1,
     path_length=1,  # cm
     broadening_method="fft",  # @ dev: Doesn't work with 'voigt'
-    databank="exomol",  # uses the recommended database. Use ('exomol', "EBJT") for a specific database ("EBJT")
+    databank=("exomol", "EBJT"),  # Simply use 'exomol' for the recommended database
 )
 s.apply_slit(1, "cm-1")  # simulate an experimental slit
 s.plot("radiance")
+
+
+#%% See line data:
+from radis.io.exomol import fetch_exomol
+
+df = fetch_exomol("SiO", database="EBJT", isotope="1", load_wavenum_max=5000)
+print(df)
+
+
+#%% See the list of recommended databases for the 1st isotope of SiO :
+from radis.io.exomol import get_exomol_database_list, get_exomol_full_isotope_name
+
+databases, recommended = get_exomol_database_list(
+    "SiO", get_exomol_full_isotope_name("SiO", 1)
+)
+print("Databases for SiO: ", databases)
+print("Database recommended by ExoMol: ", recommended)