### Define and Create Advanced Atmospheric Models in POSEIDON

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/reflection_emission_advanced_shiny_gray_decks.ipynb

This example defines various atmospheric parameters including PT profile, gas mixing ratios, and different cloud/aerosol configurations. It then uses `make_atmosphere` to construct several distinct atmospheric models for simulation.

```python
from POSEIDON.core import make_atmosphere

R_p_ref               =   1.12 * R_J

# PT Params
log_kappa_IR          =   -4.82 
log_gamma             =   -0.42 
T_int                 =   255.7 
T_equ                 =   1159.5

# Gas Phase Mixing Ratios
log_CO                =   -5.30 
log_CO2               =   -3.07 
log_H2O               =   -5.33  
log_Na                =   -6.59 

# Gray Deck Parameters
log_P_cloud = -1
albedo_deck = 0.5 # We will assume this albedo of the gray deck for all models 
phi_cloud = 1

# Fuzzy Deck Parameters
log_P_top_deck_SiO2 = 0     # <---- Top of the opaque deck is at 1 bar (extends from 100 to 1 bar)
log_r_m_SiO2        = -2    # <---- Mean particle size of the SiO2 aerosols is 1e-2 microns
log_n_max_SiO2      = 10    # <---- The number density of SiO2 at the top of the opaque deck (at 1 bar)
f_SiO2              = 0.3   # <---- The fuzziness of aerosols (how number density evolves above the cloud deck)

# Slab Parameters
log_P_top_slab_SiO2 = -5  # <---- The top of the slab in pressure space (at 1e-5 bars)
Delta_log_P_SiO2    = 2   # <---- Extend of the slab in pressure space (extends down to 1e-3 bars)
log_r_m_SiO2        = -2  # <---- Mean particle size of the SiO2 aerosols is 1e-2 microns
log_X_SiO2          = -12 # <---- Volume mixing ratio of aerosol in the slab (1e-5 to 1e-3 bars)


PT_params = np.array([log_kappa_IR, log_gamma, T_int, T_equ])
log_X_params = np.array([log_CO, log_CO2, log_H2O, log_Na,])

# Shiny Opaque Deck 
cloud_params = np.array([log_P_cloud, albedo_deck, phi_cloud])

atmosphere_shiny_deck = make_atmosphere(planet, model_shiny_deck, P, P_ref, R_p_ref, 
                                              PT_params, log_X_params, cloud_params)

# Shiny SiO2 Fuzzy Deck 
cloud_params = np.array([albedo_deck, log_P_top_deck_SiO2, log_r_m_SiO2, log_n_max_SiO2, f_SiO2])

atmosphere_shiny_fuzzy_deck = make_atmosphere(planet, model_shiny_fuzzy_deck, P, P_ref, R_p_ref, 
                                              PT_params, log_X_params, cloud_params)

# Shiny Opaque Deck + SiO2 Uniform Mixing Ratio Aerosols
cloud_params = np.array([albedo_deck, log_P_cloud, log_r_m_SiO2 ,log_X_SiO2])

atmosphere_shiny_deck_ux = make_atmosphere(planet, model_shiny_deck_UX, P, P_ref, R_p_ref, 
                                              PT_params, log_X_params, cloud_params)

# Shiny Opaque Deck + SiO2 Slab
cloud_params = np.array([albedo_deck,log_P_cloud, log_P_top_slab_SiO2, Delta_log_P_SiO2, log_r_m_SiO2 ,log_X_SiO2])

atmosphere_shiny_deck_slab = make_atmosphere(planet, model_shiny_deck_slab, P, P_ref, R_p_ref, 
                                              PT_params, log_X_params, cloud_params)
```

--------------------------------

### Plot Retrieved Spectra (TiO Comparison Setup)

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/retrieval_basic.ipynb

Sets up the data structures for plotting retrieved spectrum confidence regions, comparing a reference model with a model lacking TiO. This is a preparatory step for visualizing the spectral fit differences.

