Cross Section using DITΒΆ

Update: July 11/2021, Hajime Kawahara

We demonstarte the Discrete Integral Transform (DIT) method proposed by D.C.M van den Bekeroma and E.Pannier. DIT takes advantage especially for the case that the number of the molecular line is large (typically > 10,000). We here compare the results by DIT with the direct computation (LPF).

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import jax.numpy as jnp
plt.style.use('bmh')

from exojax.spec.hitran import SijT, doppler_sigma, gamma_hitran, gamma_natural
from exojax.spec import moldb

# Setting wavenumber bins and loading HITRAN database
nus=np.linspace(1900.0,2300.0,350000,dtype=np.float64)
mdbCO=moldb.MdbHit('05_hit12.par',nus)

# set T, P and partition function
Mmol=28.01 # 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.
qt=mdbCO.Qr_layer_HAPI([Tfix])[0]

# compute Sij, gamma_L, sigmaD
Sij=SijT(Tfix,mdbCO.logsij0,mdbCO.nu_lines,mdbCO.elower,qt)
gammaL = gamma_hitran(Pfix,Tfix, Ppart, mdbCO.n_air, \
                      mdbCO.gamma_air, mdbCO.gamma_self) \
+ gamma_natural(mdbCO.A)
sigmaD=doppler_sigma(mdbCO.nu_lines,Tfix,Mmol)

DIT uses a grid of sigmaD, gammaL, and wavenumber. set_ditgrid makes a 1D grid for sigmaD and gamma.

from exojax.spec.dit import set_ditgrid
sigmaD_grid=set_ditgrid(sigmaD)
gammaL_grid=set_ditgrid(gammaL)

# we can change the resolution using res option
#sigmaD_grid=set_ditgrid(sigmaD,res=0.1)
#gammaL_grid=set_ditgrid(gammaL,res=0.1)

#show the grids
plt.plot(sigmaD,gammaL,".")
for i in sigmaD_grid:
    plt.axvline(i,lw=1,alpha=0.5,color="C1")
for i in gammaL_grid:
    plt.axhline(i,lw=1,alpha=0.5,color="C1")
../_images/output_6_0.png

We need to precompute the contribution for wavenumber and pmarray. These can be computed using init_dit:

from exojax.spec import initspec
cnu,indexnu,pmarray=initspec.init_dit(mdbCO.nu_lines,nus)

Then, you can compute the cross section as

from exojax.spec.dit import xsvector
xs=xsvector(cnu,indexnu,pmarray,sigmaD,gammaL,Sij,nus,sigmaD_grid,gammaL_grid)