`SpectralModel`

Written by Doriann Blain.

The goal of this notebook is to give an introduction to the possibilites offered by `SpectralModel`

.

`SpectralModel`

objects are just like `Radtrans`

objects (the former is actually a child of the latter). However, with `SpectralModel`

, all of the steps required to build a forward model can be taken care of internally by built-in functions. This includes, for example, calculating the temperature profile, or the mean molar mass. Then, the spectrum is calculated via the relevant `Radtrans`

functions. The output can then be modified by, for example, applying a convolution, or shifting the
wavelengths.

`SpectralModel`

**built-in functions:** these built-in functions, starting with `compute_`

, are **modular**, and can be modified. For example, you can change the way the temperature profile is computed, and even the order in which the spectral parameters are calculated. More on this at the end of this notebook.

Below is a flowchart of the `SpectralModel`

workflow. The “Spectral function” label refers to `calcualate_emission_spectrum`

or `calculate_transit_radii`

.

We make some useful imports below.

```
[1]:
```

```
import os
import matplotlib.pyplot as plt
import numpy as np
import petitRADTRANS.physical_constants as cst
```

## Basic usage

In this section we will cover the most basic usage of `SpectralModel`

, that is, using it as a regular `Radtrans`

object.

`SpectralModel`

is imported as follows:

```
[2]:
```

```
from petitRADTRANS.spectral_model import SpectralModel
```

Let’s start out by creating a `SpectralModel`

object.

This will load the requested opacities and create an object ready to calculate spectra.

```
[3]:
```

```
spectral_model = SpectralModel(
pressures=np.logspace(-6, 2, 100),
line_species=[
'H2O',
'CO-NatAbund',
'CH4',
'CO2',
'Na',
'K'
],
rayleigh_species=['H2', 'He'],
gas_continuum_contributors=['H2-H2', 'H2-He'],
wavelength_boundaries=[0.3, 15]
)
```

```
Loading Radtrans opacities...
Loading line opacities of species 'H2O' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/H2O/1H2-16O/1H2-16O__POKAZATEL.R1000_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO-NatAbund' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CO/C-O-NatAbund/C-O-NatAbund__HITEMP.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CH4' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CH4/12C-1H4/12C-1H4__HITEMP.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO2' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CO2/12C-16O2/12C-16O2__UCL-4000.R1000_0.3-50mu.ktable.petitRADTRANS.h5.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'Na' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/Na/23Na/23Na__Allard.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'K' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/K/39K/39K__Allard.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Successfully loaded all line opacities
Loading CIA opacities for H2-H2 from file '/Users/molliere/Desktop/input_data_v3/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/molliere/Desktop/input_data_v3/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
Successfully loaded all opacities
```

As you can see, this `SpectralModel`

was initialised in the exact same way as a `Radtrans`

object. However, it is more interesting to use the extra features of `SpectralModel`

that we will touch below.

## Intermediate usages

In this section we will cover the concept of model parameters and introduce most of them.

### Calculating a transmission spectrum: a simple example

Here we will generate the same transmission spectrum as in the “Getting Started” section, but making full use of the `SpectralModel`

features. The new arguments will be introduced below.

```
[4]:
```

```
spectral_model = SpectralModel(
# Radtrans parameters
pressures=np.logspace(-6, 2, 100),
line_species=[
'H2O',
'CO-NatAbund',
'CH4',
'CO2',
'Na',
'K'
],
rayleigh_species=['H2', 'He'],
gas_continuum_contributors=['H2-H2', 'H2-He'],
wavelength_boundaries=[0.3, 15],
# SpectralModel additional parameters
# Planet parameters
planet_radius=1 * cst.r_jup_mean,
reference_gravity=10 ** 3.5,
reference_pressure=1e-2,
# Temperature profile parameters (isothermal by default)
temperature=1200,
# Mass fractions
imposed_mass_fractions={ # these can also be arrays of the same size as pressures
'H2O': 1e-3,
'CO-NatAbund': 1e-2,
'CH4': 1e-5,
'CO2': 1e-4,
'Na': 1e-4,
'K': 1e-6
},
filling_species={ # automatically fill the atmosphere with H2 and He, such that the sum of MMRs is equal to 1 and H2/He = 37/12
'H2': 37,
'He': 12
}
)
```

```
Loading Radtrans opacities...
Loading line opacities of species 'H2O' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/H2O/1H2-16O/1H2-16O__POKAZATEL.R1000_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO-NatAbund' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CO/C-O-NatAbund/C-O-NatAbund__HITEMP.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CH4' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CH4/12C-1H4/12C-1H4__HITEMP.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO2' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CO2/12C-16O2/12C-16O2__UCL-4000.R1000_0.3-50mu.ktable.petitRADTRANS.h5.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'Na' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/Na/23Na/23Na__Allard.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'K' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/K/39K/39K__Allard.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Successfully loaded all line opacities
Loading CIA opacities for H2-H2 from file '/Users/molliere/Desktop/input_data_v3/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/molliere/Desktop/input_data_v3/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
Successfully loaded all opacities
```

Arguments like `temperature`

, `imposed_mass_fractions`

and `filling_species`

are **model parameters** (`SpectralModel.model_parameters`

). They are used as inputs to the `SpectralModel`

**built-in functions**. Note for example that we did not give the mean molar mass of the atmosphere. This is because `SpectralModel`

used a built-in function to calculate it from the mass fractions.

**Model parameters:** the list of available default `SpectralModel`

model parameters can be accessed by calling the `get_default_parameters`

function of a `SpectralModel`

instance, e.g. `spectral_model.get_default_parameters()`

.

There is no documentation yet for the meaning of these parameters.

**Argument ``filling_species``:** this argument is used to automatically fill the atmosphere with the given species, with the requested weights, without changing the `imposed_mass_fractions`

, and such that the sum of mass fractions is equal to 1. The filling species can be modified at will (e.g. `{'N2': 78, 'O2': 21, 'CO2': 0.5, 'H2': 0.5, 'Xe': 25}`

, priority is always given to the `imposed_mass_fractions`

). The weights are proportional to the relative mass fractions between the filling
species. They do not have to be normalized.

Here, we set the H2/He ratio to 37/12, which is approximately the ratio found in Jupiter. Alternatively, `filling_species`

can be set to `None`

