Getting Started

Welcome to petitRADTRANS! This series of tutorials will outline the basic features of pRT, from how to generate emission and transmission spectra to running atmospheric retrievals. If you’d like to run the tutorial notebooks yourself, they’re available on gitlab.

Preparing a radiative transfer object

After installation has finished, let’s try to calculate a first spectrum! First we load required external packages though. We’ll use the jax.numpy package instead of the standard numpy, just to familiarize ourselves with it, but it’s almost exactly the same. However, by default JAX uses 32-bit floats, so we’ll set it to 64-bit mode to ensure that we have full numerical precision. Transmission spectra and emission spectra without scattering can be calculated single precision, but scattering in emission requires setting the double precision mode.

JAX floating point precision must be set before importing any petitRADTRANS modules!

from jax import config
config.update("jax_enable_x64", True)
import jax.numpy as jnp
import matplotlib.pyplot as plt

Then we’ll load pRT by importing:

import petitRADTRANS
from petitRADTRANS.radtrans import Radtrans
print(f"petitRADTRANS version: {petitRADTRANS.__version__}")

petitRADTRANS version: 4.0.0a35

We also need the physical constants module:

from petitRADTRANS import physical_constants as cst

The Radtrans object

Let’s start out by creating a radiative transfer object using the Radtrans object.

This will load the requested opacities and create an object ready to calculate spectra once the atmospheric parameters (temperature, abundances, mean molar mass etc.) are set. We will explain the meaning of the function’s inputs below.

radtrans = Radtrans(
    pressures=jnp.logspace(-6, 2, 100),
    line_species=("H2O.R500", "CO-NatAbund.R500", "CH4.R500", "CO2.R500", "Na.R500", "K.R500"),
    rayleigh_species=("H2", "He"),
    gas_continuum_contributors=("H2--H2", "H2--He"),
    wavelength_boundaries=(0.3, 15.0),
)

Binned-down file not found, starting rebinning...
Searching for the file of species 'H2O' with default spectral sampling (R1000)...
 File found.

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[50], line 1
----> 1 radtrans = Radtrans(
      2     pressures=jnp.logspace(-6, 2, 100),
      3     line_species=("H2O.R500", "CO-NatAbund.R500", "CH4.R500", "CO2.R500", "Na.R500", "K.R500"),
      4     rayleigh_species=("H2", "He"),
      5     gas_continuum_contributors=("H2--H2", "H2--He"),
      6     wavelength_boundaries=(0.3, 15.0),
      7 )

File ~/python-packages/petitRADTRANS/petitRADTRANS/radtrans.py:241, in Radtrans.__init__(self, pressures, wavelength_boundaries, line_species, gas_continuum_contributors, rayleigh_species, cloud_species, line_opacity_mode, line_by_line_opacity_sampling, scattering_in_emission, emission_angle_grid, anisotropic_cloud_scattering, retain_line_opacities_for_plotting, path_input_data)
    233 self.__absorber_present = not (
    234     len(self._line_species) == 0
    235     and len(self._gas_continuum_contributors) == 0
    236     and len(self._rayleigh_species) == 0
    237     and len(self._cloud_species) == 0
    238 )
    240 # Initialize line parameters
--> 241 self._frequencies, self._frequency_bins_edges = self._init_frequency_grid()
    243 # Initialize loaded line opacities variables
    244 self._lines_loaded_opacities = LockedDict.build_and_lock(
    245     {
    246         "temperature_pressure_grid": {},
   (...)
    264     }
    265 )

File ~/python-packages/petitRADTRANS/petitRADTRANS/radtrans.py:1764, in Radtrans._init_frequency_grid(self)
   1748 """Initialize the Radtrans frequency grid, used to calculate the spectra.
   1749 The frequency grid comes from the requested opacity files, in the following priority order:
   1750     1. lines,
   (...)
   1761         for correlated-k only, number of points used to sample the g-space (1 in the case lbl is used)
   1762 """
   1763 if len(self._line_species) > 0:
-> 1764     frequencies, frequency_bins_edges = self._init_frequency_grid_from_lines()
   1765 elif len(self._gas_continuum_contributors) > 0:
   1766     hdf5_file = cia.get_cia_aliases(self._gas_continuum_contributors[0])

