High-resolution spectra#
See also SpectralModel
(tutorial here and retrieval tutorial here), which comes with many useful features to work with high-resolution data.
Let’s set up a Radtrans
object as in the “Getting Started” example, but this time for high-resolution spectra:
[1]:
import matplotlib.pyplot as plt
import numpy as np
from petitRADTRANS.radtrans import Radtrans
atmosphere = Radtrans(
pressures=np.logspace(-10,2,130),
line_species=[
'1H2-16O',
'12C-16O'
],
rayleigh_species=['H2', 'He'],
gas_continuum_contributors=['H2-H2', 'H2-He'],
wavelength_boundaries=[2.2, 2.4],
line_opacity_mode='lbl'
)
Loading Radtrans opacities...
Loading line opacities of species '1H2-16O' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/line_by_line/H2O/1H2-16O/1H2-16O__HITRAN.R1e6_0.3-28mu.xsec.petitRADTRANS.h5'... Done.
Loading line opacities of species '12C-16O' from file '/home/dblain/petitRADTRANS/input_data/opacities/lines/line_by_line/CO/12C-16O/12C-16O__HITEMP.R1e6_0.3-28mu.xsec.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
Successfully loaded all opacities
Here we included only water and carbon monoxide, as only those play a role for the abundances and wavelength range chosen for our example here.
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.
Here we directly addressed which isotopologues we wanted to consider, when loading the opacities (see “Available opacity species” for a full list of available species). Note that we also could have written line_species = ['H2O', 'CO']
, in this case pRT loads the main isotopologues, which is also what we requested specifically above. In addition, we have invoked the high resolution mode by setting the keyword argument line_opacity_mode
to 'lbl'
, which
stands for line-by-line. We also loaded a smaller wavelength range, because the high-resolution opacities are pretty big…
We define the planetary radius and gravity at reference pressure \(P_0\), as well as the temperatures, abundances and mean molecular weight like before, see “Getting Started” for more information:
[2]:
from petitRADTRANS import physical_constants as cst
from petitRADTRANS.physics import temperature_profile_function_guillot_global
planet_radius = 1.0 * cst.r_jup_mean
reference_gravity = 10 ** 3.5
reference_pressure = 0.01
pressures = atmosphere.pressures*1e-6 # cgs to bar
infrared_mean_opacity = 0.01
gamma = 0.4
intrinsic_temperature = 200
equilibrium_temperature = 1500
temperatures = temperature_profile_function_guillot_global(
pressures=pressures,
infrared_mean_opacity=infrared_mean_opacity,
gamma=gamma,
gravities=reference_gravity,
intrinsic_temperature=intrinsic_temperature,
equilibrium_temperature=equilibrium_temperature
)
mass_fractions = {
'H2': 0.74 * np.ones_like(temperatures),
'He': 0.24 * np.ones_like(temperatures),
'1H2-16O': 1e-3 * np.ones_like(temperatures),
'12C-16O': 1e-2 * np.ones_like(temperatures)
}
mean_molar_masses = 2.33 * np.ones_like(temperatures) # 2.33 is a typical value for H2-He dominated atmospheres
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()
.
Transmission spectrum#
Now, let’s calculate and plot the transmission spectrum:
[3]:
wavelengths, transit_radii, _ = atmosphere.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
)
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(wavelengths * 1e4, transit_radii / cst.r_jup_mean)
ax.set_xlabel('Wavelength [micron]')
ax.set_ylabel(r'Transit radius [$\rm R_{Jup}$]')
[3]:
Text(0, 0.5, 'Transit radius [$\\rm R_{Jup}$]')
Let’s zoom-in a bit, to see individual lines
[4]:
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(wavelengths * 1e4, transit_radii / cst.r_jup_mean)
ax.set_xlabel('Wavelength [micron]')
ax.set_ylabel(r'Transit radius [$\rm R_{Jup}$]')
ax.set_xlim([2.3,2.3025])
[4]:
(2.3, 2.3025)
Emission Spectrum#
As before, the flux can be calculated like this
[5]:
frequencies, flux, _ = atmosphere.calculate_flux(
temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses = mean_molar_masses,
reference_gravity = reference_gravity,
frequencies_to_wavelengths=False
)
fig, ax = plt.subplots(figsize=(10,6))
ax.plot(cst.c / frequencies * 1e4, flux)
ax.set_xlabel('Wavelength [micron]')
ax.set_ylabel(r'Planet flux $F_\nu$ (10$^{-6}$ erg cm$^{-2}$ s$^{-1}$ Hz$^{-1}$)')
[5]:
Text(0, 0.5, 'Planet flux $F_\\nu$ (10$^{-6}$ erg cm$^{-2}$ s$^{-1}$ Hz$^{-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()
. This is also what we did in the example above.
Scattering and petitRADTRANS: remember that 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 how to do this). We neglect the scattering here.
Zooming in here as well:
[6]:
fig, ax = plt.subplots(figsize=(10,6))
ax.plot(cst.c/frequencies/1e-4, flux)
ax.set_xlabel('Wavelength [micron]')
ax.set_ylabel(r'Planet flux $F_\nu$ (10$^{-6}$ erg cm$^{-2}$ s$^{-1}$ Hz$^{-1}$)')
ax.set_xlim([2.3,2.3025])
[6]:
(2.3, 2.3025)