(its default value), and the `H2`

and `He`

mass fractions can be added manually to `imposed_mass_fractions`

.

Let’s calculate the spectrum.

```
[5]:
```

```
wavelengths, transit_radii = spectral_model.calculate_spectrum(
mode='transmission',
update_parameters=True # this will build notably the temperature and mass fractions profile
)
```

Let’s now plot the transit radius

```
[6]:
```

```
fig, ax = plt.subplots(figsize = (10,6))
ax.plot(wavelengths[0] * 1e4, transit_radii[0] / cst.r_jup_mean)
ax.set_xscale('log')
ax.set_xlabel('Wavelength [microns]')
ax.set_ylabel(r'Transit radius [$\rm R_{Jup}$]')
```

```
[6]:
```

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

### Calculating an emission spectrum: a more complex example

Emission spectra usually require more complex temperature profiles and mass fraction profiles than transmission spectra to accurately fit the data. Here we will introduce another default `SpectralModel`

temperature profile (the Guillot 2010 temperature profile introduced in the Getting Started section), as well as `SpectralModel`

’s behaviour with mass fractions. We will also introduce more model parameters, relevant for emission
spectra modelling.

In order to show the default behaviour of `SpectralModel`

regarding mass fractions, we will use equilibrium chemistry and filling species, and manually set an unphysical \(\rm H_2O\) mass fraction profile (with abundances > 1).

#### Flux

Let’s first generate our `SpectralModel`

.

```
[7]:
```

```
spectral_model = SpectralModel(
# Radtrans parameters
pressures=np.logspace(-6, 2, 100),
line_species=[
'H2O',
'CO-NatAbund',
'CH4',
'CO2',
'Na',
'K'
],
rayleigh_species=['H2', 'He'],
gas_continuum_contributors=['H2-H2', 'H2-He'],
wavelength_boundaries=[0.3, 15],
scattering_in_emission=True, # replace do_scat_emis from pRT2
# SpectralModel parameters
# Planet parameters
planet_radius=1 * cst.r_jup_mean,
reference_gravity=10 ** 3.5,
reference_pressure=1e-2,
# Star, system, orbit
is_observed=True, # return the flux observed at system_distance
is_around_star=True, # if True, calculate a PHOENIX stellar spectrum and add it to the emission spectrum
system_distance=10 * cst.s_cst.light_year * 1e2, # m to cm, used to scale the spectrum
star_effective_temperature=5500, # used to get the PHOENIX stellar spectrum model
star_radius=1 * cst.r_sun, # used to get the star flux irradiating the planet
orbit_semi_major_axis=0.01 * cst.au, # used to get the star flux irradiating the planet
# Temperature profile parameters
temperature_profile_mode='guillot',
temperature=1500,
intrinsic_temperature=200,
guillot_temperature_profile_gamma=0.4,
guillot_temperature_profile_kappa_ir_z0=0.01,
# Mass fractions
use_equilibrium_chemistry=True,
metallicity=3, # times solar
co_ratio=0.1,
imposed_mass_fractions={
'H2O': np.logspace(1, -12, 100) # we use a H2O mass fraction > 1 to demonstrate how SpectralModel deals with it
},
filling_species={ # automatically fill the atmosphere with H2 and He, such that the sum of MMRs is equal to 1 and H2/He = 37/12, H2/Ne = 37/0.06
'H2': 37,
'He': 12,
'Ne': 0.06
}
)
```

```
Loading Radtrans opacities...
Loading line opacities of species 'H2O' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/H2O/1H2-16O/1H2-16O__POKAZATEL.R1000_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO-NatAbund' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CO/C-O-NatAbund/C-O-NatAbund__HITEMP.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CH4' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CH4/12C-1H4/12C-1H4__HITEMP.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO2' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CO2/12C-16O2/12C-16O2__UCL-4000.R1000_0.3-50mu.ktable.petitRADTRANS.h5.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'Na' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/Na/23Na/23Na__Allard.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'K' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/K/39K/39K__Allard.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Successfully loaded all line opacities
Loading CIA opacities for H2-H2 from file '/Users/molliere/Desktop/input_data_v3/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/molliere/Desktop/input_data_v3/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
Successfully loaded all opacities
```

There are 2 parameters in particular that may require further explanation: - `is_observed`

: if `True`

, return `flux * (planet_radius / system_distance) ** 2`

instead of just `flux`

. This is only used in emission mode. - `is_around_star`

: if `True`

, calculate a PHOENIX stellar spectrum and use it to compute the emission spectrum (which will then contain a contribution of the scattered stellar light if scattering is turned on). This is used only in emission mode. This has the same
effect as adding stellar parameters to `calculate_flux()`

when using a `Radtrans`

object directly, see “Scattering for Emission Spectra”.

Let’s now calculate the spectrum.

```
[8]:
```

```
wavelengths, flux = spectral_model.calculate_spectrum(
mode='emission',
update_parameters=True # this will build notably the temperature and mass fractions profile
)
```

```
Loading chemical equilibrium chemistry table from file '/Users/molliere/Desktop/input_data_v3/input_data/pre_calculated_chemistry/equilibrium_chemistry/equilibrium_chemistry.chemtable.petitRADTRANS.h5'... Done.
Loading PHOENIX star table in file '/Users/molliere/Desktop/input_data_v3/input_data/stellar_spectra/phoenix/phoenix.startable.petitRADTRANS.h5'... Done.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
```

Now let’s plot the spectrum. **Note that it will look very weird because of the unphysical** \(\rm H_2O\) **abundance we specified.**

```
[9]:
```

```
fig, ax = plt.subplots(figsize = (10,6))
ax.plot(wavelengths[0] * 1e4, flux[0])
ax.set_xscale('log')
ax.set_xlabel('Wavelength [microns]')
ax.set_ylabel(r'Observerved planet flux, $F_{\lambda}$ [erg cm$^{-2}$ s$^{-1}$ cm$^{-1}$]')
```

```
[9]:
```

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

#### Temperature profile

The temperature profile calculation is controlled by the function `compute_temperature_profile`

.

We can plot the temperature profile. This is the same Guillot temperature profile than in the getting started section.

```
[10]:
```