```python
# Read retrieved spectrum confidence regions
wl, spec_low2, spec_low1, spec_median, \
spec_high1, spec_high2 = read_retrieved_spectrum(planet_name, model_ref['model_name'])

# Create composite spectra objects for plotting
spectra_median = plot_collection(spec_median, wl, collection = [])
spectra_low1 = plot_collection(spec_low1, wl, collection = []) 
spectra_low2 = plot_collection(spec_low2, wl, collection = []) 
spectra_high1 = plot_collection(spec_high1, wl, collection = []) 
spectra_high2 = plot_collection(spec_high2, wl, collection = [])

# Read retrieved spectrum confidence regions
wl, spec_low2, spec_low1, spec_median, \
spec_high1, spec_high2 = read_retrieved_spectrum(planet_name, model_no_TiO['model_name'])

# Create composite spectra objects for plotting
spectra_median = plot_collection(spec_median, wl, collection = spectra_median)
spectra_low1 = plot_collection(spec_low1, wl, collection = spectra_low1) 
spectra_low2 = plot_collection(spec_low2, wl, collection = spectra_low2) 
spectra_high1 = plot_collection(spec_high1, wl, collection = spectra_high1) 
spectra_high2 = plot_collection(spec_high2, wl, collection = spectra_high2)
```

--------------------------------

### Set Up Atmosphere with Cloud Parameters

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/transmission_clouds.ipynb

Use `make_atmosphere` to create an atmospheric model with specified PT and X parameters, along with cloud parameters for haze and cloud deck opacity.

```python
from POSEIDON.core import make_atmosphere

# PT and X parameters
T = 1200             # Temperature
log_H2O = -4         # H2O mixing ratio

# Cloud Parameters 
log_a       = 1.7 # <---- Rayleigh enhancement factor of the power-law haze
gamma       = -8  # <---- Scattering slope of the power-law haze
log_P_cloud = -2  # <---- log-pressure of the top of the infinite opacity deck (bar)

PT_params = np.array([T])
log_X_params = np.array([log_H2O])
cloud_params = np.array([log_a, gamma, log_P_cloud])

# Make atmosphere 
atmosphere_deck_haze = make_atmosphere(planet, model_deck_haze, P, P_ref, R_p_ref, PT_params, log_X_params, cloud_params)
```

--------------------------------

### Define Planetary Properties for 2D Model

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/transmission_2D.ipynb

Defines the properties of a planet, such as radius, gravity, and equilibrium temperature, for use in atmospheric modeling. This setup is for a warm Neptune example.

```python
#***** Define planet properties *****#

planet_name = 'HAT-P-26b'  # Planet name used for plots, output files etc.

R_p = 0.63*R_J      # Planetary radius (m)
g_p = 4.3712        # Gravitational field of planet (m/s^2)
T_eq = 1043.8       # Equilibrium temperature (K)

# Create the planet object
planet = create_planet(planet_name, R_p, gravity = g_p, T_eq = T_eq)
```

--------------------------------

### Initialize Wavelength Grid and Load Refractive Indices

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/aerosol_database.ipynb

Sets up the wavelength grid and loads real and imaginary refractive indices from a specified file. Ensure the file path is correct.

```python
wl_min = 0.2      # Minimum wavelength (um)
wl_max = 13.0     # Maximum wavelength (um)
R = 10000         # Spectral resolution of grid (R = wl/dwl)

wl = wl_grid_constant_R(wl_min, wl_max, R)

# Wavelength grid that Mie properties will be computed on 
wl_Mie = wl_grid_constant_R(wl[0], wl[-1], 1000)

# Preload the refractive indices from the file 
refractive_index_path = '../../../POSEIDON/reference_data/refractive_indices_txt_files/aerosol_database/WS15/'
file_name = refractive_index_path + 'H2O_complex.txt'

r_i_real, r_i_complex = load_refractive_indices_from_file(wl, file_name) # <---- Load in the real and imaginary indices
```

--------------------------------

### Install PyMultiNest

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/installation.md

Installs both MultiNest and its Python wrapper, PyMultiNest, using Conda-forge. This is the recommended and simplest installation method.

```bash
conda install -c conda-forge pymultinest
```

--------------------------------

### Set Up Atmospheric Grid and Parameters

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/aerosol_advanced_directionality.ipynb

Define nominal values for atmospheric parameters and create a log-pressure grid for simulations. Specify reference radius and pressure.

```python
# Some nominal values for W17b's atmosphere 
R_p_ref             =   1.69 * R_J
T                   =   1271.9
log_H2O             =   -2.96 
log_P_top_slab      =   -6.60
Delta_log_P         =   1.96 
log_r_m             =   -1.85

# Atmospheric pressure grid
P_min = 1.0e-8   # 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 = 10.0   # Retrieved R_p_ref parameter will be the radius at 10 bar
```

