Get Started
Last update: July 2nd (2023) Hajime Kawahara for v1.4
First, we recommend 64-bit unless you can think precision seriously. Use jax.config to set 64-bit:
from jax.config import config
config.update("jax_enable_x64", True)
The following schematic figure explains how ExoJAX works; (1) loading
databases (*db
), (2) calculating opacity (opa
), (3) running
atmospheric radiative transfer (art
), (4) applying operations on the
spectrum (sop
)
from IPython.display import Image
Image("../exojax.png")
1. Loading a molecular database using mdb.
ExoJAX has an API for molecular databases, called mdb
(or adb
for atomic datbases). Prior to loading the database, define the
wavenumber range first.
from exojax.utils.grids import wavenumber_grid
nu_grid, wav, resolution = wavenumber_grid(1900.,
2300.,
100000,
unit="cm-1",
xsmode="premodit")
xsmode = premodit
xsmode assumes ESLOG in wavenumber space: mode=premodit
/home/kawahara/exojax/src/exojax/utils/grids.py:126: UserWarning: Resolution may be too small. R=523403.606697253
warnings.warn('Resolution may be too small. R=' + str(resolution),
Then, let’s load the molecular database. We here use Carbon monooxide in
Exomol. CO/12C-16O/Li2015
means
Carbon monooxide/ isotopes = 12C + 16O / database name
. You can
check the database name in the ExoMol website (https://www.exomol.com/).
from exojax.spec.api import MdbExomol
mdb = MdbExomol(".database/CO/12C-16O/Li2015", nurange=nu_grid)
/home/kawahara/exojax/src/exojax/utils/molname.py:133: FutureWarning: e2s will be replaced to exact_molname_exomol_to_simple_molname.
warnings.warn(
HITRAN exact name= (12C)(16O)
Background atmosphere: H2
Reading .database/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
2. Computation of the Cross Section using opa
ExoJAX has various opacity calculator classes, so-called opa
. Here,
we use a memory-saved opa, OpaPremodit
. We assume the robust
tempreature range we will use is 500-1500K.
from exojax.spec.opacalc import OpaPremodit
opa = OpaPremodit(mdb, nu_grid, auto_trange=[500.0, 1500.0])
OpaPremodit: params automatically set.
Robust range: 484.50562701065246 - 1804.6009417674848 K
Tref changed: 296.0K->521.067611616332K
Tref_broadening is set to 866.0254037844389 K
# of reference width grid : 4
# of temperature exponent grid : 2
uniqidx: 100%|██████████| 2/2 [00:00<00:00, 5174.96it/s]
Premodit: Twt= 1153.8856089961712 K Tref= 521.067611616332 K
Making LSD:|####################| 100%
Making LSD:|####################| 100%
Making LSD:|####################| 100%
Then let’s compute cross section for two different temperature 500 and 1500 K for P=1.0 bar. opa.xsvector can do that!
P = 1.0 #bar
T_1 = 500.0 #K
xsv_1 = opa.xsvector(T_1, P) #cm2
T_2 = 1500.0 #K
xsv_2 = opa.xsvector(T_2, P) #cm2
Plot them. It can be seen that different lines are stronger at different temperatures.
import matplotlib.pyplot as plt
plt.plot(nu_grid,xsv_1,label=str(T_1)+"K") #cm2
plt.plot(nu_grid,xsv_2,alpha=0.5,label=str(T_2)+"K") #cm2
plt.legend()
plt.xlabel("wavenumber (cm-1)")
plt.ylabel("cross section (cm2)")
plt.show()
You can also plot the line strengths at T=1500K. We can first change the
mdb
reference temperature and then plot the line intensity.
mdb.change_reference_temperature(T_2)
plt.plot(mdb.nu_lines,mdb.line_strength_ref,".")
plt.xlabel("wavenumber (cm-1)")
plt.ylabel("line strength (cm)")
plt.yscale("log")
plt.show()
Tref changed: 521.067611616332K->1500.0K
3. Atmospheric Radiative Transfer
ExoJAX can solve the radiative transfer and derive the emission
spectrum. To do so, ExoJAX has art
class. ArtEmisPure
means
Atomospheric Radiative Transfer for Emission with Pure absorption. So,
ArtEmisPure
does not include scattering. We set the number of the
atmospheric layer to 100 (nlayer) and the pressure at bottom and top
atmosphere to 100 and 1.e-8 bar.
from exojax.spec.atmrt import ArtEmisPure
art = ArtEmisPure(nu_grid=nu_grid, pressure_btm=1.e2, pressure_top=1.e-8, nlayer=100)
/home/kawahara/exojax/src/exojax/spec/dtau_mmwl.py:14: FutureWarning: dtau_mmwl might be removed in future.
warnings.warn("dtau_mmwl might be removed in future.", FutureWarning)
Let’s assume the power law temperature model, within 500 - 1500 K.
\(T = T_0 P^\alpha\)
where \(T_0=1200\) K and \(\alpha=0.1\).
art.change_temperature_range(500.0, 1500.0)
Tarr = art.powerlaw_temperature(1200.0,0.1)
Also, the mass mixing ratio of CO (MMR) should be defined.
mmr_profile = art.constant_mmr_profile(0.01)
Surface gravity is also important quantity of the atmospheric model, which is a function of planetary radius and mass. Here we assume 1 RJ and 10 MJ.
from exojax.utils.astrofunc import gravity_jupiter
gravity = gravity_jupiter(1.0,10.0)
In addition to the CO cross section, we would consider collisional
induced
absorption
(CIA) as a continuum opacity. cdb
class can be used.
from exojax.spec.contdb import CdbCIA
from exojax.spec.opacont import OpaCIA
cdb = CdbCIA(".database/H2-H2_2011.cia",nurange=nu_grid)
opacia = OpaCIA(cdb, nu_grid=nu_grid)
H2-H2
Before running the radiative transfer, we need cross sections for
layers, called xsmatrix
for CO and logacia_matrix
for CIA
(strictly speaking, the latter is not cross section but coefficient
because CIA intensity is proportional density square).
xsmatrix = opa.xsmatrix(Tarr, art.pressure)
logacia_matrix = opacia.logacia_matrix(Tarr)
Convert them to opacity
dtau_CO = art.opacity_profile_lines(xsmatrix, mmr_profile, mdb.molmass, gravity)
vmrH2 = 0.855 #VMR of H2
mmw = 2.33 # mean molecular weight of the atmosphere
dtaucia = art.opacity_profile_cia(logacia_matrix, Tarr, vmrH2, vmrH2, mmw, gravity)
Add two opacities.
dtau = dtau_CO + dtaucia
Then, run the radiative transfer
F = art.run(dtau, Tarr)
fig=plt.figure(figsize=(15,4))
plt.plot(nu_grid,F)
plt.xlabel("wavenumber (cm-1)")
plt.ylabel("flux (erg/s/cm2/cm-1)")
plt.show()
You can check the contribution function too! You should check if the
dominant contribution is within the layer. If not, you need to change
pressure_top
and pressure_btm
in ArtEmisPure
from exojax.plot.atmplot import plotcf
cf=plotcf(nu_grid, dtau, Tarr,art.pressure, art.dParr)
Spectral Operators: rotational broadening, instrumental profile, Doppler velocity shift and so on, any operation on spectra.
The above spectrum is called “raw spectrum” in ExoJAX. The effects
applied to the raw spectrum is handled in ExoJAX by the spectral
operator (sop
). First, we apply the spin rotational broadening of a
planet.
from exojax.spec.specop import SopRotation
sop_rot = SopRotation(nu_grid, resolution, vsini_max=100.0)
vsini = 50.0
u1=0.0
u2=0.0
Frot = sop_rot.rigid_rotation(F, vsini, u1, u2)
/home/kawahara/exojax/src/exojax/utils/grids.py:126: UserWarning: Resolution may be too small. R=523403.606697253
warnings.warn('Resolution may be too small. R=' + str(resolution),
fig=plt.figure(figsize=(15,4))
plt.plot(nu_grid,F, label="raw spectrum")
plt.plot(nu_grid,Frot, label="rotated")
plt.xlabel("wavenumber (cm-1)")
plt.ylabel("flux (erg/s/cm2/cm-1)")
plt.legend()
plt.show()
Then, the instrumental profile with relative radial velocity shift is
applied. Also, we need to match the computed spectrum to the data grid.
This process is called sampling
(but just interpolation though).
from exojax.spec.specop import SopInstProfile
from exojax.utils.instfunc import resolution_to_gaussian_std
sop_inst = SopInstProfile(nu_grid, resolution, vrmax=1000.0)
RV=40.0 #km/s
resolution_inst = 3000.0
beta_inst = resolution_to_gaussian_std(resolution_inst)
Finst = sop_inst.ipgauss(Frot, beta_inst)
nu_obs = nu_grid[::50]
Fobs = sop_inst.sampling(Finst, RV, nu_obs)
42.43671169022172
/home/kawahara/exojax/src/exojax/utils/grids.py:126: UserWarning: Resolution may be too small. R=523403.606697253
warnings.warn('Resolution may be too small. R=' + str(resolution),
fig=plt.figure(figsize=(15,4))
plt.plot(nu_grid,Frot, label="rotated")
plt.plot(nu_grid,Finst, label="rotated+IP")
plt.plot(nu_obs,Fobs, ".", label="rotated+IP (sampling)")
plt.xlabel("wavenumber (cm-1)")
plt.ylabel("flux (erg/s/cm2/cm-1)")
plt.legend()
plt.show()
That’s it.