```
fig, ax = plt.subplots(figsize = (10,6))
ax.plot(spectral_model.temperatures, spectral_model.pressures * 1e-6)
ax.set_yscale('log')
ax.set_ylim([1e2, 1e-6])
ax.set_xlabel('T [K]')
ax.set_ylabel('P [bar]')
```

```
[10]:
```

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

#### Mass fractions

The mass fractions calculation is controlled by the function `compute_mass_fractions`

.

To calculate the mean molar masses, the function used is `compute_mean_molar_masses`

.

We can plot the mass fractions and see `SpectralModel`

’s default behaviour.

```
[11]:
```

```
fig, ax = plt.subplots(figsize = (10,6))
for species, mass_fraction in spectral_model.mass_fractions.items():
if species in spectral_model.line_species:
ax.loglog(mass_fraction, spectral_model.pressures * 1e-6, label=species)
for species, mass_fraction in spectral_model.mass_fractions.items():
if species in spectral_model.model_parameters['filling_species']:
ax.loglog(mass_fraction, spectral_model.pressures * 1e-6, label=species, ls=':')
ax.loglog(np.sum(list(spectral_model.mass_fractions.values()), axis=0), spectral_model.pressures * 1e-6, label=r'$\sum$ MMR', color='k')
ax.set_ylim([1e2, 1e-6])
ax.set_xlabel('MMR')
ax.set_ylabel('P [bar]')
ax.legend()
```

```
[11]:
```

```
<matplotlib.legend.Legend at 0x125e23b90>
```

Here we can see the priorities of `SpectralModel`

, by default: 1. Always ensure that the sum of MMR is 1. 2. Ensure that the imposed mass fractions are strictly respected, unless their sum is > 1. 3. Ensure that the equilibrium chemistry mass fractions are respected, unless they conflict with the imposed mass fractions. 4. Ensure that the filling species ratios are respected, unless they conflict with the above.

While \(\rm H_2O\) is by default in the `PreCalculatedEquilibriumChemistryTable`

(that we requested with `use_equilibrium_chemistry=True`

), we overrode its mass fraction by using `imposed_mass_fractions`

.

The species Ne is not in the `PreCalculatedEquilibriumChemistryTable`

, in contrast with H2 and He. Hence, where \(\rm H_2O\) has become too abundant, Ne has been removed from the atmosphere to maintain the equilibrium chemistry relative mass fractions.

Note that while the key `'CO-NatAbund'`

is not in the `PreCalculatedEquilibriumChemistryTable`

, `SpectralModel`

automatically linked it to the `'CO'`

key.

We can also plot the mean molar masses (MMW).

```
[12]:
```

```
fig, ax = plt.subplots(figsize = (10,6))
ax.semilogy(spectral_model.mean_molar_masses, spectral_model.pressures * 1e-6, label=species)
ax.set_ylim([1e2, 1e-6])
ax.set_xlabel('Mean molar masses')
ax.set_ylabel('P [bar]')
```

```
[12]:
```

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

Starting with the expected MMW for a jupiter-like \(\rm H_2\)/He atmosphere, as \(\rm H_2O\) becomes the dominant species, the MMW increases until it takes on the \(\rm H_2O\) molar mass value where the \(\rm H_2O\) mass fraction is 1.

## Calculating a time-varying high-resolution spectrum

Let’s now introduce more advanced `SpectralModel`

features: the spectral modification parameters. We will do this by going through the calculation of a time-varying high-resolution spectrum.

### Introduction: Doppler-shifting and relative velocity

`SpectralModel`

can generate time-varying models. This includes most importantly changes in the relative velocity between the observer and the planet.

During the observation, the planet rotates around the star, so its radial velocity relative to the observer will change, Doppler-shifting the spectral lines. At low resolution this effect is usually negligible as the spectral lines are not resolved. At high resolution, taking this effect into account is crucial.

But that’s not all. We also need to take the motion of the star relative to the observer into account. In `SpectralModel`

, this parameter is called `system_observer_radial_velocities`

. It is usually divided into 2 components:

The radial velocity between the planet’s system barycenter and the barycenter of the solar system (

`star_radial_velocity`

, \(V_\textrm{sys}\))The relative velocity between the observer and the barycenter of the solar system (

`barycentric_velocities`

, \(V_\textrm{bary}\))

Often, we also add a correction term `rest_frame_velocity_shift`

(\(V_\textrm{rest}\)) to the total velocity. This can be used as a proxy to model the effect of atmospheric winds, for example.

The `relative_velocities`

(\(V\)) in `SpectalModel`

are calculated following this equation:

\begin{equation} V(t) = \sqrt{\frac{G M_\ast}{a_p}} \sin(i_p) \sin(2 \pi \phi(t)) + V_\textrm{sys}(t) + V_\textrm{bary}(t) + V_\textrm{rest} \end{equation}

The first term is the planet radial velocity around its star relative to the observer. It is composed of the planet radial velocity semi-amplitude, often called \(K_p\), multiplied by the \(\sin\) of the orbital longitude (i.e. \(2\pi\) time the **orbital phase**, \(\phi\)). In `SpectralModel`

, \(K_p\) can be given (`radial_velocity_semi_amplitude`

), or calculated using:

`star_mass`

(\(M_\ast\)),`orbit_semi_major_axis`

(\(a_p\)),`orbital_inclination`

(\(i_p\), optional, \(90^\circ\) by default).

### Loading a `Planet`

It can be convenient to use the Planet object to get the parameters we need. Note that this is not required.

```
[13]:
```

```
from petitRADTRANS.planet import Planet
planet = Planet.get('HD 189733 b')
```

### Initializing a time-varying transmission spectrum

We will generate mock data wavelengths and times below. We will assume that we observed 19 exposures. To keep things simple, we will assume that `star_radial_velocity`

is fixed, and `barycentric_velocities`

varies linearly. We will fix `rest_frame_velocity_shift`

to -5 km.s-1, and let `SpectralModel`

calculate \(K_p\) for us.

```
[14]:
```

```
from petitRADTRANS.math import resolving_space
n_exposures = 19
data_wavelengths = resolving_space(1.519, 1.522, 2e5) * 1e-4 # (cm) generate wavelengths at a constant resolving power
times = 0.85 * planet.transit_duration * (np.linspace(0, 1, n_exposures) - 0.5) # covering 85% of the transit duration
orbital_phases = times / planet.orbital_period
mid_transit_time = 0 # (s)
star_radial_velocity = planet.star_radial_velocity # (cm.s-1) V_sys
barycentric_velocities = np.linspace(-13.25e5, -13.55e5, times.size) # (cm.s-1) V_bary
# Uncertainties assuming a S/N of 2000
data_uncertainties = 5e-4 * np.ones((1, n_exposures, data_wavelengths.size))
```

