Comparing HITEMP and ExoMol

from exojax.spec.lpf import auto_xsection
from exojax.spec.hitran import line_strength, doppler_sigma, gamma_hitran, gamma_natural
from exojax.spec.exomol import gamma_exomol
from exojax.spec import api
import numpy as np
import matplotlib.pyplot as plt

2023-03-15 08:37:58.525388: E external/org_tensorflow/tensorflow/stream_executor/cuda/cuda_driver.cc:265] failed call to cuInit: CUDA_ERROR_UNKNOWN: unknown error
WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)

First of all, set a wavenumber bin in the unit of wavenumber (cm-1). Here we set the wavenumber range as $1000 \le \nu \le 10000$ (1/cm) with the resolution of 0.01 (1/cm).

We call moldb instance with the path of par file. If the par file does not exist, moldb will try to download it from HITRAN website.

# Setting wavenumber bins and loading HITEMP database
wav=np.linspace(22930.0,23000.0,4000,dtype=np.float64) #AA
nus=1.e8/wav[::-1] #cm-1

mdbCO_HITEMP=api.MdbHitemp('CO',nus, isotope=1, gpu_transfer=True) # we use istope=1 for comparison

emf='CO/12C-16O/Li2015' #this is isotope=1 12C-16O
mdbCO_Li2015=api.MdbExomol(emf,nus, gpu_transfer = True)

HITRAN exact name= (12C)(16O)
Background atmosphere:  H2
Reading CO/12C-16O/Li2015/12C-16O__Li2015.trans.bz2
.broad is used.
Broadening code level= a0
default broadening parameters are used for  71  J lower states in  152  states

/home/kawahara/exojax/src/exojax/utils/molname.py:133: FutureWarning: e2s will be replaced to exact_molname_exomol_to_simple_molname.
  warnings.warn(

Define molecular weight of CO ($\sim 12+16=28$), temperature (K), and pressure (bar). Also, we here assume the 100 % CO atmosphere, i.e. the partial pressure = pressure.

from exojax.spec import molinfo
molecular_mass=molinfo.molmass("CO") # molecular weight
Tfix=1300.0 # we assume T=1300K
Pfix=0.99 # we compute P=1 bar=0.99+0.1
Ppart=0.01 #partial pressure of CO. here we assume a 1% CO atmosphere (very few).

partition function ratio $q(T)$ is defined by

$q(T) = Q(T)/Q(T_{ref})$; $T_{ref}$=296 K

Here, we use the partition function from HAPI

#mdbCO_HITEMP.ExomolQT(emf) #use Q(T) from Exomol/Li2015
qt_HITEMP=mdbCO_HITEMP.qr_interp(1,Tfix)
qt_Li2015=mdbCO_Li2015.qr_interp(Tfix)

Let us compute the line strength S(T) at temperature of Tfix.

$S (T;s_0,\nu_0,E_l,q(T)) = S_0 \frac{Q(T_{ref})}{Q(T)} \frac{e^{- h c E_l /k_B T}}{e^{- h c E_l /k_B T_{ref}}} \frac{1- e^{- h c \nu /k_B T}}{1-e^{- h c \nu /k_B T_{ref}}}= q_r(T)^{-1} e^{ s_0 - c_2 E_l (T^{-1} - T_{ref}^{-1})} \frac{1- e^{- c_2 \nu_0/ T}}{1-e^{- c_2 \nu_0/T_{ref}}}$

$s_0=\log_{e} S_0$ : logsij0

$\nu_0$: nu_lines

$E_l$ : elower

Why the input is $s_0 = \log_{e} S_0$ instead of $S_0$ in SijT? This is because the direct value of $S_0$ is quite small and we need to use float32 for jax.

Sij_HITEMP=line_strength(Tfix,mdbCO_HITEMP.logsij0,mdbCO_HITEMP.nu_lines,\
         mdbCO_HITEMP.elower,qt_HITEMP)
Sij_Li2015=line_strength(Tfix,mdbCO_Li2015.logsij0,mdbCO_Li2015.nu_lines,\
                mdbCO_Li2015.elower,qt_Li2015)

Then, compute the Lorentz gamma factor (pressure+natural broadening)

$\gamma_L = \gamma^p_L + \gamma^n_L$

where the pressure broadning (HITEMP)

$\gamma^p_L = (T/296K)^{-n_{air}} [ \alpha_{air} ( P - P_{part})/P_{atm} + \alpha_{self} P_{part}/P_{atm}]$

$P_{atm}$: 1 atm in the unit of bar (i.e. = 1.01325)

or

the pressure broadning (ExoMol)

$:raw-latex:gamma`^p_L = :raw-latex:alpha`{ref} ( T/T{ref} )^{-n_{texp}} ( P/P_{ref}), $

and the natural broadening

$\gamma^n_L = \frac{A}{4 \pi c}$

gammaL_HITEMP = gamma_hitran(Pfix,Tfix, Ppart, mdbCO_HITEMP.n_air, \
                      mdbCO_HITEMP.gamma_air, mdbCO_HITEMP.gamma_self) \
+ gamma_natural(mdbCO_HITEMP.A)

gammaL_Li2015 = gamma_exomol(Pfix,Tfix,mdbCO_Li2015.n_Texp,mdbCO_Li2015.alpha_ref)\
+ gamma_natural(mdbCO_Li2015.A)

Thermal broadening

$\sigma_D^{t} = \sqrt{\frac{k_B T}{M m_u}} \frac{\nu_0}{c}$

# thermal doppler sigma
sigmaD_HITEMP=doppler_sigma(mdbCO_HITEMP.nu_lines,Tfix,molecular_mass)
sigmaD_Li2015=doppler_sigma(mdbCO_Li2015.nu_lines,Tfix,molecular_mass)

Then, the line center…

In HITRAN database, a slight pressure shift can be included using $\delta_{air}$: $\nu_0(P) = \nu_0 + \delta_{air} P$. But this shift is quite a bit.

#line center
nu0_HITEMP=mdbCO_HITEMP.nu_lines
nu0_Li2015=mdbCO_Li2015.nu_lines

We use Direct LFP.

from exojax.spec.initspec import init_lpf
from exojax.spec.lpf import xsvector

numatrix_HITEMP = init_lpf(mdbCO_HITEMP.nu_lines, nus)
xsv_HITEMP=xsvector(numatrix_HITEMP, sigmaD_HITEMP, gammaL_HITEMP, Sij_HITEMP)

numatrix_Li2015 = init_lpf(mdbCO_Li2015.nu_lines, nus)
xsv_Li2015=xsvector(numatrix_Li2015, sigmaD_Li2015, gammaL_Li2015, Sij_Li2015)

fig=plt.figure(figsize=(10,4))
ax=fig.add_subplot(111)
plt.plot(wav[::-1],xsv_HITEMP,lw=2,label="HITEMP2019")
plt.plot(wav[::-1],xsv_Li2015,lw=2,ls="dashed",label="Exomol w/ .broad")
plt.xlim(22970,22976)
plt.xlabel("wavelength ($\AA$)",fontsize=14)
plt.ylabel("cross section ($cm^{2}$)",fontsize=14)
plt.legend(loc="upper left",fontsize=14)
plt.tick_params(labelsize=12)
plt.savefig("co_comparison.pdf", bbox_inches="tight", pad_inches=0.0)
plt.savefig("co_comparison.png", bbox_inches="tight", pad_inches=0.0)
plt.title("T=1300K,P=1bar")
plt.show()