--------------------------------

### Install mpi4py

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/installation.md

Installs the mpi4py package, a dependency for PyMultiNest, using Conda-forge. This command should be run within your activated Conda environment.

```bash
conda install -c conda-forge mpi4py
```

--------------------------------

### Define High-Resolution Wavelength Grid and Read Opacities

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/high_res.ipynb

Initializes POSEIDON's opacity database at native resolution for line-by-line calculations and defines the wavelength grid for optical wavelengths. Imports necessary modules from POSEIDON.

```python
from POSEIDON.core import wl_grid_line_by_line
from POSEIDON.core import read_opacities

#***** Wavelength grid *****#
 
wl_min = 0.3      # Minimum wavelength (um)
wl_max = 1.0      # Maximum wavelength (um)    

wl_high_res = wl_grid_line_by_line(wl_min, wl_max)

#***** Read opacity data *****#

opacity_treatment = 'line_by_line'

# Prepare high-resolution opacity data
opac_high_res = read_opacities(model, wl_high_res, opacity_treatment)
```

--------------------------------

### Install POSEIDON Package

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/installation.md

Navigates into the cloned POSEIDON directory and installs the package in editable mode using pip. This command should be run after cloning the repository.

```bash
cd POSEIDON
pip install -e .
```

--------------------------------

### Set Data Directory

Source: https://github.com/martiancolonist/poseidon/blob/main/tests/notebook_tests.ipynb

Define the directory where data files are stored. This is a common setup step for data loading and processing.

```python
data_dir = './data/WASP-121b'
```

--------------------------------

### Initialize PT and Log X Parameters

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/transmission_clouds.ipynb

Sets up the temperature (T) and log of water abundance (log_H2O) for atmospheric modeling. These parameters are crucial for defining the atmospheric state.

```python
from POSEIDON.core import make_atmosphere

T = 1200
log_H2O = -4

PT_params = np.array([T])
log_X_params = np.array([log_H2O])
```

--------------------------------

### Clone POSEIDON Repository

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/installation.md

Downloads the POSEIDON source code from its GitHub repository. This command should be run in the directory where you want to install POSEIDON.

```bash
git clone https://github.com/MartianColonist/POSEIDON.git
```

--------------------------------

### Initialize POSEIDON Model Components

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/reflection_emission_advanced_aerosols.ipynb

Sets up the wavelength grid, planet, and star objects for a POSEIDON model. Ensure the aerosol database is updated for v1.3.1.

```python
from POSEIDON.constants import R_Sun, R_J, M_J
from POSEIDON.core import create_star, create_planet, load_data, define_model, \
                          wl_grid_constant_R, set_priors, read_opacities
from POSEIDON.visuals import plot_data, plot_spectra_retrieved, plot_PT_retrieved
from POSEIDON.retrieval import run_retrieval
from POSEIDON.utility import read_retrieved_spectrum, read_retrieved_PT, \
                             read_retrieved_log_X, plot_collection
from POSEIDON.corner import generate_cornerplot

import numpy as np
from scipy.constants import au
from scipy.constants import parsec as pc

#***** Model wavelength grid *****#

wl_min = 0.2     # Minimum wavelength (um)
wl_max = 13     # Maximum wavelength (um)
R = 10000         # Spectral resolution of grid

# We need to provide a model wavelength grid to initialise instrument properties
wl = wl_grid_constant_R(wl_min, wl_max, R)

#***** Define planet properties *****#

planet_name = 'HD 189733b'  # Planet name used for plots, output files etc.

R_p = 1.13*R_J    # Planetary radius (m)
M_p = 1.129*M_J
d = 19.7638*pc
a_p = 0.03142*au

# Create the planet object
planet = create_planet(planet_name, R_p, mass = M_p, a_p = a_p)


#***** Define stellar properties *****#

R_s = 0.78*R_Sun      # Stellar radius (m)
T_s = 5014         # Stellar effective temperature (K)
Met_s = 0.13          # Stellar metallicity [log10(Fe/H_star / Fe/H_solar)] <--- note: for PHOENIX, only the solar metallicity models are used  
log_g_s = 4.58       # 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, stellar_grid = 'phoenix')

#***** Define models *****#

model_name_MgSiO3_SiO2_overlapping = 'HD189-Emission-and-Reflection-MgSiO3-SiO2-Overlapping'

bulk_species = ['H2', 'He']
param_species = ['CO', 'CO2', 'H2O', 'K','Na',]
aerosol_species = ['MgSiO3','SiO2']

```

