Note

You can download below the full example code. and run it with 🔬 Radis-Lab,

Compare OH(A-X) electronic spectra: RADIS vs Specair

Auto-download the MoLLIST-OH line database from ExoMol and compute the OH A-X electronic band (300-340 nm) under equilibrium and non-equilibrium conditions. Compare with Specair reference data.

Note

Electronic spectra computation in RADIS is a work in progress. RADIS currently shows ~10-20% discrepancy in integrated radiance compared to Specair, likely due to database differences. See radis.test.validation.test_validation_vs_specair_OH_AX for the full validation test suite.

from os.path import join

import numpy as np
import pandas as pd

from radis import Spectrum, SpectrumFactory, plot_diff
from radis.test.utils import getValidationCase


def _load_specair_csv(csv_path):
    """Load a Specair reference CSV and return a Spectrum object."""
    data = pd.read_csv(csv_path)
    header = ",".join(data.columns)
    w = data.iloc[:, 0].values
    I = data.iloc[:, 1].values

    mask = ~np.isnan(I)
    w, I = w[mask], I[mask]
    idx = np.argsort(w)
    w, I = w[idx], I[idx]

    if "W_cm2_nm_sr" in header:
        I = I * w**2 / 1e4

    return Spectrum.from_array(
        w, I, "radiance", wunit="nm", unit="mW/cm2/sr/cm-1", name="Specair"
    )


def _integrate_radiance(s, wunit="nm", Iunit="mW/cm2/sr/cm-1"):
    """Get integrated radiance, filtering NaN and sorting."""
    w, I = s.get("radiance", wunit=wunit, Iunit=Iunit)
    mask = ~np.isnan(I)
    w, I = w[mask], I[mask]
    idx = np.argsort(w)
    return np.trapezoid(I[idx], w[idx])

sf = SpectrumFactory(
    wavelength_min=300,
    wavelength_max=340,
    molecule="OH",
    isotope="1",
    pressure=1,
    path_length=1,
    mole_fraction=1e-3,
    wstep=0.001,
    verbose=False,
    self_absorption=True,
)
sf.warnings["ElectronicSpectraWarning"] = "ignore"
sf.fetch_databank("exomol", "MoLLIST-OH", load_energies=True)