Now let’s initalize our `SpectralModel`

. We will add an opaque cloud layer at 100 mbar. Our instrument will have a resolving power of \(R = 80\,000\).

**Argument ``wavelength_boundaries``:** if you provide parameters such as the wavelengths of your data (`rebinning_wavelengths`

), you no longer need to manually set `wavelength_boundaries`

: `SpectralModel`

will automatically calculate the optimal wavelengths boundaries for your data, taking into account the effect of time-varying Doppler-shift.

We will use `planet`

to get some of the required parameters.

We will introduce a few new model parameters, that will be described in more detail in the next sections.

```
[15]:
```

```
spectral_model = SpectralModel(
# Radtrans parameters
pressures=np.logspace(-6, 2, 100),
line_species=[
'CO-NatAbund',
'H2O'
],
rayleigh_species=['H2', 'He'],
gas_continuum_contributors=['H2-H2', 'H2-He'],
line_opacity_mode='lbl',
line_by_line_opacity_sampling=4,
# SpectralModel parameters
# Planet parameters
planet_radius=planet.radius,
reference_gravity=planet.reference_gravity,
reference_pressure=1e-2,
star_radius=planet.star_radius,
transit_duration=planet.transit_duration,
orbital_period=planet.orbital_period,
# Velocity paramters
star_mass=planet.star_mass,
orbit_semi_major_axis=planet.orbit_semi_major_axis,
orbital_inclination=planet.orbital_inclination,
rest_frame_velocity_shift=-5e5, # (cm.s-1) V_rest
system_observer_radial_velocities=star_radial_velocity - barycentric_velocities,
# Temperature profile parameters
temperature_profile_mode='isothermal',
temperature=planet.equilibrium_temperature,
# Cloud parameters
opaque_cloud_top_pressure=1e-2, # added opaque cloud
# Mass fractions
use_equilibrium_chemistry=False,
imposed_mass_fractions={
'CO-NatAbund': 1e-2,
'H2O': 1e-3,
},
filling_species={
'H2': 37,
'He': 12
},
# Observation parameters
rebinned_wavelengths=data_wavelengths, # (cm) used for the rebinning, and also to set the wavelengths boundaries
rebin_range_margin_power=4, # used to set the wavelengths boundaries, adding a margin of ~1 Angstrom (1e-4 * ~1 µm)
convolve_resolving_power=8e4, # used for the convolution
mid_transit_time=mid_transit_time,
times=times,
# Preparation parameters
tellurics_mask_threshold=0.8, # mask the fitted transmittances if it is below this value
polynomial_fit_degree=2, # degree of the polynomial fit
uncertainties=data_uncertainties
)
```

```
Loading Radtrans opacities...
Loading line opacities of species 'CO-NatAbund' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/line_by_line/CO/C-O-NatAbund/C-O-NatAbund__HITEMP.R1e6_0.3-28mu.xsec.petitRADTRANS.h5'... Done.
Loading line opacities of species 'H2O' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/line_by_line/H2O/1H2-16O/1H2-16O__HITEMP.R1e6_0.3-28mu.xsec.petitRADTRANS.h5'... Done.
Successfully loaded all line opacities
Loading CIA opacities for H2-H2 from file '/Users/molliere/Desktop/input_data_v3/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/molliere/Desktop/input_data_v3/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
Successfully loaded all opacities
```

To generate an exploitable time-varying model, we need to do several operations, in that order:

**Scale**the transit radii to the planet’s star radius.Doppler-

**shift**the spectra at each exposure, in order to take the relative velocity (that needs to be calculated as well) between the planet and the observer into account.Take into account the

**transit light loss**due to the planet’s ingress and egress.**Convolve**the spectrum to the instrument’s LSF.**Rebin**the spectrum to the instrument’s wavelengths.

With `SpectralModel`

, these steps can be done easily, by adding some arguments to the `calculate_spectrum`

function:

```
[16]:
```

```
wavelengths_rebinned, transit_radii = spectral_model.calculate_spectrum(
mode='transmission',
update_parameters=True,
scale=True, # scale the spectrum
shift=True, # Doppler-shift the spectrum to the planet's radial velocities relative to the observer
use_transit_light_loss=True, # apply the effect of transit ingress and egress
convolve=True, # convolve to the instrument's resolving power
rebin=True # rebin to the instrument's wavelengths
)
```

Let’s take a look at the outputs dimensions.

```
[17]:
```

```
print(f"The shape of the rebinned wavelengths is {wavelengths_rebinned.shape}")
print(f"The shape of the rebinned transit_radii is {transit_radii.shape}")
```

```
The shape of the rebinned wavelengths is (1, 395)
The shape of the rebinned transit_radii is (1, 19, 395)
```

The `transit_radii`

has one dimension to store respectively, the spectral **orders**, the **exposures** and the **wavelengths**. The wavelengths `wavelengths_rebinned`

has one dimension to store respectively, the spectral **orders**, the **wavelengths**. In general, the instrument’s wavelengths of reduced data is the same across exposures.

The `data_wavelengths`

had 1 dimension, so `SpectralModel`

assumed that we had 1 order. To add more orders, `data_wavelengths`

should have 2 dimensions, the first one corresponding to the spectral orders, and the second to the wavelengths.

Let’s now plot the spectrum. Since we have one order, we can plot `transit_radii[0]`

, which varies along exposures and wavelengths.

```
[18]:
```

```
fig, ax = plt.subplots(figsize = (10,6))
ax.pcolormesh(wavelengths_rebinned[0] * 1e4, orbital_phases, transit_radii[0])
ax.set_xlabel('Wavelength [microns]')
ax.set_ylabel('Orbital phase')
```

```
[18]:
```

```
Text(0, 0.5, 'Orbital phase')
```

Note that you can use the modifications arguments (`convolve`

, `rebin`

, etc.) in any combination (e.g., only `rebin`

, only `scale`

and `convolve`

, etc.)