--------------------------------

### Initialize Atmosphere and Wavelength Grid

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/brown_dwarf.ipynb

Imports necessary functions for creating an atmospheric model and defining a wavelength grid with constant spectral resolution.

```python
from POSEIDON.core import make_atmosphere
from POSEIDON.core import wl_grid_constant_R
```

--------------------------------

### List Available Aerosol Species

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/aerosol_database.ipynb

Use `aerosol_supported_species` to get a list of all aerosol species available in the database. This is useful for checking compatibility before querying.

```python
from POSEIDON.supported_chemicals import aerosol_supported_species

print(aerosol_supported_species)
```

--------------------------------

### Configure Forward Model Wavelength Grid and Planet/Star Properties

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/PT_profiles.ipynb

Sets up the wavelength grid, planet, and star properties required for configuring a forward model in POSEIDON. Ensure correct units and values for planetary/stellar parameters.

```python
from POSEIDON.constants import R_Sun, R_J, M_J
from POSEIDON.core import create_star, create_planet, define_model, wl_grid_constant_R
import numpy as np

#***** Model wavelength grid *****#

wl_min = 0.2    # Minimum wavelength (um)
wl_max = 30     # Maximum wavelength (um)
R = 5000        # Spectral resolution of grid

# We need to provide a model wavelength grid to initialise instrument properties
wl_1 = np.linspace(wl_min, 1.0, 1000)[:-1]
wl_2 = wl_grid_constant_R(1.0, wl_max, R)
wl_2 = wl_2[1:]    # Indexing to avoid 1.0 um being repeated twice
wl = np.concatenate((wl_1, wl_2))

#***** Define planet properties *****#

planet_name = 'HD 189733b'  # Planet name used for plots, output files etc.

R_p = 1.13*R_J    # Planetary radius (m)
M_p = 1.129*M_J   # Planetary mass (kg)
T_eq = 1200       # Equilibrium temperature (K)

planet = create_planet(planet_name, R_p, mass = M_p, T_eq = T_eq)

#***** Define stellar properties *****#

R_s = 0.78*R_Sun   # Stellar radius (m)
T_s = 5014         # Stellar effective temperature (K)
Met_s = 0.13       # Stellar metallicity [log10(Fe/H_star / Fe/H_solar)] <--- note: for PHOENIX, only the solar metallicity models are used  
log_g_s = 4.58     # 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, stellar_grid = 'phoenix')
```

--------------------------------

### Get Retrieved Atmosphere

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/advanced_visuals.ipynb

Generates the atmosphere object for the best-fitting set of retrieved parameters using the specified planet, model, and reference pressure.

```python
from POSEIDON.retrieval import get_retrieved_atmosphere

atmosphere_median = get_retrieved_atmosphere(planet, model_2, P, P_ref_set = P_ref)
```

--------------------------------

### Define RV and CCF Parameters

Source: https://context7.com/martiancolonist/poseidon/llms.txt

Defines the range for radial velocity (RV) and cross-correlation function (CCF) parameters. This setup is necessary before performing cross-correlation.

```python
Kp_range   = np.linspace(0, 250, 126)   # km/s
Vsys_range = np.linspace(-100, 100, 201)  # km/s
phi        = data_hr['phi']              # orbital phases

RV_range = get_RV_range(Kp_range, Vsys_range, phi)
```

--------------------------------

### Set Up Fine Grids for Opacity and Atmosphere

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/reflection_emission_surfaces.ipynb

Configures fine temperature and pressure grids for opacity calculations and specifies atmospheric pressure grid settings. This is crucial for accurate modeling, especially with thin atmospheres.

```python
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 

```

--------------------------------

### Adjust Aspect Ratio with subplot_mosaic

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/advanced_visuals.ipynb

Modify the aspect ratio of a subplot by repeating symbols in the subplot_mosaic layout string. This example stretches the 'A' subplot to be taller.