File ~/python-packages/petitRADTRANS/petitRADTRANS/radtrans.py:1874, in Radtrans._init_frequency_grid_from_lines(self)
   1869 import h5py
   1871 if self._line_opacity_mode == "c-k":  # correlated-k
   1872     # Get dimensions of molecular opacity arrays for a given P-T point, they define the resolution
   1873     # Use the first entry of self.line_species for this, if given
-> 1874     opacities_file = load.get_line_opacity_file(
   1875         path_input_data=self._path_input_data,
   1876         species=self._line_species[0],
   1877         category="correlated_k_opacities",
   1878     )
   1880     with h5py.File(opacities_file, "r") as f:
   1881         frequency_bins_edges = cst.c * f["bin_edges"][:][::-1]

File ~/python-packages/petitRADTRANS/petitRADTRANS/opacities/load.py:282, in get_line_opacity_file(path_input_data, species, category)
    271 default_hdf5_file = CorrelatedKOpacity.find(
    272     path_input_data=path_input_data,
    273     category=category,
   (...)
    276     search_online=True,
    277 )
    278 print(
    279     " File found.",
    280 )
--> 282 state = rebin_ck_line_opacities(
    283     input_file=default_hdf5_file,
    284     target_resolving_power=target_resolving_power,
    285     wavenumber_grid=None,
    286     rewrite=False,
    287 )
    289 if state == -1:
    290     raise RuntimeError(
    291         "unable to perform binning down, please install exo_k"
    292     )

File ~/python-packages/petitRADTRANS/petitRADTRANS/opacities/load.py:423, in rebin_ck_line_opacities(input_file, target_resolving_power, wavenumber_grid, rewrite)
    420 _, comm, rank = _get_mpi_context()
    422 try:
--> 423     import exo_k
    424 except ImportError:
    425     # Only raise a warning to give a chance to download the binned-down file
    426     warnings.warn(
    427         "binning down of opacities requires exo_k to be installed, no binning down has been performed"
    428     )

File /opt/anaconda3/envs/jaxprt/lib/python3.12/site-packages/exo_k/__init__.py:8
      2 """
      3 @author: jeremy leconte
      4 __init__ module to load the exo_k library
      5 """
      6 from .__version__ import __version__
----> 8 from exo_k.util.user_func import *
      9 from exo_k.two_stream import two_stream_crisp, two_stream_toon, two_stream_lmdz
     10 from exo_k.util.radiation import *

File /opt/anaconda3/envs/jaxprt/lib/python3.12/site-packages/exo_k/util/__init__.py:6
      1 # -*- coding: utf-8 -*-
      2 """
      3 @author: jeremy leconte
      4 __init__ module to load util subpackage
      5 """
----> 6 from .radiation import *
      7 from .cst import *
      8 from .spectrum import Spectrum

File /opt/anaconda3/envs/jaxprt/lib/python3.12/site-packages/exo_k/util/radiation.py:11
      9 import os
     10 import numpy as np
---> 11 import numba
     12 import scipy.integrate as integrate
     13 from exo_k.util.cst import PI,PLANCK,C_LUM,KBOLTZ,SIG_SB

File /opt/anaconda3/envs/jaxprt/lib/python3.12/site-packages/numba/__init__.py:77
     74 from numba.core import types, errors
     76 # Re-export typeof
---> 77 from numba.misc.special import (
     78     typeof, prange, pndindex, gdb, gdb_breakpoint, gdb_init,
     79     literally, literal_unroll,
     80 )
     82 # Re-export error classes
     83 from numba.core.errors import *

File /opt/anaconda3/envs/jaxprt/lib/python3.12/site-packages/numba/misc/special.py:3
      1 import numpy as np
----> 3 from numba.core.typing.typeof import typeof
      4 from numba.core.typing.asnumbatype import as_numba_type
      7 def pndindex(*args):

File /opt/anaconda3/envs/jaxprt/lib/python3.12/site-packages/numba/core/typing/__init__.py:1
----> 1 from .context import BaseContext, Context
      2 from .templates import (signature, make_concrete_template, Signature,
      3                         fold_arguments)

File <frozen importlib._bootstrap>:1360, in _find_and_load(name, import_)

File <frozen importlib._bootstrap>:1331, in _find_and_load_unlocked(name, import_)