/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/latest/radis/io/exomol.py:263: AccuracyWarning: ExoMol cache file /home/docs/.radisdb/exomol/OH/16O-1H/MoLLIST-OH/16O-1H__MoLLIST-OH.trans.h5 is missing columns ['gupper', 'glower']. This may be from an older version of RADIS. Regenerate the cache with `cache='regen'` to fix this.
  warnings.warn(
/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/latest/radis/api/exomolapi.py:1737: UserWarning: Could not load `16O-1H__self.broad`. The default broadening parameters are used.

  warnings.warn(
The default broadening parameters are used.
7.77s - Loaded database

s_specair_eq = _load_specair_csv(
    getValidationCase(
        join("test_validation_vs_specair_OH_AX_data", "radiance_spectrum.csv")
    )
)

s_radis_eq = sf.non_eq_spectrum(
    Ttrans=400, Trot=400, Tvib=400, Telec=10000, name="RADIS"
)
s_radis_eq.apply_slit(0.1)

fig, [ax0, ax1] = plot_diff(
    s_radis_eq,
    s_specair_eq,
    "radiance",
    wunit="nm",
    Iunit="mW/cm2/sr/cm-1",
    method="diff",
)
ratio = _integrate_radiance(s_radis_eq) / _integrate_radiance(s_specair_eq)
ax0.set_title(f"Equilibrium 400 K: RADIS vs Specair (ratio={ratio:.2f})")
Equilibrium 400 K: RADIS vs Specair (ratio=1.23)
/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/latest/radis/misc/warning.py:443: MissingPressureShiftWarning: Pressure-shift coefficient not given in database: assumed 0 pressure shift
  warnings.warn(WarningType(message))
/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/latest/radis/misc/curve.py:267: UserWarning: First spectrum has more resolution than 2nd. Reverse your spectra in interpolation/comparison for a better accuracy
  warnings.warn(

Text(0.5, 1.0, 'Equilibrium 400 K: RADIS vs Specair (ratio=1.23)')

s_eq = sf.eq_spectrum(Tgas=400, name="RADIS-eq")
s_noneq = sf.non_eq_spectrum(Tvib=400, Trot=400, Ttrans=400, name="RADIS-noneq")

fig2, [ax2, ax3] = plot_diff(s_eq, s_noneq, "radiance_noslit")
ratio2 = s_eq.get_integral("radiance_noslit") / s_noneq.get_integral("radiance_noslit")
ax2.set_title(f"Self-consistency: eq vs noneq at 400 K (ratio={ratio2:.4f})")
Self-consistency: eq vs noneq at 400 K (ratio=0.9978)
/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/latest/radis/misc/warning.py:443: MissingPressureShiftWarning: Pressure-shift coefficient not given in database: assumed 0 pressure shift
  warnings.warn(WarningType(message))
/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/latest/radis/misc/warning.py:443: MissingPressureShiftWarning: Pressure-shift coefficient not given in database: assumed 0 pressure shift
  warnings.warn(WarningType(message))

Text(0.5, 1.0915471193184172, 'Self-consistency: eq vs noneq at 400 K (ratio=0.9978)')

s_specair_c2 = _load_specair_csv(
    getValidationCase(
        join(
            "test_validation_vs_specair_OH_AX_data",
            "radiance_spectrum_Tgas300_Trot500_Tvib2000_Telec10000.csv",
        )
    )
)

s_radis_c2 = sf.non_eq_spectrum(
    Ttrans=300, Trot=500, Tvib=2000, Telec=10000, name="RADIS"
)
s_radis_c2.apply_slit(0.1)

fig3, [ax4, ax5] = plot_diff(
    s_radis_c2,
    s_specair_c2,
    "radiance",
    wunit="nm",
    Iunit="mW/cm2/sr/cm-1",
    method="diff",
)
ratio3 = _integrate_radiance(s_radis_c2) / _integrate_radiance(s_specair_c2)
ax4.set_title(f"Non-eq Trot=500 Tvib=2000 Telec=10000 (ratio={ratio3:.2f})")
Non-eq Trot=500 Tvib=2000 Telec=10000 (ratio=1.22)
/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/latest/radis/misc/warning.py:443: MissingPressureShiftWarning: Pressure-shift coefficient not given in database: assumed 0 pressure shift
  warnings.warn(WarningType(message))
/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/latest/radis/misc/curve.py:267: UserWarning: First spectrum has more resolution than 2nd. Reverse your spectra in interpolation/comparison for a better accuracy
  warnings.warn(

Text(0.5, 1.0, 'Non-eq Trot=500 Tvib=2000 Telec=10000 (ratio=1.22)')

s_specair_c3 = _load_specair_csv(
    getValidationCase(
        join(
            "test_validation_vs_specair_OH_AX_data",
            "radiance_spectrum_Tgas300_Trot2000_Tvib500_Telec15000.csv",
        )
    )
)

s_radis_c3 = sf.non_eq_spectrum(
    Ttrans=300, Trot=2000, Tvib=500, Telec=15000, name="RADIS"
)
s_radis_c3.apply_slit(0.1)

fig4, [ax6, ax7] = plot_diff(
    s_radis_c3,
    s_specair_c3,
    "radiance",
    wunit="nm",
    Iunit="mW/cm2/sr/cm-1",
    method="diff",
)
ratio4 = _integrate_radiance(s_radis_c3) / _integrate_radiance(s_specair_c3)
ax6.set_title(f"Non-eq Trot=2000 Tvib=500 Telec=15000 (ratio={ratio4:.2f})")
Non-eq Trot=2000 Tvib=500 Telec=15000 (ratio=1.12)
/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/latest/radis/misc/warning.py:443: MissingPressureShiftWarning: Pressure-shift coefficient not given in database: assumed 0 pressure shift
  warnings.warn(WarningType(message))
/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/latest/radis/misc/curve.py:267: UserWarning: First spectrum has more resolution than 2nd. Reverse your spectra in interpolation/comparison for a better accuracy
  warnings.warn(

Text(0.5, 1.0, 'Non-eq Trot=2000 Tvib=500 Telec=15000 (ratio=1.12)')

Total running time of the script: (0 minutes 44.545 seconds)