documents/tutorials/Transmission_beta.rst
Computes Transmission Spectra
=============================
We compute a high-resolution transmission spectrum using HITRAN CO and
compare it with the results by the different method (by Yui Kawashima).
Note that ``ArtTransPure`` is ``art`` for the transmission spectrum.
.. code:: ipython3
from jax import config
import pandas as pd
import numpy as np
from exojax.utils.grids import wavenumber_grid
from exojax.spec.opacalc import OpaPremodit
from exojax.spec.atmrt import ArtTransPure
from exojax.utils.constants import RJ, Rs
from exojax.spec.api import MdbHitran
from exojax.utils.astrofunc import gravity_jupiter
config.update("jax_enable_x64", True)
.. parsed-literal::
/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)
To load a reference spectrum, we make the following method.
.. code:: ipython3
def read_kawashima_data(filename):
dat = pd.read_csv(filename, delimiter=" ")
wav = dat["Wavelength[um]"]
mask = (wav > 2.25) & (wav < 2.6)
return wav[mask], dat["Rp/Rs"][mask]
The ref file should be in ExoJAX repo.
.. code:: ipython3
filename = "/home/kawahara/exojax/tests/integration/comparison/transmission/spectrum/CO100percent_500K.dat"
Here is the core code. We define ``art`` using ``ArtTransPure``.
``art.constant_mmr_profile`` sets the mass mixing ratio layer profile to
a constant value. Here we assume 100% CO everywhere. For a transmission
spectrum, the constant gravity assumption is not good approximation. So,
we need to compute the gravity profile. ``art.gravity_profile`` does it.
After computing opacity, we can use ``art.run`` to generate the
transmission spectrum.
.. code:: ipython3
def compare_with_kawashima_code():
mu_fid = 28.00863 # mean molecular weight
T_fid = 500.0
Nx = 50000
nu_grid, wav, res = wavenumber_grid(
22900.0, 24000.0, Nx, unit="AA", xsmode="premodit"
)
art = ArtTransPure(pressure_top=1.0e-15, pressure_btm=1.0e1, nlayer=100)
art.change_temperature_range(490.0, 510.0)
Tarr = T_fid * np.ones_like(art.pressure)
mmr_arr = art.constant_mmr_profile(1.0)
gravity_btm = gravity_jupiter(1.0, 1.0)
radius_btm = RJ
mdb_iso1 = MdbHitran("CO", nu_grid, gpu_transfer=False, isotope=1)
mdb_iso2 = MdbHitran("CO", nu_grid, gpu_transfer=False, isotope=2)
mdb_iso4 = MdbHitran("CO", nu_grid, gpu_transfer=False, isotope=4)
mmw = mu_fid * np.ones_like(art.pressure)
gravity = art.gravity_profile(Tarr, mmw, radius_btm, gravity_btm)
# HITRAN CO has not so many lines. So we use MODIT here instead of PreMODIT
opa_iso1 = OpaPremodit(
mdb=mdb_iso1, nu_grid=nu_grid, auto_trange=[T_fid - 10.0, T_fid + 10.0]
)
opa_iso2 = OpaPremodit(
mdb=mdb_iso2, nu_grid=nu_grid, auto_trange=[T_fid - 10.0, T_fid + 10.0]
)
opa_iso4 = OpaPremodit(
mdb=mdb_iso4, nu_grid=nu_grid, auto_trange=[T_fid - 10.0, T_fid + 10.0]
)
xsmatrix_iso1 = opa_iso1.xsmatrix(Tarr, art.pressure)
xsmatrix_iso2 = opa_iso2.xsmatrix(Tarr, art.pressure)
xsmatrix_iso4 = opa_iso4.xsmatrix(Tarr, art.pressure)
xsmatrix = xsmatrix_iso1 + xsmatrix_iso2 + xsmatrix_iso4
dtau = art.opacity_profile_xs(xsmatrix, mmr_arr, opa_iso1.mdb.molmass, gravity)
Rp2 = art.run(dtau, Tarr, mmw, radius_btm, gravity_btm)
return nu_grid, np.sqrt(Rp2) * radius_btm / Rs
Compares the results. Not so bad. Note that the Kawashima-san’s
reference includes the Rayleigh scattering, but we do not. That makes
this small difference.
.. code:: ipython3
import matplotlib.pyplot as plt
wav, rprs = read_kawashima_data(filename)
diffmode = 1
nus_hitran, Rp_hitran = compare_with_kawashima_code()
from exojax.spec.unitconvert import nu2wav
wav_exojax = nu2wav(nus_hitran, unit="um", wavelength_order="ascending")
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(wav, rprs * Rs / RJ, label="Kawashima")
ax.plot(wav_exojax[::-1], Rp_hitran * Rs / RJ, label="ExoJAX",ls="dashed")
plt.legend()
plt.xlabel("wavenumber cm-1")
plt.legend()
plt.ylabel("Rp (RJ)")
plt.xlim(2.29,2.35)
plt.savefig("comparison_iso.png")
plt.show()
.. parsed-literal::
/tmp/ipykernel_9525/4154571916.py:2: ParserWarning: Falling back to the 'python' engine because the 'c' engine does not support regex separators (separators > 1 char and different from '\s+' are interpreted as regex); you can avoid this warning by specifying engine='python'.
dat = pd.read_csv(filename, delimiter=" ")
/home/kawahara/exojax/src/exojax/spec/unitconvert.py:62: UserWarning: Both input wavelength and output wavenumber are in ascending order.
warnings.warn(
/home/kawahara/exojax/src/exojax/spec/atmrt.py:53: UserWarning: nu_grid is not given. specify nu_grid when using 'run'
warnings.warn(
.. parsed-literal::
xsmode = premodit
xsmode assumes ESLOG in wavenumber space: xsmode=premodit
======================================================================
The wavenumber grid should be in ascending order.
The users can specify the order of the wavelength grid by themselves.
Your wavelength grid is in *** descending *** order
======================================================================
integration: simpson
Simpson integration, uses the chord optical depth at the lower boundary and midppoint of the layers.
radis engine = vaex
radis engine = vaex
radis engine = vaex
OpaPremodit: params automatically set.
default elower grid trange (degt) file version: 2
Robust range: 485.78039920454563 - 680.34913472759 K
OpaPremodit: Tref_broadening is set to 499.89998999799934 K
OpaPremodit: gamma_air and n_air are used. gamma_ref = gamma_air/Patm
# of reference width grid : 7
# of temperature exponent grid : 2
.. parsed-literal::
uniqidx: 100%|██████████| 5/5 [00:00<00:00, 8852.48it/s]
.. parsed-literal::
Premodit: Twt= 637.741375149819 K Tref= 510.72954778036876 K
Making LSD:|####################| 100%
OpaPremodit: params automatically set.
default elower grid trange (degt) file version: 2
.. parsed-literal::
.. parsed-literal::
Robust range: 485.78039920454563 - 680.34913472759 K
OpaPremodit: Tref_broadening is set to 499.89998999799934 K
OpaPremodit: gamma_air and n_air are used. gamma_ref = gamma_air/Patm
# of reference width grid : 7
# of temperature exponent grid : 2
.. parsed-literal::
uniqidx: 100%|██████████| 5/5 [00:00<00:00, 31679.03it/s]
.. parsed-literal::
Premodit: Twt= 637.741375149819 K Tref= 510.72954778036876 K
Making LSD:|####################| 100%
OpaPremodit: params automatically set.
default elower grid trange (degt) file version: 2
Robust range: 485.78039920454563 - 680.34913472759 K
OpaPremodit: Tref_broadening is set to 499.89998999799934 K
OpaPremodit: gamma_air and n_air are used. gamma_ref = gamma_air/Patm
# of reference width grid : 6
# of temperature exponent grid : 2
.. parsed-literal::
uniqidx: 100%|██████████| 4/4 [00:00<00:00, 24745.16it/s]
.. parsed-literal::
Premodit: Twt= 637.741375149819 K Tref= 510.72954778036876 K
Making LSD:|####################| 100%
.. parsed-literal::
/home/kawahara/exojax/src/exojax/spec/unitconvert.py:62: UserWarning: Both input wavelength and output wavenumber are in ascending order.
warnings.warn(
.. image:: Transmission_beta_files/Transmission_beta_10_11.png