File <frozen importlib._bootstrap>:935, in _load_unlocked(spec)

File <frozen importlib._bootstrap_external>:991, in exec_module(self, module)

File <frozen importlib._bootstrap_external>:1087, in get_code(self, fullname)

File <frozen importlib._bootstrap_external>:1186, in get_data(self, path)

KeyboardInterrupt:

The first time you request an opacity, you may be prompted with a message asking you to select an opacity file among a numbered list. In that case, type the integer corresponding to the file you want to download, then press “Enter”. More information about loading the opacities in the next section.

First we define the pressure array of the atmosphere. Note that the pressures must always be sorted in increasing order, and be equidistant in log-space. The pressure is given to pRT in units of bar, although the internal units in petitRADTRANS are in cgs. Typically we recommend using around 100 layers in your computations. However, be careful, internally pRT converts the pressure to cgs units!

Next we define the line absorbers to be used for the calculation (line_species = ...). If you have a fresh pRT installation, your opacity folder will be empty. pRT will automatically download the species which are missing. Here we asked for the following molecular or atomic line absorbers to be included: H2O, CO-NatAbund (all isotopologues at naturally occurring isotope abundances on Earth), CH4, CO2, Na, and K. More information on the data loading process in the next section.

Additionally, the Rayleigh scattering cross-sections for H2 and He were loaded, as well as the collision induced absorption (CIA) cross-sections for the H2--H2 and H2--He pairs, called gas_continuum_contributors above.

The wavelength range that was loaded for the opacities here extends from 0.3 to 15 microns.

And we are done! This step only needs to be carried out once, after this we are free to play with this radiative transfer (Radtrans) object, for example change the temperature, abundances, etc. Only if additional opacity species need to be loaded is it necessary to create a new Radtrans instance.

Units in petitRADTRANS: all units inside petitRADTRANS are in cgs. However, when interfacing with the code, you are expected to provide pressures in bars (more intuitive). Pressures will be converted to cgs units within the code.

Loading opacities

Opacities can have different sources. You can either select one of them by adding tags such as __HITEMP or __Allard to the absorber names (e.g., H2O__HITEMP or H2O__Exomol) or by just writing H2O and setting a default line list, that will be saved for future uses (see information in the next section for more details). pRT will do this by automatically asking you the first time you request a species (e.g., H2O) for which multiple options are available.

If you do not want to run at pRT’s standard resolution for correlated-k opacities (\(\lambda/\Delta\lambda = 1,000\)) you can request another spectral resolution with .R300, so H2O.R300, for example. This will bin down the opacities accordingly, and add them to the input data folder in this new resolution. This is described more in “Rebinning opacities”.

More information on the file naming convention and opacity request are available in the opacities section.

Configuring the input_data folder

One of the most important things to setup pRT is to give the path to the input_data folder, where critical data such as opacities are stored.

This path is managed by the pRT config file. You may have seen that it has been automatically generated after importing a Radtrans object for the first time.

By default, the input_data folder is placed into your home directory, under <home>/petitRADTRANS/input_data. To set the input_data path manually, the simplest way is to execute the following, giving the path corresponding to your setup:

from petitRADTRANS.config import petitradtrans_config_parser
# petitradtrans_config_parser.set_input_data_path(r'~/petitRADTRANS/input_data')  # replace this path with the actual path to your input_data directory  # noqa: E501

petitRADTRANS config file: petitradtrans_config_parser is an object derived from the ConfigParser Python standard object. It is a way to easily manage customizable software configuration. When instantiated for the first time, petitradtrans_config_parser will create a petitradtrans_config_file.ini file in a .petitradtrans directory in your home directory. You can access this file and modify it manually if you want to. The default opacities to load (in case there are multiple options, e.g., H2O__HITEMP or H2O__Exomol for water), as well as the link to download opacities, are stored in this file too.

Locating the input_data folder

The input_data folder can be located as follows:

petitradtrans_config_parser.get_input_data_path()

'/Users/nasedkin/python-packages/petitRADTRANS/input_data'

Calculating a transmission spectrum

Let’s calculate a spectrum for a hot Jupiter, first in transmission. Let’s assume an isothermal temperature profile, for now. We start by specifying the temperature array, and abundance dictionary, containing the mass fractions of the atmospheric material, at every layer, and the mean molar masses, at every layer, in units of g:

# note that radtrans.pressures is in cgs units now, multiply by 1e-6 to get bars
temperatures = 1200 * jnp.ones_like(radtrans.pressures)

mass_fractions = {
    "H2": 0.74 * jnp.ones(temperatures.size),
    "He": 0.24 * jnp.ones(temperatures.size),
    "H2O": 1e-3 * jnp.ones(temperatures.size),
    "CO-NatAbund": 1e-2 * jnp.ones(temperatures.size),
    "CO2": 1e-4 * jnp.ones(temperatures.size),
    "CH4": 1e-5 * jnp.ones(temperatures.size),
    "Na": 1e-4 * jnp.ones(temperatures.size),
    "K": 1e-6 * jnp.ones(temperatures.size),
}

#  2.33 is a typical value for H2-He dominated atmospheres
mean_molar_masses = 2.33 * jnp.ones(temperatures.size)

Note the P-T profile of a real planet will likely be different (e.g., not isothermal). Also the abundances are made up for this simple example here (see “Interpolating chemical equilibrium abundances” for how to improve this.)

Abundances in petitRADTRANS: abundances in pRT are in units of mass fractions, not number fractions (aka volume mixing ratio, VMR). One can convert between mass fractions and VMRs by using \begin{equation} X_i = \frac{\mu_i}{\mu}n_i, \end{equation} where \(X_i\) is the mass fraction of species \(i\), \(\mu_i\) the molar mass of a molecule/atom/ion/… of species \(i\), \(\mu\) is the atmospheric mean molar mass, and \(n_i\) is the VMR of species \(i\). This is implemented in petitRADTRANS.chemistry.utils.mass_fractions2volume_mixing_ratios() and petitRADTRANS.chemistry.utils.volume_mixing_ratios2mass_fractions().

Next, let’s assume a planetary radius and gravity at a given reference pressure (in bar!) in the atmosphere, and we are all set!

planet_radius = 1.0 * cst.r_jup_mean
reference_gravity = 10**3.5
reference_pressure = 0.01

Now we simply generate a transmission spectrum. The calculate_transit_radii function returns wavelengths in cm and transit radius in cm, as well as a dictionary of additional optional outputs, which is empty in our example, since we did not request anything in addition.

wavelengths, transit_radii, _ = radtrans.calculate_transit_radii(
    temperatures=temperatures,
    mass_fractions=mass_fractions,
    mean_molar_masses=mean_molar_masses,
    reference_gravity=reference_gravity,
    planet_radius=planet_radius,
    reference_pressure=reference_pressure,
)

Let’s plot the transit radius!

fig, ax = plt.subplots(figsize=(10, 6))

ax.plot(wavelengths * 1e4, transit_radii / cst.r_jup_mean)
ax.set_xscale("log")
ax.set_xlabel("Wavelength [microns]")
ax.set_ylabel(r"Transit radius [$\rm R_{Jup}$]")

Text(0, 0.5, 'Transit radius [$\\rm R_{Jup}$]')

from getting_started_transmission_explorer import display_getting_started_transmission_explorer
display_getting_started_transmission_explorer(height=1040)

Calculating an emission spectrum

Scattering and petitRADTRANS: scattering is included for emission spectra in petitRADTRANS only if requested specifically when generating the Radtrans object, as it increases the runtime (see “Scattering for Emission Spectra” for an example on how to do this). We neglect the scattering here.

Let’s calculate an emission spectrum. Using the isothermal temperature structure from above will simply result in a black body spectrum. pRT returns all spectra in cgs units at the planetary “surface” (i.e., the top of the atmosphere).

wavelengths, flux, _ = radtrans.calculate_flux(
    temperatures=temperatures,
    mass_fractions=mass_fractions,
    mean_molar_masses=mean_molar_masses,
    reference_gravity=reference_gravity,
)

fig, ax = plt.subplots(figsize=(10, 6))

ax.plot(wavelengths * 1e4, flux)

plt.xscale("log")
ax.set_xscale("log")
ax.set_xlabel("Wavelength [microns]")
ax.set_ylabel(r"Planet flux, $F_{\lambda}$ [erg cm$^{-2}$ s$^{-1}$ cm$^{-1}$]")

