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

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

from exojax.spec.lpf import auto_xsection
from exojax.spec import SijT, doppler_sigma,  gamma_natural
from exojax.spec.exomol import gamma_exomol
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 exomole files.

# Setting wavenumber bins and loading HITRAN database
nus=np.linspace(1000.0,10000.0,900000,dtype=np.float64) #cm-1
emf='CO/12C-16O/Li2015'
mdbCO=api.MdbExomol(emf,nus,gpu_transfer=True)

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

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.010446441149536 # molecular weight
Tfix=1000.0 # we assume T=1000K
Pfix=1.e-3 # we compute P=1.e-3 bar

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 the interpolation of partition function

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

$: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 = gamma_exomol(Pfix,Tfix,mdbCO.n_Texp,mdbCO.alpha_ref)\
+ gamma_natural(mdbCO.A)

gamma_exomol(Pfix,Tfix,mdbCO.n_Texp,mdbCO.alpha_ref)

DeviceArray([3.1183732e-05, 3.8084123e-05, 3.8084123e-05, ...,
             3.8084123e-05, 3.1183732e-05, 3.1455678e-05], dtype=float32)

fig=plt.figure()
fig.add_subplot(211)
plt.plot(mdbCO.jlower,mdbCO.n_Texp,".")
fig.add_subplot(212)
plt.plot(mdbCO.jlower,mdbCO.alpha_ref,".")

[<matplotlib.lines.Line2D at 0x7efbb8248ca0>]
../_images/opacity_exomol_14_1.png

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…

#line center
nu0=mdbCO.nu_lines

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%|██████████| 8257/8257 [05:14<00:00, 26.22it/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_exomol.pdf", bbox_inches="tight", pad_inches=0.0)
plt.show()
../_images/opacity_exomol_22_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_exomol.pdf", bbox_inches="tight", pad_inches=0.0)
plt.show()
../_images/opacity_exomol_23_0.png

Important Note

Use float64 for wavenumber bin and line center.

Below, we see the difference of opacity between float64 case and float 32.

xsv_32=auto_xsection(np.float32(nus),np.float32(nu0),sigmaD,gammaL,Sij,memory_size=30)

100%|██████████| 8257/8257 [02:51<00:00, 48.08it/s]
/home/kawahara/exojax/src/exojax/spec/lpf.py:363: UserWarning: The wavenumber grid is not np.float64 but float32
  warnings.warn('The wavenumber grid is not np.float64 but '+str(nu.dtype),UserWarning)
/home/kawahara/exojax/src/exojax/spec/lpf.py:365: UserWarning: The line centers (nu_lines) are not np.float64 but float32
  warnings.warn('The line centers (nu_lines) are not np.float64 but '+str(nu.dtype),UserWarning)

fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(211)
plt.plot(1.e8/nus,xsv,".",lw=1,label="64",markersize=1)
plt.plot(1.e8/nus,xsv_32,".",lw=1,label="32",markersize=1)
plt.xlim(22985.,23025)
plt.yscale("log")
plt.ylabel("xsv $cm^{2}$")
ax=fig.add_subplot(212)
plt.plot(1.e8/nus,(xsv_32-xsv)/xsv,lw=1,label="difference")
plt.xlabel("wavelength ($\AA$)")
plt.ylabel("Difference")
plt.xlim(22985.,23025)
plt.legend(loc="upper left")
plt.show()
../_images/opacity_exomol_27_0.png

We found ~ 10 % error when using float32 as an wavenumber and line center