Retrievals: Exploring the built-in models#

Written by Evert Nasedkin. Please cite pRT’s retrieval package (Nasedkin et al. 2024) in addition to pRT (Mollière et al. 2019) if you make use of the retrieval package for your work.

In this tutorial we’ll demonstrate how to use the built in set of retrieval models.

For simplicity, we’ll use the model directly. In the context of a retrieval, you’ll need to add the corresponding parameters to RetrievalConfig for the chosen forward model and set up their priors or fixed values. Then you need to associate the model function with the data that you’re performing the retrieval on.

We recommend using the built in models for most standard use cases, or to use SpectralModels (see the SpectralModel retrieval notebook).

[1]:

# Let's start by importing everything we need
import os

import matplotlib.pyplot as plt
import numpy as np

import petitRADTRANS as prt
from petitRADTRANS import physical_constants as cst
from petitRADTRANS.retrieval.parameter import Parameter
from petitRADTRANS.plotlib.style import set_petitradtrans_plot_style

set_petitradtrans_plot_style()

Using pRT Plotting style!

Emission spectrum models#

Let’s start with the emission models. For every model we need to supply two of the mass, R_pl or log_g in order to compute the surface gravity of the object. All of the models can combine various cloud implementations, see the cloud tutorial for more information. Most of the models can use either (dis)equilibrium or free chemistry, or some combination of the two. For the specific details of the required parameters, see the API documentation.

emission_model_diseq: This model is based on Mollière et al. (2020). It uses an adiabatic profile deep in the atmosphere, an Eddington profile in the mid regions and a spline profile in the upper atmosphere. It requires the use of (dis)equilibrium chemistry, as the temperature profile depends on both the atmospheric metallicity Fe/H (actually meaning [Fe/H]) and the C/O ratio.
emission_model_diseq_patchy_clouds: As above, this model is based on Mollière et al. (2020). However, here we’ll need to supply parameters for an additional temperature profile to describe the clear atmosphere region, as well as a patchiness parameter that sets the fraction of the planet covered in clouds.
emission_model_diseq_simple_patchy_clouds: Similarly, this is a patchy cloud model, but the entire planet uses the same temperature profile.
guillot_emission: This model is based on the temperature profile of Guillot (2010). It can be used with either free or disequilibrium chemistry.
guillot_patchy_emission: Also using the Guillot profile, but allowing for patchy clouds.
interpolated_profile_emission: This model uses a spline profile throughout the atmosphere. We can choose the number of nodes between which to interpolate, and whether to use a linear or cubic spline.
gradient_profile_emission: This model is based of the temperature gradient profile described in Zhang et al. (2023). This model uses tight priors on the temperature gradient at fixed pressure levels to determine the overall temperature profile.

[2]:

from petitRADTRANS.retrieval.models import (
    emission_model_diseq,
    emission_model_diseq_patchy_clouds,
    emission_model_diseq_simple_patchy_clouds,
    guillot_emission,
    guillot_patchy_emission,
    interpolated_profile_emission,
    gradient_profile_emission
)

We’ll only set up one model here, but the idea is the same for the rest of the models. We need to define all of the parameters necessary, setup a Radtrans object, and then we can calculate the spectrum. Here we’ll set up an emission spectrum for a planet similar to those in the HR 8799 system. Again, remember that this is different from the retrieval workflow, where the Retrieval object makes the Radtrans object; here we just want to calculate a model in isolation.

[3]:

parameters = {
    # Planet parameters
    'D_pl': Parameter('D_pl', False, value=10.0 * cst.pc),
    'mass': Parameter('mass', False, value=5 * cst.m_jup),
    'planet_radius': Parameter('planet_radius', False, value=1 * cst.r_jup_mean),
    # Temperature parameters
    'T_int': Parameter('T_int', False, value=1500.0),
    'T1': Parameter('T1', False, value=0.5),
    'T2': Parameter('T2', False, value=0.4),
    'T3': Parameter('T3', False, value=0.8),
    'log_delta': Parameter('log_delta', False, value=0.65),
    'alpha': Parameter('alpha', False, value=1.70),
    # Chemical parameters
    'Fe/H': Parameter('[Fe/H]', False, value=1.0),
    'C/O': Parameter('C/O', False, value=0.7),
    'log_pquench': Parameter('log_pquench', False, value=2.5),
    # Cloud parameters
    'sigma_lnorm': Parameter('sigma_lnorm', False, value=1.25),
    'fsed_MgSiO3(s)': Parameter('fsed_MgSiO3(s)', False, value=0.05),
    'fsed_Fe(s)': Parameter('fsed_Fe(s)', False, value=7.5),
    'log_kzz': Parameter('log_kzz', False, value=5),
    'eq_scaling_Fe(s)': Parameter('eq_scaling_Fe(s)', False, value=-0.3),
    'eq_scaling_MgSiO3(s)': Parameter('eq_scaling_MgSiO3(s)', False, value=-1.00)
}

[4]:

line_species = [
    'H2O__POKAZATEL.R300',
    'CO-NatAbund.R300',
    'CH4.R300',
    'CO2.R300',
    'HCN.R300',
    'FeH.R300',
    'H2S.R300',
    'NH3.R300',
    'PH3.R300',
    'Na.R300',
    'K.R300',
    'TiO.R300',
    'VO.R300',
    'SiO.R300'
]

rayleigh_species = ['H2', 'He']
continuum_opacities = ['H2-H2', 'H2-He']
cloud_species = ['MgSiO3(s)_crystalline__DHS', 'Fe(s)_crystalline__DHS']
pressures = np.logspace(-6, 2, 100)

atmosphere = prt.radtrans.Radtrans(
    pressures=pressures,
    line_species=line_species,
    rayleigh_species=rayleigh_species,
    gas_continuum_contributors=continuum_opacities,
    cloud_species=cloud_species,
    wavelength_boundaries=[1, 5]
)

Loading Radtrans opacities...
Loading line opacities of species 'H2O__POKAZATEL.R300' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/correlated_k/H2O/1H2-16O/1H2-16O__POKAZATEL.R300_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO-NatAbund.R300' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/correlated_k/CO/C-O-NatAbund/C-O-NatAbund__HITEMP.R300_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CH4.R300' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/correlated_k/CH4/12C-1H4/12C-1H4__YT34to10.R300_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO2.R300' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/correlated_k/CO2/12C-16O2/12C-16O2__UCL-4000.R300_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'HCN.R300' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/correlated_k/HCN/1H-12C-14N/1H-12C-14N__Harris.R300_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'FeH.R300' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/correlated_k/FeH/56Fe-1H/56Fe-1H__MoLLIST.R300_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'H2S.R300' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/correlated_k/H2S/1H2-32S/1H2-32S__AYT2.R300_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'NH3.R300' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/correlated_k/NH3/14N-1H3/14N-1H3__CoYuTe.R300_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'PH3.R300' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/correlated_k/PH3/31P-1H3/31P-1H3__SAlTY.R300_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'Na.R300' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/correlated_k/Na/23Na/23Na__Allard.R300_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'K.R300' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/correlated_k/K/39K/39K__Allard.R300_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'TiO.R300' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/correlated_k/TiO/48Ti-16O/48Ti-16O__McKemmish.R300_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'VO.R300' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/correlated_k/VO/51V-16O/51V-16O__VOMYT.R300_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'SiO.R300' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/correlated_k/SiO/28Si-16O/28Si-16O__SiOUVenIR.R300_0.1-50mu.ktable.petitRADTRANS.h5'... Done.
Successfully loaded all line opacities
Loading CIA opacities for H2-H2 from file '/home/dblain/petitRADTRANS/input_data/opacities/continuum/collision_induced_absorptions/H2--H2/H2--H2-NatAbund/H2--H2-NatAbund__BoRi.R831_0.6-250mu.ciatable.petitRADTRANS.h5'... Done.
Loading CIA opacities for H2-He from file '/home/dblain/petitRADTRANS/input_data/opacities/continuum/collision_induced_absorptions/H2--He/H2--He-NatAbund/H2--He-NatAbund__BoRi.DeltaWavenumber2_0.5-500mu.ciatable.petitRADTRANS.h5'... Done.
Successfully loaded all CIA opacities
Loading opacities of cloud species 'MgSiO3(s)_crystalline__DHS' from file '/home/dblain/petitRADTRANS/input_data/opacities/continuum/clouds/MgSiO3(s)_crystalline_000/Mg-Si-O3-NatAbund(s)_crystalline_000/Mg-Si-O3-NatAbund(s)_crystalline_000__DHS.R39_0.1-250mu.cotable.petitRADTRANS.h5'(crystalline_000, using DHS scattering)... Done.
Loading opacities of cloud species 'Fe(s)_crystalline__DHS' from file '/home/dblain/petitRADTRANS/input_data/opacities/continuum/clouds/Fe(s)_crystalline_000/Fe-NatAbund(s)_crystalline_000/Fe-NatAbund(s)_crystalline_000__DHS.R39_0.1-250mu.cotable.petitRADTRANS.h5'(crystalline_000, using DHS scattering)... Done.
Successfully loaded all clouds opacities
Successfully loaded all opacities

[5]:

wavelength, model = emission_model_diseq(atmosphere, parameters, amr=False, pt_plot_mode=False)
pressure,temperature = emission_model_diseq(atmosphere, parameters, amr=False, pt_plot_mode=True)

Loading chemical equilibrium chemistry table from file '/home/dblain/petitRADTRANS/input_data/pre_calculated_chemistry/equilibrium_chemistry/equilibrium_chemistry.chemtable.petitRADTRANS.h5'... Done.

[6]:

fig, ax = plt.subplots()

ax.plot(wavelength, model, linewidth=3)
ax.set_xlabel("Wavelength [micron]")
ax.set_ylabel(r"F$_{\lambda}$ [W/m$^{2}/\mu$m]")

