Utility functions
pRT contains some useful utility functions such as a spectral library, pre-implemented pressure-temperature profiles, etc. The use of some of them is shown below. pRT also comes with a package of physical constants, most of which are defined by importing astropy and scipy constants, however.
Constants
Can be accessed by importing from petitRADTRANS import physical_constants as cst
. All units are in cgs.
cst.c
: speed of lightcst.h
: Planck constantcst.kB
: Boltzman constantcst.nA
: Avogadro constantcst.e
: electron chargecst.G
: Gravitational constantcst.m_elec
: electron masscst.sigma
: Stefan-Boltzman constantcst.L0
: Loschmidt constantcst.R
: universal gas constantcst.bar
: 1 bar in cgscst.atm
: 1 atmosphere in cgscst.au
: Astronomical unitcst.pc
: parseccst.light_year
: light yearcst.amu
: atomic mass unit in gcst.r_sun
: solar radiuscst.r_jup
: Jupiter equatorial radiuscst.r_jup_mean
: Jupiter mean radiuscst.r_earth
: Earth radiuscst.m_sun
: Solar masscst.m_jup
: Jupiter masscst.m_earth
: Earth masscst.s_earth
: solar_irradiance
Planck function
The planck function \(B_\nu(T)\), in units of erg/cm\(^2\)/s/Hz/steradian, for a given frequency array, can be generated like this:
[1]:
from petitRADTRANS import physical_constants as cst
from petitRADTRANS import physics as phys
import numpy as np
# Define wavelength array, in cm
lamb = np.logspace(-5,-2,100)
# Convert to frequencies
freq = cst.c / lamb
# Calculate Planck function at 5750 K
planck = phys.planck_function_hz(5750., freq)
Let’s plot the Planck function:
[2]:
# Plot Planck function
import pylab as plt
plt.rcParams['figure.figsize'] = (10, 6)
plt.plot(lamb/1e-4, planck)
plt.xscale('log')
plt.xlabel('Wavelength (micron)')
plt.ylabel('Intensity (erg/cm$^2$/s/Hz/sterad)')
plt.show()
plt.clf()
<Figure size 1000x600 with 0 Axes>
PHOENIX and ATLAS9 stellar model spectra
Within petitRADTRANS the PHOENIX and ATLAS9 stellar spectra can be used, as described in Appendix A of van Boekel et al. (2012). The PHOENIX model refrence, for stellar effective temperatures < 10,000 K is Husser et al. (2013). The ATLAS9 model references for effective temperatures > 10,000 K are Kurucz (1979, 1992, 1994).
The models can be acessed like this, this is for a 5750 K effective temperature star on the main sequence:
[3]:
from petitRADTRANS.stellar_spectra.phoenix import PhoenixStarTable
star = PhoenixStarTable()
stellar_spec, __ = star.compute_spectrum(5750)
wlen_in_cm = stellar_spec[:,0]
flux_star = stellar_spec[:,1]
Loading PHOENIX star table in file '/Users/molliere/Desktop/input_data_v3/input_data/stellar_spectra/phoenix/phoenix.startable.petitRADTRANS.h5'... Done.
Let’s plot the spectrum, and also overplot the black body flux from the previous section (note the required factor of \(\pi\) to convert the black body intensity to flux):
[4]:
import pylab as plt
plt.plot(wlen_in_cm/1e-4, flux_star, label = 'PHOENIX model')
plt.plot(lamb/1e-4, np.pi*planck, label = 'black body flux')
plt.title(r'$T_{\rm eff}=5750$ K')
plt.xscale('log')
plt.xlabel('Wavelength (micron)')
plt.ylabel('Surface flux (erg/cm$^2$/s/Hz/sterad)')
plt.legend(loc = 'best', frameon = False)
plt.show()
plt.clf()
<Figure size 1000x600 with 0 Axes>
Guillot temperature model
In petitRADTRANS, one can use analytical atmospheric P-T 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}\) the atmospheric opacity in the IR wavelengths (i.e. the cross-section per unit mass), \(g\) the atmospheric surface gravity, \(\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.
Let’s define an example, all units are cgs units, except for the pressure, which is in bars:
[5]:
from petitRADTRANS.physics import temperature_profile_function_guillot_global
pressures_bar = np.logspace(-6, 2, 100)
reference_gravity = 10**3.5
kappa_IR = 0.01
gamma = 0.4
T_int = 200.
T_equ = 1500.
temperatures = temperature_profile_function_guillot_global(pressures_bar,
kappa_IR,
gamma,
reference_gravity,
T_int,
T_equ)
Let’s plot the P-T profile:
[6]:
plt.plot(temperatures, pressures_bar)
plt.yscale('log')
plt.ylim([1e2, 1e-6])
plt.xlabel('T (K)')
plt.ylabel('P (bar)')
plt.show()
plt.clf()
<Figure size 1000x600 with 0 Axes>