# -*- coding: utf-8 -*- """ Summary ------- :class:`~radis.spectrum.spectrum.Spectrum` class holder Contains both the emission and absorption features calculated by a SpectrumFactory, or generated by another spectral code or an experiment and imported in RADIS. Allows use post-processing functions such as applying an instrumental slit, saving to disk, plot all spectral quantities with arbitrary units. See the :ref:`The Spectrum object ` for the list of all post-processing functions, how to save a Spectrum, or how to generate a Spectrum from text. Examples -------- Typical use:: from radis calculated_spectrum s = calculated_spectrum(w, I, conditions={'case':'previously calculated by ##'}) s.plot('radiance_noslit') s.apply_slit(0.5, shape='triangular') s.plot('radiance') Spectrum objects can be modified, stored, resampled, rescaled, or retrieved after they have been created, with :py:meth:`~radis.spectrum.spectrum.Spectrum.store`, :py:meth:`~radis.spectrum.spectrum.Spectrum.rescale_path_length`, :py:meth:`~radis.spectrum.spectrum.Spectrum.rescale_mole_fraction`, :py:meth:`~radis.spectrum.spectrum.Spectrum.resample`, :py:meth:`~radis.spectrum.spectrum.Spectrum.store`, :py:func:`~radis.tools.database.load_spec` :: from radis import load_spec s = load_spec('co_calculation.spec') s.rescale_path_length(0.5) # calculate for new path_length s.rescale_mole_fraction(0.02) # calculate for new mole fraction s.resample(w_new) # resample on new wavespace s.store('co_calculation2.spec') More in :ref:`The Spectrum object `. """ from copy import deepcopy from os.path import basename from warnings import warn import numpy as np import pandas as pd from numpy import abs, diff from radis.db.references import doi # from radis.lbl.base import print_conditions from radis.misc.arrays import ( anynan, count_nans, evenly_distributed, first_nonnan_index, last_nonnan_index, nantrapz, ) from radis.misc.config import get_config from radis.misc.debug import printdbg from radis.misc.plot import split_and_plot_by_parts from radis.misc.profiler import Profiler from radis.misc.signal import resample, resample_even from radis.misc.warning import GPUInitWarning from radis.phys.air import air2vacuum, vacuum2air from radis.phys.convert import cm2nm, conv2, nm2cm from radis.phys.units import Unit, convert_universal from radis.spectrum.rescale import rescale_mole_fraction, rescale_path_length, update from radis.spectrum.utils import ( CONVOLUTED_QUANTITIES, WAVELEN_UNITS, WAVELENVAC_UNITS, WAVENUM_UNITS, cast_waveunit, format_xlabel, make_up, make_up_unit, print_conditions, ) from radis.tools.track_ref import RefTracker # %% Spectrum class to hold results def _is_running_in_notebook(): """Check if the code is running in a Jupyter Notebook.""" try: shell = get_ipython().__class__.__name__ if shell == "ZMQInteractiveShell": return True # Jupyter notebook or qtconsole elif shell == "TerminalInteractiveShell": return False # Terminal running IPython else: return False # Spyder except NameError: return False # Probably standard Python interpreter class Spectrum(object): """This class holds results calculated with the :py:class:`~radis.lbl.factory.SpectrumFactory` calculation, with other radiative codes, or experimental data. It can be used to plot different quantities a posteriori, or manipulate output units (for instance convert a spectral radiance per wavelength units to a spectral radiance per wavenumber). See more information on how to generate, edit or combine Spectrum objects on :ref:`the Spectrum object guide `. Parameters ---------- quantities: dict of tuples ``{'quantity':(w, a)}`` or dict ``{'wavelength/wavenumber': w, quantity': a}`` where quantities are spectral quantities (absorbance, radiance, etc.) and wavenum is in :math:`cm^{-1}` or :math:`nm` (see ``waveunit``) Example:: # w, k, I are numpy arrays for wavenumbers, absorption coefficient, and radiance. from radis import Spectrum s = Spectrum({"wavenumber":w, "abscoeff":k, "radiance_noslit":I}, wunit='cm-1', units={"radiance_noslit":"mW/cm2/sr/nm", "abscoeff":"cm-1"}) Or:: s = Spectrum({"abscoeff":(w,k), "radiance_noslit":(w,I)}, wunit="cm-1" units={"radiance_noslit":"mW/cm2/sr/nm", "abscoeff":"cm-1"}) See also: :py:meth:`~radis.spectrum.spectrum.Spectrum.from_array` and :py:meth:`~radis.spectrum.spectrum.Spectrum.from_txt` units: dict units for quantities Other Parameters ---------------- conditions: dict physical conditions and calculation parameters wunit: ``'nm'``, ``'cm-1'``, ``'nm_vac'`` or ``None`` wavelength in air (``'nm'``), wavenumber (``'cm-1'``), or wavelength in vacuum (``'nm_vac'``). If ``None``, ``'wavespace'`` must be defined in ``conditions``. Quantities should be evenly distributed along this space for fast convolution with the slit function cond_units: dict units for conditions Other Parameters ---------------- name: str, or None Give a name to this Spectrum object (automatically used in plots; useful for multislab configurations). Default ``None`` populations: dict a dictionary of all species, and levels. Should be compatible with other radiative codes such as Specair output. Suggested format: {molecules: {isotopes: {elec state: rovib levels}}} e.g:: {'CO2':{1: 'X': df}} # with df a Pandas Dataframe lines: pandas Dataframe all lines in databank (necessary for using :meth:`~radis.spectrum.spectrum.Spectrum.line_survey`). Warning if you want to play with the lines content: The signification of columns in `lines` may be specific to a database format. Plus, some additional columns may have been added by the calculation (e.g: `Ei` and `S` for emission integral and linestrength in SpectrumFactory). Refer to the code to know what they mean (and their units) gpu_app: gpuApp specifies the gpuApp object that was initialized. this value is auto-generated by SpectrumFactory; don't provide a manual input. Only not None when spectrum object is produced with one of the GPU methods :meth:`~radis.lbl.factory.SpectrumFactory.eq_spectrum_gpu` or :meth:`~radis.lbl.factory.SpectrumFactory.eq_spectrum_gpu_interactive`. references: dict a dict of ``doi`` of references used to compute this object. Automatically returned with the full bibtex entry by :py:meth:`~radis.spectrum.spectrum.Spectrum.cite()` It can also be set a posteriori. Example :: s = Spectrum() s.references = {"10.1016/j.jqsrt.2010.05.001": "HITEMP-2010 database", "10.1016/j.jqsrt.2018.09.027":["calculation", "post-processing"], # RADIS main paper. Automatically added "10.1016/j.jqsrt.2020.107476":"DIT algorithm"} ) s.cite() .. raw:: html
Returns : .. code-block:: python Used for DIT algorithm ---------------------- @article{van_den_Bekerom_2021, doi = {10.1016/j.jqsrt.2020.107476}, url = {https://doi.org/10.1016%2Fj.jqsrt.2020.107476}, year = 2021, month = {mar}, publisher = {Elsevier {BV}}, volume = {261}, pages = {107476}, author = {D.C.M. van den Bekerom and E. Pannier}, title = {A discrete integral transform for rapid spectral synthesis}, journal = {Journal of Quantitative Spectroscopy and Radiative Transfer} } Used for HITEMP-2010 database ----------------------------- @article{Rothman_2010, doi = {10.1016/j.jqsrt.2010.05.001}, url = {https://doi.org/10.1016%2Fj.jqsrt.2010.05.001}, year = 2010, month = {oct}, publisher = {Elsevier {BV}}, volume = {111}, number = {15}, pages = {2139--2150}, author = {L.S. Rothman and I.E. Gordon and R.J. Barber and H. Dothe and R.R. Gamache and A. Goldman and V.I. Perevalov and S.A. Tashkun and J. Tennyson}, title = {{HITEMP}, the high-temperature molecular spectroscopic database}, journal = {Journal of Quantitative Spectroscopy and Radiative Transfer} } Used for calculation, post-processing ------------------------------------- @article{Pannier_2019, doi = {10.1016/j.jqsrt.2018.09.027}, url = {https://doi.org/10.1016%2Fj.jqsrt.2018.09.027}, year = 2019, month = {jan}, publisher = {Elsevier {BV}}, volume = {222-223}, pages = {12--25}, author = {Erwan Pannier and Christophe O. Laux}, title = {{RADIS}: A nonequilibrium line-by-line radiative code for {CO}2 and {HITRAN}-like database species}, journal = {Journal of Quantitative Spectroscopy and Radiative Transfer} } .. raw:: html
warnings: boolean if ``True``, test if inputs are valid, e.g, spectra are evenly distributed in wavelength, and raise a warning if not. Note that this take ~ 3.5 ms for a 20k points spectrum, when the rest of the creation process is only ~ 1.8ms (makes it 3 times longer, and can be a problem if hundreds of spectra are created in a row). Default ``True`` Examples -------- Manipulate a Spectrum calculated by RADIS:: s = calc_spectrum(2125, 2300, Tgas=2000, databank='CDSD') s.print_conditions() s.plot('absorbance') s.line_survey(overlay='absorbance') s.plot('radiance_noslit', wunits='cm-1', Iunits='W/m2/sr/cm-1') s.apply_slit(5) s.plot('radiance') w, t = s.get('transmittance_noslit') # for use in multi-slabs configs Any tuple of numpy arrays (w, I) can also be converted into a Spectrum object from the :class:`~radis.spectrum.spectrum.Spectrum` class directly, or using the :func:`~radis.spectrum.models.calculated_spectrum` function. All the following methods are equivalent:: from radis import Spectrum, calculated_spectrum s1 = calculated_spectrum(w, I, wunit='nm', Iunit='mW/cm2/sr/nm') s2 = Spectrum.from_array(w, I, 'radiance_noslit', wunit='nm', unit='mW/cm2/sr/nm') s3 = Spectrum({'radiance_noslit': (w, I)}, units={'radiance_noslit':'mW/cm2/sr/nm'}, wunit='nm') See more examples in the [Spectrum]_ page. Spectrum objects can be stored, retrieved, rescaled, resampled:: from radis import load_spec s = load_spec('co_calculation.spec') s.rescale_path_length(0.5) # calculate for new path_length s.rescale_mole_fraction(0.02) # calculate for new mole fraction s.resample(w_new) # resample on new wavespace s.store('co_calculation2.spec') .. minigallery:: radis.spectrum.spectrum.Spectrum Notes ----- Implementation: quantities are stored in the ``self._q`` dictionary. They are better accessed with the :meth:`~radis.spectrum.spectrum.Spectrum.get` method that deals with units and wavespace Wavebase: quantities are stored either in wavenum or wavelength base, but this doesnt matter as they are retrieved / plotted with the :meth:`~radis.spectrum.spectrum.Spectrum.get` and :meth:`~radis.spectrum.spectrum.Spectrum.plot` methods which have units as input arguments Attributes ---------- conditions : dict Stores computation / measurement conditions c: dict convenience wrapper to ``conditions``:: s.c["calculation_time"] is s.conditions["calculation_time"] >> True populations: dict Stores molecules, isotopes, electronic states and vibrational or rovibrational populations See Also -------- :func:`~radis.spectrum.models.calculated_spectrum`, :func:`~radis.spectrum.models.transmittance_spectrum`, :func:`~radis.spectrum.models.experimental_spectrum` :meth:`~radis.spectrum.spectrum.Spectrum.from_array` :meth:`~radis.spectrum.spectrum.Spectrum.from_txt` :func:`~radis.tools.database.load_spec` References ---------- .. [Spectrum] See the :ref:`Spectrum object page ` """ # hardcode attribute names, but can save a lot of memory if hundreds of spectra __slots__ = [ "_q", "units", "c", "conditions", "cond_units", "gpu_app", "populations", "references", # doi entries "lines", "name", "_slit", "file", "profiler", ] def __init__( self, quantities, units=None, conditions=None, cond_units=None, gpu_app=None, populations=None, lines=None, wunit=None, name=None, references={}, check_wavespace=True, **kwargs, ): # Check inputs # --------------- # ... Replace None attributes with dictionaries if conditions is None: conditions = {} if units is None: units = {} if cond_units is None: cond_units = {} if populations is None: populations = {} if references is None: references = {} self._init_annotations() # typing hints for get() and plot() if "warnings" in kwargs: check_wavespace = kwargs.pop("warnings") warn( "Spectrum(warning=) renamed Spectrum(check_wavespace=) in 0.9.30", DeprecationWarning, ) # ... waveunit in Spectrum renamed wunit if "waveunit" in kwargs: assert wunit is None # not given twice warn( "waveunit parameter in Spectrum(waveunit=) is now named `wunit`", DeprecationWarning, ) wunit = kwargs.pop("waveunit") waveunit = wunit # (but not renamed within the code or conditions yet) if len(kwargs) > 0: raise ValueError(f"Unexpected input: {list(kwargs.keys())}") # Deal with deprecated inputs # ... wavespace renamed waveunit if "wavespace" in conditions: warn( "wavespace key in conditions is now named waveunit", DeprecationWarning ) conditions["waveunit"] = conditions["wavespace"] del conditions["wavespace"] # ... make sure waveunit is given. Either implicitly, or explicitly: if waveunit is None and not "waveunit" in conditions: if "wavelength" in quantities: raise AssertionError( "waveunit ('nm', 'nm_vac' ?) has to be defined in `conditions`" + "or with `wunit='nm'` or `wunit='nm_vac'`" ) elif "wavenumber" in quantities: raise AssertionError( "waveunit ('cm-1'?) has to be defined in `conditions`" + "or with `wunit='cm-1'`" ) else: raise AssertionError( "waveunit ('nm', 'cm-1'?) has to be defined in `conditions`" + "or with `wunit=`" ) # ... Cast in standard waveunit format if waveunit is not None: waveunit = cast_waveunit(waveunit) if "waveunit" in conditions: conditions["waveunit"] = cast_waveunit(conditions["waveunit"]) # ... Make sure unit match if "waveunit" in conditions: if waveunit is not None and conditions["waveunit"] != waveunit: raise ValueError( f"waveunit defined in conditions ({conditions['waveunit']}) and directly ({waveunit}) dont match" ) elif waveunit is not None: # ... or define them in dictionary conditions["waveunit"] = waveunit if "wavelength" in quantities: assert "wavenumber" not in quantities and "wavespace" not in quantities assert conditions["waveunit"] in WAVELEN_UNITS + WAVELENVAC_UNITS if "wavenumber" in quantities: assert "wavelength" not in quantities and "wavespace" not in quantities assert conditions["waveunit"] in WAVENUM_UNITS if "wavespace" in quantities: assert "wavenumber" not in quantities and "wavelength" not in quantities # Check quantities format if len(quantities) == 0: raise AssertionError( "Spectrum is created with no quantities. Add " + "`radiance`, `transmittance` etc... with dict " + "format. e.g: {`wavenumber`:w, `radiance`: I}" ) # Create the arrays # ----------------- self._q = {} #: dict: stores spectral arrays self._slit = {} #: dict: hold slit function # infer format: tuple_format = ( "wavelength" not in quantities and "wavenumber" not in quantities and "wavespace" not in quantities ) k0 = list(quantities.keys())[0] v0 = list(quantities.values())[0] for k, v in quantities.items(): if len(v) != len(v0): raise AssertionError( f"Input arrays should have the same length. Got : {k} : len {len(v)}, {k0} : len {len(v0)}" ) if tuple_format and len(v) != 2: raise AssertionError( "Input arrays should have format `{'quantity':(wavespace, array)} or {'wavelength/wavenumber':wavespace, 'quantity':array}` but not both. Got :" + f"{quantities}", ) if tuple_format: for k, (w, I) in quantities.items(): # creates and stores a copy of w,I self._add_quantity(k, w, I, check_wavespace=check_wavespace) else: for k, I in quantities.items(): if k in ["wavelength", "wavenumber", "wavespace"]: continue # add: if not "wavespace" in self._q: # adding wavenumber space if "wavespace" in quantities.keys(): w = quantities["wavespace"] elif "wavelength" in quantities.keys(): w = quantities["wavelength"] elif "wavenumber" in quantities.keys(): w = quantities["wavenumber"] # creates and stores a copy of w,I : self._add_quantity(k, w, I, check_wavespace=True) else: # this also checks that all arrays have same length self._add_quantity(k, w, I, check_wavespace=False) # Check condition units : for k, v in conditions.items(): if hasattr(v, "unit"): if k in cond_units: # assert units match assert v.unit == cond_units[k] else: # add in cond units cond_units[k] = v.unit # Finally, add our attributes self.gpu_app = gpu_app self.conditions = self.c = conditions """ dict: computation conditions, or experimental parameters, or any metadata you need to store with the Spectrum object. """ self.populations = populations """ populations of rovibrational levels used to compute partition functions of the spectrum See also :py:meth:`~radis.spectrum.spectrum.Spectrum.get_populations`, :py:meth:`~radis.spectrum.spectrum.Spectrum.plot_populations` .. minigallery:: radis.spectrum.spectrum.Spectrum.populations """ self.lines = lines """ informations on emitting or absorbing lines that contribute to the spectrum. See also :py:meth:`~radis.spectrum.spectrum.Spectrum.line_survey` .. minigallery:: radis.spectrum.spectrum.Spectrum.lines """ self.units = units """ dict: units for spectral quantities. """ self.cond_units = cond_units self.name = name self.file = None # used to store filename when loaded from a file self.profiler = None # Add references self.references = RefTracker(**references) if not doi["RADIS-2018"] in self.references: self.references.add( doi["RADIS-2018"], "post-processing" ) # Radis main paper (TODO: replace with Community paper when published ) # %% Constructors @classmethod def from_array( self, w, I, quantity, wunit=None, Iunit=None, waveunit=None, unit=None, *args, **kwargs, ): r"""Construct Spectrum from 2 arrays. Parameters ---------- w, I: array waverange and vector quantity: str spectral quantity name wunit: ``'nm'``, ``'cm-1'``, ``'nm_vac'`` unit of waverange: wavelength in air (``'nm'`` or ``'nm_air'``), wavenumber (``'cm-1'``), or wavelength in vacuum (``'nm_vac'``). If ``None``, then ``w`` must be a dimensioned array. Iunit: str spectral quantity unit (arbitrary). Ex: ``'mW/cm2/sr/nm'`` for radiance_noslit If ``None``, then ``I`` must be a dimensioned array. *args, **kwargs see :class:`~radis.spectrum.spectrum.Spectrum` doc Other Parameters ---------------- conditions: dict physical conditions and calculation parameters cond_units: dict units for conditions populations: dict a dictionary of all species, and levels. Should be compatible with other radiative codes such as Specair output. Suggested format: {molecules: {isotopes: {elec state: rovib levels}}} e.g:: {'CO2':{1: 'X': df}} # with df a Pandas Dataframe lines: pandas Dataframe all lines in databank (necessary for using :meth:`~radis.spectrum.spectrum.Spectrum.line_survey`). Warning if you want to play with the lines content: The signification of columns in `lines` may be specific to a database format. Plus, some additional columns may have been added by the calculation (e.g: `Ei` and `S` for emission integral and linestrength in SpectrumFactory). Refer to the code to know what they mean (and their units) Returns ------- Spectrum creates a :class:`~radis.spectrum.spectrum.Spectrum` object Examples -------- Create a spectrum:: from radis import Spectrum s = Spectrum.from_array(w, I, 'radiance_noslit', wunit='nm', unit='mW/cm2/sr/nm') Dimensioned arrays can also be used directly :: 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') To create a spectrum with absorption and emission components (e.g: ``radiance_noslit`` and ``transmittance_noslit``, or ``emisscoeff`` and ``abscoeff``) call the :class:`~radis.spectrum.spectrum.Spectrum` class directly. Ex:: from radis import Spectrum s = Spectrum({'abscoeff': (w, A), 'emisscoeff': (w, E)}, units={'abscoeff': 'cm-1', 'emisscoeff':'W/cm2/sr/nm'}, wunit='nm') .. minigallery:: radis.spectrum.spectrum.Spectrum.from_array See Also -------- :class:`~radis.spectrum.spectrum.Spectrum`, :func:`~radis.spectrum.models.calculated_spectrum`, :func:`~radis.spectrum.models.transmittance_spectrum`, :func:`~radis.spectrum.models.experimental_spectrum` :meth:`~radis.spectrum.spectrum.Spectrum.from_txt` :func:`~radis.tools.database.load_spec`, :ref:`the Spectrum page ` """ # Deprecated inputs if waveunit is not None: warn( "`waveunit=` parameter in from_array is now named `wunit=`", DeprecationWarning, ) wunit = waveunit if unit is not None: warn( "`unit=` parameter in from_array is now named `Iunit=`", DeprecationWarning, ) Iunit = unit # Check if dimensioned arrays try: I.unit except AttributeError: # dimensionless if Iunit is None: raise ValueError( "I must be a dimensioned array, or ``Iunit=`` must be given." ) else: if Iunit is not None: raise ValueError( f"I is a dimensioned array (in {Iunit.unit.to_string()}), therefore ``Iunit`` cannot be given too (got {Iunit}). Set ``Iunit=None``" ) Iunit = I.unit.to_string() I = I.value try: w.unit except AttributeError: # dimensionless if wunit is None: raise ValueError( "``w`` must be a dimensioned array, or ``wunit=`` must be given." ) else: wunit_input = wunit wunit = w.unit.to_string() if wunit not in WAVELEN_UNITS + WAVENUM_UNITS: # Convert to something we know, for instance 'nm' : if w.unit.is_equivalent("nm"): w = w.to("nm") wunit = "nm" elif w.unit.is_equivalent("1 / cm"): w = w.to("1 / cm") wunit = "cm-1" # else, an error will be raised anyway on Spectrum creation. w = w.value # Check units. # In particular, deal with confusions arising from "nm" being either "nm_vac" or "nm_air" # TODO: generalize to "µm" too if wunit_input is None: if wunit == "nm": warn( "Input wunit not given, risk of confusion : RADIS assumes wavelengths are given as seen in air, not vacuum. To remove this warning be explicit by writing `Spectrum.from_array(..., wunit='nm_air')`" ) else: pass else: # Check units are the same if wunit == "nm" and ( wunit_input == "nm_vac" or wunit_input == "nm_air" ): # extra precision given by the user, not an error. Use the input value wunit = wunit_input elif wunit != wunit_input: raise ValueError( f"``w`` is a dimensioned array (in {wunit}), therefore ``wunit`` should not be given ( `wunit=None`), or be the same (got `wunit={wunit_input}`). Set ``wunit=None``" ) else: pass quantities = {quantity: (w, I)} units = {quantity: Iunit} # Update Spectrum conditions conditions = kwargs.pop("conditions", {}) for k in ["path_length"]: # usual conditions if k in kwargs: conditions[k] = kwargs.pop(k) return self( quantities, units, wunit=wunit, conditions=conditions, *args, **kwargs ) @classmethod def from_txt(self, file, quantity, wunit, unit, waveunit=None, *args, **kwargs): r"""Construct Spectrum from txt file. Parameters ---------- file: str file name quantity: str spectral quantity name wunit: ``'nm'``, ``'cm-1'``, ``'nm_vac'`` unit of waverange: wavelength in air (``'nm'``), wavenumber (``'cm-1'``), or wavelength in vacuum (``'nm_vac'``). unit: str spectral quantity unit *args, **kwargs the following inputs are forwarded to loadtxt: ``'delimiter'``, ``'skiprows'`` The rest if forwarded to Spectrum. see :class:`~radis.spectrum.spectrum.Spectrum` doc Other Parameters ---------------- delimiter: ``','``, etc. see :py:func:`numpy.loadtxt` skiprows: int see :py:func:`numpy.loadtxt` argsort: bool sorts the arrays in ``file`` by wavespace. Convenient way to load a file where points have been manually added at the end. Default ``False``. **Optional Spectrum parameters:** conditions: dict physical conditions and calculation parameters cond_units: dict units for conditions populations: dict a dictionary of all species, and levels. Should be compatible with other radiative codes such as Specair output. Suggested format: {molecules: {isotopes: {elec state: rovib levels}}} e.g:: {'CO2':{1: 'X': df}} # with df a Pandas Dataframe lines: pandas Dataframe all lines in databank (necessary for using :meth:`~radis.spectrum.spectrum.Spectrum.line_survey`). Warning if you want to play with the lines content: The signification of columns in `lines` may be specific to a database format. Plus, some additional columns may have been added by the calculation (e.g: `Ei` and `S` for emission integral and linestrength in SpectrumFactory). Refer to the code to know what they mean (and their units) Returns ------- s: Spectrum creates a :class:`~radis.spectrum.spectrum.Spectrum` object Examples -------- Generate an experimental spectrum from txt. In that example the ``delimiter`` key is forwarded to :py:func:`~numpy.loadtxt`:: from radis import Spectrum s = Spectrum.from_txt('spectrum.csv', 'radiance', wunit='nm', unit='W/cm2/sr/nm', delimiter=',') To create a spectrum with absorption and emission components (e.g: ``radiance_noslit`` and ``transmittance_noslit``, or ``emisscoeff`` and ``abscoeff``) call the :class:`~radis.spectrum.spectrum.Spectrum` class directly. Ex:: from radis import Spectrum s = Spectrum({'abscoeff': (w, A), 'emisscoeff': (w, E)}, units={'abscoeff': 'cm-1', 'emisscoeff':'W/cm2/sr/nm'}, wunit='nm') .. minigallery:: radis.spectrum.spectrum.Spectrum.from_txt Notes ----- Internally, the numpy :py:func:`~numpy.loadtxt` function is used and transposed:: w, I = np.loadtxt(file).T You can use ``'delimiter'`` and '``skiprows'`` as arguments. See Also -------- :func:`~radis.spectrum.models.calculated_spectrum`, :func:`~radis.spectrum.models.transmittance_spectrum`, :func:`~radis.spectrum.models.experimental_spectrum`, :meth:`~radis.spectrum.spectrum.Spectrum.from_array`, :func:`~radis.tools.database.load_spec`, :ref:`the Spectrum page ` """ # Deprecated inputs if waveunit is not None: warn( "`waveunit=` parameter in from_array is now named `wunit=`", DeprecationWarning, ) wunit = waveunit # Get input for loadtxt kwloadtxt = {} for k in ["delimiter", "skiprows"]: if k in kwargs: kwloadtxt[k] = kwargs.pop(k) argsort = kwargs.pop("argsort", False) w, I = np.loadtxt(file, **kwloadtxt).T if argsort: b = np.argsort(w) w, I = w[b], I[b] quantities = {quantity: (w, I)} units = {quantity: unit} # Update Spectrum conditions conditions = kwargs.pop("conditions", {}) for k in ["path_length"]: # usual conditions if k in kwargs: conditions[k] = kwargs.pop(k) s = self(quantities, units, wunit=wunit, conditions=conditions, *args, **kwargs) # Store filename s.file = file return s @classmethod def from_hdf5( self, file, wmin=None, wmax=None, wunit=None, columns=None, engine="pytables" ): r"""Generates a Spectrum from an HDF5 file. Uses :py:func:`~radis.io.spec_hdf.hdf2spec` Other Parameters ---------------- wmin, wmax: float range of wmin, wmax to load , using ``wunit`` (if None, load everything) columns: list of str spectral arrays to load (if None, load everything) Examples -------- :: Spectrum.from_hdf5("rad_hdf.h5", wmin=2100, wmax=2200, columns=['abscoeff', 'emisscoeff']) .. minigallery:: radis.Spectrum.from_hdf5 See Also -------- :py:func:`~radis.spectrum.spectrum.Spectrum.to_hdf5` """ from radis.io.spec_hdf import hdf2spec s = hdf2spec( file, wmin=wmin, wmax=wmax, wunit=wunit, columns=columns, engine=engine ) # Store filename s.file = file return s @classmethod def from_spec(self, file): r"""Generates a Spectrum from a .spec [json] file. Uses :py:func:`~radis.tools.database.load_spec`""" from radis.tools.database import load_spec return load_spec(file) @classmethod def from_specutils(self, spectrum, var="radiance"): r"""Convert a :py:mod:`specutils` :py:class:`~specutils.spectra.spectrum1d.Spectrum1D` to a ``radis`` :py:class:`~radis.spectrum.spectrum.Spectrum` object. Parameters ---------- spectrum: a ``specutils`` :py:class:`specutils.spectra.spectrum1d.Spectrum1D` var: str spectral array, default ``"radiance"`` Examples -------- Taken from the Specutils website (https://specutils.readthedocs.io/en/stable/#getting-started-with-specutils) :: from astropy.io import fits from astropy import units as u from specutils import Spectrum1D f = fits.open('https://data.sdss.org/sas/dr16/sdss/spectro/redux/26/spectra/1323/spec-1323-52797-0012.fits') # The spectrum is in the second HDU of this file. specdata = f[1].data lamb = 10**specdata['loglam'] * u.AA flux = specdata['flux'] * 10**-17 * u.Unit('erg cm-2 s-1 AA-1') spec = Spectrum1D(spectral_axis=lamb, flux=flux) from radis import Spectrum s = Spectrum.from_specutils(spec) s.plot(wunit='nm') .. minigallery:: radis.Spectrum.from_specutils See Also -------- :py:func:`~radis.spectrum.spectrum.Spectrum.to_specutils` """ try: import specutils except ModuleNotFoundError as err: raise ModuleNotFoundError( "Specutils is required to use this function." "Install it with:\n" " conda install -c conda-forge specutils" "\nor\n" " pip install specutils" ) from err assert isinstance(spectrum, specutils.Spectrum1D) if "waveunit" in spectrum.meta and spectrum.meta["waveunit"] in WAVELEN_UNITS: raise ValueError( f"Specutils only handles wavelengths in vacuum. Expected `spectrum.meta['waveunit']` to be one of {WAVELENVAC_UNITS}, got {spectrum.meta['waveunit']}`" ) conditions = spectrum.meta.copy() if spectrum.spectral_axis.unit == "nm": # note : in Specutils wavelengths are given as wavelength in vacuum (1e7/X conversion # to wavenumbers) waverange = spectrum.wavelength waveunit = "nm_vac" if "waveunit" in conditions: conditions["waveunit"] = waveunit else: waverange = spectrum.spectral_axis waveunit = None # use default # Convert to radis Spectrum return Spectrum.from_array( waverange, spectrum.flux, quantity=var, wunit=waveunit, conditions=conditions, ) @classmethod def from_mat( self, file, quantity, wunit, unit, data_key=None, w_name="nu", I_name=None, index=None, *args, **kwargs, ): r"""Construct Spectrum from Matlab ``.mat`` file. Parameters ---------- file: str file name quantity: str spectral quantity name wunit: ``'nm'``, ``'cm-1'``, ``'nm_vac'`` unit of waverange: wavelength in air (``'nm'``), wavenumber (``'cm-1'``), or wavelength in vacuum (``'nm_vac'``). unit: str spectral quantity unit index: int index within Matlab ``.mat`` array. data_key: str data key to use within Matlab ``.mat`` array. If ``None``, guess. w_name, I_name: str key to use to parse the ``data[data_key]`` array and return waverange and quantity. *args, **kwargs the following inputs are forwarded to :py:func:`~scipy.io.loadmat`: ``'simplify_cells'``, ``'skiprows'`` The rest if forwarded to Spectrum and will be registered as a Spectrum condition. see :class:`~radis.spectrum.spectrum.Spectrum` doc Returns ------- s: Spectrum creates a :class:`~radis.spectrum.spectrum.Spectrum` object Examples -------- .. minigallery:: radis.Spectrum.from_mat Notes ----- Internally, the scipy :py:func:`~scipy.io.loadmat` function is used and corresponds to :: data = scipy.io.loadmat(file, **kwloadmat) data = data[data_key] w, I = data[w_name], data[I_name] Special keywords can be given to ``kwloadmat``. See docs of ``kwargs``. See Also -------- :meth:`~radis.spectrum.spectrum.Spectrum.from_array`, :func:`~radis.tools.database.load_spec` """ from scipy.io import loadmat # default loadmat reading arguments kwloadmat = {"simplify_cells": True} # override arguments for k in ["simplify_cells"]: if k in kwargs: kwloadmat[k] = kwargs.pop(k) data = loadmat(file, **kwloadmat) if data_key is None: guessed_key = [k for k in data.keys() if not k.startswith("_")] if len(guessed_key) == 1: data_key = guessed_key[0] else: raise ValueError( f"data_key not given and could not be guessed. Available keys : {guessed_key}. Choose one and set `data_key=`" ) data_array = data[data_key] try: w = data_array[w_name] except KeyError: raise KeyError( f"Could not find key {w_name} to parse waverange in data. Available keys : [{data_array.keys()}]" ) if I_name is None: guessed_I_name = [k for k in data_array.keys() if k != w_name] if len(guessed_I_name) == 1: I_name = guessed_I_name[0] else: raise ValueError( f"I_name not given and could not be guessed. Available keys : {guessed_I_name}. Choose one and set `I_name=`" ) try: I = data_array[I_name] except KeyError: raise KeyError( f"Could not find key {I_name} to parse quantity in data. Available keys : [{data_array.keys()}]" ) if index is not None: I = I[:, index] quantities = {quantity: (w, I)} units = {quantity: unit} # Update Spectrum conditions conditions = kwargs.copy() s = self(quantities, units, wunit=wunit, conditions=conditions) # Store filename s.file = file return s @classmethod def from_xsc( self, file, conditions={}, ): """Generates a Spectrum from a manually downloaded HITRAN cross-section file under `.xsc` format. Uses :py:func:`~radis.api.hitran.hitranapi.hitranxsc` Parameters ---------- file: str file name conditions: dictionary Conditions metadata added to the Spectrum object. see :py:class:`~radis.spectrum.spectrum.Spectrum` doc. Note that ``Tgas`` and ``pressure`` are automatically read from the cross-section file. """ from radis.api.hitranapi import hitranxsc data = hitranxsc(file) from radis.phys.units import Unit as u conditions.update( { "thermal_equilibrium": True, "pressure": data["P"] * u("Torr"), "Tgas": data["T"] * u("K"), "molecule": data["molec"], } ) return Spectrum.from_array( data["wavenumber"], data["spectrum"], name=data["name"], quantity="xsection", # cm2/molecule wunit="cm-1", Iunit="cm2/molecule", conditions=conditions, ) # Public functions # %% ====================================================================== # ---------------- # XXX ===================================================================== def get( self, var, wunit="default", Iunit="default", copy=True, trim_nan=False, return_units=False, ): r"""Retrieve a spectral quantity from a Spectrum object. You can select wavespace unit, intensity unit, or propagation medium. Parameters ---------- var: variable (``'absorbance'``, ``'transmittance'``, ``'xsection'`` etc.) Should be a defined quantity among :data:`~radis.spectrum.utils.CONVOLUTED_QUANTITIES` or :data:`~radis.spectrum.utils.NON_CONVOLUTED_QUANTITIES`. To get the full list of spectral arrays defined in this Spectrum object use the :meth:`~radis.spectrum.spectrum.Spectrum.get_vars` method. wunit: ``'nm'``, ``'cm'``, ``'nm_vac'``. wavespace unit: wavelength in air (``'nm'``), wavenumber (``'cm-1'``), or wavelength in vacuum (``'nm_vac'``). if ``"default"``, default unit for waveunit is used. See :py:meth:`~radis.spectrum.spectrum.Spectrum.get_waveunit`. Iunit: unit for variable ``var`` if ``"default"``, default unit for quantity `var` is used. See the :py:attr:`~radis.spectrum.spectrum.Spectrum.units` attribute. For ``var="radiance"``, one can use per wavelength (~ 'W/m2/sr/nm') or per wavenumber (~ 'W/m2/sr/cm-1') units .. note:: if using ``default``, the returned unit may be different from Spectrum.units[var]. e.g, for a Spectrum with radiance stored in 'mW/cm2/sr/nm', getting it with `wunit='cm-1', Iunit='default'` will return Iunit in 'mW/cm2/sr/cm-1' for consistency. Use ``return_units`` to be sure about which units are used. Other Parameters ---------------- copy: bool if ``True``, returns a copy of the stored quantity (modifying it wont change the Spectrum object). Default ``True``. trim_nan: bool if ``True``, removes ``nan`` on the sides of the spectral array (and corresponding wavespace). Default ``False``. return_units: bool if ``True``, return dimensioned Astropy arrays. Is using ``return_units = 'as_str'``, return `wunit` and `Iunit` as extra variables, i.e, ``w, I, wunit, Iunit = s.get(..., return_units='as_str')`` Default ``False`` Returns ------- w, I: array-like wavespace, quantity (ex: wavelength, radiance_noslit). For numpy users, note that these are copies (values) of the Spectrum quantity and not a view (reference): if you modify them the Spectrum is not changed Examples -------- Get transmittance in cm-1:: w, I = s.get('transmittance_noslit', wunit='cm-1') Get radiance (in wavelength in air):: _, R = s.get('radiance_noslit', wunit='nm', Iunit='W/cm2/sr/nm') Use with `return_units` to get dimensioned Astropy Quantities :: w, R = s.get('radiance_noslit', return_units=True) See Also -------- :meth:`~radis.spectrum.spectrum.Spectrum.get_radiance`, :meth:`~radis.spectrum.spectrum.Spectrum.get_radiance_noslit`, :ref:`the Spectrum page ` """ # check input if not var in self.get_vars(): if var + "_noslit" in self.get_vars(): raise KeyError( (f"`{var}` doesnt exist but `{var + '_noslit'}` does") + ". Have you used .apply_slit()?" ) else: # Try to compute it automatically : try: self.update(var, verbose=False) except ValueError as err: raise ValueError( f"{var} not in Spectrum arrays {self.get_vars()}. An error occurred while trying to recompute it from the available arrays and conditions. See above" ) from err # Get quantity I = self._q[var] # Copy if copy: I = I.copy() # Get wavespace (in correct unit, and correct medium) if wunit == "default": if "default_output_unit" in self.conditions: wunit = self.c["default_output_unit"] else: wunit = self.get_waveunit() # Handle explicit nm_air before cast_waveunit collapses it to "nm" force_air = wunit in ("nm_air",) wunit = cast_waveunit(wunit) if wunit == "cm-1": w = self.get_wavenumber(copy=copy) elif wunit == "nm": if force_air: medium = "air" else: medium = self.conditions.get("medium", "air") w = self.get_wavelength(medium=medium, copy=copy) elif wunit == "nm_vac": w = self.get_wavelength(medium="vacuum", copy=copy) else: raise ValueError(wunit) # Trim nan if needed if trim_nan: m = first_nonnan_index(I) M = last_nonnan_index(I) if m is None: raise ValueError("All values are nan. Check your data?") if m > 0 or M < len(I) - 1: # else, no change to be made w = w[m : M + 1] I = I[m : M + 1] # Convert y unit if necessary Iunit0 = self.units[var] if Iunit == "default": # Add unit consistent with the namespace # https://github.com/radis/radis/issues/456 if wunit in ["nm", "nm_air", "nm_vac"]: if var in [ "radiance", "radiance_noslit", "emisscoeff", "emisscoeff_continuum", ]: Iunit = Iunit0.replace("/cm-1", "/nm") else: Iunit = Iunit0 elif wunit in ["cm-1"]: if var in [ "radiance", "radiance_noslit", "emisscoeff", "emisscoeff_continuum", ]: Iunit = Iunit0.replace("/nm", "/cm-1") else: Iunit = Iunit0 else: raise ValueError(wunit) # Retrieve data (with correct unit) if Iunit != "default" and Iunit != Iunit0: if var in ["radiance", "radiance_noslit"]: # deal with the case where we want to get a radiance in per # wavelength unit (~ W/sr/cm2/nm) in wavenumber units (~ W/sr/cm2/cm-1), # or the other way round w_cm = self.get_wavenumber(copy=False) if trim_nan: w_cm = w_cm[m : M + 1] I = convert_universal( I, Iunit0, Iunit, w_cm, per_nm_is_like="mW/cm2/sr/nm", per_cm_is_like="mW/cm2/sr/cm-1", ) elif var in ["emisscoeff"]: # idem for emisscoeff in (~ W/sr/cm3/nm) or (~ W/sr/cm3/cm-1) w_cm = self.get_wavenumber(copy=False) if trim_nan: w_cm = w_cm[m : M + 1] I = convert_universal( I, Iunit0, Iunit, w_cm, per_nm_is_like="mW/cm3/sr/nm", per_cm_is_like="mW/cm3/sr/cm-1", ) elif var in ["absorbance"]: # no unit assert Iunit in ["", "1"] # dont change the variable: I has no dimension elif var in ["transmittance"]: # no unit assert Iunit in ["", "1"] # dont change the variable: I has no dimension else: I = conv2(I, Iunit0, Iunit) else: pass if return_units == "as_str": # only used internally return w, I, wunit, Iunit elif return_units: if wunit == "nm_vac": wunit = "nm" # fix for Astropy units return w * Unit(wunit), I * Unit(Iunit) else: return w, I def _get_wavespace(self, copy=True): r"""Return wavespace (if the same for all quantities) Other Parameters ---------------- copy: boolean if ``True``, returns a copy of the stored waverange (modifying it wont change the Spectrum object). Default ``True``. Returns ------- w: array_like (a copy of) spectrum wavespace for convoluted or non convoluted quantities """ w = self._q["wavespace"] if copy: w = w.copy() return w def get_wavelength(self, medium="air", copy=True): r"""Return wavelength in defined medium. Parameters ---------- medium: ``'air'``, ``'vacuum'`` returns wavelength as seen in air, or vacuum. Default ``'air'``. See :func:`~radis.phys.air.vacuum2air`, :func:`~radis.phys.air.air2vacuum` Other Parameters ---------------- copy: boolean if ``True``, returns a copy of the stored waverange (modifying it wont change the Spectrum object). Default ``True``. Returns ------- w: array_like (a copy of) spectrum wavelength for convoluted or non convoluted quantities See Also -------- :ref:`the Spectrum page ` """ # Check input if not medium in ["air", "vacuum"]: raise NotImplementedError(f"Unknown propagating medium: {medium}") # Now convert stored wavespace to the output unit w = self._get_wavespace(copy=copy) if self.get_waveunit() == "cm-1": w = cm2nm(w) # vacuum wavelength # Correct for propagation medium (air, vacuum) if medium == "air": w = vacuum2air(w) else: # medium == 'vacuum' pass # no change needed elif self.get_waveunit() == "nm": # nm air if medium == "air": pass # no change needed else: w = air2vacuum(w) elif self.get_waveunit() == "nm_vac": # nm vacuum if medium == "air": w = vacuum2air(w) else: pass # no change needed else: raise ValueError(self.get_waveunit()) return w def get_wavenumber(self, copy=True): r"""Return wavenumber (if the same for all quantities) Other Parameters ---------------- copy: boolean if ``True``, returns a copy of the stored waverange (modifying it wont change the Spectrum object). Default ``True``. Returns ------- w: array_like (a copy of) spectrum wavenumber for convoluted or non convoluted quantities """ w = self._get_wavespace(copy=copy) if self.get_waveunit() == "cm-1": # pass elif self.get_waveunit() == "nm": # wavelength air w = air2vacuum(w) w = nm2cm(w) elif self.get_waveunit() == "nm_vac": # wavelength vacuum w = nm2cm(w) else: raise ValueError(self.get_waveunit()) return w def get_radiance(self, Iunit="mW/cm2/sr/nm", copy=True): r"""Return radiance in whatever unit, and can even convert from ~1/nm to ~1/cm-1 (and the other way round) Other Parameters ---------------- copy: boolean if ``True``, returns a copy of the stored waverange (modifying it wont change the Spectrum object). Default ``True``. See Also -------- :meth:`~radis.spectrum.spectrum.Spectrum.get`, :meth:`~radis.spectrum.spectrum.Spectrum.get_radiance_noslit`, :ref:`the Spectrum page ` """ return self.get("radiance", Iunit=Iunit, copy=copy)[1] def get_radiance_noslit(self, Iunit="mW/cm2/sr/nm", copy=True): r"""Return radiance (non convoluted) in whatever unit, and can even convert from ~1/nm to ~1/cm-1 (and the other way round) Other Parameters ---------------- copy: boolean if ``True``, returns a copy of the stored waverange (modifying it wont change the Spectrum object). Default ``True``. See Also -------- :meth:`~radis.spectrum.spectrum.Spectrum.get`, :meth:`~radis.spectrum.spectrum.Spectrum.get_radiance`, :ref:`the Spectrum page ` """ return self.get("radiance_noslit", Iunit=Iunit, copy=copy)[1] def get_name(self): r"""Return Spectrum name. If not defined, returns either the :attr:`~radis.spectrum.spectrum.Spectrum.file` name if Spectrum was loaded from a file, or the ``'spectrum{id}'`` with the Python ``id`` object """ try: self.name except AttributeError: warn( DeprecationWarning( "Spectrum has no .name attribute and is probably outdated. Update!" ) ) self.name = None if self.name is not None: name = self.name elif self.file is not None: name = f"{basename(self.file)}" else: name_params = [] for key, unit in [("species", ""), ("dbformat", ""), ("Tgas", "K")]: if key in self.conditions: name_params.append(f"{self.conditions[key]}{unit}") name = "-".join(name_params) + "-#" + str(id(self))[-4:] return name def savetxt(self, filename, var, wunit="default", Iunit="default"): r"""Export spectral quantity var to filename. (note that by doing this you will loose additional information, such as the calculation conditions or the units. You better save a Spectrum object under a .spec file with :py:meth:`~radis.spectrum.spectrum.Spectrum.store` and load it afterwards with :py:func:`~radis.tools.database.load_spec`) Parameters ---------- filename: str file name var: str which spectral variable ot export Other Parameters ---------------- wunit, Iunit, medium: str see :meth:`~radis.spectrum.spectrum.Spectrum.get` for more information Notes ----- Export variable as:: np.savetxt(filename, np.vstack(self.get(var, wunit=wunit, Iunit=Iunit, medium=medium)).T, header=header) See Also -------- :meth:`~radis.spectrum.spectrum.Spectrum.store`, :meth:`~radis.spectrum.spectrum.Spectrum.save`, :ref:`the Spectrum page ` """ # Get units to export if wunit == "default": wunit = self.get_waveunit() wunit = cast_waveunit(wunit) if wunit == "cm-1": xlabel = "Wavenumber (cm-1)" elif wunit == "nm": xlabel = "Wavelength [air] (nm)" elif wunit == "nm_vac": xlabel = "Wavelength [vacuum] (nm)" else: raise ValueError(wunit) if Iunit == "default": try: yunit = self.units[var] except KeyError: # unit not defined in dictionary yunit = "a.u" else: yunit = Iunit header = f"{xlabel}\t{var} ({yunit})" np.savetxt( filename, np.vstack(self.get(var, wunit=wunit, Iunit=Iunit)).T, header=header, ) def update(self, quantity="all", optically_thin=None, verbose=True): r"""Calculate missing quantities: ex: if path_length and emisscoeff are given, recalculate radiance_noslit. Parameters ---------- quantity: str name of the spectral quantity to recompute. If 'same', only the quantities in the Spectrum are recomputed. If 'all', then all quantities that can be derived are recomputed. Default 'all'. Examples -------- Initialize a spectrum from the absorption coefficient, retrieve the transmittance or the emission coefficient : :: s = Spectrum.from_array([1900. , 1900.01, 1900.02, 1900.03, 1900.04, 1900.05, 1900.06, 1900.07, 1900.08, 1900.09], [2.71414065e-06, 2.88341489e-06, 3.06942277e-06, 3.27445689e-06, 3.50121831e-06, 3.75290756e-06, 4.03334037e-06, 4.34709612e-06, 4.69971017e-06, 5.09792551e-06], 'abscoeff', wunit='cm-1', Iunit='cm-1', conditions={'path_length':1, # cm 'thermal_equilibrium':True, 'Tgas':700, # K }) s.update('transmittance_noslit') s.update('emisscoeff') See Also -------- :ref:`the Spectrum page ` Notes ----- To compute radiance in the optically thin approximation (without self-absorption), set ``self.conditions['self_absorption'] = False`` before calling this method. """ return update( self, quantity=quantity, optically_thin=optically_thin, verbose=verbose ) # Rescale functions def rescale_path_length( self, new_path_length, old_path_length=None, inplace=True, force=False ): r"""Rescale spectrum to new path length. Starts from absorption coefficient and emission coefficient, and solves the RTE again for the new path length Convoluted values (with slit) are dropped in the process. Parameters ---------- new_path_length: float new path length old_path_length: float, or None if None, current path length (conditions['path_length']) is used Other Parameters ---------------- inplace: boolean if ``True``, modifies the Spectrum object directly. Else, returns a copy. Default ``True``. force: boolean if False, won't allow rescaling to 0 (not to loose information). Default ``False`` Returns ------- s: Spectrum Cropped Spectrum. If ``inplace=True``, Spectrum has been updated directly anyway. Allows :ref:`chaining ` Examples -------- :: for path in [0.1, 10, 100]: s.rescale_path_length(10, inplace=False).plot(nfig='same') Additionally, you can also use astropy units in the input arguments, for example: :: # preparing a test spectrum : import radis s = radis.spectrum_test() # rescaling : import astropy.units as u s.rescale_path_length(1 * u.km).plot() .. minigallery:: radis.Spectrum.rescale_path_length Notes ----- Implementation: To deal with all the input cases, we first make a list of what has to be recomputed, and what has to be recalculated See Also -------- :ref:`the Spectrum page ` :ref:`the provided example ` """ return rescale_path_length( self, new_path_length=new_path_length, old_path_length=old_path_length, inplace=inplace, force=force, ) def rescale_mole_fraction( self, new_mole_fraction, old_mole_fraction=None, inplace=True, ignore_warnings=False, force=False, verbose=True, ): r"""Update spectrum with new molar fraction Convoluted values (with slit) are dropped in the process. Parameters ---------- new_mole_fraction: float new mole fraction old_mole_fraction: float, or None if None, current mole fraction (conditions['mole_fraction']) is used Other Parameters ---------------- inplace: boolean if ``True``, modifies the Spectrum object directly. Else, returns a copy. Default ``True``. force: boolean if False, won't allow rescaling to 0 (not to loose information). Default ``False`` Returns ------- s: Spectrum Cropped Spectrum. If ``inplace=True``, Spectrum has been updated directly anyway. Allows :ref:`chaining ` Examples -------- :: s.rescale_mole_fraction(0.2) .. minigallery:: radis.Spectrum.rescale_mole_fraction Notes ----- Implementation: similar to rescale_path_length() but we have to scale abscoeff & emisscoeff Note that this is valid only for small changes in mole fractions. Then, the change in line broadening becomes significant See Also -------- :ref:`the Spectrum page ` """ return rescale_mole_fraction( self, new_mole_fraction=new_mole_fraction, old_mole_fraction=old_mole_fraction, inplace=inplace, ignore_warnings=ignore_warnings, force=force, verbose=verbose, ) def crop(self, wmin=None, wmax=None, wunit="default", inplace=True): r"""Crop spectrum to ``wmin-wmax`` range in ``wunit`` (inplace) Parameters ---------- wmin, wmax: float, or None boundaries of spectral range (in ``wunit``) wunit: ``'nm'``, ``'cm-1'``, ``'nm_vac'`` which waveunit to use for ``wmin, wmax``. If ``default``: use the default Spectrum wavespace defined with :meth:`~radis.spectrum.spectrum.Spectrum.get_waveunit`. Other Parameters ---------------- inplace: bool if ``True``, modifies the Spectrum object directly. Else, returns a copy. Default ``True``. Returns ------- s: Spectrum Cropped Spectrum. If ``inplace=True``, Spectrum has been updated directly anyway. Allows :ref:`chaining ` Examples -------- Crop to experimental Spectrum, and compare:: from radis import calc_spectrum, load_spec, plot_diff s = calc_spectrum(...) s_exp = load_spec('typical_result.spec') s.crop(s_exp.get_wavelength.min(), s_exp.get_wavelength.max(), 'nm') plot_diff(s_exp, s) .. minigallery:: radis.Spectrum.crop See Also -------- :func:`radis.spectrum.operations.crop` """ from radis.spectrum.operations import crop if wunit == "default": wunit = self.get_waveunit() return crop(self, wmin=wmin, wmax=wmax, wunit=wunit, inplace=inplace) def trim(self, inplace=True): r"""Remove :py:attr:`~numpy.nan` common to all arrays on each side of the Spectrum. Returns a smaller Spectrum (inplace or not). Returns ------- s: Spectrum : trimmed Spectrum. If ``inplace=True``, Spectrum has been updated directly anyway. Allows :ref:`chaining `. """ if inplace: s = self else: s = self.copy() trim_left = [] trim_right = [] for k in s.get_vars(): m = first_nonnan_index(s._q[k]) M = last_nonnan_index(s._q[k]) if m is None: raise ValueError( f"All values are nan in {k}: Spectrum would be trimmed entirely. Are you sure?" ) trim_left.append(m) trim_right.append(M) m = np.max(trim_left) M = np.min(trim_right) if m > 0 or M < len(s) - 1: # else, no change to be made for k, v in s._q.items(): s._q[k] = v[m : M + 1] return s def sort(self, inplace=True): r"""Sort the Spectrum by wavelength / wavenumber. Parameters ---------- inplace : bool, optional if ``True``, modifies the Spectrum object directly. The default is ``True``. Returns ------- s: Spectrum sorted Spectrum. If ``inplace=True``, Spectrum has been updated directly anyway. Allows :ref:`chaining `. Examples -------- .. minigallery:: radis.spectrum.spectrum.Spectrum.sort """ if inplace: s = self else: s = self.copy() b = np.argsort(s._q["wavespace"]) for k, v in s._q.items(): s._q[k] = v[b] return s def offset(self, offset, unit, inplace=True): # type: (Spectrum, float, str) -> Spectrum """Offset the spectrum by a wavelength or wavenumber (inplace) Parameters ---------- offset: float Constant to add to all quantities in the Spectrum. unit: 'nm' or 'cm-1' unit for ``offset`` Other Parameters ---------------- inplace: bool if ``True``, modifies the Spectrum object directly. Else, returns a copy. Default ``True``. Returns ------- s: Spectrum Offset Spectrum. If ``inplace=True``, Spectrum has been updated directly anyway. Allows :ref:`chaining `. Examples -------- :: s.offset(5, 'nm') .. minigallery:: radis.Spectrum.offset See Also -------- :func:`radis.spectrum.operations.offset`, :ref:`the Spectrum page ` """ from radis.spectrum.operations import offset as offset_func return offset_func(self, offset, unit, inplace=inplace) def get_integral(self, var, wunit="default", Iunit="default", **kwargs): r"""Returns integral of variable 'var' over waverange. Parameters ---------- var: str spectral quantity to integrate wunit: str over which waverange to integrated. If ``default``, use the default Spectrum wavespace defined with :meth:`~radis.spectrum.spectrum.Spectrum.get_waveunit`. Iunit: str default ``'default'`` .. warning:: this is the unit of the quantity, not the unit of the integral. Don't forget to multiply by ``wunit`` Other Parameters ---------------- kwargs: **dict forwarded to :meth:`~radis.spectrum.spectrum.Spectrum.get` Returns ------- integral: float integral in [Iunit]*[wunit] See Also -------- :meth:`~radis.spectrum.spectrum.Spectrum.get_power`, :ref:`the Spectrum page ` """ w, I = self.get(var, wunit=wunit, Iunit=Iunit, **kwargs) return abs(nantrapz(I, w)) def get_power(self, unit="mW/cm2/sr"): r"""Returns integrated radiance (no slit) power density. Parameters ---------- Iunit: str power unit. Returns ------- P: float radiated power in ``unit`` Examples -------- :: s.get_power('W/cm2/sr') .. minigallery:: radis.spectrum.spectrum.Spectrum.offset See Also -------- :meth:`~radis.spectrum.spectrum.Spectrum.get_integral`, :ref:`the Spectrum page ` """ # New calculation as suggested in bug #460: P = self.get_integral("radiance_noslit", wunit="cm-1", Iunit="mW/cm2/sr/cm-1") # Old calculation before bug #460. # Mixes wavelengths in vacuum and air, which creates an error of 1.00027. # P = self.get_integral("radiance_noslit", wunit="nm", Iunit="mW/cm2/sr/nm") # P is in mW/cm2/sr/cm-1 * cm-1 return conv2(P, "mW/cm2/sr", unit) def has_nan(self, ignore_wavespace=True) -> bool: r""" Parameters ---------- s : Spectrum radis Spectrum. Returns ------- b : bool returns whether Spectrum has ``nan`` Note ---- ``print(s)`` will also show which spectral quantities have ````nan. """ for k, v in self._q.items(): if k == "wavespace" and ignore_wavespace: continue if anynan(v): return True return False # %% Plotting routines def get_vars(self): r"""Returns all spectral quantities stored in this object (convoluted or non convoluted) """ # remove wavespace varlist = [k for k in self._q.keys() if k != "wavespace"] return varlist def get_quantities(self): r"""Returns all spectral quantities stored in this object (convoluted or non convoluted). Wrapper to :py:meth:`~radis.spectrum.spectrum.get_vars` """ return self.get_vars() def _get_items(self) -> dict: r"""Return a dictionary of tuples, e.g:: {'radiance':(w,I), 'transmittance_noslit':(w_ns,T)} In the general case, users should use the :meth:`~radis.spectrum.spectrum.Spectrum.get` method .. warning:: real quantities are returned, not copies. """ items = { k: (self._q["wavespace"], v) for k, v in self._q.items() if k != "wavespace" } return items def _plot_matplotlib( self, var=None, wunit="default", Iunit="default", show_points=False, nfig=None, yscale="linear", show_medium="vacuum_only", normalize=False, force=False, plot_by_parts=False, show=True, show_ruler=False, **kwargs, ): r"""Plot a :py:class:`~radis.spectrum.spectrum.Spectrum` object. Plots the spectrum data using Matplotlib. .. note:: default plotting library and templates can be edited in :py:attr:`radis.config` ["plot"] Parameters ---------- var: variable (`absorbance`, `transmittance`, `transmittance_noslit`, `xsection`, etc.) For full list see :py:meth:`~radis.spectrum.spectrum.Spectrum.get_vars()`. If ``None``, plot the first thing in the Spectrum. Default ``None``. wunit: ``'default'``, ``'nm'``, ``'cm-1'``, ``'nm_vac'``, wavelength air, wavenumber, or wavelength vacuum. If ``'default'``, Spectrum :py:meth:`~radis.spectrum.spectrum.Spectrum.get_waveunit` is used. Iunit: unit for variable if `default`, default unit for quantity `var` is used. for radiance, one can use per wavelength (~ `W/m2/sr/nm`) or per wavenumber (~ `W/m2/sr/cm-1`) units Other Parameters ---------------- show_points: boolean show calculated points. Default ``True``. nfig: int, None, or 'same' plot on a particular figure. 'same' plots on current figure. For instance:: s1._plot_matplotlib() s2._plot_matplotlib(nfig='same') show_medium: bool, ``'vacuum_only'`` if ``True`` and ``wunit`` are wavelengths, plot the propagation medium in the xaxis label (``[air]`` or ``[vacuum]``). If ``'vacuum_only'``, plot only if ``wunit=='nm_vac'``. Default ``'vacuum_only'`` (prevents from inadvertently plotting spectra with different propagation medium on the same graph). yscale: 'linear', 'log' plot yscale normalize: boolean, or tuple. option to normalize quantity to 1 (ex: for radiance). Default ``False`` plot_by_parts: bool if ``True``, look for discontinuities in the wavespace and plot the different parts without connecting lines. Useful for experimental spectra produced by overlapping step-and-glue. Additional parameters read from ``kwargs`` : ``split_threshold`` and ``cutwings``. See more in :py:func:`~radis.misc.plot.split_and_plot_by_parts`. force: bool plotting on an existing figure is forbidden if labels are not the same. Use ``force=True`` to ignore that. show: bool show figure. Default ``False``. Will still show the figure in interactive mode, e.g, `%matplotlib inline` in a Notebook. show_ruler: bool if `True`, add a ruler tool to the Matplotlib toolbar. Convenient to measure distances between peaks, etc. .. warning:: still experimental ! Try it, feedback welcome ! **kwargs: **dict kwargs forwarded as argument to plot (e.g: lineshape attributes: `lw=3, color='r'`) Returns ------- line: line plot """ import matplotlib.pyplot as plt from radis.misc.plot import fix_style, set_style # Get variable x, y, wunit, Iunit = self.get( var, wunit=wunit, Iunit=Iunit, return_units="as_str" ) # Get labels xlabel = format_xlabel(wunit, show_medium) # cosmetic changes ylabel = f"{make_up(var)} ({make_up_unit(Iunit, var)})" # Plot # ------- if normalize: if isinstance(normalize, tuple): from radis.misc.arrays import norm_on wmin, wmax = normalize y = norm_on(y, x, wmin=wmin, wmax=wmax) else: # y /= y.max() y /= np.nanmax(y) Iunit = "norm" set_style() if nfig == "same": nfig = plt.gcf().number fig = plt.figure(nfig) # If figure exist, ensures xlabel and ylabel are the same (prevents some # users errors if plotting difference units!)... Note that since # 'radiance' and 'radiance_noslit' are now plotted under the same name, # they cannot be differenced. But at least this allows user to plot # both on the same figure if they want to compare [and have the same unit] def clean_error_msg(string): string = string.replace(r"$^\mathregular{", "^") string = string.replace(r"}$", "") return string if not force and (fig.gca().get_xlabel().lower() not in ["", xlabel.lower()]): raise ValueError( f"Error while plotting {var}. Cannot plot " + f"on a same figure with different xlabel: Current: {clean_error_msg(fig.gca().get_xlabel())}, New: {clean_error_msg(xlabel)}" + ". If it's a unit problem, change unit with `s.plot(..., wunit=...)`. Use force=True if you really want to plot as is." ) label1 = clean_error_msg(fig.gca().get_ylabel().lower()) label2 = clean_error_msg(ylabel.lower()) if not force and (label1 not in ["", label2]): raise ValueError( f"Error while plotting {var}. Cannot plot " + f"on a same figure with different ylabel: \nCurrent: {clean_error_msg(fig.gca().get_ylabel())}\nNew: {clean_error_msg(ylabel)}" + "\nIf it's a unit problem, you can change the unit with `s.plot(..., Iunit=)`. Use `force=True` if you really want to plot with these units." ) # Add extra plotting parameters if "lw" not in kwargs and "linewidth" not in kwargs: kwargs["lw"] = 0.5 # Add a label. Not shown by default but User can set it if using plt.legend() # (useful when plotting multiple plots on same figure) label = kwargs.pop("label", self.get_name()) ax = fig.gca() # Actual plot : # ... note: '-k' by default with style origin for first plot if not plot_by_parts: (line,) = ax.plot(x, y, label=label, **kwargs) else: (line,) = split_and_plot_by_parts(x, y, ax=fig.gca(), label=label, **kwargs) # note: split_and_plot_by_parts pops 'cutwing' & 'split_threshold' from kwargs if show_points: ax.plot(x, y, "o", color="lightgrey", **kwargs) # Labels ax.ticklabel_format(useOffset=False, axis="x") ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) ax.set_yscale(yscale) if "label" in kwargs: ax.legend() fix_style() from radis.phys.convert import cm2nm, cm2nm_air, div_safe, nm2cm, nm_air2cm if ax.child_axes != []: pass elif wunit == "cm-1": secx = ax.secondary_xaxis( "top", functions=(div_safe(cm2nm), div_safe(nm2cm)) ) secx.set_xlabel("Wavelength (nm)") elif wunit == "nm": medium = self.conditions.get("medium", "air") if medium == "vacuum": secx = ax.secondary_xaxis( "top", functions=(div_safe(nm2cm), div_safe(cm2nm)) ) else: secx = ax.secondary_xaxis( "top", functions=(div_safe(nm_air2cm), div_safe(cm2nm_air)) ) secx.set_xlabel("wavenumber (cm⁻¹)") elif wunit == "nm_air": secx = ax.secondary_xaxis( "top", functions=(div_safe(nm_air2cm), div_safe(cm2nm_air)) ) secx.set_xlabel("wavenumber (cm⁻¹)") elif wunit == "nm_vac": secx = ax.secondary_xaxis( "top", functions=(div_safe(nm2cm), div_safe(cm2nm)) ) secx.set_xlabel("wavenumber (cm⁻¹)") # Add plotting tools # ... Add cursor try: fig.cursor # if already exist, do not add again except AttributeError: from matplotlib.widgets import Cursor fig.cursor = Cursor(fig.gca(), useblit=True, color="r", lw=1, alpha=0.2) # ... Add Ruler if show_ruler: from radis.tools.plot_tools import add_ruler add_ruler(fig, wunit=wunit, Iunit=Iunit) plt.tight_layout() if show: plt.show() return line def _plot_plotly( self, var=None, wunit="default", Iunit="default", show_points=False, nfig=None, yscale="linear", show_medium="vacuum_only", normalize=False, force=False, plot_by_parts=False, show=True, show_ruler=False, template="plotly_dark", **kwargs, ): r"""Plot a :py:class:`~radis.spectrum.spectrum.Spectrum` object. Plots the spectrum data using Plotly. .. note:: default plotting library and templates can be edited in :py:attr:`radis.config` ["plot"] Parameters ---------- var: variable (`absorbance`, `transmittance`, `transmittance_noslit`, `xsection`, etc.) For full list see :py:meth:`~radis.spectrum.spectrum.Spectrum.get_vars()`. If ``None``, plot the first thing in the Spectrum. Default ``None``. wunit: ``'default'``, ``'nm'``, ``'cm-1'``, ``'nm_vac'``, wavelength air, wavenumber, or wavelength vacuum. If ``'default'``, Spectrum :py:meth:`~radis.spectrum.spectrum.Spectrum.get_waveunit` is used. Iunit: unit for variable if `default`, default unit for quantity `var` is used. for radiance, one can use per wavelength (~ `W/m2/sr/nm`) or per wavenumber (~ `W/m2/sr/cm-1`) units Other Parameters ---------------- show_points: boolean show calculated points. Default ``True``. nfig: int, None, or 'same' plot on a particular figure. 'same' plots on current figure. For instance:: s1._plot_plotly() show_medium: bool, ``'vacuum_only'`` if ``True`` and ``wunit`` are wavelengths, plot the propagation medium in the xaxis label (``[air]`` or ``[vacuum]``). If ``'vacuum_only'``, plot only if ``wunit=='nm_vac'``. Default ``'vacuum_only'`` (prevents from inadvertently plotting spectra with different propagation medium on the same graph). yscale: 'linear', 'log' plot yscale normalize: boolean, or tuple. option to normalize quantity to 1 (ex: for radiance). Default ``False`` plot_by_parts: bool if ``True``, look for discontinuities in the wavespace and plot the different parts without connecting lines. Useful for experimental spectra produced by overlapping step-and-glue. Additional parameters read from ``kwargs`` : ``split_threshold`` and ``cutwings``. See more in :py:func:`~radis.misc.plot.split_and_plot_by_parts`. force: bool plotting on an existing figure is forbidden if labels are not the same. Use ``force=True`` to ignore that. show: bool show figure. Default ``False``. Will still show the figure in interactive mode, e.g, `%matplotlib inline` in a Notebook. show_ruler: bool if `True`, add a ruler tool to the Matplotlib toolbar. Convenient to measure distances between peaks, etc. .. warning:: still experimental ! Try it, feedback welcome ! template : str, optional The Plotly template to use for the plot. Default is 'plotly_dark'. **kwargs: **dict kwargs forwarded as argument to plot (e.g: lineshape attributes: `lw=3, color='r'`) Returns ------- line: line plot """ # Get variable x, y, wunit, Iunit = self.get( var, wunit=wunit, Iunit=Iunit, return_units="as_str" ) # Get labels xlabel = format_xlabel(wunit, show_medium) ylabel = f"{make_up(var)} ({make_up_unit(Iunit, var)})" import plotly.express as px fig = px.line(x=x, y=y, template=template) fig.update_layout( xaxis_title=xlabel, yaxis_title=ylabel, yaxis=dict(tickmode="linear", dtick=0.01), ) fig.show() def plot( self, var=None, wunit="default", Iunit="default", show_points=False, nfig=None, yscale="linear", show_medium="vacuum_only", normalize=False, force=False, plot_by_parts=False, show=True, show_ruler=False, plotting_library="default", **kwargs, ): r"""Plot a :py:class:`~radis.spectrum.spectrum.Spectrum` object. .. note:: default plotting library and templates can be edited in :py:attr:`radis.config` ["plot"] Parameters ---------- var: variable (`absorbance`, `transmittance`, `transmittance_noslit`, `xsection`, etc.) For full list see :py:meth:`~radis.spectrum.spectrum.Spectrum.get_vars()`. If ``None``, plot the first thing in the Spectrum. Default ``None``. wunit: ``'default'``, ``'nm'``, ``'cm-1'``, ``'nm_vac'``, wavelength air, wavenumber, or wavelength vacuum. If ``'default'``, Spectrum :py:meth:`~radis.spectrum.spectrum.Spectrum.get_waveunit` is used. Iunit: unit for variable if `default`, default unit for quantity `var` is used. for radiance, one can use per wavelength (~ `W/m2/sr/nm`) or per wavenumber (~ `W/m2/sr/cm-1`) units Other Parameters ---------------- show_points: boolean show calculated points. Default ``True``. nfig: int, None, or 'same' plot on a particular figure. 'same' plots on current figure. For instance:: s1.plot() s2.plot(nfig='same') show_medium: bool, ``'vacuum_only'`` if ``True`` and ``wunit`` are wavelengths, plot the propagation medium in the xaxis label (``[air]`` or ``[vacuum]``). If ``'vacuum_only'``, plot only if ``wunit=='nm_vac'``. Default ``'vacuum_only'`` (prevents from inadvertently plotting spectra with different propagation medium on the same graph). yscale: 'linear', 'log' plot yscale normalize: boolean, or tuple. option to normalize quantity to 1 (ex: for radiance). Default ``False`` plot_by_parts: bool if ``True``, look for discontinuities in the wavespace and plot the different parts without connecting lines. Useful for experimental spectra produced by overlapping step-and-glue. Additional parameters read from ``kwargs`` : ``split_threshold`` and ``cutwings``. See more in :py:func:`~radis.misc.plot.split_and_plot_by_parts`. force: bool plotting on an existing figure is forbidden if labels are not the same. Use ``force=True`` to ignore that. show: bool show figure. Default ``False``. Will still show the figure in interactive mode, e.g, `%matplotlib inline` in a Notebook. show_ruler: bool if `True`, add a ruler tool to the Matplotlib toolbar. Convenient to measure distances between peaks, etc. .. warning:: still experimental ! Try it, feedback welcome ! plotting_library (str): Specifies the plotting library to use for rendering the plot. Options: - "auto" (default): Automatically selects the plotting backend based on the environment. - If running in a Jupyter notebook, it defaults to Plotly for interactive plots. - If running in a standard Python environment, it defaults to Matplotlib for static plots. - "plotly": Forces the use of Plotly for interactive plots. - "matplotlib": Forces the use of Matplotlib for static plots. **kwargs: **dict kwargs forwarded as argument to plot (e.g: lineshape attributes: `lw=3, color='r'`) Returns ------- line: line plot Examples -------- Plot an :py:func:`~radis.spectrum.models.experimental_spectrum` in arbitrary units:: s = experimental_spectrum(..., Iunit='mW/cm2/sr/nm') s.plot(Iunit='W/cm2/sr/cm-1') See more examples in :ref:`the plot Spectral quantities page `. .. minigallery:: radis.Spectrum.plot See Also -------- :py:func:`~radis.spectrum.compare.plot_diff`, :ref:`the Spectrum page ` """ # Deprecated if "plot_medium" in kwargs: show_medium = kwargs.pop("plot_medium") warn(DeprecationWarning("`plot_medium` was renamed to `show_medium`")) # Check inputs, get defaults # ------ if var in ["intensity", "intensity_noslit"]: raise ValueError("`intensity` not defined. Use `radiance` instead") if var is None: # if nothing is defined, try these first: params = self.get_vars() if "radiance" in params: var = "radiance" elif "radiance_noslit" in params: var = "radiance_noslit" elif "transmittance" in params: var = "transmittance" elif "transmittance_noslit" in params: var = "transmittance_noslit" else: # or plot the first variable we find var = list(params)[0] if var.replace("_noslit", "") in params: # favour convolved quantities var = var.replace("_noslit", "") use_plotly = False if plotting_library == "default": config = get_config() plotting_library = config["plotting_library"] if plotting_library == "auto": use_plotly = _is_running_in_notebook() elif plotting_library == "plotly": use_plotly = True elif plotting_library == "matplotlib": use_plotly = False else: raise ValueError( "plotting_library must be 'auto', 'matplotlib' or 'plotly'" ) if use_plotly: return self._plot_plotly( var=var, wunit=wunit, Iunit=Iunit, show_points=show_points, nfig=nfig, yscale=yscale, show_medium=show_medium, normalize=normalize, force=force, plot_by_parts=plot_by_parts, show=show, show_ruler=show_ruler, **kwargs, ) else: return self._plot_matplotlib( var=var, wunit=wunit, Iunit=Iunit, show_points=show_points, nfig=nfig, yscale=yscale, show_medium=show_medium, normalize=normalize, force=force, plot_by_parts=plot_by_parts, show=show, show_ruler=show_ruler, **kwargs, ) def get_populations( self, molecule=None, isotope=None, electronic_state=None, show_warning=True ): r"""Return populations that are featured in the spectrum, either as upper or lower levels. Parameters ---------- molecule: str, or None if None, only one molecule must be defined. Else, an error is raised isotope: int, or None isotope number. if None, only one isotope must be defined. Else, an error is raised electronic_state: str if None, only one electronic state must be defined. Else, an error is raised show_warning: bool if False, turns off warning about meaning of populations, see Notes and discussion on https://github.com/radis/radis/issues/508. Returns ------- pandas dataframe of levels, where levels are the index, and 'Evib' and 'nvib' are featured Notes ----- Structure:: {molecule: {isotope: {electronic_state: {'vib': pandas Dataframe, # (copy of) vib levels 'rovib': pandas Dataframe, # (copy of) rovib levels 'Ia': float # isotopic abundance }}}} (If Spectrum generated with RADIS, structure should match that of SpectrumFactory.get_populations()) Examples -------- An example on how different are populations used for partition function and spectrum calculations :: #%% CO2 example # For instance, plot populations of a given vibrational level, v1,v2,v3=(0,1,0) import radis s = radis.spectrum_test(molecule="CO2", Tvib=3000, Trot=1000, export_lines=True, export_populations="rovib", isotope=1) pops = s.get_populations("CO2")["rovib"] import matplotlib.pyplot as plt pops.query("v1==0 & v2==1 & v3==0").plot("j", "n", label="pops. used to compute partition functions") s.lines.query("v1l==0 & v2l==1 & v3l==0").plot("jl", "nl", ax=plt.gca(), kind="scatter", color="r", label="pops. of visible absorbing lines") plt.legend() plt.xlim((0,70)) .. minigallery:: radis.spectrum.spectrum.Spectrum.get_populations See Also -------- :py:meth:`~radis.spectrum.spectrum.Spectrum.plot_populations`, :py:meth:`~radis.spectrum.spectrum.Spectrum.line_survey` """ # Warn user on the meaning of these population warn( UserWarning( """Populations valid for partition function calculation but sometimes NOT for spectra calculations, e.g. CO2. See help on how to use 's.lines.query' instead, for instance in https://github.com/radis/radis/issues/508. Turn off warning with 'show_warning=False' """ ) ) # Check inputs, get default values populations = self.populations if populations is None or populations == {}: raise ValueError( "Populations not defined. Use export_populations='vib' or 'rovib' in calc_spectrum" ) if type(populations) != dict: raise TypeError("Method defined for populations as dictionary") if molecule is None: if len(list(populations.keys())) != 1: raise ValueError( f"Please choose which molecule among: {list(populations.keys())}" ) molecule = list(populations.keys())[0] if isotope is None: if len(list(populations[molecule].keys())) != 1: raise ValueError( f"Please choose which isotope among: {list(populations[molecule].keys())}" ) isotope = list(populations[molecule].keys())[0] if electronic_state is None: if len(list(populations[molecule][isotope].keys())) != 1: raise ValueError( f"Please choose which electronic state among: {list(populations[molecule][isotope].keys())}" ) electronic_state = list(populations[molecule][isotope].keys())[0] if __debug__: printdbg( f"get vib populations for {molecule}({electronic_state})[iso{isotope}]" ) # Return return populations[molecule][isotope][electronic_state] def get_vib_levels( self, molecule=None, isotope=None, electronic_state=None, first=None ): r"""Return vibrational levels in the spectrum (energies, populations) Parameters ---------- molecule: str, or None if None, only one molecule must be defined. Else, an error is raised isotope: int, or None isotope number. if None, only one isotope must be defined. Else, an error is raised electronic_state: str if None, only one electronic state must be defined. Else, an error is raised first: int, or 'all' or None only show the first N levels. If None or 'all', all levels are shown Returns ------- out: pandas DataFrame pandas dataframe of levels, where levels are the index, and 'Evib' and 'nvib' are featured """ pops = self.get_populations( molecule=molecule, isotope=isotope, electronic_state=electronic_state ) try: vib_pops = pops["vib"] except KeyError: raise KeyError( "Vibrational levels not defined in this Spectrum object. " + "If using RADIS, make sure you chose export_populations='vib'" ) if first is not None: if not "nvib" in vib_pops: raise KeyError( "Vibrational populations (nvib) not calculated in this " + "Spectrum object. Cant get first most populated levels. " + "If using RADIS, make sure you used a non_eq_spectrum " + "calculation" ) if first == "all": first = None out = vib_pops.sort_values(by="nvib", ascending=False)[:first] else: out = vib_pops # Return return out def get_rovib_levels( self, molecule=None, isotope=None, electronic_state=None, first=None ): r"""Return rovibrational levels calculated in the spectrum (energies, populations) Parameters ---------- molecule: str, or None if None, only one molecule must be defined. Else, an error is raised isotope: int, or None isotope number. if None, only one isotope must be defined. Else, an error is raised electronic_state: str if None, only one electronic state must be defined. Else, an error is raised first: int, or 'all' or None only show the first N levels. If None or 'all', all levels are shown Returns ------- out: pandas DataFrame pandas dataframe of levels, where levels are the index, and 'Evib' and 'nvib' are featured """ pops = self.get_populations( molecule=molecule, isotope=isotope, electronic_state=electronic_state ) try: rovib_pops = pops["rovib"] except KeyError: raise KeyError( "Vibrational levels not defined in this Spectrum object. " + "If using RADIS, make sure you chose export_populations='vib'" ) if first is not None: if not "n" in rovib_pops: raise KeyError( "Rovibrational populations (n) not calculated in this " + "Spectrum object. Cant get first most populated levels. " + "If using RADIS, make sure you used a non_eq_spectrum " + "calculation" ) if first == "all": first = None out = rovib_pops.sort_values(by="n", ascending=False)[:first] else: out = rovib_pops # Return return out def plot_populations( self, what=None, nunit="", correct_for_abundance=False, **kwargs ): r"""Plots vib populations if given and format is valid. Parameters ---------- what: 'vib', 'rovib', None if None plot everything nunit: '', 'cm-3' plot either in a fraction of vibrational levels, or a molecule number in in cm-3 correct_for_abundance: boolean if ``True``, multiplies each population by the isotopic abundance (as it is done during the calculation of emission integral) kwargs: **dict are forwarded to the plot Examples -------- .. minigallery:: radis.spectrum.spectrum.Spectrum.plot_populations See Also -------- :py:meth:`~radis.spectrum.spectrum.Spectrum.get_populations`, :py:meth:`~radis.spectrum.spectrum.Spectrum.line_survey`, """ import matplotlib.pyplot as plt from radis.misc.plot import fix_style, set_style # Check input, get defaults pops = self.populations if not isinstance(pops, dict): raise ValueError("Populations not defined in given Spectrum") # Plot function def _plot_elec_state(what, df, state_name, Ia, fig): if what == "vib": E, n, g = df["Evib"], df["nvib"], df["gvib"] ylabel = "n/g$_{vib}$" title = "Vibrational populations" elif what == "rovib": E, n, g = df["E"], df["n"], df["gj"] ylabel = "n/(2J+1)" title = "Rovibrational populations" if correct_for_abundance: n = n * Ia if nunit == "cm-3": from radis.phys.constants import k_b try: P = self.conditions["pressure"] # bar T = self.conditions["Tgas"] mfrac = self.conditions["mole_fraction"] except KeyError: raise KeyError( "P (pressure), T (Tgas) and n (mole_fraction) " + "are needed to calculate total number density in (cm-3)" ) N = P * 1e5 / k_b / T * mfrac * 1e-6 n = n * N unitlabel = " [cm-3]" elif nunit == "": unitlabel = " [fraction]" else: raise ValueError(f"Unknown unit: {nunit}") # Plot if fig is None: fig = plt.figure() plt.xlabel("Energy (cm-1)") plt.ylabel(ylabel + unitlabel) plt.yscale("log") plt.title(title) ax = fig.gca() ax.plot(E, n / g, "o", label=state_name, **kwargs) return fig # Initialize figures, styles fig_vib = None fig_rovib = None set_style() # Loop over all molecules, all isotopes, all electronic states # Note that the below works for both dict and pandas dataframe for molecule, isotopes in pops.items(): for isotope, elec_states in isotopes.items(): for elec_state, content in elec_states.items(): state_name = f"{molecule}({elec_state})(iso{isotope})" Ia = None if correct_for_abundance: if "Ia" in list(content.keys()): Ia = content["Ia"] else: raise KeyError( "Ia: isotopic abundance not defined in " + "Spectrum populations" ) for k in content.keys(): if k == "rovib" and (what == k or what is None): df = self.get_rovib_levels( molecule, isotope, elec_state, first="all" ) # sort by n + check is defined) fig_rovib = _plot_elec_state( "rovib", df, state_name, Ia, fig_rovib ) if k == "vib" and (what == k or what is None): df = self.get_vib_levels( molecule, isotope, elec_state, first="all" ) # sort by n + check is defined) fig_vib = _plot_elec_state( "vib", df, state_name, Ia, fig_vib ) # Update for fig in [fig_vib, fig_rovib]: if fig is not None: ax = fig.gca() ax.legend() fix_style(ax=ax) # %% ------------------ Instrumental Slit Function --------------------- def _apply_slit_nm_on_cm1( self, varlist, slit_function, norm_by, shape, center_wavespace, mode, plot_slit, store, auto_recenter_crop, assert_evenly_spaced, verbose, *args, **kwargs, ): from scipy.interpolate import interp1d from radis.tools.slit import convolve_with_slit, get_slit_function w_cm = self._q["wavespace"] w_nm = cm2nm(w_cm[::-1]) # increasing order w_uniform = np.linspace(w_nm.min(), w_nm.max(), len(w_nm)) I_conv_slices = {} for qns in varlist: q = qns[:-7] I = self._q[qns][::-1] interp_to_nm = interp1d(w_nm, I, kind="linear", fill_value="extrapolate") I_interp = interp_to_nm(w_uniform) wslit_nm, Islit_nm = get_slit_function( slit_function, unit="nm", norm_by=norm_by, shape=shape, center_wavespace=center_wavespace, return_unit="nm", wstep=np.diff(w_uniform).min(), auto_recenter_crop=auto_recenter_crop, verbose=verbose, plot=plot_slit, *args, **kwargs, ) _, I_conv_nm = convolve_with_slit( w_uniform, I_interp, wslit_nm, Islit_nm, mode=mode, wunit="nm", verbose=verbose, assert_evenly_spaced=assert_evenly_spaced, ) interp_back = interp1d( w_uniform, I_conv_nm, kind="linear", fill_value="extrapolate" ) I_final = interp_back(w_nm)[::-1] I_conv_slices[q] = [I_final] new_units = {} for q in I_conv_slices.keys(): qns = q + "_noslit" if norm_by == "area": new_units[q] = self.units[qns] elif norm_by == "max": new_units[q] = (Unit("nm") * Unit(self.units[qns])).to_string() else: raise ValueError(f"Unknown normalization type: {norm_by}") for q in I_conv_slices.keys(): I_conv = np.hstack(I_conv_slices[q]) self._q[q] = I_conv self.units[q] = new_units[q] if store: self._slit["wavespace"] = wslit_nm self._slit["intensity"] = Islit_nm self.conditions.update( { "slit_function": slit_function, "slit_unit": "nm", "slit_dispersion": None, "slit_dispersion_threshold": None, "slit_shape": shape, "norm_by": norm_by, } ) return self def apply_slit( self, slit_function, unit="nm", shape="triangular", center_wavespace=None, norm_by="area", mode="valid", plot_slit=False, store=True, slit_dispersion=None, slit_dispersion_threshold=0.01, auto_recenter_crop=True, assert_evenly_spaced=True, verbose=True, inplace=True, *args, **kwargs, ): r"""Apply an instrumental slit function to all quantities in Spectrum. Slit function can be generated with usual shapes (see ``shape=``) or imported from an experimental slit function (path to a text file or numpy array of shape n*2). Convoluted spectra are cut on the edge compared to non-convoluted spectra, to remove side effects. See ``mode=`` to change this behaviour. Warning with units: read about ``'unit'`` and ``'return_unit'`` parameters. Parameters ---------- slit_function: float or str or array If ``float``: generate slit function with FWHM of slit function (in nm or cm-1 depending on ``unit=``). A (top, base) tuple of (``float``, ``float``) is required when asking for a trapezoidal slit function. If ``.txt``: import experimental slit function from .txt file: format must be 2-columns with wavelengths and intensity (doesn't have to be normalized) If ``array``: format must be 2-columns with wavelengths and intensity (doesn't have to be normalized) It is recommended to truncate the input slit function to its minimum useful spectral extension (see Notes of :func:`~radis.tools.slit.convolve_with_slit`). unit: ``'nm'`` or ``'cm-1'`` unit of slit_function (FWHM, or imported file) shape: ``'triangular'``, ``'trapezoidal'``, ``'gaussian'``, or any of :data:`~radis.tools.slit.SLIT_SHAPES` which shape to use when generating a slit. Will call, respectively, :func:`~radis.tools.slit.triangular_slit`, :func:`~radis.tools.slit.trapezoidal_slit`, :func:`~radis.tools.slit.gaussian_slit`. Default 'triangular' center_wavespace: float, or ``None`` center of slit when generated (in unit). Not used if slit is imported. norm_by: ``'area'``, ``'max'`` normalisation type: - ``'area'`` normalizes the slit function to an area of 1. It conserves energy, and keeps the same units. - ``'max'`` normalizes the slit function to a maximum of 1. The convoluted spectrum units change (they are multiplied by the spectrum waveunit, e.g: a radiance non convoluted in mW/cm2/sr/nm on a wavelength (nm). range will yield a convoluted radiance in mW/cm2/sr. Note that the slit is set to 1 in the Spectrum wavespace (i.e: a Spectrum calculated in cm-1 will have a slit set to 1 in cm-1). Default ``'area'`` mode: ``'valid'``, ``'same'`` ``'same'`` returns output of same length as initial spectra, but boundary effects are still visible. ``'valid'`` returns output of length len(spectra) - len(slit) + 1, for which lines outside of the calculated range have no impact. Default ``'valid'``. Other Parameters ---------------- assert_evenly_spaced: boolean, or ``'resample'`` for the convolution to be accurate, ``w`` should be evenly spaced. If ``assert_evenly_spaced=True``, then we check this is the case, and raise an error if arrays is not evenly spaced. If ``'resample'``, then we resample ``w`` and ``I`` if needed. Recommended, but it takes some time. auto_recenter_crop: bool if ``True``, recenter slit and crop zeros on the side when importing an experimental slit. Default ``True``. See :func:`~radis.tools.slit.recenter_slit`, :func:`~radis.tools.slit.crop_slit` plot_slit: boolean if ``True``, plot slit store: boolean if ``True``, store slit in the Spectrum object so it can be retrieved with :meth:`~radis.spectrum.spectrum.Spectrum.get_slit` and plot with :meth:`~radis.spectrum.spectrum.Spectrum.plot_slit`. Default ``True`` slit_dispersion: func of (lambda, in ``'nm'``), or ``None`` spectrometer reciprocal function : dλ/dx(λ) (in ``nm``) If not ``None``, then the slit_dispersion function is used to correct the slit function for the whole range. Can be important if slit function was measured far from the measured spectrum (e.g: a slit function measured at 632.8 nm will look broader at 350 nm because the spectrometer dispersion is higher at 350 nm. Therefore it should be corrected) Default ``None`` .. warning:: slit dispersion function is assumed to be given in ``nm`` if your spectrum is stored in ``cm-1`` the wavenumbers are converted to wavelengths before being forwarded to the dispersion function See :func:`~radis.test.tools.test_slit.test_auto_correct_dispersion` for an example of the slit dispersion effect. A Python implementation of the slit dispersion: >>> def f(lbd): >>> return w/(2*f)*(tan(Φ)+sqrt((2*d/m/(w*1e-9)*cos(Φ))^2-1)) Theoretical / References: >>> dλ/dx ~ d/mf # at first order >>> dλ/dx = w/(2*f)*(tan(Φ)+sqrt((2*d/m/(w)*cos(Φ))^2-1)) # cf with: - Φ: spectrometer angle (°) - f: focal length (mm) - m: order of dispersion - d: grooves spacing (mm) = 1/gr with gr in (gr/mm) See Laux 1999 "Experimental study and modeling of infrared air plasma radiation" for more information slit_dispersion_warning_threshold: float if not ``None``, check that slit dispersion is about constant (< ``threshold`` change) on the calculated range. Default 0.01 (1%). See :func:`~radis.tools.slit.offset_dilate_slit_function` inplace: bool if ``True``, adds convolved arrays directly in the Spectrum. If ``False``, returns a new Spectrum with only the convolved arrays. Note: if you want a new Spectrum with both the convolved and non convolved quantities, use :: s.copy().apply_slit() *args, **kwargs are forwarded to slit generation or import function verbose: bool print stuff energy_threshold: float tolerance fraction when resampling. Default ``1e-3`` (0.1%) If areas before and after resampling differ by more than that an error is raised. Returns ------- Spectrum : same Spectrum, with new spectral arrays. Allows :ref:`chaining `. If ``inplace=False``, return a new Spectrum with the new spectral arrays only. Notes ----- Units: the slit function is first converted to the wavespace (wavelength/wavenumber) that the Spectrum is stored in, and applied to the spectral quantities in their native wavespace. Implementation: :func:`~radis.tools.slit.convolve_with_slit` is applied to all quantities in :meth:`~radis.spectrum.spectrum.Spectrum.get_vars` that ends with _noslit. Generate a triangular instrumental slit function (or any other shape depending of shape=) with base ``slit_function_base`` (Uses the central wavelength of the spectrum for the slit function generation) We deal with several special cases (which makes the code a little heavy, but the method very versatile): - when slit unit and spectrum unit arent the same - when spectrum is not evenly spaced Examples -------- :: s.apply_slit(1.2, 'nm') This applies the instrumental function to all available spectral arrays. To manually apply the slit to a particular spectral array, use :py:meth:`~radis.spectrum.spectrum.Spectrum.take` :: s.take('transmittance_noslit').apply_slit(1.2, 'nm') See :func:`~radis.tools.slit.convolve_with_slit` for more details on Units and Normalization The slit is made considering the "center wavelength" which is the mean wavelength of the full spectrum you are applying it to. .. minigallery:: radis.Spectrum.apply_slit :add-heading: See Also -------- :func:`~radis.tools.slit.get_slit_function`, :func:`~radis.tools.slit.convolve_with_slit`, :ref:`the Spectrum page ` """ # TODO: add warning if FWHM >= wstep(spectrum)/5 from radis.tools.slit import ( cast_waveunit, convolve_with_slit, get_slit_function, normalize_slit, offset_dilate_slit_function, ) # Check inputs # --------- if "slit_unit" in kwargs: unit = kwargs.pop("slit_unit") warn(DeprecationWarning("slit_unit was renamed unit")) unit = cast_waveunit(unit) varlist = [k for k in self.get_vars() if k.endswith("_noslit")] if len(varlist) == 0: raise AssertionError( "No spectral arrays to apply slit on. Spectral arrays " + "to be convolved should end with _noslit. " + f"Available arrays in Spectrum : {self.get_vars()}" ) # Forward relevant inputs to convolution instead of slit function generation kwargsconvolve = {} for kw in ["slit_dispersion", "verbose"]: if kw in kwargs: kwargsconvolve.update({kw: kwargs.pop(kw)}) # For non evenly distributed cases we take the minimum wstep among the # spectral range (and resample in convolve_with_slit) # Note: by construction all variables now have the same wavespace w = self._q["wavespace"] # non convoluted wavespace wstep = abs(diff(w)).min() assert wstep > 0 waveunit = self.get_waveunit() if __debug__: printdbg( f"apply_slit: {slit_function} in {unit}, center `{center_wavespace}`{unit}, applied in waveunit {waveunit}" ) if center_wavespace is None: # center_wavespace should be ~ unit center_wavespace = w[len(w) // 2] # w ~ waveunit if waveunit == "cm-1" and unit == "nm": center_wavespace = cm2nm(center_wavespace) # wavenum > wavelen elif waveunit == "nm" and unit == "cm-1": center_wavespace = nm2cm(center_wavespace) # wavelen > wavenum # Get slit once and for all (and convert the slit unit # to the Spectrum `waveunit` if wavespaces are different) # ------- wslit0, Islit0 = get_slit_function( slit_function, unit=unit, norm_by=norm_by, shape=shape, center_wavespace=center_wavespace, return_unit=waveunit, wstep=wstep, auto_recenter_crop=auto_recenter_crop, verbose=verbose, plot=plot_slit, *args, **kwargs, ) # Check if dispersion is specified if slit_dispersion is not None: if waveunit == "nm": w_nm = w wslit0_nm = wslit0 else: w_nm = cm2nm(w) wslit0_nm = cm2nm(wslit0) slice_windows = _cut_slices( w_nm, wslit0_nm, slit_dispersion, slit_dispersion_threshold ) else: slice_windows = [np.ones_like(w, dtype=bool)] # Create dictionary to store convolved I_conv_slices = {} for qns in varlist: # Convolve and store the output in a new variable name (quantity name minus `_noslit`) # Create if required q = qns[:-7] # new name (minus '_noslit') w_conv_slices = [] I_conv_slices[q] = [] if unit == "nm" and waveunit == "cm-1": # special handling for wavelength-space convolution # for constant nm slit return self._apply_slit_nm_on_cm1( varlist, slit_function, norm_by, shape, center_wavespace, mode, plot_slit, store, auto_recenter_crop, assert_evenly_spaced, verbose, *args, **kwargs, ) # Loop over all waverange slices (needed if slit changes over the spectral range) for islice, slice_window in enumerate(slice_windows): print("Slice {0}/{1}".format(islice + 1, len(slice_windows))) # Scale slit if slit_dispersion is not None: # apply spectrometer linear dispersion function. # dont forget it has to be added in nm and not cm-1 wslit, Islit = offset_dilate_slit_function( wslit0_nm, Islit0, w_nm[slice_window], slit_dispersion, threshold=slit_dispersion_threshold, verbose=verbose, ) # Convert it back if needed if waveunit == "cm-1": wslit = nm2cm(wslit) # We need to renormalize now that Islit has changed wslit, Islit = normalize_slit(wslit, Islit, norm_by=norm_by) else: wslit = wslit0 Islit = Islit0 # no need to renormalize it # Apply to all variables # --------- for i, q in enumerate(I_conv_slices.keys()): # Convolve and store the output in a new variable name (quantity name minus `_noslit`) # Create if required qns = q + "_noslit" w_window = w[slice_window] I_window = self._q[qns][slice_window] # Apply convolution w_conv_window, I_conv_window = convolve_with_slit( w_window, I_window, wslit, Islit, mode=mode, wunit=waveunit, verbose=verbose, assert_evenly_spaced=assert_evenly_spaced, # assumes Spectrum is correct by construction **kwargsconvolve, ) if i == 0: w_conv_slices.append(w_conv_window) I_conv_slices[q].append(I_conv_window) # Get units new_units = {} for q in I_conv_slices.keys(): qns = q + "_noslit" # Get units if norm_by == "area": new_units[q] = self.units[qns] elif norm_by == "max": new_unit = (Unit(unit) * Unit(self.units[qns])).to_string() # because it's like if we multiplied by slit FWHM in the wavespace # it was generated new_units[q] = new_unit # Note: there was another mode called 'max2' where, unlike 'max', # unit was multiplied by [unit] not [return_unit] # Removed for simplification. You should stay with norm_by='area' anyway else: raise ValueError(f"Unknown normalization type: {norm_by}") # Merge and store all variables # --------- w_conv = np.hstack(w_conv_slices) if inplace: if len(self._q["wavespace"]) != len(w_conv) or not np.allclose( self._q["wavespace"], w_conv ): raise AssertionError( "Wavespace of convolved arrays are different and they cannot be " "stored in the same Spectrum object. You can use " "Spectrum.apply_slit(inplace=False) to return a new spectrum " "with only the convolved arrays." ) for q in I_conv_slices.keys(): # Merge all slices I_conv = np.hstack(I_conv_slices[q]) # Store self._q[q] = I_conv # Get units self.units[q] = new_units[q] s_out = self else: # make new spectrum quantities = { "wavespace": w_conv } # a copy will be created in Spectrum creation for q in I_conv_slices.keys(): # Merge all slices quantities[q] = np.hstack(I_conv_slices[q]) s_out = Spectrum( quantities, units=new_units, conditions=self.conditions.copy(), cond_units=self.cond_units.copy(), wunit=waveunit, name=self.name, check_wavespace=False, ) # Store slit in Spectrum, in the Spectrum unit if store: s_out._slit["wavespace"] = wslit0 # in 'waveunit' s_out._slit["intensity"] = Islit0 # Update conditions s_out.conditions["slit_function"] = slit_function s_out.conditions["slit_unit"] = unit # input slit unit s_out.conditions["slit_dispersion"] = slit_dispersion s_out.conditions["slit_dispersion_threshold"] = slit_dispersion_threshold s_out.conditions["slit_shape"] = shape # TODO: probably removed after Spectrum is stored. s_out.conditions["norm_by"] = norm_by return s_out # to be able to chain: s.apply_slit().plot() def get_slit(self, wunit="same"): r"""Get slit function that was applied to the Spectrum. Returns ------- wslit, Islit: array slit function with wslit in Spectrum ``waveunit``. See :meth:`~radis.spectrum.spectrum.Spectrum.get_waveunit` """ if not wunit in ["same", self.get_waveunit()]: raise NotImplementedError( f"Unit must be Spectrum waveunit: {self.get_waveunit()}" ) # Make sure that slit is stored already try: wslit = self._slit["wavespace"] # stored in Spectrum waveunit Islit = self._slit["intensity"] except KeyError: raise KeyError( "Slit function not found in Spectrum " + "conditions. Have you used Spectrum.apply_slit " + "with store=True?" ) return wslit, Islit def plot_slit(self, wunit=None, waveunit=None): r"""Plot slit function that was applied to the Spectrum. If dispersion was used (see :meth:`~radis.spectrum.spectrum.Spectrum.apply_slit`) the different slits are built again and plotted too (dotted). Parameters ---------- wunit: ``'nm'``, ``'cm-1'``, or ``None`` plot slit in wavelength or wavenumber. If ``None``, use the unit the slit in which the slit function was given. Default ``None`` Returns ------- fix, ax: matplotlib objects figure and ax Examples -------- :: s = radis.spectrum_test() s.apply_slit((0.4, 0.6), "nm") # add trapezoidal slit function s.plot_slit() # plot the slit function that was applied .. minigallery:: radis.spectrum.spectrum.Spectrum.plot_slit See Also -------- :ref:`the Spectrum page ` """ # Deprecated inputs if waveunit is not None: warn( "`waveunit=` parameter in from_array is now named `wunit=`", DeprecationWarning, ) wunit = waveunit from radis.tools.slit import ( normalize_slit, offset_dilate_slit_function, plot_slit, ) # Check inputs assert wunit in ["nm", "cm-1", "nm_vac", None] # @dev: note: wunit in 'nm_vac' also implemented for consistency, # although not mentioned in docs. if wunit is None: wunit = self.conditions["slit_unit"] # Get slit arrays (in Spectrum.waveunit) wslit0, Islit0 = self.get_slit() # as imported # Get slit unit norm_by = self.conditions["norm_by"] waveunit = self.get_waveunit() if norm_by == "area": Iunit = f"1/{waveunit}" elif norm_by == "max": # set maximum to 1 Iunit = "" elif norm_by is None: Iunit = None else: raise ValueError(f"Unknown normalization type: `norm_by` = {norm_by}") # Plot in correct unit (plot_slit deals with the conversion if needed) fig, ax = plot_slit( wslit0, Islit0, wunit=waveunit, plot_unit=wunit, Iunit=Iunit ) # Plot other slit functions if dispersion was applied: if ( "slit_dispersion" in self.conditions and "slit_dispersion_threshold" in self.conditions ): slit_dispersion = self.conditions["slit_dispersion"] slit_dispersion_threshold = self.conditions["slit_dispersion_threshold"] if slit_dispersion is not None: waveunit = self.get_waveunit() # Get slit in air wavelength: if waveunit == "nm": wslit0_nm = wslit0 elif waveunit == "nm_vac": wslit0_nm = vacuum2air(wslit0) else: wslit0_nm = cm2nm(wslit0) w_nm = self.get_wavelength(medium="air") slice_windows = _cut_slices( w_nm, wslit0, slit_dispersion, slit_dispersion_threshold ) # Loop over all waverange slices (needed if slit changes over the spectral range) for islice, slice_window in enumerate(slice_windows): w_nm_sliced = w_nm[slice_window] w_min = w_nm_sliced.min() w_max = w_nm_sliced.max() # apply spectrometer linear dispersion function. # dont forget it has to be added in nm and not cm-1 wslit, Islit = offset_dilate_slit_function( wslit0_nm, Islit0, w_nm[slice_window], slit_dispersion, threshold=slit_dispersion_threshold, verbose=False, ) # Convert it back if needed if waveunit == "cm-1": wslit = nm2cm(wslit) elif waveunit == "nm_vac": wslit = air2vacuum(wslit) # We need to renormalize now that Islit has changed wslit, Islit = normalize_slit(wslit, Islit, norm_by=norm_by) plot_slit( wslit, Islit, wunit=waveunit, plot_unit=wunit, Iunit=Iunit, ls="--", title=f"Slit used on range {w_min:.2f}-{w_max:.2f} nm", ) return fig, ax def line_survey( self, overlay=None, wunit="cm-1", Iunit="hitran", medium="air", cutoff=None, plot="S", lineinfo=["int", "A", "El"], barwidth="hwhm_voigt", # 0.01, yscale="log", writefile=None, *args, **kwargs, ): r"""Plot Line Survey (all linestrengths used for calculation) Output in Plotly (html) Parameters ---------- overlay: tuple (w, I, [name], [units]), or list or tuples plot (w, I) on a secondary axis. Useful to compare linestrength with calculated / measured data:: LineSurvey(overlay='abscoeff') wunit: ``'nm'``, ``'cm-1'`` wavelength / wavenumber units Iunit: ``'hitran'``, ``'splot'`` Linestrength output units: - ``'hitran'``: (cm-1/(molecule/cm-2)) - ``'splot'`` : (cm-1/atm) (Spectraplot units [2]_) Note: if not None, cutoff criteria is applied in this unit. Not used if plot is not 'S' medium: ``'air'``, ``'vacuum'`` show wavelength either in air or vacuum. Default ``'air'`` plot: str what to plot. Default ``'S'`` (scaled line intensity). But it can be any key in the lines, such as population (``'nu'``), or Einstein coefficient (``'A'``) lineinfo: list, or ``'all'`` extra line information to plot. Should be a column name in the databank (s.lines). For instance: ``'int'``, ``'selbrd'``, etc... Default [``'int'``] Other Parameters ---------------- writefile: str if not ``None``, a valid filename to save the plot under .html format. If ``None``, use the ``fig`` object returned to show the plot. yscale: ``'log'``, ``'linear'`` Default ``'log'`` barwidth: float or str if float, width of bars, in ``wunit``, as a fraction of full-range; i.e. :: barwidth=0.01 makes bars span 1% of the full range. if ``str``, uses the column as width. Example :: barwidth = 'hwhm_voigt' Returns ------- fig: a Plotly figure. If using a Jupyter Notebook, the plot will appear. Else, use ``writefile`` to export to an html file. Plot in Plotly. See Output in [1]_ Examples -------- An example using the :class:`~radis.lbl.factory.SpectrumFactory` to generate a spectrum:: from radis import SpectrumFactory sf = SpectrumFactory( wavenum_min=2380, wavenum_max=2400, mole_fraction=400e-6, path_length=100, # cm isotope=[1], export_lines=True, # required for LineSurvey! db_use_cached=True) sf.load_databank('HITRAN-CO2-TEST') s = sf.eq_spectrum(Tgas=1500) s.apply_slit(0.5) s.line_survey(overlay='radiance_noslit', barwidth=0.01, lineinfo="all") # or barwidth='hwhm_voigt' See the output in :ref:`Examples ` .. minigallery:: radis.spectrum.spectrum.Spectrum.line_survey References ---------- .. [1] `RADIS Online Documentation (LineSurvey) `__ .. [2] `SpectraPlot `__ See Also -------- :func:`~radis.tools.line_survey.LineSurvey`, :ref:`the Spectrum page ` """ from radis.tools.line_survey import LineSurvey # Check inputs if wunit == "default": wunit = self.get_waveunit() def get_overlay(overlay): """Overlay line survey with a spectral quantity (like radiance or transmittance)""" if isinstance(overlay, str): # either get it from the Spectrum w, I, _, units = self.get(overlay, wunit=wunit, return_units="as_str") name = overlay return (w, I, name, units) else: # Or use a given tuple or arrays try: w, I = overlay except (TypeError, ValueError): raise ValueError( "Overlay has to be string, or (w,I) tuple of " + "arrays" ) return (w, I, "", "") if overlay is not None: if type(overlay) is not list: overlay = [overlay] overlay = [get_overlay(ov) for ov in overlay] return LineSurvey( self, overlay=overlay, wunit=wunit, Iunit=Iunit, medium=medium, cutoff=cutoff, plot=plot, lineinfo=lineinfo, barwidth=barwidth, yscale=yscale, writefile=writefile, *args, **kwargs, ) def get_conditions(self): r"""Get all physical / computational parameters. .. minigallery:: radis.spectrum.spectrum.Spectrum.get_conditions :add-heading: See Also -------- :py:meth:`~radis.spectrum.spectrum.Spectrum.get_conditions`, :py:meth:`~radis.spectrum.spectrum.Spectrum.print_conditions`, :ref:`the Spectrum page ` """ return self.conditions def get_condition(self, condition, unit="default", return_unit=False): """Get condition in arbitrary unit If condition was stored as a dimensioned value, return a dimensioend value (converted to the right unit) If condition was stored as a non dimensioned value, return a non dimensioned value (still converted to the right unit) Parameters ---------- condition: str condition name unit: str unit name return_unit: bool if True, return value and unit separately (value is not dimensioned in this case) Examples -------- Get pressure (as float, non dimensioned) in Pascal:: s = radis.spectrum_test() pressure_Pa, _ = s.get_condition("pressure", "Pa", return_unit=True) Get pressure and convert it to a dimensioned value:: from radis.phys.units import Unit as u s = radis.spectrum_test() pressure_value, pressure_unit = s.get_condition("pressure", return_unit=True) pressure = pressure_value * u(pressure_unit) """ cond = self.conditions[condition] if unit == "default": if return_unit: if hasattr(cond, "unit"): return cond.value, cond.unit else: if condition not in self.cond_units: raise ValueError( f"Unit of condition {condition} not defined in s.cond_units" ) return cond, self.cond_units[condition] else: return cond else: import astropy.units as u if return_unit: if hasattr(cond, "unit"): return cond.to(unit).value, unit else: if condition not in self.cond_units: raise ValueError( f"Unit of condition {condition} not defined in s.cond_units" ) return (cond * u.Unit(self.cond_units[condition])).to( unit ).value, unit else: # If was dimensioned, return a dimensioend value (converted to the right unit) # If not dimensioned, return a non dimensioned value (still converted to the right unit) if hasattr(cond, "unit"): return cond.to(unit) else: return (cond * u.Unit(self.cond_units[condition])).to(unit).value return def print_conditions(self, **kwargs): r"""Prints all physical / computational parameters. Parameters ---------- kwargs: dict refer to :py:func:`~radis.spectrum.utils.print_conditions` Examples -------- :: s = radis.spectrum_test() s.print_conditions() You can also simply print the Spectrum object directly:: s = radis.spectrum_test() print(s) Both syntaxes above will return (in radis==0.15):: # output >> Spectrum Name: CO-hitran-700K-#3680 Spectral Quantities ---------------------------------------- abscoeff [cm-1] (40,002 points) absorbance (40,002 points) emissivity_noslit (40,002 points) transmittance_noslit (40,002 points) radiance_noslit [mW/cm2/sr/cm-1] (40,002 points) Physical Conditions ---------------------------------------- Tgas 700 K Trot 700 K Tvib 700 K isotope 1,2,3 mole_fraction 0.1 molecule CO overpopulation None path_length 1 cm pressure 1.01325 bar rot_distribution boltzmann self_absorption True state X thermal_equilibrium True vib_distribution boltzmann wavenum_max 2300.0000 cm-1 wavenum_min 1900.0000 cm-1 Computation Parameters ---------------------------------------- NwG 3 NwL 5 Tref 296 K broadening_method voigt_poly cutoff 1e-27 cm-1/(#.cm-2) dbformat hitran dbpath C:\Users\erwan\.radisdb\hitran\CO.hdf5 default_output_unit cm-1 diluents {'air': 0.9} folding_thresh 1e-06 include_neighbouring_lines True memory_mapping_engine auto neighbour_lines 0 cm-1 optimization simple parsum_mode full summation profiler {'spectrum_calculation': {'check_line_databank': ... pseudo_continuum_threshold 0 radis_version 0.14 sparse_ldm True spectral_points 40000.0 truncation 50 cm-1 waveunit cm-1 wstep 0.01 cm-1 zero_padding 40002 Config parameters ---------------------------------------- DEFAULT_DOWNLOAD_PATH ~/.radisdb GRIDPOINTS_PER_LINEWIDTH_ERROR_THRESHOLD 1 GRIDPOINTS_PER_LINEWIDTH_WARN_THRESHOLD 3 SPARSE_WAVERANGE auto Information ---------------------------------------- calculation_time 0.11549490000000162 s chunksize None db_use_cached True dxG 0.1375350788016573 dxL 0.20180288881201608 export_lines False export_populations None export_rovib_fraction True levelsfmt None lines_calculated 742 lines_cutoff 0 lines_in_continuum 0 load_energies False lvl_use_cached True parfuncpath None total_lines 742 warning_broadening_threshold 0.01 warning_linestrength_cutoff 0.01 wavenum_max_calc 2300.0000 cm-1 wavenum_min_calc 1900.0000 cm-1 ---------------------------------------- .. minigallery:: radis.spectrum.spectrum.Spectrum.print_conditions :add-heading: See Also -------- :py:meth:`~radis.spectrum.spectrum.Spectrum.get_conditions`, :py:func:`~radis.spectrum.utils.print_conditions`, :ref:`the Spectrum page ` """ return print_conditions(self.get_conditions(), self.cond_units, **kwargs) def store( self, path, discard=["lines", "populations"], compress=True, add_info=None, add_date=None, if_exists_then="error", verbose=True, ): r"""Save a Spectrum object in JSON format. Object can be recovered with :func:`~radis.tools.database.load_spec`. If many Spectrum are saved in a same folder you can view their properties with the :class:`~radis.tools.database.SpecDatabase` structure. Parameters ---------- path: str or file-like object If a string: path to folder (database) or file. If a folder, file is saved to database and name is generated automatically. If a file name, then Spectrum is saved to this file and the later formatting options dont apply. If a file-like object (e.g., ``io.BytesIO``, ``io.StringIO``): writes directly to the object without any file system operations. Use ``BytesIO`` when ``compress=True`` (binary output), use ``StringIO`` when ``compress=False`` (text output). When using file-like objects, the ``add_date``, ``add_info``, and ``if_exists_then`` parameters are ignored. compress: boolean if ``False``, save under text format, readable with any editor. if ``True``, saves under binary format. Faster and takes less space. If ``2``, removes all quantities that can be regenerated with s.update(), e.g, transmittance if abscoeff and path length are given, radiance if emisscoeff and abscoeff are given in non-optically thin case, etc. Default ``True``. add_info: list append these parameters and their values if they are in conditions example:: add_info = ['Tvib', 'Trot'] discard: list of str parameters to exclude, for instance to save some memory for instance Default [`lines`, `populations`]: retrieved Spectrum looses the :meth:`~radis.spectrum.spectrum.Spectrum.line_survey` ability, and :meth:`~radis.spectrum.spectrum.Spectrum.plot_populations` (but it saves tons of memory!) if_exists_then: ``'increment'``, ``'replace'``, ``'error'``, ``'ignore'`` what to do if file already exists. If increment an incremental digit is added. If replace file is replaced (yeah). If ``'ignore'`` no file is created. If error (or anything else) an error is raised. Default `error` Returns ------- str or file-like object If path was a string: filename used (may be different from given path as new info or incremental identifiers are added). If path was a file-like object: returns the same object. Notes ----- If many spectra are stored in a folder, it may be time to set up a :class:`~radis.tools.database.SpecDatabase` structure to easily see all Spectrum conditions and get Spectrum that suits specific parameters. Examples -------- Store a spectrum in compressed mode, regenerate quantities after loading:: from radis import load_spec s.store('test.spec', compress=True) # s is a Spectrum s2 = load_spec('test.spec') s2.update() # regenerate missing quantities Store a spectrum to a BytesIO buffer (no disk write):: from io import BytesIO buffer = BytesIO() s.store(buffer, compress=True) buffer.seek(0) # Reset position for reading # buffer.getvalue() contains the serialized spectrum .. minigallery:: radis.spectrum.spectrum.Spectrum.store :add-heading: See Also -------- :class:`~radis.tools.database.SpecDatabase`, :func:`~radis.tools.database.load_spec`, :meth:`~radis.spectrum.spectrum.Spectrum.save`, :meth:`~radis.spectrum.spectrum.Spectrum.savetxt` """ # TODO # # - in case a spectrometer linear dispersion function is used in # :meth:`~radis.spectrum.spectrum.Spectrum.apply_slit`, it probably isn't # stored with the current code. Find a workaround? from radis.tools.database import save if isinstance(discard, str): discard = [discard] return save( self, path, discard=discard, compress=compress, add_info=add_info, add_date=add_date, if_exists_then=if_exists_then, verbose=verbose, ) def save(self, *args, **kwargs): r"""Alias to Spectrum.store. See Spectrum.store for documentation """ return self.store(*args, **kwargs) def to_json(self, *args, **kwargs): r"""Alias to Spectrum.store(compress=False). See Spectrum. :py:meth:`~radis.spectrum.spectrum.Spectrum.store` for documentation """ return self.store(compress=False, *args, **kwargs) def to_hdf5(self, file, engine="pytables"): r"""Convert a ``radis`` :py:class:`~radis.spectrum.spectrum.Spectrum` object to HDF5 format. Uses :py:func:`~radis.io.spec_hdf.spec2hdf` Parameters ---------- file: str path to save the HDF5 file engine : str, optional Which HDF5 library to use for reading/writing data. Options are: - 'pytables' (default) : Use PyTables library Returns ------- The spectrum is saved to disk in HDF5 format .. note:: The HDF5 file will contain all spectral arrays in an 'arrays' group with their units as metadata, the lines database in a 'lines' group if present, populations in metadata if present, and conditions and references in metadata. Use :py:meth:`~radis.spectrum.spectrum.Spectrum.from_hdf5` to load the spectrum back. Examples -------- :: s.to_hdf5('spectrum.h5') Load back with:: from radis import Spectrum s2 = Spectrum.from_hdf5('spectrum.h5') .. minigallery:: radis.Spectrum.to_hdf5 See Also -------- :py:meth:`~radis.spectrum.spectrum.Spectrum.from_hdf5` :py:func:`~radis.io.spec_hdf.spec2hdf` :py:meth:`~radis.spectrum.spectrum.Spectrum.store` """ from radis.io.spec_hdf import spec2hdf return spec2hdf(self, file, engine=engine) def to_pandas(self, copy=True): r"""Convert a Spectrum to a Pandas DataFrame Returns ------- pd.DataFrame : pandas DataFrame where columns are spectral arrays, and units are stored in attributes ``df.attrs`` Notes ----- Pandas does not support units yet. pint-pandas is an advanced project but not fully working. See discussion in https://github.com/pandas-dev/pandas/issues/10349 For the moment, we store units as metadata""" df = pd.DataFrame(self._q, copy=copy) df.attrs = self.units return df def to_specutils(self, var=None, wunit="default", Iunit="default"): r"""Convert a ``radis`` :py:class:`~radis.spectrum.spectrum.Spectrum` object to ``specutils`` :py:class:`specutils.spectra.spectrum1d.Spectrum1D` Parameters ---------- var: 'radiance', 'radiance_noslit', etc. which spectral array to convert. If ``None`` and only one spectral array is defined, use it. wunit: ``'nm'``, ``'cm'``, ``'nm_vac'``. wavespace unit: wavelength in air (``'nm'``), wavenumber (``'cm-1'``), or wavelength in vacuum (``'nm_vac'``). if ``"default"``, default unit for waveunit is used. See :py:meth:`~radis.spectrum.spectrum.Spectrum.get_waveunit`. .. note:: ``specutils`` handles wavelengths in vacuum only. If using ``'nm'`` wavelengths will be converted to wavelengths in vacuum. We recommend using ``'nm_air'`` or ``'nm_vac'`` to avoid any confusion. Iunit: unit for variable ``var`` if ``"default"``, default unit for quantity `var` is used. See the :py:attr:`~radis.spectrum.spectrum.Spectrum.units` attribute. For ``var="radiance"``, one can use per wavelength (~ 'W/m2/sr/nm') or per wavenumber (~ 'W/m2/sr/cm-1') units .. note:: `nan` that may have been added on the wings of the spectra if a slit has been applied are removed using ``trim_nan`` parameter of :py:meth:`~radis.spectrum.spectrum.Spectrum.get` . The waverange in the 1DSpectrum object may therefore be cropped. Examples -------- :: spectrum = s.to_specutils() Add uncertainties by reading a spectrum noise region with :py:class:`~specutils.SpectralRegion` and :py:func:`~specutils.manipulation.noise_region_uncertainty` :: from specutils import SpectralRegion from specutils.manipulation import noise_region_uncertainty noise_region = SpectralRegion(2012 / u.cm, 2009 / u.cm) spectrum = noise_region_uncertainty(spectrum, noise_region) Find lines out of the noise using :py:func:`~specutils.fitting.find_lines_threshold` :: from specutils.fitting import find_lines_threshold lines = find_lines_threshold(spectrum) .. minigallery:: radis.Spectrum.to_specutils See Also -------- :py:func:`~radis.spectrum.spectrum.Spectrum.from_specutils` """ try: from specutils.spectra import Spectrum1D except ModuleNotFoundError as err: raise ModuleNotFoundError( "Specutils is required to use this function." "Install it with:\n" " conda install -c conda-forge specutils" "\nor\n" " pip install specutils" ) from err meta = self.get_conditions().copy() if var is None: var = self._get_unique_var(operation_name="to_specutils") if wunit == "default": wunit = self.get_waveunit() if wunit in WAVELEN_UNITS: # wavelength units in air # AFAIK, specutil's Spectrum1D will only handle wavelengths as seen in vacuum # so below we request wavelengths in vac : wunit = "nm_vac" if wunit in WAVELEN_UNITS: # wavelength units in air raise ValueError( f"specutil's Spectrum1D will only handle wavelengths as seen in vacuum. Use one of `s.to_specutils(..., wunit={WAVELENVAC_UNITS}`)" ) # Update waveunit stored in conditions (in particular, update if wavelengths are in vacuum or air) if "waveunit" in meta: meta["waveunit"] = wunit if Iunit == "default": Iunit = self.units[var] w, I = self.get( var, wunit=wunit, Iunit=Iunit, copy=True, trim_nan=True, return_units=True ) # w, I are dimensioned arrays return Spectrum1D( flux=I, spectral_axis=w, meta=meta, ) def resample( self, w_new, unit="same", out_of_bounds="nan", energy_threshold="default", print_conservation=False, inplace=True, **kwargs, ): r"""Resample spectrum over a new wavelength/wavenumber range. .. warning:: This may result in information loss. Resampling is done with oversampling and spline interpolation. These parameters can be adjusted, and energy conservation ensured with the appropriate parameters. To minimize information loss, always resample the high-resolution spectrum over the low-resolution spectrum, i.e. :: s_highres.resample(s_lowres) Fills with ``'nan'`` or transparent medium (transmittance 1, radiance 0) when out of bound (see ``out_of_bounds``) Parameters ---------- w_new: array, or Spectrum new wavespace to resample the spectrum on. Must be inclosed in the current wavespace (we won't extrapolate) One can also give a Spectrum directly:: s1.resample(s2.get_wavenumber()) s1.resample(s2) # also valid unit: ``'same'``, ``'nm'``, ``'cm-1'``, ``'nm_vac'`` unit of new wavespace. It ``'same'`` it is assumed to be the current waveunit. Default ``'same'``. The spectrum waveunit is changed to this unit after resampling (i.e: a spectrum calculated and stored in `cm-1` but resampled in `nm` will be stored in `nm` from now on). If ``'nm'``, wavelength in air. If ``'nm_vac'``, wavelength in vacuum. out_of_bounds: ``'transparent'``, ``'nan'``, ``'error'`` what to do if resampling is out of bounds. ``'transparent'``: fills with transparent medium. 'nan': fill with nan. ``'error'``: raises an error. Default ``'nan'`` medium: ``'air'``, ``'vacuum'``, or ``'default'`` in which medium is the new waverange is calculated if it is given in 'nm'. Ignored if unit='cm-1' Other Parameters ---------------- energy_threshold: float or ``None`` or ``'default' if energy conservation (integrals on the intersecting range) is above this threshold, raise an error. If ``None``, dont check for energy conservation. If ``'default'``, look up the value in :py:attr:`radis.config` ["RESAMPLING_TOLERANCE_THRESHOLD"] Default ``'default'`` print_conservation: boolean if ``True``, prints energy conservation. Default ``False``. inplace: boolean if ``True``, modifies the Spectrum object directly. Else, returns a copy. Default ``True``. **kwargs: **dict all other arguments are sent to :func:`radis.misc.signal.resample` Returns ------- Spectrum : resampled Spectrum object. If using ``inplace=True``, the Spectrum object has been modified anyway. Examples -------- Convert a Spectrum in wavenumber to wavelengths:: s_nm = s.resample(s.get_wavelength(), "nm", inplace=False) .. minigallery:: radis.spectrum.spectrum.Spectrum.resample See Also -------- :func:`radis.misc.signal.resample`, :py:meth:`radis.spectrum.spectrum.Spectrum.resample_even` """ if inplace: s = self else: s = self.copy() # Get wavespace units stored_waveunit = s.get_waveunit() # spectrum unit if unit == "same": # resampled unit unit = stored_waveunit else: unit = cast_waveunit(unit) # Get output wavespace (it's the w_new array, unless w_new is a Spectrum) if isinstance(w_new, Spectrum): wstep_new = w_new.conditions.get("wstep", None) if unit == "nm": w_new = w_new.get_wavelength(medium="air") elif unit == "nm_vac": w_new = w_new.get_wavelength(medium="vacuum") elif unit == "cm-1": w_new = w_new.get_wavenumber() else: raise ValueError(unit) else: # wavespace already given as array: w_new = w_new diff = np.diff(w_new) if not np.all(diff == diff[0]): from radis.misc.warning import default_warning_status, warn warn( "When resampling the spectrum, the new waverange had unequal spacing.", "UnevenWaverangeWarning", status=default_warning_status, ) wstep_new = f"Unequal spacing, average ={np.mean(diff):.2f}" else: wstep_new = diff[0] # Get current waverange in output unit -> w if unit == "nm": w = s.get_wavelength(medium="air") elif unit == "nm_vac": w = s.get_wavelength(medium="vacuum") elif unit == "cm-1": w = s.get_wavenumber() else: raise ValueError(f"Unknown unit: {unit}") # Update stored_waveunit to new unit if unit != stored_waveunit: s.conditions["waveunit"] = unit if wstep_new is not None: s.conditions["wstep"] = wstep_new # Now let's resample def get_filling(variable): """Get out of bounds values for spectral quantity `variable`""" if out_of_bounds == "transparent": # Fill with optically transparent medium if variable in ["transmittance", "transmittance_noslit"]: fill_with = 1 else: fill_with = 0 elif out_of_bounds == "nan": fill_with = "nan" elif out_of_bounds == "error": fill_with = "error" else: raise ValueError(f"Unexpected value for out_of_bound: {out_of_bounds}") return fill_with # There are different cases depending on the unit of w_new # ... Note @devs: we're looping over dictionaries directly rather than # ... using the (safer) .get() function because it's much faster (the # ... air2vacuum conversion in particular is quite slow, but has been # ... done once for all with get_wavelength() above ) for k, I in s._q.items(): if k == "wavespace": continue fill_with = get_filling(k) Inew = resample( w, I, w_new, ext=fill_with, energy_threshold=energy_threshold, print_conservation=False, **kwargs, ) s._q[k] = Inew # update wavespace s._q["wavespace"] = w_new return s def resample_even( self, energy_threshold=5e-3, print_conservation=False, inplace=True ): r"""Resample spectrum over the same waverange, but evenly spaced. .. warning:: This may result in information loss. Resampling is done with oversampling and spline interpolation. These parameters can be adjusted, and energy conservation ensured with the appropriate parameters. Parameters ---------- energy_threshold: float or ``None`` if energy conservation (integrals on the intersecting range) is above this threshold, raise an error. If ``None``, dont check for energy conservation Default 5e-3 (0.5%) print_conservation: boolean if ``True``, prints energy conservation. Default ``False``. inplace: boolean if ``True``, modifies the Spectrum object directly. Else, returns a copy. Default ``True``. **kwargs: **dict all other arguments are sent to :func:`radis.misc.signal.resample` Returns ------- Spectrum : resampled Spectrum object. If using ``inplace=True``, the Spectrum object has been modified anyway. Examples -------- .. minigallery:: radis.spectrum.spectrum.Spectrum.resample See Also -------- :func:`radis.misc.signal.resample_even`, :py:meth:`radis.spectrum.spectrum.Spectrum.resample` """ if inplace: s = self else: s = self.copy() w = s._q["wavespace"] # wavelength or wavenumber range for k, I in s._q.items(): if k == "wavespace": continue w_new, Inew = resample_even( w, I, resfactor=2, print_conservation=print_conservation ) s._q[k] = Inew # update wavespace s._q["wavespace"] = w_new return s # %% Fonctions to fit a lineshape model # def fit_model( self, model, plot=False, confidence=0.9545, verbose=False, debug=False ): """Fit a lineshape model to the spectrum. The model can be a simple lineshape model among :py:class:`~astropy.modeling.models.Gaussian1D`, :py:class:`~astropy.modeling.models.Lorentz1D`, :py:class:`~astropy.modeling.models.Voigt1D`, or a combination of multiple models. Parameters ---------- model: astropy.modeling.Model or a list of models model to fit to the spectrum. If a list is given, a sum of models is fitted plot: bool if True, plot the difference between the model and the spectrum Other Parameters ---------------- confidence: 0.6827, 0.9545, or 0.9973 confidence interval to use. Returns ------- g_fit the fitted model y_err uncertainty on the fitted parameters calculated as the square root of the diagonal of the covariance matrix Examples -------- :: from astropy.modeling import models s = radis.spectrum_test().crop(2201.7, 2205.1) g_fit, y_err = s.fit_model(models.Lorentz1D(), plot=True) Example with 6 Voigt lines: :: from astropy.modeling import models s = radis.spectrum_test().crop(2201.7, 2225.1) g_fit_list, y_err = s.fit_model([models.Voigt1D() for _ in range(6)], plot=True) Other Examples -------------- .. minigallery:: radis.Spectrum.fit_model """ # TODO : return error on fitted parameters by default from astropy.modeling.fitting import LevMarLSQFitter from specutils.fitting import fit_lines # make the model a list of models : if not (isinstance(model, list) or isinstance(model, tuple)): model = [model] # from astropy.modeling import models # from astropy.modeling.core import Fittable1DModel # assert isinstance(model, Fittable1DModel) fitter = LevMarLSQFitter(calc_uncertainties=True) w_fit = self.get(self.get_vars()[0], return_units=True)[0] # Initialize model # ---------------- # Fit the spectrum and calculate the fitted flux values (``y_fit``) # ... make some reasonable first assumption on amplitude, center, fwhm (if available) # If multiple models, initialize them sequentially and remove the fitted # values from the spectrum before fitting the next model # fitter0 = LevMarLSQFitter() if len(model) > 1: s0 = self.copy(copy_lines=False, copy_arrays=True) else: s0 = self for i, mod in enumerate(model): # Note from Astropy Modeling: # > ("The .value property on parameters should be set" # " to unitless values, not Quantity objects. To set" # "a parameter to a quantity simply set the " # "parameter directly without using .value") # If values are set to the default, then make a safe guess : if "amplitude" in mod.param_names: if ( mod.__getattribute__("amplitude").value == mod.__getattribute__("amplitude").default ): mod.__setattr__("amplitude", s0.max()) # dimensioned if "amplitude_L" in mod.param_names: # Lorentzian amplitude (?) in Voigt mode if ( mod.__getattribute__("amplitude_L").value == mod.__getattribute__("amplitude_L").default ): mod.__setattr__("amplitude_L", s0.max()) # dimensioned if "x_0" in mod.param_names: if ( mod.__getattribute__("x_0").value == mod.__getattribute__("x_0").default ): mod.__setattr__("x_0", s0.argmax()) # dimensioned def get_fwhm_fast(s0): """a very rough 1st-guess for HWHM: we measure the distance before the intensity drops below half the maximum """ w, I = s0.get(s0.get_vars()[0]) mask_below_half = (I - I.min()) < (I.max() - I.min()) / 2 i_fwhm_left = np.argmin( ((w[mask_below_half] - s0.argmax(value_only=True))) ** 2 ) if i_fwhm_left + 1 > len(w[mask_below_half]): fwhm = abs(w[-1] - w[0]) else: fwhm = abs( w[mask_below_half][i_fwhm_left + 1] - w[mask_below_half][i_fwhm_left] ) return fwhm if "fwhm" in mod.param_names: if ( mod.__getattribute__("fwhm").value == mod.__getattribute__("fwhm").default ): mod.__setattr__("fwhm", get_fwhm_fast(s0)) # NOT dimensioned TODO if "fwhm_L" in mod.param_names and "fwhm_G" in mod.param_names: # Lorentzian and Gaussian width in Voigt profile if ( mod.__getattribute__("fwhm_L").value == mod.__getattribute__("fwhm_L").default and mod.__getattribute__("fwhm_G").value == mod.__getattribute__("fwhm_G").default ): # distribute FWHM on the two profiles: fwhm = get_fwhm_fast(s0) mod.__setattr__("fwhm_L", fwhm / 2) # NOT dimensioned TODO mod.__setattr__("fwhm_G", fwhm / 2) # NOT dimensioned TODO if len(model) > 1: # Remove fitted model from spectrum, and initialize the next model g_one_mod = fit_lines(s0.to_specutils(), mod, fitter=fitter) y_one_mod = g_one_mod(w_fit) s0 -= y_one_mod if debug: s0.plot() import matplotlib.pyplot as plt plt.title(f"s0 after iteration {i} of guessing initial parameters") plt.tight_layout() plt.plot(w_fit, y_one_mod, "--r") # Final fit # --------- # This time we fit the full spectrum ('self') with all models (model) from astropy.modeling.fitting import DEFAULT_MAXITER maxiter = len(model) * DEFAULT_MAXITER # arbitrary g_fit_list = fit_lines( self.to_specutils(), model, fitter=fitter, maxiter=maxiter ) y_fit_list = [g_fit(w_fit) for g_fit in g_fit_list] if plot: import matplotlib.pyplot as plt plt_line = self.plot( lw=5, color="grey", show=False, # needed when using `inline` ploting (e.g. default in Spyder) ) ax = plt_line.figure.axes[0] for i, y_fit in enumerate(y_fit_list): g_fit = g_fit_list[i] label = " ".join( [ f"{g_fit.param_names[k]}"[:6] + f"={v:.2f}" for k, v in enumerate(g_fit.parameters) ] ) ax.plot(w_fit, y_fit, label=label) # plot legend with small font size: ax.legend(fontsize=14) plt.show() # needed if `plt.isinteractive() = False` if verbose: print(fitter.fit_info["ierr"], fitter.fit_info["message"]) # Get uncertainties if fitter.fit_info["param_cov"] is not None: cov_matrix_diagonal = fitter.fit_info["param_cov"].diagonal() std = np.sqrt(cov_matrix_diagonal) if confidence == 0.6827: y_err = 1 * std elif confidence == 0.9545: y_err = 2 * std elif confidence == 0.9973: y_err = 3 * std else: raise ValueError("Use confidence as one of 0.6827, 0.9545, or 0.9973") # Add units y_err = [ ( mod.__getattribute__(param).unit * y_err[i] if mod.__getattribute__(param).unit else y_err[i] ) for mod in model for i, param in enumerate(mod.param_names) ] # TODO refactor there is probably a better way to write it else: y_err = None # Compute area under the curve of models.Voigt1D or models.Gaussian1D or models.Lorentz1D from astropy.modeling import models from scipy.integrate import trapezoid for index, line in enumerate(g_fit_list): if isinstance(line, models.Voigt1D): area = np.abs(trapezoid(y_fit_list[index], w_fit)) # slow, but accurate elif isinstance(line, models.Gaussian1D): area = line.amplitude * line.stddev * np.sqrt(2 * np.pi) elif isinstance(line, models.Lorentz1D): area = line.amplitude * np.pi * line.fwhm / 2 else: raise ValueError( f"Unexpected model type: {type(line).__name__}. \nExpected `models.Voigt1D` or `Gaussian1D` or `Lorentz1D`" ) # return area g_fit_list[index].area = area return g_fit_list, y_err # %% ====================================================================== # Semi public functions # ---------------- # Access object properties (but don't manipulate the spectrum itself) # XXX ===================================================================== def get_waveunit(self): r"""Returns whether this spectrum is defined in wavelength (nm) or wavenumber (cm-1)""" return self.conditions["waveunit"] def is_at_equilibrium(self, check="warn", verbose=False): r"""Returns whether this spectrum is at (thermal) equilibrium. Reads the ``thermal_equilibrium`` key in Spectrum conditions. It does not imply chemical equilibrium (mole fractions are still arbitrary) If they are defined, also check that the following assertions are True: Tvib = Trot = Tgas self_absorption = True overpopulation = None If they are not, still trust the value in Spectrum conditions, but raise a warning. Other Parameters ---------------- check: ``'warn'``, ``'error'``, ``'ignore'`` what to do if Spectrum conditions dont match the given equilibrium state: raise a warning, raise an error, or just ignore and dont even check. Default ``'warn'``. verbose: bool if ``True``, print why is the spectrum is not at equilibrium, if applicable. """ conditions = self.conditions # Get value try: equilibrium = conditions["thermal_equilibrium"] except KeyError: raise KeyError( "We need to know if Spectrum is at equilibrium, but " + "`thermal_equilibrium` is not defined in conditions. Please add the " + "value manually with s.conditions['thermal_equilibrium']=..." ) if check == "ignore": return equilibrium # Check output match the rest of the spectrum conditions try: assert conditions["Tgas"] != r"N/A" if "Tvib" in conditions: assert conditions["Tvib"] == conditions["Tgas"] if "Trot" in conditions: assert conditions["Trot"] == conditions["Tgas"] if "overpopulation" in conditions: assert conditions["overpopulation"] is None if "self_absorption" in conditions: assert conditions["self_absorption"] # is True except AssertionError: guess = False if verbose: # Print which equilibrium test failed print("Spectrum not at equilibrium because the following test failed:") import sys import traceback _, _, tb = sys.exc_info() tb_info = traceback.extract_tb(sys.exc_info()[2]) print(tb_info[-1][-1]) # @dev: see https://stackoverflow.com/a/11587247/5622825 except KeyError as err: warn(f"Condition missing to know if spectrum is at equilibrium: {err}") guess = not equilibrium else: guess = True if equilibrium != guess: msg = ( f"Declared value of equilibrium ({equilibrium}) does not match the inferred one ({guess})" + ". Update your Spectrum conditions" ) if check == "warn": warn(msg) elif check == "error": raise AssertionError(msg) return equilibrium def is_optically_thin(self): r"""Returns whether the spectrum is optically thin, based on the value on the self_absorption key in conditions. If not given, raises an error """ try: return not self.conditions["self_absorption"] except KeyError: raise KeyError( "We need to know if Spectrum is optically thin, but " + "`self_absorption` is not defined in conditions. Please add the " + "value manually with s.conditions['self_absorption']=..." ) def copy(self, copy_lines=True, quantity="all", copy_arrays=True): r"""Returns a copy of this Spectrum object (performs a smart deepcopy) Parameters ---------- copy_lines: bool default ``True`` quantity: 'all', or one of 'radiance_noslit', 'absorbance', etc. if not 'all', copy only one quantity. Default ``'all'`` copy_arrays: bool if ``False``, returned array's quantity is a pointer to the original Spectrum. Faster, but warning, changing them will then change the original Spectrum. Default ``True`` Examples -------- .. minigallery:: radis.spectrum.spectrum.Spectrum.copy :add-heading: """ try: return self.__copy__( copy_lines=copy_lines, quantity=quantity, copy_arrays=copy_arrays ) except MemoryError: raise MemoryError( "during copy of Spectrum. If you don't need them, " + "dropping lines before copying may save a lot of space: " + "del s.lines ; or, use copy_lines=False" ) def __copy__(self, copy_lines=True, quantity="all", copy_arrays=True): r"""Generate a new spectrum object. Note: using deepcopy would work but then the Spectrum object would be pickled and unpickled again. It's a little faster here Parameters ---------- copy_lines: bool default ``True`` quantity: 'all', or one of 'radiance_noslit', 'absorbance', etc. if not 'all', copy only one quantity. Default ``'all'`` copy_arrays: bool if ``False``, returned array's quantity is a pointer to the original Spectrum. Faster, but warning, changing them will then change the original Spectrum. Default ``True`` Notes ----- Performance: deepcopy: 3.32 ms initial : 10 ms no check test, optimised: 1.8 ms ... asymptote: without evenly spaced check, without copies: 1.84 ms """ # quantities = {s:(v[0].copy(), v[1].copy()) for (s,v) in self.items()} #╪ 1.8 ms # quantities = dict(self.items()) # 912 ns, not a copy but no need as # # Spectrum() recreates a copy anyway if quantity == "all": quantities = dict(self._get_items()) elif quantity in self.get_vars(): quantities = { quantity: self.get( quantity, wunit=self.get_waveunit(), Iunit=self.units[quantity], copy=copy_arrays, ) } else: self.update(quantity) # needed to get self.units[quantity] quantities = { quantity: self.get( quantity, wunit=self.get_waveunit(), Iunit=self.units[quantity], copy=copy_arrays, ) # copy=True still needed as the waverange is the one of the original array } del self._q[quantity] # no need to keep it in this Spectrum # del self.units[quantity] # we can keep the units (+ they'll be copied below) try: units = deepcopy(self.units) except AttributeError: units = None conditions = deepcopy(self.conditions) # 39 us try: cond_units = deepcopy(self.cond_units) except AttributeError: cond_units = None lines = None if copy_lines: try: # 143 us (2 ms with deepcopy(lines)) lines = self.lines.copy(deep=True) except AttributeError: pass try: populations = self.populations except AttributeError: populations = None references = self.references.copy() waveunit = self.get_waveunit() # 163 ns name = self.name # Generate copied Spectrum s = Spectrum( # 1.51 ms quantities=quantities, conditions=conditions, cond_units=cond_units, populations=populations, lines=lines, units=units, wunit=waveunit, references=references, name=name, check_wavespace=False, # no need, as Spectrum was normally already properly defined # saves about 3.5 ms on the Performance test object ) # Add extra information # ... file name (if exists) s.file = self.file # ... slit information try: wslit, Islit = self.get_slit() s._slit["wavespace"] = wslit # in 'waveunit' s._slit["intensity"] = Islit except KeyError: # no slit to apply pass return s def compare_with( self, other, spectra_only=False, plot=True, wunit="default", verbose=True, rtol=1e-5, ignore_nan=False, ignore_outliers=False, normalize=False, **kwargs, ): r"""Compare Spectrum with another Spectrum object. Parameters ---------- other: type Spectrum another Spectrum to compare with spectra_only: boolean, or str if ``True``, only compares spectral quantities (in the same waveunit) and not lines or conditions. If str, compare a particular quantity name. If False, compare everything (including lines and conditions and populations). Default ``False`` plot: boolean if ``True``, use plot_diff to plot all quantities for the 2 spectra and the difference between them. Default ``True``. wunit: ``"nm"``, ``"cm-1"``, ``"default"`` in which wavespace to compare (and plot). If ``"default"``, natural wavespace of first Spectrum is taken. rtol: float relative difference to use for spectral quantities comparison ignore_nan: boolean if ``True``, nans are ignored when comparing spectral quantities ignore_outliers: boolean, or float if not False, outliers are discarded. i.e, output is determined by:: out = (~np.isclose(I, Ie, rtol=rtol, atol=0)).sum()/len(I) < ignore_outliers normalize: bool Normalize the spectra to be plotted Other Parameters ---------------- kwargs: dict arguments are forwarded to :func:`~radis.spectrum.compare.plot_diff` Returns ------- equals: boolean return True if spectra are equal (respective to tolerance defined by rtol and other input conditions) Examples -------- Compare two Spectrum objects, or specifically the transmittance:: s1.compare_with(s2) s1.compare_with(s2, 'transmittance') Note that you can also simply use `s1 == s2`, that uses :meth:`~radis.spectrum.spectrum.Spectrum.compare_with` internally:: s1 == s2 # will return True or False .. minigallery:: radis.spectrum.spectrum.Spectrum.compare_with :add-heading: See Also -------- :func:`~radis.spectrum.compare.compare_spectra` """ from radis.spectrum.compare import compare_spectra return compare_spectra( self, other, spectra_only=spectra_only, plot=plot, wunit=wunit, verbose=verbose, rtol=rtol, ignore_nan=ignore_nan, ignore_outliers=ignore_outliers, normalize=normalize, **kwargs, ) def cite(self, format="bibentry"): r"""Prints bibliographic references used to compute this spectrum, as stored in the :py:attr:`~radis.spectrum.spectrum.Spectrum.references` dictionary. Default references known to RADIS are listed in :py:data:`radis.db.references.doi`. Parameters ---------- format: default ``'bibentry'``. See more in :py:func:`habanero.content_negotiation` Examples -------- :: from radis import calc_spectrum s = calc_spectrum( 1900, 2300, # cm-1 molecule="CO", isotope="1,2,3", pressure=1.01325, # bar Tvib=2000, # Trot=300, mole_fraction=0.1, path_length=1, # cm databank="hitran", ) s.cite() .. raw:: html
Returns : .. code-block:: python Used for algorithm ------------------ @article{van_den_Bekerom_2021, doi = {10.1016/j.jqsrt.2020.107476}, url = {https://doi.org/10.1016%2Fj.jqsrt.2020.107476}, year = 2021, month = {mar}, publisher = {Elsevier {BV}}, volume = {261}, pages = {107476}, author = {D.C.M. van den Bekerom and E. Pannier}, title = {A discrete integral transform for rapid spectral synthesis}, journal = {Journal of Quantitative Spectroscopy and Radiative Transfer} } Used for calculation, rovibrational energies -------------------------------------------- @article{Pannier_2019, doi = {10.1016/j.jqsrt.2018.09.027}, url = {https://doi.org/10.1016%2Fj.jqsrt.2018.09.027}, year = 2019, month = {jan}, publisher = {Elsevier {BV}}, volume = {222-223}, pages = {12--25}, author = {Erwan Pannier and Christophe O. Laux}, title = {{RADIS}: A nonequilibrium line-by-line radiative code for {CO}2 and {HITRAN}-like database species}, journal = {Journal of Quantitative Spectroscopy and Radiative Transfer} } Used for data retrieval ----------------------- @article{Ginsburg_2019, doi = {10.3847/1538-3881/aafc33}, url = {https://doi.org/10.3847%2F1538-3881%2Faafc33}, year = 2019, month = {feb}, publisher = {American Astronomical Society}, volume = {157}, number = {3}, pages = {98}, author = {Adam Ginsburg and Brigitta M. Sip{\H{o}}cz and C. E. Brasseur and Philip S. Cowperthwaite and Matthew W. Craig and Christoph Deil and James Guillochon and Giannina Guzman and Simon Liedtke and Pey Lian Lim and Kelly E. Lockhart and Michael Mommert and Brett M. Morris and Henrik Norman and Madhura Parikh and Magnus V. Persson and Thomas P. Robitaille and Juan-Carlos Segovia and Leo P. Singer and Erik J. Tollerud and Miguel de Val-Borro and Ivan Valtchanov and Julien Woillez and}, title = {astroquery: An Astronomical Web-querying Package in Python}, journal = {The Astronomical Journal} } Used for line database ---------------------- @article{Gordon_2017, doi = {10.1016/j.jqsrt.2017.06.038}, url = {https://doi.org/10.1016%2Fj.jqsrt.2017.06.038}, year = 2017, month = {dec}, publisher = {Elsevier {BV}}, volume = {203}, pages = {3--69}, author = {I.E. Gordon and L.S. Rothman and C. Hill and R.V. Kochanov and Y. Tan and P.F. Bernath and M. Birk and V. Boudon and A. Campargue and K.V. Chance and B.J. Drouin and J.-M. Flaud and R.R. Gamache and J.T. Hodges and D. Jacquemart and V.I. Perevalov and A. Perrin and K.P. Shine and M.-A.H. Smith and J. Tennyson and G.C. Toon and H. Tran and V.G. Tyuterev and A. Barbe and A.G. Cs{\'{a}}sz{\'{a}}r and V.M. Devi and T. Furtenbacher and J.J. Harrison and J.-M. Hartmann and A. Jolly and T.J. Johnson and T. Karman and I. Kleiner and A.A. Kyuberis and J. Loos and O.M. Lyulin and S.T. Massie and S.N. Mikhailenko and N. Moazzen-Ahmadi and H.S.P. Müller and O.V. Naumenko and A.V. Nikitin and O.L. Polyansky and M. Rey and M. Rotger and S.W. Sharpe and K. Sung and E. Starikova and S.A. Tashkun and J. Vander Auwera and G. Wagner and J. Wilzewski and P. Wcis{\l}o and S. Yu and E.J. Zak}, title = {The {HITRAN}2016 molecular spectroscopic database}, journal = {Journal of Quantitative Spectroscopy and Radiative Transfer} } Used for partition function --------------------------- @article{Gamache_2021, doi = {10.1016/j.jqsrt.2021.107713}, url = {https://doi.org/10.1016%2Fj.jqsrt.2021.107713}, year = 2021, month = {sep}, publisher = {Elsevier {BV}}, volume = {271}, pages = {107713}, author = {Robert R. Gamache and Bastien Vispoel and Michaël Rey and Andrei Nikitin and Vladimir Tyuterev and Oleg Egorov and Iouli E. Gordon and Vincent Boudon}, title = {Total internal partition sums for the {HITRAN}2020 database}, journal = {Journal of Quantitative Spectroscopy and Radiative Transfer} } @article{Kochanov_2016, doi = {10.1016/j.jqsrt.2016.03.005}, url = {https://doi.org/10.1016%2Fj.jqsrt.2016.03.005}, year = 2016, month = {jul}, publisher = {Elsevier {BV}}, volume = {177}, pages = {15--30}, author = {R.V. Kochanov and I.E. Gordon and L.S. Rothman and P. Wcis{\l}o and C. Hill and J.S. Wilzewski}, title = {{HITRAN} Application Programming Interface ({HAPI}): A comprehensive approach to working with spectroscopic data}, journal = {Journal of Quantitative Spectroscopy and Radiative Transfer} } Used for spectroscopic constants -------------------------------- @article{Guelachvili_1983, doi = {10.1016/0022-2852(83)90203-5}, url = {https://doi.org/10.1016%2F0022-2852%2883%2990203-5}, year = 1983, month = {mar}, publisher = {Elsevier {BV}}, volume = {98}, number = {1}, pages = {64--79}, author = {G. Guelachvili and D. de Villeneuve and R. Farrenq and W. Urban and J. Verges}, title = {Dunham coefficients for seven isotopic species of {CO}}, journal = {Journal of Molecular Spectroscopy} } .. raw:: html
Other Examples -------------- .. minigallery:: radis.spectrum.spectrum.Spectrum.cite See Also -------- :py:class:`~radis.tools.track_ref.RefTracker`, :py:data:`~radis.db.references.doi` """ if self.references is None: raise ValueError( "You first need to register the dictionary of bibliographic references. Set `s.references= {'doi':'use in the calculation'}` " ) if not isinstance(self.references, RefTracker): self.references = RefTracker(**self.references) return self.references.cite(format=format) # %% ====================================================================== # Private functions # ---------------- # XXX ===================================================================== def _init_annotations(self): r"""Annotations are used to give typing hints for get() and plot() functions, based on what quantities are available in the Spectrum object.""" from radis.tools.slit import SLIT_SHAPES try: # Python >3.6 only self.get.__annotations__["var"] = [] self.plot.__annotations__["var"] = [] self.apply_slit.__annotations__["shape"] = SLIT_SHAPES except AttributeError: pass # old Python version def _add_quantity(self, name, w, I, check_wavespace=True): r"""Add quantity. Note: creates a copy of the input array If check_wavespace, check that array matches existing wavespace, and that new wavespace is evenly spaced (required for the slit) Note: this check takes a lot of time! (few ms) """ assert len(w) == len(I) # Add wavespace if "wavespace" in self._q: if check_wavespace: # Check new wavespace match the existing one if not np.allclose(w, self._q["wavespace"]): raise ValueError( f"wavespace for {name} doesnt correspond to existing wavespace" + " for non convoluted quantities" ) else: pass # nothing to be done else: if name in CONVOLUTED_QUANTITIES: self._q["wavespace"] = np.array(w) # copy # no need to check if wavespace is evenly spaced: we won't # apply the slit function again else: # Check Wavelength/wavenumber is evenly spaced if check_wavespace: if not evenly_distributed(w, atolerance=1e-5): warn( f"Wavespace is not evenly spaced ({np.abs(np.diff(w)).max() / w.mean() * 100:.3f}%) for {name}." + " This may create problems if later convolving with slit function (`s.apply_slit()`). You can use `s.resample_even()`" ) self._q["wavespace"] = np.array(w) # copy # Add quantity itself self._q[name] = np.array(I) # copy # also make the quantity accessible with s.[name] like Pandas dataframes (removed eventually) # setattr(self, name, quantity) # Warning this makes another copy of it (it's a tuple!) # add to annotations (Python >3.6) try: self.get.__annotations__["var"].append(name) self.plot.__annotations__["var"].append(name) except AttributeError: # old Python version pass def __eq__(self, other): r"""Override the default Equals behavior.""" return self.compare_with(other, spectra_only=True, verbose=False, plot=False) def __ne__(self, other): r"""Define a non-equality test.""" return not self.__eq__(other) def __dir__(self): r"""Names shown with tab completion: remove certain attributes to simplify the use of this class (@minou).""" # attrs = super(Spectrum, self).__dir__() attrs = dir(type(self)) # Python 2 and 3 compatible exclude = [ "clear", "fromkeys", "items", "pop", "popitem", "setdefault", "values", ] return [k for k in attrs if not k in exclude] def __str__(self): r"""Print all Spectrum attributes.""" # Print name print("Spectrum Name: ", self.get_name()) # Print spectral quantities print("Spectral Quantities") print("-" * 40) for k, v in self._q.items(): if k == "wavespace": continue # print number of points with a comma separator unit_str = ( f"\t[{self.units[k]}]" if (k in self.units and self.units[k] != "") else "" ) nan_str = f", {count_nans(v)} nans" if anynan(v) else "" print(f" {k}{unit_str}\t({len(v):,d} points{nan_str})") # Print populations if self.populations: print("Populations Stored") print("-" * 40) try: for k, v in self.populations.items(): print(" " * 2, k, "\t\t", list(v.keys())) except Exception: pass # Print conditions if self.conditions: self.print_conditions() return "" # self.print_conditions() def take(self, var, copy_lines=False, copy_arrays=True): r""" Parameters ---------- var : str spectral quantity (``'absorbance'``, ``'transmittance'``, ``'xsection'`` etc.) Other Parameters ---------------- copy_lines: bool if ``True``, export ``s.lines``. Default ``False`` copy_arrays: bool if ``False``, returned array's quantity is a pointer to the original Spectrum. Faster, but warning, changing them will then change the original Spectrum. Default ``True`` Returns ------- s: Spectrum same Spectrum with only the `var` spectral quantity Examples -------- Use ``take`` to chain other commands :: s.take('radiance').normalize().plot() .. minigallery:: radis.spectrum.spectrum.Spectrum.take :add-heading: """ return self.copy(quantity=var, copy_lines=copy_lines, copy_arrays=copy_arrays) # %% Add min, max, normalize operations def _get_unique_var(self, operation_name="algebraic"): quantities = self.get_vars() if len(quantities) > 1: raise ValueError( f"There is an ambiguity with the Spectrum `{operation_name}()` operation. " + f"There are currently multiple spectral arrays in the Spectrum {self.get_name()} ({ ', '.join(self.get_vars())})\n\n" + f"Use `s.take('transmittance').{operation_name}()` or `s.take('radiance').{operation_name}()`" f", etc. to apply `{operation_name}()` to the one spectral array you want." ) elif len(quantities) == 0: raise ValueError(f"No spectral array defined in Spectrum {self.get_name()}") else: var = quantities[0] return var def max(self, value_only=False, return_wmax=False): r"""Maximum of the Spectrum, if only one spectral quantity is available:: s.max() Returns a dimensioned quantity by default, with the default unit of the spectrum. If ``value_only=False``, a float is returned, without dimensions. If there are multiple arrays in your Spectrum, use :py:meth:`~radis.spectrum.spectrum.Spectrum.take`, e.g. :: s.take('radiance').max() Parameters ---------- value_only: bool return_position: bool if True, returns wavelength or wavenumber of maximum, in the Spectrum unit. Examples -------- Normalize a spectrum (we could also have used :py:meth:`~radis.spectrum.spectrum.Spectrum.normalize` directly) :: s /= s.max() Below is a slightly more advanced use case, convenient to process non-calibrated experimental data. Imagine we want to remove a baseline that we have identified to originate from a particular species, modelled by spectrum ``s``, but whose intensity (and units) are not the same as that of the experimental spectrum. We substract exactly the intensity that correspond to the intensity on a given range ``w1, w2`` of the experimental spectrum ``s_exp``, using a combination of :py:meth:`~radis.spectrum.spectrum.Spectrum.take`, :py:meth:`~radis.spectrum.spectrum.Spectrum.normalize`, :py:meth:`~radis.spectrum.spectrum.Spectrum.crop` and :py:meth:`~radis.spectrum.spectrum.Spectrum.max`: :: s_exp -= s.take('radiance').normalize() * s_exp.crop((w1, w2)), inplace=False).max() Other Examples -------------- Get max position ``wmax`` :: s = radis.spectrum_test() max, wmax = s.max(return_wmax=True) Note that the later can also be achieved with :py:meth:`~radis.spectrum.spectrum.Spectrum.argmax` :: wmax = s.argmax() .. minigallery:: radis.Spectrum.max See Also -------- :py:meth:`~radis.spectrum.spectrum.Spectrum.argmax` """ # TODO refactor: value_only can be replaced by ~return_units, and used in get() var = self._get_unique_var(operation_name="max") w, I = self.get( var, wunit=self.get_waveunit(), Iunit=self.units[var], copy=False ) nan_mask = ~np.isnan(I) if return_wmax: imax = I[nan_mask].argmax() Imax = I[nan_mask][imax] wmax = w[nan_mask][imax] if value_only: return Imax, wmax else: return Imax * Unit(self.units[var]), wmax * Unit(self.get_waveunit()) else: Imax = I[nan_mask].max() if value_only: return Imax else: return Imax * Unit(self.units[var]) def argmax(self, value_only=False): r"""Return wave position of maximum of the Spectrum Equivalent to the following use of the function :py:meth:`~radis.spectrum.spectrum.Spectrum.max`:: s.max(value_only=value_only, return_wmax=True)[1] .. minigallery:: radis.Spectrum.argmax See Also -------- :py:meth:`~radis.spectrum.spectrum.Spectrum.max` """ return self.max(value_only=value_only, return_wmax=True)[1] def min(self, value_only=True, return_wmin=False): r"""Minimum of the Spectrum, if only one spectral quantity is available :: s.min() Returns a dimensioned quantity by default, with the default unit of the spectrum. If ``value_only=False``, a float is returned, without dimensions. If there are multiple arrays in your Spectrum, use :py:meth:`~radis.spectrum.spectrum.Spectrum.take`, e.g. :: s.take('radiance').min() Other Examples -------------- Get min position ``wmin`` :: s = radis.spectrum_test() min, wmin = s.min(return_wmin=True) Note that the later can also be achieved with :py:meth:`~radis.spectrum.spectrum.Spectrum.argmin` :: wmax = s.argmin() .. minigallery:: radis.Spectrum.min See Also -------- :py:meth:`~radis.spectrum.spectrum.Spectrum.argmin` """ var = self._get_unique_var(operation_name="min") w, I = self.get( var, wunit=self.get_waveunit(), Iunit=self.units[var], copy=False ) nan_mask = ~np.isnan(I) if return_wmin: imin = I[nan_mask].argmin() Imin = I[nan_mask][imin] wmin = w[nan_mask][imin] if value_only: return Imin, wmin else: return Imin * Unit(self.units[var]), wmin * Unit(self.get_waveunit()) else: Imin = I[nan_mask].min() if value_only: return Imin.min() else: return Imin.min() * Unit(self.units[var]) def argmin(self, value_only=False): r"""Return wave position of minimum of the Spectrum Equivalent to the following use of the function :py:meth:`~radis.spectrum.spectrum.Spectrum.min`:: s.min(value_only=value_only, return_wmin=True)[1] .. minigallery:: radis.Spectrum.argmin See Also -------- :py:meth:`~radis.spectrum.spectrum.Spectrum.min` """ return self.min(value_only=value_only, return_wmin=True)[1] def normalize( self, normalize_how="max", wrange=(), wunit=None, inplace=False, force=False, verbose=False, return_norm=False, ): r"""Normalise the Spectrum, if only one spectral quantity is available. Parameters ---------- normalize_how: ``'max'``, ``'area'``, ``'mean'`` how to normalize. ``'max'`` is the default but may not be suited for very noisy experimental spectra. ``'area'`` will normalize the integral to 1. ``'mean'`` will normalize by the mean amplitude value wrange: tuple of floats, float normalize on this range if a tuple, else normalise at that point. If an empty tuple, use the entire range. wunit: ``"nm"``, ``"cm-1"``, ``"nm_vac"`` unit of the normalisation range above. If ``None``, use the spectrum default waveunit. inplace: bool if ``True``, changes the Spectrum. Other Parameters ---------------- force: boolean By default, normalizing some parameters such as transmittance is forbidden because considered non-physical. Use force=True if you really want to. Examples -------- :: s.normalize("max", (4200, 4800), inplace=True).plot() If there are multiple arrays in your Spectrum, use :py:meth:`~radis.spectrum.spectrum.Spectrum.take`, e.g. :: s.take('radiance').normalize(wrange=(4200, 4800)).plot() .. minigallery:: radis.spectrum.spectrum.Spectrum.normalize """ from radis.spectrum.operations import multiply var = self._get_unique_var(operation_name="normalize") if var in ["transmittance", "transmittance_noslit"] and not force: raise ValueError( f"Cannot normalize {var}. Use force=True if you really want." ) s = self if wunit is None: wunit = s.get_waveunit() if isinstance(wrange, (int, float)): # Use user provided point. Only makes sense to normalize using # normalize_how='max', otherwise raise an error. w, I = s.get( var, wunit=wunit, Iunit=s.units[var], copy=False ) # (faster not to copy) if normalize_how == "max": # Find the value of I closest to the given w norm = I[np.argmin(np.abs(w - wrange))] norm_unit = s.units[var] elif normalize_how == "mean" or normalize_how == "area": raise ValueError( f"Unexpected `normalize_how`: {normalize_how}. " "For a single value wrange only `normalize_how='max'` is " "defined." ) else: raise ValueError(f"Unexpected `normalize_how`: {normalize_how}") out = multiply(s, 1 / (norm * Unit(norm_unit)), inplace=inplace) elif len(wrange) == 0: # Use entire wavelength range if normalize_how == "max": norm = np.nanmax( s.get(var, wunit=wunit, Iunit=s.units[var], copy=False)[1] ) norm_unit = s.units[var] elif normalize_how == "mean": norm = np.nanmean( s.get(var, wunit=wunit, Iunit=s.units[var], copy=False)[1] ) norm_unit = s.units[var] elif normalize_how == "area": w, I = s.get(var, wunit=wunit, Iunit=s.units[var], copy=False) norm = nantrapz(I, w) norm_unit = Unit(s.units[var]) * Unit(wunit) else: raise ValueError(f"Unexpected `normalize_how`: {normalize_how}") # Ensure we use the same unit system! out = multiply(s, 1 / (norm * Unit(norm_unit)), inplace=inplace) elif len(wrange) == 2: # Use the user provided range wmin, wmax = wrange w, I = s.get( var, wunit=wunit, Iunit=s.units[var], copy=False ) # (faster not to copy) b = (w > wmin) & (w < wmax) if normalize_how == "max": norm = np.nanmax(I[b]) norm_unit = s.units[var] elif normalize_how == "mean": norm = np.nanmean(I[b]) norm_unit = s.units[var] elif normalize_how == "area": norm = np.abs(nantrapz(I[b], w[b])) norm_unit = Unit(s.units[var]) * Unit(wunit) else: raise ValueError(f"Unexpected `normalize_how`: {normalize_how}") out = multiply(s, 1 / (norm * Unit(norm_unit)), inplace=inplace) else: raise ValueError(f"Unexpected `wrange`: {wrange}") if verbose: print(f"Normalization factor : {norm}") if return_norm: return out, norm * Unit(norm_unit) return out def get_baseline(self, algorithm="als", **kwargs): r""" Calculate and returns a baseline Parameters ---------- s: Spectrum Spectrum which needs a baseline var: str on which spectral quantity to read the baseline. Default ``'radiance'``. See :py:data:`~radis.spectrum.utils.SPECTRAL_QUANTITIES` algorithm: 'als', 'polynomial' Asymmetric least square or Polynomial fit **kwargs: dict additional parameters to send to the algorithm. By default, for 'polynomial': **kwargs = **{"deg": 1, "max_it":500} for 'als': **kwargs = {"asymmetry_param": 0.05, "smoothness_param": 1e6} Returns ------- baseline: Spectrum Spectrum object where intensity is the baseline of s Examples -------- .. minigallery:: radis.Spectrum.get_baseline See also -------- :py:func:`~radis.spectrum.operations.get_baseline` """ from radis.spectrum.operations import get_baseline var = self._get_unique_var(operation_name="get_baseline") return get_baseline(self, var, **kwargs) # %% Performance / profilers def print_perf_profile(self, number_format="{:.3f}", precision=16): r"""Prints Profiler output dictionary in a structured manner. Example ------- :: Spectrum.print_perf_profile() # output >> spectrum_calculation 0.189s ████████████████ check_line_databank 0.000s check_non_eq_param 0.042s ███ fetch_energy_5 0.015s █ calc_weight_trans 0.008s reinitialize 0.002s copy_database 0.000s memory_usage_warning 0.002s reset_population 0.000s calc_noneq_population 0.041s ███ part_function 0.035s ██ population 0.006s scaled_non_eq_linestrength 0.005s map_part_func 0.001s corrected_population_se 0.003s calc_emission_integral 0.006s applied_linestrength_cutoff 0.002s calc_lineshift 0.001s calc_hwhm 0.007s generate_wavenumber_arrays 0.001s calc_line_broadening 0.074s ██████ precompute_LDM_lineshapes 0.012s LDM_Initialized_vectors 0.000s LDM_closest_matching_line 0.001s LDM_Distribute_lines 0.001s LDM_convolve 0.060s █████ others 0.001s calc_other_spectral_quan 0.003s generate_spectrum_obj 0.000s others -0.016s Other Parameters ---------------- precision: int, optional total number of blocks. Default 16. See Also -------- :py:meth:`~radis.spectrum.spectrum.Spectrum.generate_perf_profile` """ from radis.spectrum.utils import print_perf_profile profiler = getattr(self, "profiler", None) if profiler is None: warn("No profiler attached to this Spectrum instance.") return None total_time = profiler["spectrum_calculation"]["value"] return print_perf_profile( profiler, total_time, number_format=number_format, precision=precision, first_line=self.get_name() + " profiler :", ) def generate_perf_profile(self): r"""Generate a visual/interactive performance profile diagram using ``tuna`` .. note:: requires a profiler attached to ``Spectrum.profiler`` .. warning:: deprecated in favor of :py:meth:`~radis.spectrum.spectrum.Spectrum.print_perf_profile` Examples -------- :: s = calc_spectrum(...) s.generate_perf_profile() See typical output in https://github.com/radis/radis/pull/325 .. image:: https://user-images.githubusercontent.com/16088743/128018032-6049be72-1881-46ac-9d7c-1ed89f9c4f42.png :alt: https://user-images.githubusercontent.com/16088743/128018032-6049be72-1881-46ac-9d7c-1ed89f9c4f42.png :target: https://user-images.githubusercontent.com/16088743/128018032-6049be72-1881-46ac-9d7c-1ed89f9c4f42.png .. note:: You can also profile with `tuna` directly:: python -m cProfile -o program.prof your_radis_script.py tuna your_radis_script.py See Also -------- :py:meth:`~radis.spectrum.spectrum.Spectrum.print_perf_profile` """ from radis.spectrum.utils import generate_perf_profile profiler_all = getattr(self, "profiler", None) if profiler_all is None: warn("No profiler attached to this Spectrum instance.") return None profiler = profiler_all.get("spectrum_calculation", {}).copy() # Add total calculation time: profiler.update({"value": self.conditions.get("calculation_time")}) return generate_perf_profile(profiler) # %% Define Spectrum Algebra # +, -, *, ^ operators # Notes: # s + s = 2*s # valid only if optically thin # Other possibility: # use s1>s2 for SerialSlabs # and s1//s2 for MergeSlabs # Or: # use s1*s2 for Serial Slabs # and s1+s2 for MergeSlabs # For the moment the first version is implemented # Plus def __add__(self, other): r"""Override '+' behavior Add is defined as : - for numeric values: add a baseline (returns a copy) - for 2 Spectra: not defined (not physical) """ if isinstance(other, float) or isinstance(other, int): from radis.spectrum.operations import add_constant return add_constant(self, other, inplace=False) elif isinstance(other, np.ndarray): from radis.spectrum.operations import add_array return add_array(self, other, inplace=False) elif isinstance(other, Spectrum): from radis.spectrum.operations import add_spectra return add_spectra(self, other) else: # warn("You shouldn't use the '+'. See '//' or '>' for more details", Warning) raise NotImplementedError( f"+ not implemented for a Spectrum and a {type(other)} object" ) def __radd__(self, other): r"""Right side addition.""" return self.__add__(other) def __iadd__(self, other): r"""Override '+=' behavior Add is defined as : - for numeric values: add a baseline (inplace) - for 2 Spectra: not defined (not physical) """ if isinstance(other, float) or isinstance(other, int): from radis.spectrum.operations import add_constant return add_constant(self, other, inplace=True) elif isinstance(other, np.ndarray): from radis.spectrum.operations import add_array return add_array(self, other, inplace=True) else: warn("You shouldn't use the '+'. See '//' or '>' for more details", Warning) raise NotImplementedError( f"+ not implemented for a Spectrum and a {type(other)} object" ) # Minus def __sub__(self, other): r"""Override '-' behavior Add is defined as : - for numeric values: substract a baseline (returns a copy) - for 2 Spectra: defined only for baseline substraction """ if isinstance(other, float) or isinstance(other, int): from radis.spectrum.operations import add_constant return add_constant(self, -other, inplace=False) elif isinstance(other, np.ndarray): from radis.spectrum.operations import add_array return add_array(self, -other, inplace=False) elif isinstance(other, Spectrum): from radis.spectrum.operations import substract_spectra return substract_spectra(self, other) else: raise NotImplementedError( f"- not implemented for a Spectrum and a {type(other)} object" ) def __rsub__(self, other): r"""Right side substraction.""" raise NotImplementedError( "right substraction (-) not implemented for Spectrum objects" ) def __isub__(self, other): r"""Override '-=' behavior Add is defined as : - for numeric values: substract a baseline (inplace) - for 2 Spectra: defined only for baseline substraction """ if isinstance(other, float) or isinstance(other, int): from radis.spectrum.operations import add_constant return add_constant(self, -other, inplace=True) elif isinstance(other, np.ndarray): from radis.spectrum.operations import add_array return add_array(self, -other, inplace=True) elif isinstance(other, Spectrum): from radis.spectrum.operations import substract_spectra return substract_spectra(self, other) else: raise NotImplementedError( f"-= not implemented for a Spectrum and a {type(other)} object" ) # Times def __mul__(self, other): r"""Override '*' behavior Multiply is defined as : - for numeric values: multiply (equivalent to optically thin scaling) (only if in front, i.e: 2*s works but s*2 is not implemented) (returns a copy) - for 2 Spectra: not defined """ import astropy.units as u if ( isinstance(other, float) or isinstance(other, int) or isinstance(other, u.quantity.Quantity) ): from radis.spectrum.operations import multiply return multiply(self, other, inplace=False) elif isinstance(other, Spectrum): raise NotImplementedError( "* not implemented for 2 Spectrum objects. Use > to combine them along the line of sight, as in SerialSlabs" ) elif isinstance(other, np.array): raise NotImplementedError( f"* not implemented for a Spectrum and a {type(other)} object. Use a dimensioned (with units) array" ) else: raise NotImplementedError( f"* not implemented for a Spectrum and a {type(other)} object" ) def __rmul__(self, other): r"""Right side multiplication.""" import astropy.units as u if ( isinstance(other, float) or isinstance(other, int) or isinstance(other, u.quantity.Quantity) ): from radis.spectrum.operations import multiply return multiply(self, other, inplace=False) elif isinstance(other, Spectrum): raise NotImplementedError( "* not implemented for 2 Spectrum objects. Use > to combine them along the line of sight, as in SerialSlabs" ) elif isinstance(other, np.array): raise NotImplementedError( f"right side * not implemented for a Spectrum and a {type(other)} object. Use a dimensioned (with units) array" ) else: raise NotImplementedError( f"right side * not implemented for a Spectrum and a {type(other)} object" ) def __imul__(self, other): r"""Override '*=' behavior Multiply is defined as : - for numeric values: multiply (equivalent to optically thin scaling) (only if in front, i.e: s *= 2) (modifies inplace) - for 2 Spectra: not defined """ import astropy.units as u if ( isinstance(other, float) or isinstance(other, int) or isinstance(other, u.quantity.Quantity) ): from radis.spectrum.operations import multiply return multiply(self, other, inplace=True) elif isinstance(other, Spectrum): raise NotImplementedError("* not implemented for 2 Spectrum objects. Use >") elif isinstance(other, np.array): raise NotImplementedError( f"*= not implemented for a Spectrum and a {type(other)} object. Use a dimensioned (with units) array" ) else: raise NotImplementedError( f"*= not implemented for a Spectrum and a {type(other)} object" ) # Divide def __truediv__(self, other): r"""Override '/' behavior Divide is defined as : - for numeric values: divide algebraically (equivalent to optically thin scaling) """ import astropy.units as u if isinstance(other, float) or isinstance(other, int): from radis.spectrum.operations import multiply return multiply(self, 1 / other, inplace=False) elif isinstance(other, u.quantity.Quantity): from radis.spectrum.operations import multiply return multiply(self, 1 / other.value, unit=1 / other.unit, inplace=False) elif isinstance(other, np.array): raise NotImplementedError( f"/ not implemented for a Spectrum and a {type(other)} object. Use a dimensioned (with units) array" ) else: raise NotImplementedError( f"/ not implemented for a Spectrum and a {type(other)} object" ) def __rtruediv__(self, other): r"""Right side division.""" raise NotImplementedError( f"right side / not implemented for a Spectrum and a {type(other)} object" ) def __itruediv__(self, other): r"""Override '/=' behavior Divide is defined as : - for numeric values: divide quantities algebraically (equivalent to optically thin scaling) """ import astropy.units as u if isinstance(other, float) or isinstance(other, int): from radis.spectrum.operations import multiply return multiply(self, 1 / other, inplace=True) elif isinstance(other, u.quantity.Quantity): from radis.spectrum.operations import multiply return multiply(self, 1 / other.value, unit=1 / other.unit, inplace=True) else: raise NotImplementedError( f"/= not implemented for a Spectrum and a {type(other)} object" ) # Line of sight operations def __bool__(self): r"""This prevents behaviors such as:: s1 > s2 > s3 which are actually interpreted by Python as "s1 > s2 and s2 > s3" and would return a wrong result (probably s2 > s3 ? ) """ raise ArithmeticError( "A Spectrum cannot be evaluated as a boolean. " + "You may have tried using syntax such as `s1>s2>s3` " + "which Python interprets as `s1>s2 and s2>s3`. " + "Use `(s1>s2)>s3)` or SerialSlabs(s1, s2, s3) instead." ) def __gt__(self, other): r"""Overloads '>' behavior no comparison: here we use > to define a ``Line of sight``. Examples -------- s_plasma is seen through s_room:: s = s_plasma > s_room Equivalent to:: s = SerialSlabs(s_plasma, s_room) """ if isinstance(other, Spectrum): from radis.los.slabs import SerialSlabs return SerialSlabs(self, other) else: raise NotImplementedError( f"> not implemented for a Spectrum and a {type(other)} object" ) # def __rgt__(self, other): # ''' Right side > ''' # # if isinstance(other, Spectrum): # from radis.los.slabs import SerialSlabs # return SerialSlabs(other, self) # else: # raise NotImplementedError('right side > not implemented for a Spectrum and a {0} object'.format( # type(other))) def __floordiv__(self, other): r"""Overloads '//' behavior not a division here: we use it to say that Slabs are ``in parallel``, i.e., as if their respective mole fractions were added in the same physical space. Won't work if they are not defined on the same waverange, but it's okay: let the user use MergeSlabs manually with the appropriate options Examples -------- s_co2 added with s_co: s = s_co2 // s_co Equivalent to:: s = MergeSlabs(s_co2, s_co) """ if isinstance(other, Spectrum): from radis.los.slabs import MergeSlabs return MergeSlabs(self, other) else: raise NotImplementedError( f"// not implemented for a Spectrum and a {type(other)} object" ) # the following is so that json_tricks.dumps and .loads can be used directly, # ie.: # >>> s == loads(s.dumps()) # Does not work yet. # see https://github.com/mverleg/pyjson_tricks#class-instances # def __json_encode__(self): # ''' Called by json_tricks.dumps ''' # from radis.tools.database import _format_to_jsondict # return _format_to_jsondict(self, discard=[], compress=[], verbose=True) # # def __json_decode__(self, **attrs): # from radis.tools.database import _json_to_spec # print(attrs) # raise # return _json_to_spec(attrs) def __len__(self): r"""Length of a Spectrum object = length of the wavespace if unique, else raises an error""" return len(self._get_wavespace(copy=False)) def recalc_gpu( self, var="abscoeff", wunit="default", Iunit="default", pressure=None, Tgas=None, mole_fraction=None, path_length=None, slit_function=None, ): """Recalculate the spectrum based on new input parameters. Can only be called for spectrum objects produced by :py:meth:`~radis.lbl.factory.SpectrumFctory.eq_spectrum_gpu`. This method is used internally by :py:meth:`~radis.lbl.factory.SpectrumFctory.eq_spectrum_gpu_interactive`. Parameters may be passed as arguments, or updated directly in :py:attr:`~radis.spectrum.spectrum.Spectrum.conditions`, after which spectrum.recalc_gpu() may be called without passing arguments. Parameters ---------- var: variable (`absorbance`, `transmittance`, `transmittance_noslit`, `xsection`, etc.) For full list see :py:meth:`~radis.spectrum.spectrum.Spectrum.get_vars()`. If ``None``, plot the first thing in the Spectrum. Default ``None``. wunit: ``'nm'``, ``'cm-1'``, ``'nm_vac'`` or ``None`` wavelength in air (``'nm'``), wavenumber (``'cm-1'``), or wavelength in vacuum (``'nm_vac'``). If ``None``, ``'wavespace'`` must be defined in ``conditions``. Iunit: unit for variable if `default`, default unit for quantity `var` is used. for radiance, one can use per wavelength (~ `W/m2/sr/nm`) or per wavenumber (~ `W/m2/sr/cm-1`) units pressure: float, optional Pressure in [bar] Tgas: float, optional Equilibrium temperature in [K] mole_fraction: float, optional Mole fraction path_length: float, optional Absroption length in [cm] slit_function: float, optional Slit FWHM in [cm-1] Returns ------- new_y: numpy.ndarray[np.fltoat32] The recalculated array specified by the ``var`` keyword. Examples -------- :: s = sf.eq_spectrum_gpu( Tgas=1100.0, # K pressure=1, # bar mole_fraction=0.8, path_length=0.2, # cm ) I = [] for T in [300.0, 1000.0, 2000.0]: I.append(s.recalc_gpu('radiance', Tgas=T) """ profiler = Profiler(False) profiler.start("spectrum_calculation", 1) profiler.start("setting_params", 2) try: self.conditions["gpu_backend"] if self.gpu_app is None: raise KeyError except KeyError: warn( "GPU not initialized, spectrum.recalc_gpu() can only be called on spectrum objects produced by sf.eq_spectrum_gpu()!", GPUInitWarning, ) # Update conditions: if Tgas is not None: self.conditions["Tgas"] = Tgas self.conditions["Tvib"] = self.conditions["Tgas"] self.conditions["Trot"] = self.conditions["Tgas"] if pressure is not None: self.conditions["pressure"] = pressure if mole_fraction is not None: self.conditions["mole_fraction"] = mole_fraction if path_length is not None: self.conditions["path_length"] = path_length if slit_function is not None: self.conditions["slit_function"] = slit_function profiler.stop("setting_params", "Input parameters set") profiler.start("gpu_iterate", 2) abscoeff, times = self.gpu_app.iterate( self.conditions["pressure"], self.conditions["Tgas"], self.conditions["mole_fraction"], verbose=0, # TODO: GPU apply_slit not supported yet # l=s.conditions["path_length"], # slit_FWHM=s.conditions["slit_function"], ) profiler.stop("gpu_iterate", "GPU calculations done") profiler.start("spectral_quantities", 2) iter_N_L, iter_N_G = self.gpu_app.get_griddims() self.conditions["NwL"] = iter_N_L self.conditions["NwG"] = iter_N_G # TODO : refactor this function and the update() mechanism. Ensure conditions are correct. for k in list(self._q.keys()): # reset all quantities if k in ["wavespace", "wavelength", "wavenumber"]: pass elif k == "abscoeff": self._q["abscoeff"] = abscoeff else: del self._q[k] _, new_y = self.get( var, copy=False, # copy = False saves some time & memory, it's a pointer/reference to the real data, which is fine here as data is just plotted wunit=wunit, Iunit=Iunit, ) profiler.stop("spectral_quantities", "Calculated spectral qunatities") profiler.stop("spectrum_calculation", "Spectrum calculated") self.conditions["calculation_time"] = profiler.final["spectrum_calculation"][ "value" ] return new_y def exit_gpu(self, event=None): if self.gpu_app is not None: self.gpu_app.__del__() self.gpu_app = None # #Not working: # def __del__(self): # self.exit_gpu() # %% Private functions # to cut def _cut_slices(w_spec_nm, w_slit_nm, slit_dispersion, slit_dispersion_threshold=0.01): """ Used to cut a waverange into slices where dispersion does not vary too much. Works whether the waverange is reversed or not. Parameters ---------- w_spec_nm : numpy array The waverange (reversed or not) to cut into slices. w_slit_nm : numpy array The waverange of the slit function. slit_dispersion : function The dispersion function of the spectrometer, must be monotone and accept arrays. slit_dispersion_threshold : float Threshold between 0 and 1, used to cut the waverange when the slit dispersion variations exceed the threshold. Returns ------- slices : list List of boolean arrays, corresponding to each slice. """ from scipy.optimize import root_scalar if slit_dispersion(w_spec_nm[0]) > slit_dispersion(w_spec_nm[-1]): w_start = w_spec_nm[0] w_end = w_spec_nm[-1] else: w_start = w_spec_nm[-1] w_end = w_spec_nm[0] # Find the wavelengths to delimit each slice w_list_slices = [w_start] slit_disp_target = slit_dispersion(w_start) * (1 - slit_dispersion_threshold) while slit_disp_target > slit_dispersion(w_end): res = root_scalar( lambda w: slit_dispersion(w) - slit_disp_target, bracket=[w_start, w_end] ) assert res.converged, "Did not converged, something went wrong... +\ Are you sure the slit dispersion function is monotone?" w_next = res.root w_list_slices.append(w_next) slit_disp_target = slit_dispersion(w_next) * (1 - slit_dispersion_threshold) w_start = w_next w_list_slices.append(w_end) if w_list_slices[0] > w_list_slices[-1]: w_list_slices[::-1] w_slices = np.array(w_list_slices) # Extend the slices accordingly to the slit function # to remove boundary effect caused by the convolution range_nm = abs(w_slit_nm[-1] - w_slit_nm[0]) / 2 slit_disp_0 = slit_dispersion(w_slit_nm[len(w_slit_nm) // 2]) w_mins, w_maxs = w_slices[:-1].copy(), w_slices[1:].copy() w_mins -= range_nm / slit_disp_0 * slit_dispersion(w_mins) w_maxs += range_nm / slit_disp_0 * slit_dispersion(w_maxs) slices = [] for w_min, w_max in zip(w_mins, w_maxs): slice_w = np.logical_and(w_min < w_spec_nm, w_spec_nm < w_max) slices.append(slice_w) if w_spec_nm[-1] < w_spec_nm[0]: slices = slices[::-1] return slices # %% ====================================================================== # Test class function # ------------------- # XXX ===================================================================== # %% Test functions if __name__ == "__main__": from radis.test.spectrum.test_spectrum import _run_testcases print("Test spectrum: ", _run_testcases(debug=False))