[6]:

Text(0, 0.5, 'F$_{\\lambda}$ [W/m$^{2}/\\mu$m]')

../../_images/content_notebooks_retrieval_models_8_1.png

[7]:

fig, ax = plt.subplots()

ax.plot(temperature, pressure, linewidth=3)
ax.set_xlabel("Temperature [K]")
ax.set_ylabel(r"Pressure [bar]")
ax.set_ylim(1e2,1e-6)
ax.set_yscale('log')

../../_images/content_notebooks_retrieval_models_9_0.png

Transmission spectrum models#

We can do the same for the transmission spectrum models:

For every model we need to supply two of the mass, R_pl or log_g in order to compute the surface gravity of the object. All of the models can combine various cloud implementations, see the cloud tutorial for more information. All of the transmission models can use either (dis)equilibrium or free chemistry, or some combination of the two. For the specific details of the required parameters, see the API documentation.

isothermal_transmission: This model assumes a constant temperature throughout the atmosphere. Clouds can be patchy or have complete coverage.
guillot_transmission: This model uses the Guillot (2010) profile,.
guillot_patchy_transmission: Also based on Guillot (2010), this model also allows for the use of patchy clouds.
madhu_seager_patchy_transmission: This model implements the Madhushudan and Seager (2009) temperature profile, and also allows for patchy clouds.

[8]:

from petitRADTRANS.retrieval.models import (
    isothermal_transmission,
    guillot_transmission,
    guillot_patchy_transmission,
    madhu_seager_patchy_transmission
)

For our transmission example, we’ll use a WASP-39b-like planet.

Again, remember that what we do below is different from the retrieval workflow, where the Retrieval object makes the Radtrans object; here we just want to calculate a model in isolation.

[9]:

parameters = {
    # Star parameters
    'stellar_radius': Parameter('stellar_radius', False, value=0.9324 * cst.r_sun),
    # Planet parameters
    'planet_radius': Parameter('planet_radius', False, value=1.3 * cst.r_jup_mean),
    'log_g': Parameter('log_g', False, value=2.75),
    # Temeprature parameters
    'T_int': Parameter('T_int', False, value=750),
    'T_equ': Parameter('T_equ', False, value=600),
    'log_kappa_IR': Parameter('log_kappa_IR', False, value=-1.0),
    'gamma': Parameter('gamma', False, value=1.0),
    # Cloud parameters
    'haze_factor': Parameter('haze_factor', False, value=0.0),
    'power_law_opacity_350nm': Parameter('power_law_opacity_350nm', False, value=1e-4),
    'power_law_opacity_coefficient': Parameter('power_law_opacity_coefficient', False, value=-1.9),
    'log_Pcloud': Parameter('log_Pcloud', False, value=0),
    # Chemical parameters
    'Fe/H': Parameter('Fe/H', False, value=1.0),
    'C/O': Parameter('C/O', False, value=0.1)
}

[10]:

line_species = [
    'H2O__POKAZATEL.R300',
    'CO-NatAbund.R300',
    'CH4.R300',
    'CO2.R300',
    'HCN.R300',
    'FeH.R300',
    'H2S.R300',
    'NH3.R300',
    'PH3.R300',
    'Na.R300',
    'K.R300',
    'TiO.R300',
    'VO.R300',
    'SiO.R300'
]

rayleigh_species = ['H2', 'He']
continuum_opacities = ['H2-H2', 'H2-He']
cloud_species = []
pressures = np.logspace(-6,2,100)

atmosphere = prt.radtrans.Radtrans(
    pressures=pressures,
    line_species=line_species,
    rayleigh_species=rayleigh_species,
    gas_continuum_contributors=continuum_opacities,
    cloud_species=cloud_species,
    wavelength_boundaries=[1, 5]
)

[11]:

wavelength, model = guillot_transmission(atmosphere, parameters, amr = False, pt_plot_mode = False)
pressure,temperature = guillot_transmission(atmosphere, parameters, amr = False, pt_plot_mode = True)

[12]:

fig, ax = plt.subplots()
ax.plot(wavelength, model*1e6, linewidth = 3)
ax.set_xlabel("Wavelength [micron]")
ax.set_ylabel(r"F$_{p}$/F$_{\rm star}$ [ppm]")

[12]:

Text(0, 0.5, 'F$_{p}$/F$_{\\rm star}$ [ppm]')

../../_images/content_notebooks_retrieval_models_16_1.png

[13]:

fig, ax = plt.subplots()
ax.plot(temperature, pressure, linewidth = 3)
ax.set_xlabel("Temperature [K]")
ax.set_ylabel(r"Pressure [bar]")
ax.set_ylim(1e2,1e-6)
ax.set_yscale('log')

../../_images/content_notebooks_retrieval_models_17_0.png

Retrievals: Exploring the built-in models

Contents

Retrievals: Exploring the built-in models#

Emission spectrum models#

Transmission spectrum models#