### Spectral modifications in details

In the sections below we will explain the effect of each new `calculate_spectrum`

argument more in-depth.

#### No modification

If we calculated the transit radii as we did in the intermediate usage section, we would generate only a single spectrum for one exposure. This is triggered by `shift = False`

, which is the default value. Hence, only a single spectrum is returned.

```
[19]:
```

```
wavelengths, transit_radii_not_modified = spectral_model.calculate_spectrum(
mode='transmission',
update_parameters=True
)
# Plot the spectrum
fig, ax = plt.subplots(figsize = (10,6))
ax.plot(wavelengths[0] * 1e4, transit_radii_not_modified[0]) # [0] is used to get the 1st and only exposure
ax.set_xlabel('Wavelength [microns]')
ax.set_ylabel(r'Transit radius [cm]')
```

```
[19]:
```

```
Text(0, 0.5, 'Transit radius [cm]')
```

#### Scaling

The scaling is controlled by the function `scale_spectrum`

.

Here “scaling” refer to modifying the spectrum such that it is expressed relative to its planet’s star’s spectum.

For transiting planets, the data collected are a combination of the star’s spectrum (\(F_\ast\)) and the planet’s spectrum (\(F_p\)).

For emission spectra, the observed spectrum before or after the (“secondary”) eclipse (\(F_\mathrm{obs}\)) can be expressed as: \begin{equation}
F_{\mathrm{obs}} = F_\ast + F_p.
\end{equation} The stellar spectrum can be acquired independently during the planet’s actual (“secondary”) eclipse, to obtain the **scaled** spectrum (after subtracting 1): \begin{equation}
F_{\mathrm{scaled}} = F_{\mathrm{obs}}/F_\ast -1 = F_p / F_\ast.
\end{equation} In `SpectralModel`

, the scaled emission spectrum is:

```
spectrum /= star_observed_spectrum
```

where `star_observed_spectrum`

is the re-binned spectrum obtained from the `compute_star_flux`

function (PHOENIX/ATLAS9 model by default), observed at `system_distance`

assuming a star radius of `star_radius`

.

In transmission (“primary” eclipse), the observed spectrum during transit (\(F_\mathrm{obs}\)) can be expressed as a function of the planet (spectral) radius (\(R_p\), aka the “transit radius”) and the star radius (\(R_\ast\)): \begin{equation}
F_{\mathrm{obs}} = F_\ast (1 - R_p / R_\ast)^2.
\end{equation} Similarly, the **scaled** spectrum can be obtained by dividing with the star spectrum: \begin{equation}
F_{\mathrm{scaled}} = 1 - (R_p / R_\ast)^2.
\end{equation} In `SpectralModel`

, the scaled transmission spectrum is:

```
spectrum = 1 - (spectrum / star_radius) ** 2
```

since spectrum is the planet’s transit radius.

To activate the scaling, set `scale=True`

:

```
[20]:
```

```
wavelengths, transit_radii_scaled = spectral_model.calculate_spectrum(
mode='transmission',
update_parameters=True,
scale=True
)
fig, ax = plt.subplots(figsize = (10,6))
ax.plot(wavelengths[0] * 1e4, transit_radii_scaled[0])
ax.set_xlabel('Wavelength [microns]')
ax.set_ylabel(r'Scaled transit spectrum [arbitrary units]')
```

```
[20]:
```

```
Text(0, 0.5, 'Scaled transit spectrum [arbitrary units]')
```

#### Shifting

The shifting is controlled by the function `shift_wavelengths`

.

Now let’s add Doppler-shifts the spectrum add the time dimension to our model. The explantions behind this step are in this section.

To activate the shifting, set `shift=True`

:

```
[21]:
```

```
wavelengths_shifted, transit_radii_shifted = spectral_model.calculate_spectrum(
mode='transmission', # can also be 'emission' to generate an emission spectrum, in that case, other model_parameters need to be provided
update_parameters=True, # the parameters that we set will be properly initialised
scale=True,
shift=True
)
```

The `transit_radii`

and `wavelengths`

shapes has been modified.

```
[22]:
```

```
print(f"The shape of the shifted wavelengths is {wavelengths_shifted.shape}")
print(f"The shape of the shifted transit_radii is {transit_radii_shifted.shape}")
```

```
The shape of the shifted wavelengths is (19, 552)
The shape of the shifted transit_radii is (19, 552)
```

Now, instead of assuming 1 exposure as before, we have 1 spectrum for each of our 19 exposures, each with a spectrum Doppler-shifted at a different **relative velocity**. The relative velocities has been calculated by the built-in function `compute_relative_velocities`

.

Now let’s plot the relative velocities `SpectralModel`

calculated.

```
[23]:
```

```
fig, ax = plt.subplots(figsize = (10,6))
ax.plot(
orbital_phases,
spectral_model.model_parameters['relative_velocities'] * 1e-5, # (cm.s-1 to km.s-1)
label='relative velocity'
)
# Plot visual indicators
ax.plot(
orbital_phases,
spectral_model.model_parameters['system_observer_radial_velocities'] * 1e-5, # cm.s-1 to km.s-1
ls='--', color='k', label='system-observer radial velocity'
)
ax.plot(
orbital_phases,
(spectral_model.model_parameters['system_observer_radial_velocities'] + spectral_model.model_parameters['rest_frame_velocity_shift']) * 1e-5,
ls='--', color='r', label='planet rest frame velocity'
)
y_lim = ax.get_ylim()
ax.vlines(0, y_lim[0], y_lim[1], color='grey', ls=':')
ax.set_xlabel('Orbital phase')
ax.set_ylabel(r'Velocity [km$\cdot$s$^{-1}$]')
ax.set_ylim(y_lim)
ax.legend()
```

```
[23]:
```

```
<matplotlib.legend.Legend at 0x1264f3dd0>
```

Now let’s plot the spectra. We will plot 3 of them for clarity.

```
[24]:
```

