[ ]:

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

import petitRADTRANS as prt
from petitRADTRANS.retrieval.models import molliere_2020_emission, guillot_transmission
from petitRADTRANS.retrieval.models import mark_legacy_model_function
from petitRADTRANS.plotlib.style import set_petitradtrans_plot_style
from petitRADTRANS.retrieval.parameter import Parameter
from petitRADTRANS import physical_constants as cst

set_petitradtrans_plot_style()

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 legacy set of retrieval models. These models are kept to maintain backwards compatibility with pRT3, but we strongly recommend using the new interface demonstrated in the retrieval model tutorial.

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).

If define a model function using the model interface, which takes in a Radtrans object and a dictionary of Parameters, you need to decorate the function with mark_legacy_model_function, which tells the retrieval package what the expected data type is for the inputs and outputs. However, built in models share the same name between the two interfaces, and the retrieval package will automatically determine which interface is correct.

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.

molliere_2020_emission: 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.
molliere_2020_emission_two_column: 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 cloud_fraction parameter that sets the fraction of the planet covered in clouds.
guillot_emission: This model is based on the temperature profile of Guillot (2010). It can be used with either free or disequilibrium chemistry.
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.

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.

[ ]:

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),
}

[ ]:

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_000__DHS", "Fe(s)_crystalline_000__DHS")
pressures = jnp.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 '/Users/nasedkin/python-packages/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 '/Users/nasedkin/python-packages/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 '/Users/nasedkin/python-packages/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 '/Users/nasedkin/python-packages/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 '/Users/nasedkin/python-packages/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 '/Users/nasedkin/python-packages/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 '/Users/nasedkin/python-packages/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 '/Users/nasedkin/python-packages/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 '/Users/nasedkin/python-packages/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 '/Users/nasedkin/python-packages/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 '/Users/nasedkin/python-packages/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 '/Users/nasedkin/python-packages/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 '/Users/nasedkin/python-packages/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 '/Users/nasedkin/python-packages/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 '/Users/nasedkin/python-packages/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 '/Users/nasedkin/python-packages/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_000__DHS' from file '/Users/nasedkin/python-packages/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_000__DHS' from file '/Users/nasedkin/python-packages/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

[ ]:

wavelength, model = molliere_2020_emission(atmosphere, parameters, adaptive_mesh_refinement=False, pt_plot_mode=False)
pressure, temperature = molliere_2020_emission(atmosphere, parameters, adaptive_mesh_refinement=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_legacy_retrieval_models_9_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_legacy_retrieval_models_10_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.

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.

[ ]:

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),
}

[3]:

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 = jnp.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_legacy_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_legacy_retrieval_models_17_0.png

Retrievals: Exploring the built-in models

Contents

Retrievals: Exploring the built-in models#

Emission spectrum models#

Transmission spectrum models#