```python
fig_combined = plt.figure(constrained_layout=True, figsize=(12, 7))

# We'll add two more rows of 'A' to alter the aspect ratio
axd = fig_combined.subplot_mosaic(
    """
    AAAAA
    AAAAA
    AAAAA
    abcde
    """
)
```

--------------------------------

### Create Atmospheric Model

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/emission_high_res_cross_correlate.ipynb

Sets up atmospheric parameters including pressure grid, reference pressure and radius, and molecular species. This is used to create a model atmosphere for spectral synthesis.

```python
P_min = 1.0e-5  # 0.1 ubar
P_max = 100     # 100 bar
N_layers = 100  # 100 layers

P = np.logspace(np.log10(P_max), np.log10(P_min), N_layers)

P_ref = 1e-2   # Reference pressure (bar)
R_p_ref = R_p  # Radius at reference pressure

log_species = [-4, -4]

PT_params = np.array(
    [0.2, 0.1, 0.17, -1.39, 1, 1500]  # a1, a2, log_P1, log_P2, log_P3, T_top
)
log_X_params = np.array([log_species])

atmosphere = make_atmosphere(planet, model, P, P_ref, R_p_ref, PT_params,
                             log_X_params, P_param_set=1e-5)
```

--------------------------------

### Create Atmosphere with Line Model

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/PT_profiles.ipynb

Sets up the atmospheric pressure grid, reference pressure, and free parameters for the 'Line' PT profile. Then, it constructs the atmosphere object using the specified planet, model, and parameters.

```python
# Atmospheric pressure grid
P_min = 1.0e-6   # 1 ubar
P_max = 100      # 100 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 = 1   # The R_p_ref parameter will be the radius at 1 bar

# Free parameters
R_p_ref = 1.12 * R_J

log_kappa_IR =   -4.69 
log_gamma    =   -0.31 
log_gamma_2  =   -1.37 
alpha_Line   =   0.11 
beta_Line    =   1.13 
T_int        =   258.8 

PT_params = np.array([log_kappa_IR, log_gamma, log_gamma_2, alpha_Line, beta_Line, T_int])
log_X_params = np.array([])
cloud_params = np.array([])

atmosphere_Line = make_atmosphere(planet, model_Line, P, P_ref, R_p_ref, 
                                  PT_params, log_X_params, cloud_params)
```

--------------------------------

### Run Atmospheric Retrieval

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/retrieval_offsets.ipynb

Execute the atmospheric retrieval process using the specified model, data, priors, and settings. This example uses the MultiNest sampling algorithm.

```python
#**** Run atmospheric retrieval *****#

run_retrieval(planet, star, model_offset_stis, opac, data_offset_stis, priors, 
              wl, P, P_ref, R = R, spectrum_type = 'transmission', 
              sampling_algorithm = 'MultiNest', N_live = 400, verbose = False)
```

--------------------------------

### MultiNest Retrieval Output

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/retrieval_basic.ipynb

Example output from the MultiNest sampling algorithm during an atmospheric retrieval. Displays progress, acceptance rates, and nested sampling log-likelihood values.