```
# Select 3 exposures to plot
exposures_to_plot = [
0,
10,
-1
]
# Set the 3 spectrum colors in the plot
colors = [
'b',
'g',
'r'
]
# Plot the un-shifted spectrum
fig, ax = plt.subplots(figsize = (10,6))
ax.plot(
wavelengths[0] * 1e4, transit_radii_scaled[0],
label='scaled spectrum', color='k', ls=':'
)
# Plot the shifted spectra
for i, exposure in enumerate(exposures_to_plot):
ax.plot(
wavelengths_shifted[exposure] * 1e4, transit_radii_shifted[exposure],
label=f"shifted spectrum (v = {spectral_model.model_parameters['relative_velocities'][exposure] * 1e-5:.2f} km.s-1)", color=colors[i]
)
ax.set_xlabel('Wavelength [microns]')
ax.set_ylabel(r'Scaled transit spectrum [arbitrary units]')
ax.legend()
```

```
[24]:
```

```
<matplotlib.legend.Legend at 0x1265974d0>
```

#### Convolving

The convolution is controlled by the function `convolve`

.

Let’s now add a convolution to the spectral modification, to better represent the instrument’s line-spread function (LSF). The LSF is the function you’d measure with the spectrograph with inifitesimally small pixels if the spectrum contained only one line of zero width (i.e., a \(\delta\) function). The LSF thus caracterizes how the spectrograph disperses light as a function of wavelength. By default, `SpectralModel`

uses a Gaussian kernel and the parameter `convolve_resolving_power`

that we provided. `convolve_resolving_power`

is defined as \(\lambda/\Delta\lambda\), where \(\lambda\) is the wavelength, and \(\Delta\lambda\) the width of the spectrograph’s resolution element. The built-in function for this operation is `convolve`

.

To activate the shifting, set `convolve=True`

:

```
[25]:
```

```
wavelengths_shifted, transit_radii_shifted_convolved = spectral_model.calculate_spectrum(
mode='transmission', # can also be 'emission' to generate an emission spectrum, in that case, other model_parameters need to be provided
update_parameters=True, # the parameters that we set will be properly initialised
scale=True,
shift=True,
convolve=True
)
```

Let’s plot the spectrum at one exposure.

```
[26]:
```

```
exposure_to_plot = 10
# Plot the spectrum
fig, ax = plt.subplots(figsize = (10,6))
ax.plot(
wavelengths_shifted[exposure_to_plot] * 1e4, transit_radii_shifted[exposure_to_plot],
label=f"shifted spectrum (v = {spectral_model.model_parameters['relative_velocities'][exposure_to_plot] * 1e-5:.2f} km.s-1)", color='k', ls=':'
)
ax.plot(
wavelengths_shifted[exposure_to_plot] * 1e4, transit_radii_shifted_convolved[exposure_to_plot],
label=f"shifted convolved spectrum (v = {spectral_model.model_parameters['relative_velocities'][exposure_to_plot] * 1e-5:.2f} km.s-1)"
)
ax.set_xlabel('Wavelength [microns]')
ax.set_ylabel(r'Scaled transit spectrum [arbitrary units]')
ax.legend()
```

```
[26]:
```

```
<matplotlib.legend.Legend at 0x1266d2b10>
```

#### Re-binning

The re-binning is controlled by the function `rebin_spectrum`

.

Let’s now re-bin the spectrum to the data wavelengths (which are ultimately defined by the pixel position and width on the detector) in order to be able to compare our model with the data properly.

To activate the shifting, set `rebin=True`

:

```
[27]:
```

```
wavelengths_rebinned, transit_radii_rebinned = spectral_model.calculate_spectrum(
mode='transmission', # can also be 'emission' to generate an emission spectrum, in that case, other model_parameters need to be provided
update_parameters=True, # the parameters that we set will be properly initialised
scale=True,
shift=True,
convolve=True,
rebin=True
)
```

The shapes of `transit_radii`

and `wavelengths`

have changed again.

```
[28]:
```

```
print(f"The shape of the rebinned wavelengths is {wavelengths_rebinned.shape}")
print(f"The shape of the rebinned transit radii is {transit_radii_rebinned.shape}")
```

```
The shape of the rebinned wavelengths is (1, 395)
The shape of the rebinned transit radii is (1, 19, 395)
```

Let’s plot the spectra.

```
[29]:
```

```
fig, ax = plt.subplots(figsize = (10,6))
ax.pcolormesh(wavelengths_rebinned[0] * 1e4, orbital_phases, transit_radii_rebinned[0])
ax.set_xlabel('Wavelength [microns]')
ax.set_ylabel('Orbital phase')
```

```
[29]:
```

```
Text(0, 0.5, 'Orbital phase')
```

#### Transit effect

The transit effect is controlled by the function `compute_transit_fractional_light_loss`

.

During a primary transit, the planet’s disk is not always fully superimposed on the star’s disk. This means that the transit depth during ingress and egress is smaller than during the full transit. In `SpectralModel`

, this is done with the built-in function `compute_transit_fractional_light_loss`

, which is based on Mandel & Agol (2002) (Eq. 1), and uses the model parameters `transit_duration`

(\(T_{14}\)) and `orbital_period`

(\(P\))..

The limb-darkening effect is currently **not** taken into account.

We can calculate the spectrum to obtain the same result as at the start of this section.

```
[30]:
```

```
wavelengths_rebinned, transit_radii_final = spectral_model.calculate_spectrum(
mode='transmission', # can also be 'emission' to generate an emission spectrum, in that case, other model_parameters need to be provided
update_parameters=True, # the parameters that we set will be properly initialised
scale=True,
shift=True,
use_transit_light_loss=True,
convolve=True,
rebin=True
)
# Plot the spectra across times at one wavelength
fig, ax = plt.subplots(figsize = (10,6))
wavelength_id = int(transit_radii_final.shape[-1] / 2)
ax.plot(orbital_phases, transit_radii_final[0, :, wavelength_id])
ax.set_xlabel('Orbital phase')
ax.set_ylabel('Scaled spectrum [arbitrary unit]')
ax.set_title(f"Spectrum at {wavelengths_rebinned[0, wavelength_id] * 1e4:.5f} microns")
```

```
[30]:
```

```
Text(0.5, 1.0, 'Spectrum at 1.52050 microns')
```

### Further modifications

This section goes into additional modification capabilities of `SpectralModel`

, that are useful for tests or data analysis.

#### Simulating data

To build simulated data, we can add telluric transmittances, setllar lines, instrumental deformations, and noise. To do this, we can make use of the following arguments: - `telluric_transmittances`

: telluric transmittances to combine with the spectrum. Those can make use of the `airmass`

