Forward Modeling of an Emission Spectrum using the PreMODIT Cross Section

Update November 23rd; Hajime Kawahara

Here, we try to compute a emission spectrum using PreMODIT. Note that PreMODIT is a beta release. Presumably, we will improve it and release better one in ExoJAX 2. We try to compute an emission spectrum in which many methane lines exist. This situation mocks a T-type brown dwarf.

from exojax.spec import rtransfer as rt

#ATMOSPHERE
NP=100
T0=800.0 #K
Parr, dParr, k=rt.pressure_layer(NP=NP)
Tarr = T0*(Parr)**0.1

A T-P profile we assume is here.

import matplotlib.pyplot as plt
plt.style.use('bmh')
plt.plot(Tarr,Parr)
plt.yscale("log")
plt.gca().invert_yaxis()
plt.show()

We set a wavenumber grid using wavenumber_grid. Specify xsmode=“modit” though it is not mandatory. MODIT uses FFT, so the (internal) wavenumber grid should be evenly spaced in log.

from exojax.utils.grids import wavenumber_grid
nus,wav,resolution=wavenumber_grid(16360,16560,10000,unit="AA",xsmode="premodit")

xsmode assumes ESLOG in wavenumber space: mode=premodit

Loading a molecular database of CH4 and CIA (H2-H2)… The point here is gpu_transfer = False. We do not need to send line information to the device.

from exojax.spec import api, contdb
mdbCH4=api.MdbExomol('.database/CH4/12C-1H4/YT10to10/',nus,gpu_transfer=False)
cdbH2H2=contdb.CdbCIA('.database/H2-H2_2011.cia',nus)

Background atmosphere:  H2
Reading .database/CH4/12C-1H4/YT10to10/12C-1H4__YT10to10__06000-06100.trans.bz2
Reading .database/CH4/12C-1H4/YT10to10/12C-1H4__YT10to10__06100-06200.trans.bz2
.broad is used.
Broadening code level= a1
default broadening parameters are used for  23  J lower states in  40  states
H2-H2

We have 80.5 million lines

len(mdbCH4.A)

80505310

In fact, this number is too large for MODIT.

from exojax.spec import molinfo
molmassCH4=molinfo.molmass("CH4")

Computing the relative partition function,

from jax import vmap
qt=vmap(mdbCH4.qr_interp)(Tarr)

PreMODIT precomputes a kind of the line shape density, called LBD (line basis density). init_premodit easily does that. It takes a while.

from exojax.spec.initspec import init_premodit

interval_contrast = 0.1
dit_grid_resolution = 0.1
Ttyp = 700.0

lbd, multi_index_uniqgrid, elower_grid, ngamma_ref_grid, n_Texp_grid, R, pmarray = init_premodit(
    mdbCH4.nu_lines,
    nus,
    mdbCH4.elower,
    mdbCH4.alpha_ref,
    mdbCH4.n_Texp,
    mdbCH4.Sij0,
    Ttyp,
    interval_contrast=interval_contrast,
    dit_grid_resolution=dit_grid_resolution,
    warning=False)

uniqidx: 100%|██████████| 2/2 [00:03<00:00,  1.84s/it]

Then, let’s compute the cross section matrix.

from exojax.spec import premodit

xsm = premodit.xsmatrix(Tarr, Parr, R, pmarray, lbd, nus, ngamma_ref_grid,
                   n_Texp_grid, multi_index_uniqgrid, elower_grid, molmassCH4, qt)

import numpy as np
fig=plt.figure(figsize=(20,4))
ax=fig.add_subplot(111)
c=plt.imshow(np.log10(xsm),cmap="bone_r",vmin=-23,vmax=-19)
plt.colorbar(c,shrink=0.8)
plt.text(50,30,"PreMODIT")
ax.set_aspect(0.2/ax.get_data_ratio())
plt.show()

/tmp/ipykernel_20422/351839861.py:4: RuntimeWarning: divide by zero encountered in log10
  c=plt.imshow(np.log10(xsm),cmap="bone_r",vmin=-23,vmax=-19)
/tmp/ipykernel_20422/351839861.py:4: RuntimeWarning: invalid value encountered in log10
  c=plt.imshow(np.log10(xsm),cmap="bone_r",vmin=-23,vmax=-19)

len(xsm[xsm<0.0]), np.min(xsm)

(42278, DeviceArray(-5.3630485e-22, dtype=float32))

This negative value is very small. For instance, jnp.abs can remove it.

import jax.numpy as jnp
xsm=jnp.abs(xsm)

computing delta tau for CH4

from exojax.spec.rtransfer import dtauM
import jax.numpy as jnp
Rp=0.88
Mp=33.2
g=2478.57730044555*Mp/Rp**2 #gravity cm/s2
MMR=0.0059 #mass mixing ratio
dtaum=dtauM(dParr,xsm,MMR*np.ones_like(Tarr),molmassCH4,g)

computing delta tau for CIA

from exojax.spec.rtransfer import dtauCIA
mmw=2.33 #mean molecular weight
mmrH2=0.74
molmassH2=molinfo.molmass("H2")
vmrH2=(mmrH2*mmw/molmassH2) #VMR
dtaucH2H2=dtauCIA(nus,Tarr,Parr,dParr,vmrH2,vmrH2,\
            mmw,g,cdbH2H2.nucia,cdbH2H2.tcia,cdbH2H2.logac)

The total delta tau is a summation of them

dtau=dtaum+dtaucH2H2

you can plot a contribution function using exojax.plot.atmplot

from exojax.plot.atmplot import plotcf
plotcf(nus,dtau,Tarr,Parr,dParr)
plt.show()

radiative transfering…

from exojax.spec import planck
from exojax.spec.rtransfer import rtrun
sourcef = planck.piBarr(Tarr,nus)
F0=rtrun(dtau,sourcef)

fig=plt.figure(figsize=(20,4))
ax=fig.add_subplot(211)
plt.plot(wav[::-1],F0,lw=1,label="PreMODIT")
plt.legend()
plt.xlabel("wavelength ($\AA$)")
plt.savefig("ch4.png")

PreMODIT uses ESLOG as the wavenunmber grid. We can directly apply the response to the raw spectrum.

#response and rotation settings
from exojax.spec.response import ipgauss_sampling
from exojax.spec.spin_rotation import convolve_rigid_rotation
from exojax.utils.grids import velocity_grid
vsini_max = 100.0
vr_array = velocity_grid(resolution, vsini_max)


from exojax.utils.constants import c
import jax.numpy as jnp

wavd=jnp.linspace(16360,16560,1500) #observational wavelength grid
nusd = 1.e8/wavd[::-1]

RV=10.0 #RV km/s
vsini=20.0 #Vsini km/s
u1=0.0 #limb darkening u1
u2=0.0 #limb darkening u2

Rinst=100000. #spectral resolution of the spectrograph
beta=c/(2.0*np.sqrt(2.0*np.log(2.0))*Rinst) #IP sigma (STD of Gaussian)
Frot = convolve_rigid_rotation(F0, vr_array, vsini, u1, u2)
F = ipgauss_sampling(nusd, nus, Frot, beta, RV)

fig=plt.figure(figsize=(20,4))
plt.plot(wav[::-1],F0,alpha=0.5)
plt.plot(wavd[::-1],F)
plt.xlabel("wavelength ($\AA$)")
plt.savefig("moditCH4.png")

Let’s save the spectrum for the retrieval.

np.savetxt("spectrum_ch4_high.txt",np.array([wavd,F]).T,delimiter=",")