```text
POSEIDON now running 'My_first_retrieval'
 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  400
 dimensionality =    3
 *****************************************************
 Starting MultiNest
 generating live points
 live points generated, starting sampling
Acceptance Rate:				0.980392
Replacements:					450
Total Samples:					459
Nested Sampling ln(Z):			-134865.155519
Acceptance Rate:				0.956023
Replacements:					500
Total Samples:					523
Nested Sampling ln(Z):			 -75434.767207
Acceptance Rate:				0.929054
Replacements:					550
Total Samples:					592
Nested Sampling ln(Z):			 -61176.237528
Acceptance Rate:				0.891530
Replacements:					600
Total Samples:					673
Nested Sampling ln(Z):			 -46455.737724
Acceptance Rate:				0.874832
Replacements:					650
Total Samples:					743
Nested Sampling ln(Z):			 -38165.846938
Acceptance Rate:				0.848485
Replacements:					700
Total Samples:					825
Nested Sampling ln(Z):			 -30172.415474
Acceptance Rate:				0.818777
Replacements:					750
Total Samples:					916
Nested Sampling ln(Z):			 -23296.754964
Acceptance Rate:				0.815494
Replacements:					800
Total Samples:					981
Nested Sampling ln(Z):			 -18983.452314
Acceptance Rate:				0.810296
Replacements:					850
Total Samples:				     1049
Nested Sampling ln(Z):			 -15336.226292
Acceptance Rate:				0.800000
Replacements:					900
Total Samples:				     1125
Nested Sampling ln(Z):			 -12064.292982
Acceptance Rate:				0.792327
Replacements:					950
Total Samples:				     1199
Nested Sampling ln(Z):			  -8998.262038
Acceptance Rate:				0.786782
Replacements:				     1000
Total Samples:				     1271
Nested Sampling ln(Z):			  -7079.140574
Acceptance Rate:				0.784167
Replacements:				     1050
Total Samples:				     1339
Nested Sampling ln(Z):			  -5655.993029
Acceptance Rate:				0.780696
Replacements:				     1100
Total Samples:				     1409
Nested Sampling ln(Z):			  -4347.502418
Acceptance Rate:				0.777027
Replacements:				     1150
Total Samples:				     1480
Nested Sampling ln(Z):			  -3140.551689
Acceptance Rate:				0.768738
Replacements:				     1200
Total Samples:				     1561
Nested Sampling ln(Z):			  -2492.241604
Acceptance Rate:				0.762660
Replacements:				     1250
Total Samples:				     1639
Nested Sampling ln(Z):			  -1893.721955
Acceptance Rate:				0.760679
Replacements:				     1300
Total Samples:				     1709
Nested Sampling ln(Z):			  -1452.120895
Acceptance Rate:				0.758853
Replacements:				     1350
Total Samples:				     1779
Nested Sampling ln(Z):			  -1064.269377
Acceptance Rate:				0.755124
Replacements:				     1400
Total Samples:				     1854
Nested Sampling ln(Z):			   -804.681885
Acceptance Rate:				0.753638
Replacements:				     1450
Total Samples:				     1924
Nested Sampling ln(Z):			   -597.203188
Acceptance Rate:				0.751127
Replacements:				     1500
Total Samples:				     1997

```

--------------------------------

### Initialize Spectra List for Plotting

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/reflection_emission_advanced_aerosols.ipynb

Initializes an empty list to store spectra data before plotting. This is a common setup step when using POSEIDON's visualization utilities.

```python
from POSEIDON.visuals import plot_spectra
from POSEIDON.utility import plot_collection

spectra = []
```

--------------------------------

### Set Up POSEIDON Atmospheric Model

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/JWST_proposal.ipynb

Configures the atmospheric model, including bulk and parameterized species, pressure-temperature profile, and cloud treatment. It also defines the atmospheric grid and opacity treatment.

```python
from POSEIDON.core import define_model, read_opacities, \
                          make_atmosphere, compute_spectrum

#***** Define model *****#

model_name = 'CH4_Search_With_JWST'

bulk_species = ['H2', 'He']
param_species = ['H2O', 'CH4']   # Only H2O and CH4 in this model

# Create the model object
model = define_model(model_name, bulk_species, param_species, 
                     PT_profile = 'isotherm', cloud_model = 'cloud-free')

# Specify the pressure grid of the atmosphere
P_min = 1.0e-7    # 1 mbar
P_max = 100       # 100 bar
N_layers = 100    # 100 layers

P = np.logspace(np.log10(P_max), np.log10(P_min), N_layers)

# Specify the reference pressure and radius
P_ref = 1.0e-2        # 10 mbar
R_p_ref = R_p         # Radius at reference pressure

# Provide a specific set of model parameters for the atmosphere
PT_params = np.array([T_eq])              # Assume terminator temperature is T_eq
log_X_params = np.array([-3.62, -7.46])   # H2O and CH4 log abundances in model

# Generate the atmosphere
atmosphere = make_atmosphere(planet, model, P, P_ref, R_p_ref, 
                             PT_params, log_X_params)

opacity_treatment = 'opacity_sampling'

# Define fine temperature grid (K)
T_fine_min = 400     # Same as prior range for T
T_fine_max = 2000    # Same as prior range for T
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 = 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)

# Pre-interpolate the opacities (note: model wavelength range was initialised in cell above)
opac = read_opacities(model, wl, opacity_treatment, T_fine, log_P_fine)
```

--------------------------------

### Define Fine Pressure Grid and Read Opacities

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/transmission_clouds.ipynb

Sets up a fine pressure grid for opacity calculations and then reads opacities, including those for SiO2 aerosols. Ensure T_fine and wl are defined prior to this step.