Text(0, 0.5, 'Planet flux, $F_{\\lambda}$ [erg cm$^{-2}$ s$^{-1}$ cm$^{-1}$]')

Standard flux units: before pRT3 flux was accessed as atmosphere.flux after running atmosphere.calc_flux(), which contained flux as \(F_\nu\), so in units of erg cm\(^{-2}\) s\(^{-1}\) Hz\(^{-1}\). pRT’s calculate_flux() method now returns wavelength and flux as \(F_\lambda\) in its standard setting, so flux in erg cm\(^{-2}\) s\(^{-1}\) cm\(^{-1}\). To return frequencies and \(F_\nu\), instead of wavelengths and \(F_\lambda\), please set the keyword frequencies_to_wavelengths=False when calling calculate_flux().

Let’s try a different temperature structure instead, namely the often-used analytical profile from Guillot (2010), his Equation 29:

\begin{equation} T^4 = \frac{3T_{\rm int}^4}{4}\left(\frac{2}{3}+\tau\right) + \frac{3T_{\rm equ}^4}{4}\left[\frac{2}{3}+\frac{1}{\gamma\sqrt{3}}+\left(\frac{\gamma}{\sqrt{3}}-\frac{1}{\gamma\sqrt{3}}\right)e^{-\gamma\tau\sqrt{3}}\right] \end{equation}

with \(\tau = P\kappa_{\rm IR}/g\). Here, \(\tau\) is the optical depth, \(P\) the pressure, \(\kappa_{\rm IR}\) is the atmospheric opacity in the IR wavelengths (i.e., the cross-section per unit mass), \(\gamma\) is the ratio between the optical and IR opacity, \(T_{\rm equ}\) the atmospheric equilibrium temperature, and \(T_{\rm int}\) is the planetary internal temperature.

from petitRADTRANS.temperature_profiles import guillot_global_temperature_profile

pressures = jnp.asarray(radtrans.pressures) * 1e-6  # cgs to bar

infrared_mean_opacity = 0.01
gamma = 0.4
intrinsic_temperature = 200
equilibrium_temperature = 1500


temperatures = guillot_global_temperature_profile(
    pressures=pressures,
    infrared_mean_opacity=infrared_mean_opacity,
    gamma=gamma,
    gravities=reference_gravity,
    intrinsic_temperature=intrinsic_temperature,
    equilibrium_temperature=equilibrium_temperature,
)

Let’s plot the P-T profile:

fig, ax = plt.subplots(figsize=(10, 6))

ax.plot(temperatures, pressures)
ax.set_yscale("log")
ax.set_ylim([1e2, 1e-6])
ax.set_xlabel("T [K]")
ax.set_ylabel("P [bar]")

Text(0, 0.5, 'P [bar]')

And recalculate the spectrum:

wavelengths, flux, _ = radtrans.calculate_flux(
    temperatures=temperatures,
    mass_fractions=mass_fractions,
    mean_molar_masses=mean_molar_masses,
    reference_gravity=reference_gravity,
)
fig, ax = plt.subplots(figsize=(10, 6))

ax.plot(wavelengths * 1e4, flux)

ax.set_xscale("log")
ax.set_xlabel("Wavelength [microns]")
ax.set_ylabel(r"Planet flux, $F_{\lambda}$ [erg cm$^{-2}$ s$^{-1}$ cm$^{-1}$]")

Text(0, 0.5, 'Planet flux, $F_{\\lambda}$ [erg cm$^{-2}$ s$^{-1}$ cm$^{-1}$]')

from getting_started_emission_explorer import display_getting_started_emission_explorer
display_getting_started_emission_explorer(height=1040)

Often when comparing models to data, we want to know what the spectrum would be as observed from earth. If we know the distance to the planet and the planet radius, we can calculate this using the function petitRADTRANS.physics.flux2irradiance().

An alternative: the SpectralModel object

The Radtrans object is pRT’s foundation. It gives a lot of control on the spectral parameterization, as shown in this notebook.

An alternative to Radtrans is the SpectralModel object, which uses Radtrans as a parent object, and automatise most of the steps required to obtain a spectrum. The SpectralModel object comes with default functions to automatically calculate the temperature, mass fractions, etc, and these functions can easily be customised. More information can be found in “SpectralModel”.