Rocky Planets with Reflecting and Emitting Surfaces
This tutorial covers how to make forward model thermal emission + reflection spectra for a terrestrial exoplanet with a reflective and emitting surface.
We both demonstrate how to generate models for rocky exoplanets with overlaying atmospheres, and rocky exoplanets without atmospheres.
For a more detailed tutorial on how to model reflection spectra of Earth-like planets, see the ‘Simulating Habitable Worlds Observatory Reflection Spectra with POSEIDON’ tutorial notebook.
Note:
As of POSEIDON v1.4, there is a surface albedo database that is now included in the inputs folder. If you already have downloaded the v1.3 version of the inputs, and don’t want to redownload them, you can access the surface_reflectivities folder here to add it manually.
Download the surface_reflectivities.zip folder, unzip it, and add it to your inputs folder.
As of 1.4 you should have four top-level folders in inputs: stellar_grids, opacity, chemistry_grids, and surface_reflectivities.
Note:
If you end up using the emission and reflection modules in POSEIDON for a publication, please consider citing relevant PICASO papers (Batalha et al. 2019, Mukherjee et al. 2021, and Mang et al. 2026) since the Toon radiative transfer used for reflection and thermal emission (w/ scattering) were originally sourced from PICASO.
Additionally, consider citing Mullens et al. 2026 (in prep), which introduced surfaces w/ wavelength dependent albedos into POSEIDON.
Surfaces w/ Wavelength-Dependent Albedo are Important!
The above figure emphasizes why including surfaces with wavelength-dependent albedos is important for generating and interpreting emergent spectra from rocky worlds. For more details on the figure above, see Mullens et al (2026)
Here, we look at Trappist 1-e and 55 Cancri e.
Rocky Exoplanets w/ Overlaying Atmospheres (TRAPPIST-1e)
Here, we build off the ‘Transmission Spectra for Terrestrial Planets’ tutorial by defining Trappist-1e with a reflecting and emitting surface.
We will explore a wavelength range that explores short wavelengths (0.3-1 um) where reflection dominates, and longer wavelengths (1-12 um) where thermal emission dominates.
[1]:
from POSEIDON.core import create_star, create_planet
from POSEIDON.constants import R_Sun, R_E, M_E
from POSEIDON.core import wl_grid_constant_R
from scipy.constants import parsec as pc
from scipy.constants import au
#***** Wavelength grid *****#
wl_min = 0.3 # Minimum wavelength (um)
wl_max = 12 # Maximum wavelength (um)
R = 1000 # Wavelength resolution
wl = wl_grid_constant_R(wl_min, wl_max, R)
#***** Define stellar properties *****#
R_s = 0.11697*R_Sun # Stellar radius (m)
T_s = 2559.0 # Stellar effective temperature (K)
Met_s = 0.04 # Stellar metallicity [log10(Fe/H_star / Fe/H_solar)]
log_g_s = 5.21 # Stellar log surface gravity (log10(cm/s^2) by convention)
# Create the stellar object
star = create_star(R_s, T_s, log_g_s, Met_s, wl = wl) # <---- Remember to provide wavelength for emission/reflection!
#***** Define planet properties *****#
planet_name = 'TRAPPIST-1e' # Planet name used for plots, output files etc.
R_p = 0.917985*R_E # Planetary radius (m)
M_p = 0.6356*M_E # Planetary mass (kg)
T_eq = 255.0 # Equilibrium temperature (K)
a_p = 0.2925 *au # Planetary distance for host star (au)
d = 12.1 * pc # System distance (pc)
# Create the planet object
planet = create_planet(planet_name, R_p, mass = M_p, T_eq = T_eq, a_p = a_p, d = d)
How to Define a Surface
In order to define a surface in our model, we set surface = True in the define_model() function. A hard surface is parameterized by the surface pressure parameter, log_P_surf.
If we want a reflecting surface, we need to set reflection = True.
If we want a thermally emitting surface, we need to set either thermal = True. If you want to use the full Toon emission functionality, which includes thermal multiple scattering, you will also need to set thermal_scattering = True.
The reflection algorithm defines the surface as having an albedo vs wavelength, whereas the thermal emission algorithm defines the surface as having an emissivity vs wavelength (where emissivity is 1 - albedo).
The temperature of the surface in models with overlaying atmospheres is equal the temperature of the pressure-temperature profile at P_surf.
Surface Specific Parameters
As of POSEIDON v1.4, there is a new parameter in define_model() called: surface_model
The default surface (where the albedo = 0 at all wavelengths), is called gray. This surface is essentially the same as an opaque gray deck (used in cloud models), and is infinitely absorbing (perfect blackbody).
The surface model constant defines a surface with a constant albedo (albedo_surf).
The surface model lab_data defines a surface that utilizes actual lab data. POSEIDON comes pre-included with a surface albedo database (see ‘Surface Albedo Database’ tab). This functionality allows for multiple surface components, which are defined in the list surface_components. Each surface component has a perecentage; NOTE: these percentages are always normalized to 100%.
There are two additional parameters for multi-component surfaces: surface_percentage_option (which is linear or log) and surface_percentage_apply_to (which is models or albedos).
For this tutorial, we will define our percentages in linear space, though defining them in log space has certain advantages for retrievals (see ‘Retrieval Tutorial: Bare Rocky Exoplanets’).
The second parameter defines how the percentages are applied: they are either applied to each model (the default; for N surfaces, N forward models are computed) or to the albedos (for N surfaces, 1 forward model is computed). For computational expediency, it is faster to allow percentages to be applied to the albedos isntead of the models.
For surfaces that utilize lab data, users can simply add their own lab data: just generate a txt file with the first column being wavelength, second column being albedo, and place it in the inputs/surface_reflectivities folder. POSEIDON will automatically detect whether or not their is a txt file with the name. For more details on the pre-included albedo database, as well as preferences on their format, see the Surface Albedo Database.
Case 1: The Default ‘Gray’ Surface (Albedo = 0 at all wavelengths)
[2]:
from POSEIDON.core import define_model
import numpy as np
#***** Define models *****#
model_name_gray = 'Gray-Surface'
bulk_species = ['N2'] # For terrestrial planets, only the single most abundant gas should be provided here
param_species = ['O2', 'H2O'] # Lets just define O2 and H2O for simplicity
# Create a model object with a gray surface
model_gray = define_model(model_name_gray, bulk_species, param_species,
PT_profile = 'file_read', X_profile = 'file_read',
radius_unit = 'R_E',
surface = True, # <----- Set surface = True
surface_model = 'gray', # <----- Default surface model
reflection = True, thermal = True, thermal_scattering = True) # <----- Set reflection and scattering to True
print(model_gray['param_names'])
['R_p_ref' 'log_O2' 'log_H2O' 'log_P_surf']
Case 2: A ‘Constant Albedo’ Surface (Albedo = albedo_surf at all wavelengths)
Here, we will show that surfaces work with and without thermal scattering turned on. Models without thermal scattering will run much faster, but will only consider absorption and emission processes, ignoring scattering effects like forward and multiple scattering (see ‘Thermal Scattering’ tutorial).
In clear models, this difference is minimal, but becomes potentially important when considering Mie scattering aerosols.
Note that reflection, by default, is a scattering process (and therefore will be relatively slow), however there is no simple replacement for reflection.
[ ]:
model_name_const = 'Constant-Albedo-Surface'
bulk_species = ['N2']
param_species = ['O2', 'H2O']
# Without thermal scattering
model_const_no_scattering = define_model(model_name_const, bulk_species, param_species,
PT_profile = 'file_read', X_profile = 'file_read',
radius_unit = 'R_E', surface = True, # <----- Set surface = True
surface_model = 'constant', # <----- Set surface_model to 'constant'
reflection = True, thermal = True, thermal_scattering = False) # <----- Set reflection and thermal to True
# With scattering
model_const_w_scattering = define_model(model_name_const, bulk_species, param_species,
PT_profile = 'file_read', X_profile = 'file_read',
radius_unit = 'R_E', surface = True, # <----- Set surface = True
surface_model = 'constant', # <----- Set surface_model to 'constant'
reflection = True, thermal = True, thermal_scattering = True) # <----- Set reflection, thermal, and thermal scattering to True
print(model_const_no_scattering['param_names'])
['R_p_ref' 'log_O2' 'log_H2O' 'log_P_surf' 'albedo_surf']
Case 3: Define a ‘Lab Data’ Surface
Here we define a multi-component surface defined by two different lab datasets. We use linear percentages (0 to 1), and test both applying the surface percentages to models and albedos.
As we will see below, there is no real difference in applying the surface percentages to either the albedos or the models. The benefit of applying the percentages to the albedos is that only one forward model is generated, which is much faster.
[4]:
model_name_lab = 'Lab-Data-Albedo-Surface'
bulk_species = ['N2']
param_species = ['O2', 'H2O']
# List surface compoenents here
# Names match the .txt files in 'surface_reflectivities' folder
# Lets use two common, expected rocky surfaces: basalt and granite
surface_components = ['Tholeiitic_basalt_H25','Granitoid_H12']
# Applying percentages to albedos (1 forward model generated)
model_lab_albedos = define_model(model_name_lab, bulk_species, param_species,
PT_profile = 'file_read', X_profile = 'file_read',
radius_unit = 'R_E', surface = True, # <----- Set surface = True
reflection = True, thermal = True, thermal_scattering = True, # <----- Set reflection and scattering to True
surface_model = 'lab_data', # <----- Set surface_model to 'lab_data'
surface_components = surface_components, # <----- Input surface_components
surface_percentage_apply_to= 'albedos', # <----- Apply percentages to albedos
surface_percentage_option = 'linear') # <----- This is the default option (percentages, 0 to 1)
# Applying percentages to models (2 forward models generated)
model_lab_models = define_model(model_name_lab, bulk_species, param_species,
PT_profile = 'file_read', X_profile = 'file_read',
radius_unit = 'R_E', surface = True, # <----- Set surface = True
reflection = True, thermal = True, thermal_scattering = True, # <----- Set reflection and scattering to True
surface_model = 'lab_data', # <----- Set surface_model to 'lab_data'
surface_components = surface_components, # <----- Input surface_components
surface_percentage_apply_to= 'models') # <----- Apply percentages to albedos
print(model_lab_albedos['param_names'])
['R_p_ref' 'log_O2' 'log_H2O' 'log_P_surf'
'Tholeiitic_basalt_H25_percentage' 'Granitoid_H12_percentage']
Load in the files that define the atmospheric state of TRAPPIST-1e (same steps as transmission_terrestrial.ipynb)
[5]:
from POSEIDON.utility import read_PT_file
#***** Create model pressure grid (same for both models) *****#
# Specify the pressure grid of the atmosphere
P_min = 1.0e-7 # 0.1 ubar
P_max = 10.0 # 10 bar (you can extend the atmosphere deeper than the surface)
N_layers = 100 # 100 layers
P = np.logspace(np.log10(P_max), np.log10(P_min), N_layers)
# Specify location of the P-T profile file
PT_file_dir = '../../../POSEIDON/reference_data/models/TRAPPIST-1e'
PT_file_name_2 = 'TRAPPIST-1e_1.0bar_100xCO2_Archean_PT.txt'
# Read the P-T profile files
T_Archean = read_PT_file(PT_file_dir, PT_file_name_2, P, skiprows = 1,
P_column = 2, T_column = 3)
from POSEIDON.utility import read_chem_file
# Specify location of the composition file
chem_file_dir = '../../../POSEIDON/reference_data/models/TRAPPIST-1e'
chem_file_name_2 = 'TRAPPIST-1e_1.0bar_100xCO2_Archean_chem.txt'
chem_species_file = ['N2', 'O2', 'O3', 'H2O', 'CH4', 'N2O', 'CO2', 'CO'] # Same ordering for both files
# Read the composition file
X_Archean = read_chem_file(chem_file_dir, chem_file_name_2, P, chem_species_file,
chem_species_in_model = model_gray['chemical_species'],
skiprows = 1)
Let’s define our atmosphere objects.
For our gray surface, all we need to do is provide log_P_surf (the pressure at which the surface top is located)
[ ]:
from POSEIDON.core import make_atmosphere
P_ref = 1 # We'll set the reference pressure at the surface
R_p_ref = R_p # Radius at reference pressure
log_P_surf = 0 # Surface pressure is 1 bar
surface_params = np.array([log_P_surf]) #<---- Put surface params into new list, surface_params
# Generate the atmosphere
atmosphere_blackbody = make_atmosphere(planet, model_gray, P, P_ref, R_p_ref,
T_input = T_Archean, X_input = X_Archean,
surface_params = surface_params) #<---- Put surface params into make_atmosphere
For our constant albedo surface, we need to adtionally provide albedo_surf (which ranges 0 to 1 where 0 is absorbing 100% and 1 is reflecting 100%)
[7]:
from POSEIDON.core import make_atmosphere
from POSEIDON.visuals import plot_geometry, plot_PT, plot_chem
P_ref = 1 # We'll set the reference pressure at the surface
R_p_ref = R_p # Radius at reference pressure
log_P_surf = 0 # Surface pressure is 1 bar
surface_albedo = 0.5 # Half reflecting, half absorbing
surface_params = np.array([log_P_surf, surface_albedo,]) #<---- Put surface params into new list, surface_params
# Generate the atmospheres
atmosphere_const_no_scattering = make_atmosphere(planet, model_const_no_scattering, P, P_ref, R_p_ref,
T_input = T_Archean, X_input = X_Archean,
surface_params = surface_params) #<---- Put surface params into make_atmosphere
atmosphere_const_w_scattering = make_atmosphere(planet, model_const_w_scattering, P, P_ref, R_p_ref,
T_input = T_Archean, X_input = X_Archean,
surface_params = surface_params) #<---- Put surface params into make_atmosphere
For lab data models, we instead need to define the surface percentages.
Note that these will always be normalized so that they equal one! (i.e. if basalt and granite both equal 1, they will be normalized down to 0.5)
[8]:
from POSEIDON.core import make_atmosphere
from POSEIDON.visuals import plot_geometry, plot_PT, plot_chem
P_ref = 1 # We'll set the reference pressure at the surface
R_p_ref = R_p # Radius at reference pressure
log_P_surf = 0 # Surface pressure is 1 bar
Basalt_percentage = 0.5 # 50% basalt
Granite_percentage = 0.5 # 50% granite
surface_params = np.array([log_P_surf,
Basalt_percentage,
Granite_percentage]) #<---- Put surface params into new list, surface_params
# Generate the atmospheres
atmosphere_lab_albedos = make_atmosphere(planet, model_lab_albedos, P, P_ref, R_p_ref,
T_input = T_Archean, X_input = X_Archean,
surface_params = surface_params) #<---- Put surface params into make_atmosphere
atmosphere_lab_models = make_atmosphere(planet, model_lab_models, P, P_ref, R_p_ref,
T_input = T_Archean, X_input = X_Archean,
surface_params = surface_params) #<---- Put surface params into make_atmosphere
[9]:
from POSEIDON.core import read_opacities, wl_grid_constant_R
#***** Read opacity data *****#
opacity_treatment = 'opacity_sampling'
# Define fine temperature grid (K)
T_fine_min = 100 # 100 K lower limit covers the TRAPPIST-1e P-T profile
T_fine_max = 300 # 300 K upper limit covers the TRAPPIST-1e P-T profile
T_fine_step = 10 # 10 K steps are a good tradeoff between accuracy and RAM
T_fine = np.arange(T_fine_min, (T_fine_max + T_fine_step), T_fine_step)
# Define fine pressure grid (log10(P/bar))
log_P_fine_min = -6.0 # 1 ubar is the lowest pressure in the opacity database
log_P_fine_max = 0.0 # 1 bar is the surface pressure, so no need to go deeper
log_P_fine_step = 0.2 # 0.2 dex steps are a good tradeoff between accuracy and RAM
log_P_fine = np.arange(log_P_fine_min, (log_P_fine_max + log_P_fine_step),
log_P_fine_step)
# Create opacity object (both models share the same molecules, so we only need one)
opac = read_opacities(model_gray, wl, opacity_treatment, T_fine, log_P_fine,
opacity_database = 'Temperate')
Reading in cross sections in opacity sampling mode...
N2-N2 done
N2-H2O done
O2-O2 done
O2-N2 done
O2 done
H2O done
Opacity pre-interpolation complete.
Comparing TRAPPIST-1e model spectra
Below, we will just compute the Fp/Fs spectrum, but note that you can get just Fp by replacing spectrum_type with direct_emission
[ ]:
from POSEIDON.core import compute_spectrum
from POSEIDON.visuals import plot_spectra
from POSEIDON.utility import plot_collection
# Compute gray surface spectrum
spectrum_gray_FpFs = compute_spectrum(planet, star, model_gray, atmosphere_blackbody, opac, wl,
spectrum_type = 'emission', use_photosphere_radius= True) # Set spectrum to emission
# Compute constant albedo surface spectra
spectrum_const_no_scattering_FpFs = compute_spectrum(planet, star, model_const_no_scattering,
atmosphere_const_no_scattering, opac, wl,
spectrum_type = 'emission', use_photosphere_radius= True)
spectrum_const_w_scattering_FpFs = compute_spectrum(planet, star, model_const_w_scattering,
atmosphere_const_w_scattering, opac, wl,
spectrum_type = 'emission', use_photosphere_radius= True)
# Compute lab data surface spectra
spectrum_lab_data_albedos_FpFs = compute_spectrum(planet, star, model_lab_albedos, atmosphere_lab_albedos,
opac, wl, spectrum_type = 'emission', use_photosphere_radius= True)
spectrum_lab_data_models_FpFs = compute_spectrum(planet, star, model_lab_models, atmosphere_lab_models,
opac, wl, spectrum_type = 'emission', use_photosphere_radius= True)
[ ]:
# Plot the gray surface
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
# Two axes, to have a zoom-in on the shorter wavelengths where reflection dominates
fig_combined = plt.figure(figsize=(16,9/2),constrained_layout=True) # Change (9,5.5) to alter the aspect ratio
# This function is the magic. Each letter corresponds to one matplotlib axis, which you can then pass to POSEIDON's plotting functions
axd = fig_combined.subplot_mosaic(
"""
AABBB
"""
)
spectra = [] # Empty plot collection
# Add one model spectra to the plot collection object
spectra = plot_collection(spectrum_gray_FpFs, wl, collection = spectra)
# Reflection Wavelengths only
fig_spec = plot_spectra(spectra, planet, R_to_bin = 1000, plot_full_res = False,
spectra_labels = ['Gray Surface'],
colour_list = ['gray'],
plt_label = ('Gray Surface, log P surf = 0\nUVIS Wavelengths'),
figure_shape = 'wide',
y_unit = 'eclipse_depth',
wl_max = 1,
y_max = 1e-8,
ax = axd['A'])
# Full Spectrum
fig_spec = plot_spectra(spectra, planet, R_to_bin = 1000, plot_full_res = False,
spectra_labels = ['Gray Surface'],
colour_list = ['gray'],
plt_label = ('Gray Surface, log P surf = 0\nFull Spectrum'),
figure_shape = 'wide',
y_unit = 'eclipse_depth',
wl_axis = 'linear',
ax = axd['B'])
<Figure size 853.36x480 with 0 Axes>
<Figure size 853.36x480 with 0 Axes>
[ ]:
# Plot the const albedo
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
# Two axes, to have a zoom-in on the shorter wavelengths where reflection dominates
fig_combined = plt.figure(figsize=(16,9/2),constrained_layout=True) # Change (9,5.5) to alter the aspect ratio
# This function is the magic. Each letter corresponds to one matplotlib axis, which you can then pass to POSEIDON's plotting functions
axd = fig_combined.subplot_mosaic(
"""
AABBB
"""
)
spectra = [] # Empty plot collection
# Add one model spectra to the plot collection object
spectra = plot_collection(spectrum_const_no_scattering_FpFs, wl, collection = spectra)
spectra = plot_collection(spectrum_const_w_scattering_FpFs, wl, collection = spectra)
# Reflection Wavelengths only
fig_spec = plot_spectra(spectra, planet, R_to_bin = 1000, plot_full_res = False,
spectra_labels = ['No Thermal Scattering', 'With Thermal Scattering'],
colour_list = ['red','yellow'],
plt_label = ('UVIS Wavelengths'),
figure_shape = 'wide',
y_unit = 'eclipse_depth',
wl_max = 1,
y_max = 1e-8,
line_width_list = [3,1],
ax = axd['A'])
# Full Spectrum
fig_spec = plot_spectra(spectra, planet, R_to_bin = 1000, plot_full_res = False,
spectra_labels = ['No Thermal Scattering', 'With Thermal Scattering'],
colour_list = ['red','yellow'],
plt_label = ('Constant Albedo (0.5) Surface, log P surf = 0\nFull Spectrum'),
figure_shape = 'wide',
y_unit = 'eclipse_depth',
wl_axis = 'linear',
line_width_list = [3,1],
ax = axd['B'])
<Figure size 853.36x480 with 0 Axes>
<Figure size 853.36x480 with 0 Axes>
[ ]:
# Plot the lab data
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
# Two axes, to have a zoom-in on the shorter wavelengths where reflection dominates
fig_combined = plt.figure(figsize=(16,9/2),constrained_layout=True) # Change (9,5.5) to alter the aspect ratio
# This function is the magic. Each letter corresponds to one matplotlib axis, which you can then pass to POSEIDON's plotting functions
axd = fig_combined.subplot_mosaic(
"""
AABBB
"""
)
spectra = [] # Empty plot collection
# Add one model spectra to the plot collection object
spectra = plot_collection(spectrum_lab_data_albedos_FpFs, wl, collection = spectra)
spectra = plot_collection(spectrum_lab_data_models_FpFs, wl, collection = spectra)
# Reflection Wavelengths only
fig_spec = plot_spectra(spectra, planet, R_to_bin = 1000, plot_full_res = False,
spectra_labels = ['Percentage on Albedos', 'Percentage on Models'],
colour_list = ['magenta','darkturquoise'],
plt_label = ('UVIS Wavelengths'),
figure_shape = 'wide',
y_unit = 'eclipse_depth',
wl_max = 1,
y_max = 1e-8,
line_width_list = [3,1],
ax = axd['A'])
# Full Spectrum
fig_spec = plot_spectra(spectra, planet, R_to_bin = 1000, plot_full_res = False,
spectra_labels = ['Percentage on Albedos', 'Percentage on Models'],
colour_list = ['magenta','darkturquoise'],
plt_label = ('Lab Data (Half Granite, Half Basalt) Surface, log P surf = 0\nFull Spectra'),
figure_shape = 'wide',
y_unit = 'eclipse_depth',
wl_axis = 'linear',
line_width_list = [3,1],
ax = axd['B'])
<Figure size 853.36x480 with 0 Axes>
<Figure size 853.36x480 with 0 Axes>
Adding Clouds
We can also add Mie scattering clouds to our models with surfaces.
In particular, we can add multiple aerosols in 1D models (non-patchy clouds), or we can add 1 species in a patchy cloud model i.e., surfaces are NOT compatible with the two-sector patchy cloud featured in ‘Advanced Aerosols in Eclipse Geometry’ tutorial.
In this portion, we try both
This part of the tutorial assumes familiarity with Mie scattering clouds in POSEIDON.
[ ]:
model_name_lab_plus_clouds = 'Lab-Data-Albedo-Surface-Plus-Clouds'
bulk_species = ['N2']
param_species = ['O2', 'H2O']
# List surface compoenents here
# Names match the .txt files in 'surface_reflectivities' folder
# Lets use two common, expected rocky surfaces: basalt and granite
surface_components = ['Tholeiitic_basalt_H25','Granitoid_H12']
# List aerosols here
# Here we include liquid h2o and ice aerosols
aerosol_species = ['H2O_ice', 'H2O']
# Applying percentages to albedos (1 forward model generated) to speed up models
# Here, we have both aerosols be part of a single slab cloud
model_lab_h2o_one_slab = define_model(model_name_lab, bulk_species, param_species,
PT_profile = 'file_read', X_profile = 'file_read',
radius_unit = 'R_E', surface = True, # <----- Set surface = True
reflection = True, thermal = True, thermal_scattering = True, # <----- Set reflection and scattering to True
cloud_model = 'Mie',cloud_type = 'one_slab', # Slab model for both aerosol species
aerosol_species = aerosol_species, # Including aerosols
surface_model = 'lab_data', # <----- Set surface_model to 'lab_data'
surface_components = surface_components, # <----- Input surface_components
surface_percentage_apply_to= 'albedos', # <----- Apply percentages to albedos
surface_percentage_option = 'linear') # <----- This is the default option (percentages, 0 to 1)
print(model_lab_h2o_one_slab['param_names'])
print()
# Here, we do patchy clouds with just ice aerosols
aerosol_species = ['H2O_ice']
model_lab_h2o_patchy_slab = define_model(model_name_lab, bulk_species, param_species,
PT_profile = 'file_read', X_profile = 'file_read',
radius_unit = 'R_E', surface = True, # <----- Set surface = True
reflection = True, thermal = True, thermal_scattering = True, # <----- Set reflection and scattering to True
cloud_model = 'Mie',cloud_type = 'slab', # Slab model for both aerosol species
aerosol_species = aerosol_species, # Including aerosols
cloud_dim = 2, # Patchy Clouds
surface_model = 'lab_data', # <----- Set surface_model to 'lab_data'
surface_components = surface_components, # <----- Input surface_components
surface_percentage_apply_to= 'albedos', # <----- Apply percentages to albedos
surface_percentage_option = 'linear') # <----- This is the default option (percentages, 0 to 1)
print(model_lab_h2o_patchy_slab['param_names'])
['R_p_ref' 'log_O2' 'log_H2O' 'log_P_top_slab' 'Delta_log_P'
'log_r_m_H2O_ice' 'log_X_H2O_ice' 'log_r_m_H2O' 'log_X_H2O' 'log_P_surf'
'Tholeiitic_basalt_H25_percentage' 'Granitoid_H12_percentage']
['R_p_ref' 'log_O2' 'log_H2O' 'f_cloud' 'log_P_top_slab_H2O_ice'
'Delta_log_P_H2O_ice' 'log_r_m_H2O_ice' 'log_X_H2O_ice' 'log_P_surf'
'Tholeiitic_basalt_H25_percentage' 'Granitoid_H12_percentage']
[15]:
# Add the aerosols above to the opac object
from POSEIDON.clouds import switch_aerosol_in_opac
opac_h2o_aerosols = switch_aerosol_in_opac(model_lab_h2o_one_slab, opac)
Reading in database for aerosol cross sections...
[ ]:
from POSEIDON.core import make_atmosphere
from POSEIDON.visuals import plot_geometry, plot_PT, plot_chem
P_ref = 1 # We'll set the reference pressure at the surface
R_p_ref = R_p # Radius at reference pressure
log_P_surf = 0 # Surface pressure is 1 bar
Basalt_percentage = 0.5 # 50% basalt
Granite_percentage = 0.5 # 50% granite
surface_params = np.array([log_P_surf,
Basalt_percentage,
Granite_percentage]) #<---- Put surface params into new list, surface_params
log_P_top_slab = -2
Delta_log_P = 1
log_r_m_H2O = -2
log_X_H2O = -12
# Assume the same r_m and X for both H2O aerosols, for simplicity
cloud_params = np.array([log_P_top_slab, Delta_log_P, log_r_m_H2O, log_X_H2O, log_r_m_H2O, log_X_H2O])
# Generate the atmospheres
atmosphere_lab_h2o_one_slab = make_atmosphere(planet, model_lab_h2o_one_slab, P, P_ref, R_p_ref,
T_input = T_Archean, X_input = X_Archean,
cloud_params = cloud_params,
surface_params = surface_params) #<---- Put surface params into make_atmosphere
# Lets assume 50% patchy cloud of just water ice
f_cloud = 0.5
cloud_params = np.array([f_cloud, log_P_top_slab, Delta_log_P, log_r_m_H2O, log_X_H2O])
atmosphere_lab_h2o_patchy_slab = make_atmosphere(planet, model_lab_h2o_patchy_slab, P, P_ref, R_p_ref,
T_input = T_Archean, X_input = X_Archean,
cloud_params = cloud_params,
surface_params = surface_params) #<---- Put surface params into make_atmosphere
[ ]:
from POSEIDON.core import compute_spectrum
from POSEIDON.visuals import plot_spectra
from POSEIDON.utility import plot_collection
# Compute spectra
# Remember to use new opac object here with the aerosols included
spectrum_lab_h2o_one_slab_FpFs = compute_spectrum(planet, star, model_lab_h2o_one_slab,
atmosphere_lab_h2o_one_slab, opac_h2o_aerosols, wl,
spectrum_type = 'emission', use_photosphere_radius= True)
spectrum_lab_h2o_patchy_slab_FpFs = compute_spectrum(planet, star, model_lab_h2o_patchy_slab,
atmosphere_lab_h2o_patchy_slab, opac_h2o_aerosols, wl,
spectrum_type = 'emission', use_photosphere_radius= True)
[ ]:
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
# Two axes, to have a zoom-in on the shorter wavelengths where reflection dominates
fig_combined = plt.figure(figsize=(16,9/2),constrained_layout=True) # Change (9,5.5) to alter the aspect ratio
# This function is the magic. Each letter corresponds to one matplotlib axis, which you can then pass to POSEIDON's plotting functions
axd = fig_combined.subplot_mosaic(
"""
AABBB
"""
)
spectra = [] # Empty plot collection
# Add the two model spectra to the plot collection object
spectra = plot_collection(spectrum_lab_h2o_one_slab_FpFs, wl, collection = spectra)
spectra = plot_collection(spectrum_lab_h2o_patchy_slab_FpFs, wl, collection = spectra)
# Reflection Wavelengths only
fig_spec = plot_spectra(spectra, planet, R_to_bin = 1000, plot_full_res = False,
spectra_labels = ['One Slab (H2O liquid + Ice)', 'Patchy Slab (50% Ice)'],
colour_list = ['navy','aqua'],
plt_label = ('UVIS Wavelengths'),
figure_shape = 'wide',
y_unit = 'eclipse_depth',
wl_max = 1,
y_max = 1e-8,
line_width_list = [3,1],
ax = axd['A'])
# Full Spectrum
fig_spec = plot_spectra(spectra, planet, R_to_bin = 1000, plot_full_res = False,
spectra_labels = ['One Slab (H2O liquid + Ice)', 'Patchy Slab (50% Ice)'],
colour_list = ['navy','aqua'],
plt_label = ('Cloudy Spectra + Lab Data Surface\nFull Spectra'),
figure_shape = 'wide',
y_unit = 'eclipse_depth',
wl_axis = 'linear',
line_width_list = [3,1],
ax = axd['B'])
<Figure size 853.36x480 with 0 Axes>
<Figure size 853.36x480 with 0 Axes>
Rocky Exoplanets w/ Overlaying Atmospheres (55 Cancri e)
Lets now look at a hotter exoplanet that may or may not have an overlaying atmosphere. We both look at the case where there is an overlaying CO\(_2\) atmosphere, and one where it is just a bare rock.
[ ]:
#***** Define planet properties *****#
planet_name_55 = '55-Cancri-e' # Planet name used for plots, output files etc.
R_p = 1.95*R_E # Planetary radius (m)
M_p = 8.8*M_E
d = 12.5855*pc
a_p = 0.01544*au
# Create the planet object
planet_55 = create_planet(planet_name_55, R_p, mass = M_p, a_p = a_p)
#***** Define stellar properties *****#
R_s = 0.964411455860500*R_Sun # Stellar radius (m)
T_s = 5317.89014944000 # Stellar effective temperature (K)
Met_s = 0.38330181143999997 # Stellar metallicity [log10(Fe/H_star / Fe/H_solar)] <--- note: for PHOENIX, only the solar metallicity models are used
log_g_s = 4.45813087040500 # Stellar log surface gravity (log10(cm/s^2) by convention)
# Create the stellar object
star_55 = create_star(R_s, T_s, log_g_s, Met_s, wl = wl, stellar_grid = 'phoenix')
#***** Specify data location and instruments *****#
model_name_55 = '55-Cancri-e'
#***** Define model *****#
bulk_species = ['CO2'] # <- Here we just assume that the entire atmosphere is made of CO2
param_species = []
surface_components = ['Tholeiitic_basalt_H25'] # <- Just one surface component , for simplicity
# Create the model object
model_55 = define_model(model_name_55, bulk_species, param_species,
PT_profile = 'gradient', # < - parameterized PT profile
radius_unit = 'R_E', surface = True,
reflection = True, thermal = True, thermal_scattering = True,
surface_model = 'lab_data',
surface_components = surface_components)
print(model_55['param_names'])
['R_p_ref' 'T_high' 'T_deep' 'log_P_surf']
We generate a new opac object, with just CO2 included in it, loaded on a much larger temperature space.
[ ]:
#***** Read opacity data *****#
opacity_treatment = 'opacity_sampling'
# Define fine temperature grid (K)
T_fine_min = 500 # Same as prior range for T
T_fine_max = 2500 # Same as prior range for T
T_fine_step = 20 # 10 K steps are a good tradeoff between accuracy and RAM
T_fine = np.arange(T_fine_min, (T_fine_max + T_fine_step), T_fine_step)
# Define fine pressure grid (log10(P/bar))
log_P_fine_min = -6.0 # 1 ubar is the lowest pressure in the opacity database
log_P_fine_max = 2.0 # 100 bar is the highest pressure in the opacity database
log_P_fine_step = 0.2 # 0.2 dex steps are a good tradeoff between accuracy and RAM
log_P_fine = np.arange(log_P_fine_min, (log_P_fine_max + log_P_fine_step),
log_P_fine_step)
#***** Specify fixed atmospheric settings *****#
# Atmospheric pressure grid
P_min = 1.0e-6 # 10 nbar
P_max = 100 # 10 bar
N_layers = 100 # 100 layers
# Let's space the layers uniformly in log-pressure
P = np.logspace(np.log10(P_max), np.log10(P_min), N_layers)
# Specify the reference pressure
P_ref = 1e-2 # Retrieved R_p_ref parameter will be the radius at 10 bar
# Pre-interpolate the opacities
opac_55 = read_opacities(model_55, wl, opacity_treatment, T_fine, log_P_fine, database_version = '1.3')
Reading in cross sections in opacity sampling mode...
CO2-CO2 done
CO2 done
Opacity pre-interpolation complete.
Here we generate spectra at four different surface pressures to see how the emergent spectra changes when the atmosphere is thick (100 and 1 bar) where CO2 features dominate, and thin (1e-2 and 1e-4 bar) where the surface features dominate.
[70]:
from POSEIDON.core import make_atmosphere
R_p_ref = 1.92 * R_E
T_high = 1668
T_deep = 2380
# Lets loop over P_surf to see the interplay between atmospheric signatures and surface signatures
spectra = []
log_P_surf_array = [2,0,-2,-4]
for log_P_surf in log_P_surf_array:
log_X_params = np.array([])
PT_params = np.array( [T_high, T_deep] )
surface_params = np.array( [log_P_surf])
atmosphere = make_atmosphere(planet_55, model_55, P, P_ref,
R_p_ref, PT_params, log_X_params = log_X_params,
surface_params = surface_params
)
spectrum = compute_spectrum(planet_55, star_55, model_55, atmosphere, opac_55, wl,
spectrum_type = 'emission')
spectra = plot_collection(spectrum, wl, collection = spectra)
[ ]:
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
# Two axes, to have a zoom-in on the shorter wavelengths where reflection dominates
fig_combined = plt.figure(figsize=(16,9/2),constrained_layout=True) # Change (9,5.5) to alter the aspect ratio
# This function is the magic. Each letter corresponds to one matplotlib axis, which you can then pass to POSEIDON's plotting functions
axd = fig_combined.subplot_mosaic(
"""
AABBB
"""
)
# Reflection Wavelengths only
fig_spec = plot_spectra(spectra, planet_55, R_to_bin = 100, plot_full_res = False,
spectra_labels = ['P surf = 100 bar', 'P surf = 1 bar', 'P surf = 1e-2 bar', 'P surf = 1e-4 bar'],
colour_list = ['gold','orange','red','crimson'],
plt_label = ('UVIS Wavelengths'),
figure_shape = 'wide',
y_unit = 'eclipse_depth',
wl_max = 1,
ax = axd['A'])
# Full Spectrum
fig_spec = plot_spectra(spectra, planet_55, R_to_bin = 100, plot_full_res = False,
spectra_labels = ['P surf = 100 bar', 'P surf = 1 bar', 'P surf = 1e-2 bar', 'P surf = 1e-4 bar'],
colour_list = ['gold','orange','red','crimson'],
plt_label = ('Changing log P surf\nFull Spectra'),
figure_shape = 'wide',
y_unit = 'eclipse_depth',
wl_axis = 'linear',
legend_location = 'lower right',
ax = axd['B'])
<Figure size 853.36x480 with 0 Axes>
<Figure size 853.36x480 with 0 Axes>
As we can see above, the basalt feature at 8.5 microns shows up in lower surface pressures (thinner atmospheres).
Bare Rocky Exoplanets w/ No Atmosphere
Finally, lets explore exoplanets without overlaying atmospheres.
We will specifically be exploring in this section how realstic, geological surfaces defined with lab data can have their own mid-infared features.
Setting a Bare Rock Model
In order to turn off the atmosphere and have a bare rocky exoplanet set surface = True, and disable_atmosphere = True in define_model(). Also, provide empty lists for the bulk and param species.
In these models, since there is not pressure-temeprature profile to assume the surface is equilibrium with, there is a new free parameter: T_surf.
This surface temperature is assumed to be homogenous for all points on the planet’s surface.
Here, we make models for a gray surface, a constant albedo surface, and a surface composed of lab data components.
[ ]:
from POSEIDON.core import define_model
import numpy as np
#***** Define models *****#
model_name_gray_bare_rock = 'Gray-Surface-Bare-Rock'
# For bare rocks, we have no bulk or param species
bulk_species = []
param_species = []
# Create a model object with a blackbody surface
model_gray_bare_rock = define_model(model_name_gray_bare_rock, bulk_species, param_species,
radius_unit = 'R_E', surface = True, # <----- Set surface = True
reflection = True, thermal_scattering = True, thermal = True,
disable_atmosphere = True) # <----- Set disable_atmosphere = True
print(model_gray_bare_rock['param_names'])
['R_p_ref' 'T_surf']
[73]:
model_name_const_bare_rock = 'Constant-Albedo-Surface-Bare-Rock'
bulk_species = []
param_species = []
model_const_bare_rock = define_model(model_name_const_bare_rock, bulk_species, param_species,
radius_unit = 'R_E', surface = True, # <----- Set surface = True
reflection = True, thermal_scattering = True, thermal = True, # <----- Set reflection and scattering to True
surface_model = 'constant', # <----- Set surface_model to 'constant'
disable_atmosphere = True) # <----- Set disable_atmosphere = True
print(model_const_bare_rock['param_names'])
['R_p_ref' 'T_surf' 'albedo_surf']
[74]:
model_name_lab_bare_rock = 'Lab-Data-Albedo-Surface-Bare-Rock'
bulk_species = []
param_species = []
surface_components = ['Tholeiitic_basalt_H25','Granitoid_H12'] # <----- List surface compoenents here
model_lab_bare_rock = define_model(model_name_lab_bare_rock, bulk_species, param_species,
radius_unit = 'R_E', surface = True, # <----- Set surface = True
reflection = True, thermal_scattering = True, thermal = True,
surface_model = 'lab_data', # <----- Set surface_model to 'lab_data'
surface_components = surface_components, # <----- Input surface_components
surface_percentage_apply_to= 'albedos', # <----- Apply percentages to albedos
surface_percentage_option = 'linear', # Default option
disable_atmosphere = True) # <----- Set disable_atmosphere = True
print(model_lab_bare_rock['param_names'])
['R_p_ref' 'T_surf' 'Tholeiitic_basalt_H25_percentage'
'Granitoid_H12_percentage']
Note that even though there is no ‘atmosphere’ anymore, we still make an atmosphere object since it is utilized by POSEIDON to store many free parameters.
[75]:
from POSEIDON.core import make_atmosphere
R_p_ref = R_p # Radius is just the planetary radius
T_surf = 300 # Surface temperature of 300K
surface_params = np.array([T_surf]) #<---- Put surface params into new list, surface_params
P = [] #<---- No atmosphere = no pressure
P_ref = []
atmosphere_gray_bare_rock = make_atmosphere(planet, model_gray_bare_rock, P, P_ref, R_p_ref,
surface_params = surface_params) #<---- Put surface params into make_atmosphere
[76]:
from POSEIDON.core import make_atmosphere
from POSEIDON.visuals import plot_geometry, plot_PT, plot_chem
R_p_ref = R_p
T_surf = 300 # Surface temperature of 300K
surface_albedo = 0.5 # Half reflecting, half absorbing
surface_params = np.array([T_surf, surface_albedo,]) #<---- Put surface params into new list, surface_params
# Generate the atmospheres
atmosphere_const_bare_rock = make_atmosphere(planet, model_const_bare_rock, P, P_ref, R_p_ref,
surface_params = surface_params) #<---- Put surface params into make_atmosphere
[77]:
from POSEIDON.core import make_atmosphere
from POSEIDON.visuals import plot_geometry, plot_PT, plot_chem
R_p_ref = R_p # Radius at reference pressure
T_surf = 300 # Surface temperature of 300K
Basalt_percentage = 0.5
Granite_percentage = 0.5
surface_params = np.array([T_surf,
Basalt_percentage,
Granite_percentage]) #<---- Put surface params into new list, surface_params
# Generate the atmospheres
atmosphere_lab_bare_rock = make_atmosphere(planet, model_lab_bare_rock, P, P_ref, R_p_ref,
surface_params = surface_params) #<---- Put surface params into make_atmosphere
For bare rocky exoplanets, there is no opac object.
Therefore, we define an empty opac list for compute_spectrum()
[78]:
from POSEIDON.core import compute_spectrum
from POSEIDON.visuals import plot_spectra
from POSEIDON.utility import plot_collection
opac_empty = [] # No opac object when there is no atmosphere
spectrum_gray_bare_rock_FpFs = compute_spectrum(planet, star, model_gray_bare_rock, atmosphere_gray_bare_rock,
opac_empty, wl,
spectrum_type = 'emission')
spectrum_const_bare_rock_FpFs = compute_spectrum(planet, star, model_const_bare_rock, atmosphere_const_bare_rock,
opac_empty, wl,
spectrum_type = 'emission')
spectrum_lab_data_bare_rock_FpFs = compute_spectrum(planet, star, model_lab_bare_rock, atmosphere_lab_bare_rock,
opac_empty, wl,
spectrum_type = 'emission')
[79]:
spectra = [] # Empty plot collection
# Add the three model spectra to the plot collection object
spectra = plot_collection(spectrum_gray_bare_rock_FpFs, wl, collection = spectra)
spectra = plot_collection(spectrum_const_bare_rock_FpFs, wl, collection = spectra)
spectra = plot_collection(spectrum_lab_data_bare_rock_FpFs, wl, collection = spectra)
# Produce figure
fig_spec = plot_spectra(spectra, planet, R_to_bin = 1000, plot_full_res = False,
spectra_labels = ['Blackbody Surface','Const Albedo (0.5)', 'Lab Data (Basalt + Granite)'],
colour_list = ['gray','darkturquoise', 'red'],
line_width_list = [10,7,3],
plt_label = ('Bare Rocky Exoplanet (T = 300K)'),
figure_shape = 'wide',
y_unit = 'eclipse_depth',
wl_axis = 'linear'
#y_max = 2e-14, y_min = 0
)
Surfaces Have Diagnostic Mid-Infrared Features
Lets explore the granite/basalt feature in the mid-infrared a bit more (and lets turn up the temperature!)
[ ]:
from POSEIDON.core import make_atmosphere
from POSEIDON.visuals import plot_geometry, plot_PT, plot_chem
R_p_ref = R_p # Radius at reference pressure
T_surf = 1500 # Surface temperature of 1500K
Basalt_percentage = 1
Granite_percentage = 0
surface_params = np.array([T_surf,
Basalt_percentage,
Granite_percentage]) #<---- Put surface params into new list, surface_params
# Generate the atmosphere
atmosphere_basalt = make_atmosphere(planet, model_lab_bare_rock, P, P_ref, R_p_ref,
surface_params = surface_params) #<---- Put surface params into make_atmosphere
Basalt_percentage = 0
Granite_percentage = 1
surface_params = np.array([T_surf,
Basalt_percentage,
Granite_percentage]) #<---- Put surface params into new list, surface_params
# Generate the atmosphere
atmosphere_granite = make_atmosphere(planet, model_lab_bare_rock, P, P_ref, R_p_ref,
surface_params = surface_params)
Basalt_percentage = 0.5
Granite_percentage = 0.5
surface_params = np.array([T_surf,
Basalt_percentage,
Granite_percentage]) #<---- Put surface params into new list, surface_params
# Generate the atmosphere
atmosphere_both = make_atmosphere(planet, model_lab_bare_rock, P, P_ref, R_p_ref,
surface_params = surface_params)
[81]:
spectrum_granite = compute_spectrum(planet, star, model_lab_bare_rock, atmosphere_granite,
opac_empty, wl,
spectrum_type = 'emission')
spectrum_basalt = compute_spectrum(planet, star, model_lab_bare_rock, atmosphere_basalt,
opac_empty, wl,
spectrum_type = 'emission')
spectrum_both = compute_spectrum(planet, star, model_lab_bare_rock, atmosphere_both,
opac_empty, wl,
spectrum_type = 'emission')
[82]:
spectra = [] # Empty plot collection
# Add the three model spectra to the plot collection object
spectra = plot_collection(spectrum_granite, wl, collection = spectra)
spectra = plot_collection(spectrum_basalt, wl, collection = spectra)
spectra = plot_collection(spectrum_both, wl, collection = spectra)
# Produce figure
fig_spec = plot_spectra(spectra, planet, R_to_bin = 1000, plot_full_res = False,
spectra_labels = ['Granite','Basalt', 'Both'],
colour_list = ['#aca69e',"#000000", "#524f4b"],
line_width_list = [3,3,3],
plt_label = ('Bare Rocky Exoplanet (T = 1500K)'),
figure_shape = 'wide',
y_unit = 'eclipse_depth',
wl_axis = 'linear',
legend_location = 'lower right'
#y_max = 2e-14, y_min = 0
)
If you were able to find evidence for both granite and basalt on the surface of an exoplanet, that would be very exciting since it could indicate the prescense of a process like plate tectonics!
Benchmark: Thin Atmosphere vs Bare Rock
Here, we compare a thin-isothermal atmosphere with a granite surface to a bare rocky planet (granite surface), to ensure the two agree.
For rocky planets with overlaying atmospheres, the pressure grid is cut at the surface pressure before radiative transfer is performed. To ensure that the thin atmosphere doesn’t cut the entire pressure grid, we generate a very fine pressure grid (N_layers = 1000).
[89]:
model_name_thin = 'Thin-Isothermal-Atmosphere'
bulk_species = ['N2']
param_species = []
# List surface compoenents here
# Names match the .txt files in 'surface_reflectivities' folder
# Lets use two common, expected rocky surfaces: basalt and granite
surface_components = ['Tholeiitic_basalt_H25','Granitoid_H12']
# Applying percentages to albedos (1 forward model generated)
model_thin = define_model(model_name_thin, bulk_species, param_species,
PT_profile = 'isotherm',
radius_unit = 'R_E', surface = True, # <----- Set surface = True
reflection = True, thermal = True, thermal_scattering = True, # <----- Set reflection and scattering to True
surface_model = 'lab_data', # <----- Set surface_model to 'lab_data'
surface_components = surface_components, # <----- Input surface_components
surface_percentage_apply_to= 'albedos', # <----- Apply percentages to albedos
surface_percentage_option = 'linear') # <----- This is the default option (percentages, 0 to 1)
[ ]:
from POSEIDON.core import read_opacities, wl_grid_constant_R
#***** Read opacity data *****#
opacity_treatment = 'opacity_sampling'
# Define fine temperature grid (K)
T_fine_min = 500 # 100 K lower limit covers the TRAPPIST-1e P-T profile
T_fine_max = 2500 # 300 K upper limit covers the TRAPPIST-1e P-T profile
T_fine_step = 20 # 10 K steps are a good tradeoff between accuracy and RAM
T_fine = np.arange(T_fine_min, (T_fine_max + T_fine_step), T_fine_step)
# Define fine pressure grid (log10(P/bar))
log_P_fine_min = -6.0 # 1 ubar is the lowest pressure in the opacity database
log_P_fine_max = 0.0 # 1 bar is the surface pressure, so no need to go deeper
log_P_fine_step = 0.2 # 0.2 dex steps are a good tradeoff between accuracy and RAM
log_P_fine = np.arange(log_P_fine_min, (log_P_fine_max + log_P_fine_step),
log_P_fine_step)
#***** Specify fixed atmospheric settings *****#
# Atmospheric pressure grid
P_min = 1.0e-6 # 10 nbar
P_max = 100 # 10 bar
N_layers = 1000 # 1000 layers, to ensure that P isn't cut off
# Let's space the layers uniformly in log-pressure
P = np.logspace(np.log10(P_max), np.log10(P_min), N_layers)
# Create opacity object (both models share the same molecules, so we only need one)
opac_thin = read_opacities(model_thin, wl, opacity_treatment, T_fine, log_P_fine,
database_version = '1.3')
Reading in cross sections in opacity sampling mode...
N2-N2 done
Opacity pre-interpolation complete.
[ ]:
from POSEIDON.core import make_atmosphere
from POSEIDON.visuals import plot_geometry, plot_PT, plot_chem
P_ref = 10**-5.85 # We'll set the reference pressure at the surface (see log P surf)
R_p_ref = R_p # Radius at reference pressure
T = 1500
PT_params = np.array([T])
log_X_params = np.array([])
log_P_surf = -5.85 # Surface pressure is very low!
Basalt_percentage = 0 # 50% basalt
Granite_percentage = 1.0 # 50% granite
surface_params = np.array([log_P_surf,
Basalt_percentage,
Granite_percentage]) #<---- Put surface params into new list, surface_params
# Generate the atmospheres
atmosphere_thin = make_atmosphere(planet, model_thin, P, P_ref, R_p_ref,
PT_params, log_X_params = log_X_params,
surface_params = surface_params) #<---- Put surface params into make_atmosphere
spectrum_thin = compute_spectrum(planet, star, model_thin, atmosphere_thin, opac_thin, wl,
spectrum_type = 'emission')
[ ]:
spectra = [] # Empty plot collection
# Add the two model spectra to the plot collection object
spectra = plot_collection(spectrum_granite, wl, collection = spectra)
spectra = plot_collection(spectrum_thin, wl, collection = spectra)
# Produce figure
fig_spec = plot_spectra(spectra, planet, R_to_bin = 1000, plot_full_res = False,
spectra_labels = ['Granite (Bare Rock)','Granite w/ Thin Isothermal Atmosphere'],
colour_list = ['#aca69e',"#000000"],
line_width_list = [10,3],
plt_label = ('Bare Rocky Exoplanet (T = 1500K)'),
figure_shape = 'wide',
y_unit = 'eclipse_depth',
wl_axis = 'linear',
legend_location = 'lower right'
#y_max = 2e-14, y_min = 0
)
Retrievals
All of the models above work with the retrieval framework in POSEIDON.
If you are more interested in bare rocky exoplanet retrievals, see ‘Retrieval Tutorial: Bare Rocky Exoplanets’