```python
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)

# Now we can pre-interpolate the sampled opacities (may take up to a minute)
opac_sio2 = read_opacities(model_fuzzy_deck, wl, opacity_treatment, T_fine, log_P_fine)
```

--------------------------------

### Plot Retrieved Transmission Spectrum

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/retrieval_equilibrium.ipynb

Reads and plots the retrieved transmission spectrum, including confidence regions. Requires prior setup of planet and model names, and data objects.

```python
from POSEIDON.utility import read_retrieved_spectrum, plot_collection
from POSEIDON.visuals import plot_spectra_retrieved
from POSEIDON.corner import generate_cornerplot

# Read retrieved spectrum confidence regions
wl, spec_low2, spec_low1, spec_median, \
spec_high1, spec_high2 = read_retrieved_spectrum(planet_name, model_name)

# Create composite spectra objects for plotting
spectra_median = plot_collection(spec_median, wl, collection = [])
spectra_low1 = plot_collection(spec_low1, wl, collection = []) 
spectra_low2 = plot_collection(spec_low2, wl, collection = []) 
spectra_high1 = plot_collection(spec_high1, wl, collection = []) 
spectra_high2 = plot_collection(spec_high2, wl, collection = [])

# Produce figure
fig_spec = plot_spectra_retrieved(spectra_median, spectra_low2, spectra_low1, 
                                  spectra_high1, spectra_high2, planet_name,
                                  data, R_to_bin = 200,
                                  data_colour_list = ['deepskyblue', 'lime', 'orange', 'crimson'],
                                  data_labels = ['NIRISS SOSS Ord2', 'NIRISS SOSS Ord1',
                                                 'NIRSpec G395H NRS1', 'NIRSpec G395H NRS2'],
                                  data_marker_list = ['o', 'o', 'o', 'o'],
                                  legend_location = 'upper right',
                                  legend_box = False,
                                  y_min = 1.26e-2, y_max = 1.54e-2,
                                  plt_label = 'Equilibrium Retrieval')
```

--------------------------------

### Set up priors for retrieval

Source: https://github.com/martiancolonist/poseidon/blob/main/tests/notebook_tests.ipynb

Prepare prior objects for atmospheric retrieval. This involves defining parameters like planet, star, model, data, and prior types/ranges.

```python
priors = set_priors(planet, star, model, data, prior_types, prior_ranges)
```

--------------------------------

### Setting up Figure and Axes for Spectral Plots

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/reflection_emission_surfaces.ipynb

Initializes a matplotlib figure and a subplot mosaic for displaying spectral data. This setup is required before calling plotting functions like plot_spectra.

```python
# 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
Full Spectra'),
                        figure_shape = 'wide',
                        y_unit = 'eclipse_depth',
                        wl_axis = 'linear',
                        line_width_list = [3,1],
                        ax = axd['B'])  
```

--------------------------------

### Define Atmospheric Parameters and Generate Atmosphere

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/transmission_basic.ipynb

Sets up pressure layers and atmospheric parameters to generate an atmospheric model. Ensure 'planet', 'star', and 'model' objects are defined prior to use.

```python
P_min = 1.0e-7    # 0.1 ubar
P_max = 100       # 100 bar
N_layers = 100    # 100 layers

# We'll space the layers uniformly in log-pressure
P = np.logspace(np.log10(P_max), np.log10(P_min), N_layers)

# Specify the reference pressure and radius
P_ref = 10.0   # Reference pressure (bar)
R_p_ref = R_p  # Radius at reference pressure

# Provide a specific set of model parameters for the atmosphere 
PT_params = np.array([1000])              # T (K)
log_X_params = np.array([-3.3, -5.0])     # log(H2O), log(CH4)

# Generate the atmosphere
atmosphere = make_atmosphere(planet, model, P, P_ref, R_p_ref, 
                             PT_params, log_X_params)
```

--------------------------------

### Parallelize Retrieval with mpirun

Source: https://github.com/martiancolonist/poseidon/blob/main/docs/content/notebooks/transmission_high_res_retrieval.ipynb

Example command to run a Python script in parallel using mpirun for faster high-resolution retrievals. This is useful for leveraging multiple CPU cores on a system.

```bash
mpirun -n 24 python -u YOUR_RETRIEVAL_SCRIPT.py
```