Computing CO cross section using HITRAN (opacity calculator = LPF)

This tutorial demonstrates how to compute the opacity of CO using HITRAN step by step.

from exojax.spec.lpf import auto_xsection
from exojax.spec.hitran import SijT, doppler_sigma, gamma_hitran, gamma_natural
from exojax.spec import api
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('bmh')

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 HITRAN database
nus = np.linspace(1000.0, 10000.0, 900000, dtype=np.float64)  #cm-1
isotope = 1
mdbCO = api.MdbHitran('CO', nus, isotope=isotope, gpu_transfer=True)

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.

Mmol=28.0 # molecular weight
Tfix=1000.0 # we assume T=1000K
Pfix=1.e-3 # we compute P=1.e-3 bar
Ppart=Pfix #partial pressure of CO. here we assume a 100% CO atmosphere.

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

qt=mdbCO.qr_interp(isotope,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=SijT(Tfix,mdbCO.logsij0,mdbCO.nu_lines,mdbCO.elower,qt)

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

\(\gamma_L = \gamma^p_L + \gamma^n_L\)

where the pressure broadning

\(\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)

and the natural broadening

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

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

Thermal broadening

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

# thermal doppler sigma
sigmaD=doppler_sigma(mdbCO.nu_lines,Tfix,Mmol)

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=mdbCO.nu_lines

#Use below if you wanna include a slight pressure line shift
#nu0=mdbCO.nu_lines+mdbCO.delta_air*Pfix

Although it depends on your GPU, you might need to devide the computation into multiple loops because of the limitation of the GPU memory. Here we assume 30MB for GPU memory (not exactly, memory size for numatrix).

xsv=auto_xsection(nus,nu0,sigmaD,gammaL,Sij,memory_size=30) #use 30MB GPU MEMORY for numax

100%|██████████| 98/98 [00:00<00:00, 179.80it/s]

Plot it!

fig=plt.figure(figsize=(10,3))
ax=fig.add_subplot(111)
plt.plot(nus,xsv,lw=0.1,label="exojax")
plt.yscale("log")
plt.xlabel("wavenumber ($cm^{-1}$)")
plt.ylabel("cross section ($cm^{2}$)")
plt.legend(loc="upper left")
plt.savefig("co_hitran.pdf", bbox_inches="tight", pad_inches=0.0)
plt.show()
../_images/opacity_20_0.png 
fig=plt.figure(figsize=(10,3))
ax=fig.add_subplot(111)
plt.plot(1.e8/nus,xsv,lw=1,label="exojax")
plt.yscale("log")
plt.xlabel("wavelength ($\AA$)")
plt.ylabel("cross section ($cm^{2}$)")
plt.xlim(22985.,23025)
plt.legend(loc="upper left")
plt.savefig("co_hitran.pdf", bbox_inches="tight", pad_inches=0.0)
plt.show()
../_images/opacity_21_0.png