model parameter to automatically build time- and wavelength-varying transmittances. A good place to download telluric transmittances is the
SKYCALC website. petitRADTRANS also has a module to directly download SKYCALC data, calculate the airmass, and even find the best day for your observations in the module `petitRADTRANS.cli.eso_skycalc_cli`

. - `telluric_transmittances_wavelengths`

: wavelengths of the provided telluric transmittances. - `instrumental_deformations`

: time and wavelength matrix of the instrumental deformations. Must be of
shape `(n_orders, n_exposures, n_wavelengths)`

. - `noise_matrix`

: time and wavelength matrix of noise. Must be of shape `(n_orders, n_exposures, n_wavelengths)`

.

The detailed `SpectralModel`

modification flowchart in that case is: 1. `scale`

, `shift`

, then add the `transit_light_loss`

, 2. Rebin the `telluric_transmittances`

to the current wavelengths, then multiply them to the spectra, 3. `convolve`

, then `rebin`

the spectra, 4. Multiply the spectra with the `instrumental_deformations`

, 5. Add the `noise_matrix`

to the spectra.

Below is an example, using simplistic `telluric_transmittances`

:

```
[31]:
```

```
# High resolution telluric transmittances (typically downloaded from SKYCALC)
# If we add airmass information, as we do below, the transmittance here needs to be given at airmass == 1.
telluric_transmittances_wavelengths = resolving_space(1.515, 1.525, 1e6) * 1e-4
telluric_transmittances = np.ones(telluric_transmittances_wavelengths.size)
# Add simplisitc lines
telluric_transmittances[2850:2870] = 0.85
telluric_transmittances[2900:2930] = 0.7
telluric_transmittances[3400:3500] = 0.1
telluric_transmittances[3900:3950] = 0.5
telluric_transmittances[4400:4475] = 0.3
# Get airmass
mid_transit_time, _, _, _ = planet.calculate_mid_transit_time(
observation_day=2458004.424877, # (JD)
day2second=False # if True, return the result in seconds since the start of the transit's day
)
airmass = planet.calculate_airmass(
time=2458004.424877 + times / cst.s_cst.day,
site_name='CAHA',
time_format='jd'
)
spectral_model.model_parameters['airmass'] = airmass
# Random wavelength-constant instrumental deformations (typically unknown on real data), use a seed for reprooducibility
instrumental_deformations = (0.4 - 0.2) * np.random.default_rng(seed=12345).random(n_exposures) + 0.2
instrumental_deformations = np.repeat(instrumental_deformations[:, np.newaxis], data_wavelengths.size, axis=-1)
# Noise assuming a S/N of 1000
data_uncertainties = 5e-4 * np.ones((1, n_exposures, data_wavelengths.size))
noise_matrix = np.random.default_rng(seed=54321).normal(loc=0, scale=data_uncertainties)
noise_factor = 20 # for visual purposes
# Simulated data
wavelengths_rebinned, noisy_data = spectral_model.calculate_spectrum(
mode='transmission',
update_parameters=True,
telluric_transmittances_wavelengths=telluric_transmittances_wavelengths,
telluric_transmittances=telluric_transmittances,
instrumental_deformations=instrumental_deformations,
noise_matrix=noise_matrix * noise_factor, # increases the noise to make it more visible in the figure
scale=True,
shift=True,
use_transit_light_loss=True,
convolve=True,
rebin=True
)
# Plot the spectra
fig, ax = plt.subplots(figsize = (10,6))
ax.pcolormesh(wavelengths_rebinned[0] * 1e4, orbital_phases, noisy_data[0])
ax.set_xlabel('Wavelength [microns]')
ax.set_ylabel('Orbital phase')
ax.set_title(f"Simulated noisy data (noise time {noise_factor})")
```

```
[31]:
```

```
Text(0.5, 1.0, 'Simulated noisy data (noise time 20)')
```

Here you can see that the planet’s spectral lines have become completely invisible, as we are dominated by the telluric lines. Note that the lines become slightly deeper as the orbital phase increase, this is the effect of the `airmass`

we added above. Also note that in the case here, where we added the airmass specifically, the transmittance needs to be given at airmass = 1, and will then be scaled by `SpectralModel`

.

#### Preparing

The preparation step is controlled by the function `preparing_pipeline`

.

The default petitRADTRANS preparing pipeline “Polyfit” uses polynomial fits on wavelengths (for each exposure), then on times (for each wavelength) to remove the `instrumental_deformations`

and `telluric_transmittances`

. See Blain et al. (2024) for more information.

We will build the same spectrum without noise so we can see the result on the lines, and use the preparing pipeline on it. This will simulate noiseless prepared data.

```
[32]:
```

```
from petitRADTRANS.retrieval.preparing import polyfit
# Simulated noiseless data (instead you would load the data here)
data_wavelengths, data = spectral_model.calculate_spectrum(
mode='transmission', # can also be 'emission' to generate an emission spectrum, in that case, other model_parameters need to be provided
update_parameters=True, # the parameters that we set will be properly initialised
telluric_transmittances_wavelengths=telluric_transmittances_wavelengths,
telluric_transmittances=telluric_transmittances,
instrumental_deformations=instrumental_deformations,
noise_matrix=None, # no noise so we can see the lines
scale=True,
shift=True,
use_transit_light_loss=True,
convolve=True,
rebin=True,
prepare=False # if the data were loaded, a preparation could not be done rigth away
)
# The data and its uncertainties must be masked arrays
data = np.ma.masked_array(data)
data_uncertainties = np.ma.masked_array(data_uncertainties)
# You can mask for example invalid pixels, here we will mask nothing
data.mask = np.zeros(data.shape, dtype=bool)
data_uncertainties.mask = np.zeros(data_uncertainties.shape, dtype=bool)
# Prepare the loaded data
prepared_data, preparation_matrix, prepared_data_uncertainties = polyfit(
spectrum=data,
uncertainties=data_uncertainties,
wavelengths=data_wavelengths,
airmass=airmass,
tellurics_mask_threshold=0.8,
polynomial_fit_degree=2,
full=True,
apply_throughput_removal=True,
apply_telluric_lines_removal=True
)
# Plot the prepared data
fig, ax = plt.subplots(figsize = (10,6))
ax.pcolormesh(data_wavelengths[0] * 1e4, orbital_phases, prepared_data[0])# Get simulated noiseless data (instead you would load the data here)
```

```
[32]:
```

```
<matplotlib.collections.QuadMesh at 0x12676a3d0>
```

We get the lines again, but they are deformed by the preparing pipeline.

To get a prepared forward model, you can use the `prepare`

argument of `calculate_spectrum`

:

```
[33]:
```

```
# Simulated prepared noiseless data
wavelengths_rebinned, prepared_model = spectral_model.calculate_spectrum(
mode='transmission', # can also be 'emission' to generate an emission spectrum, in that case, other model_parameters need to be provided
update_parameters=True, # the parameters that we set will be properly initialised
scale=True,
shift=True,
use_transit_light_loss=True,
convolve=True,
rebin=True,
prepare=True
)
# Plot the spectra
fig, ax = plt.subplots(figsize = (10,6))
ax.pcolormesh(wavelengths_rebinned[0] * 1e4, orbital_phases, prepared_model[0])
ax.set_xlabel('Wavelength [microns]')
ax.set_ylabel('Orbital phase')
ax.set_title('Prepared model')
```

```
[33]:
```

```
Text(0.5, 1.0, 'Prepared model')
```

In general, preparing pipelines performs better when the number of exposures and of wavelengths is higher, because the fit of the large trends in the data will be better.

To compare pipelines, the Bias Pipeline Metric (BPM, see Blain et al. 2024) can be used. It can be calculated via:

```
[34]:
```

```
from petitRADTRANS.retrieval.preparing import bias_pipeline_metric
bpm = bias_pipeline_metric(
prepared_true_model=prepared_model,
prepared_mock_observations=prepared_data,
mock_observations_preparation_matrix=preparation_matrix,
mock_noise=None # we are comparing with the noiseless simulated prepared data
)
print(f"BPM: {np.mean(np.abs(bpm)):.2e}")
```

```
BPM: 1.22e-06
```

The pipeline that has the mean absolute BPM the closest to 0 for a given set of simulated data can be seen as the most performant.

## Saving and loading

You can also save your `SpectralModel`

for reproducibility or latter usage.

**Caution:** currently `SpectralModel`

does not save custom functions: it is assumed that the default `SpectralModel`

has been used.

```
[35]:
```

```
# Save the model
spectral_model.save('my_model.h5')
```

If you want to reload your `SpectralModel`

with all of its `model_parameters`

, simply use:

```
[37]:
```

```
spectral_model_loaded = SpectralModel.load('my_model.h5')
```

```
Converting pressures from CGS to bar for Radtrans input...
Loading Radtrans opacities...
Loading line opacities of species 'CO-NatAbund' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/line_by_line/CO/C-O-NatAbund/C-O-NatAbund__HITEMP.R1e6_0.3-28mu.xsec.petitRADTRANS.h5'... Done.
Loading line opacities of species 'H2O' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/line_by_line/H2O/1H2-16O/1H2-16O__HITEMP.R1e6_0.3-28mu.xsec.petitRADTRANS.h5'... Done.
Successfully loaded all line opacities
Loading CIA opacities for H2-H2 from file '/Users/molliere/Desktop/input_data_v3/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/molliere/Desktop/input_data_v3/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
Successfully loaded all opacities
```

We can check that the loaded `SpectralModel`

is in the same state as the local one:

```
[38]:
```

```
loaded_wavelengths, loaded_model = spectral_model_loaded.calculate_spectrum(
**spectral_model_loaded.model_parameters['modification_parameters']
)
fig, ax = plt.subplots(figsize = (10,6))
ax.plot(loaded_wavelengths[0] * 1e4, loaded_model[0, 10], label='loaded_model')
ax.plot(wavelengths_rebinned[0] * 1e4, prepared_model[0, 10], label='saved_model')
ax.set_xlabel('Wavelength [microns]')
ax.set_ylabel('Transit spectrum (A.U.)')
```

```
[38]:
```

```
Text(0, 0.5, 'Transit spectrum (A.U.)')
```

## Custom functions

All built-in `SpectralModel`

functions can be customised.

Here is how to modify, for example the temperature profile function.

**All** the customised `SpectralModel`

functions **must** include a `**kwargs`

argument.

```
[39]:
```

```
import copy
# Copy the SpectralModel
custom_spectral_model = copy.deepcopy(spectral_model)
# Set a new temperature profile function based on some top-notch physics
def my_t_profile(pressures, my_parameter, x, **kwargs): # the **kwargs is important, even if it's not used
# Note that the "temperature" model parameter is not necessary, but in this context it could be used in place of "x"
return x * (1 + 0.3 * np.sin(np.deg2rad(my_parameter * np.log(pressures))))
# Add necessary model parameters for the new TP model
# The new parameters are added here, but can also be added during instanciation just like any other model parameter
custom_spectral_model.model_parameters['my_parameter'] = 30
custom_spectral_model.model_parameters['x'] = 1200
# Modify the temperature profile of the model
custom_spectral_model.compute_temperatures = my_t_profile
custom_spectral_model.update_model_functions_map(update_model_parameters=False)
# Get the model
wavelengths, transit_radii_new_tp = custom_spectral_model.calculate_spectrum(
mode='transmission',
update_parameters=True,
scale=True
)
# Plot the new temperature profile
fig, ax = plt.subplots(figsize = (6,6))
ax.semilogy(custom_spectral_model.temperatures, custom_spectral_model.pressures * 1e-6)
ax.set_xlabel('Temperature [K]')
ax.set_ylabel('Pressure [bar]')
ax.set_ylim(np.array([custom_spectral_model.pressures[0], custom_spectral_model.pressures[-1]]) * 1e-6)
ax.set_title('New temperature profile')
# Plot the new spectrum profile
fig, ax = plt.subplots(figsize = (10,6))
ax.plot(wavelengths[0], transit_radii_scaled[0], label='isothermal TP')
ax.plot(wavelengths[0], transit_radii_new_tp[0], label='new TP')
ax.set_xlabel('Wavelength [microns]')
ax.set_ylabel('Scaled spectrum [arbitrary units]')
ax.legend()
ax.set_title('New spectrum')
```

```
[39]:
```

```
Text(0.5, 1.0, 'New spectrum')
```