Note
Go to the end to download the full example code.
Observe a transit of a spotted star.#
This example demonstrates stellar contamination of a transit spectrum.
A transit can change the spot coverage of a star and produce a signal that is difficult to distinguish from atmospheric absorption. We aim to simulate data from Moran et al. [2023].
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
from astropy import units as u
import libpypsg
from VSPEC import ObservationModel,PhaseAnalyzer
from VSPEC import params
from VSPEC.config import MSH
SEED = 10
PRE_TRANSIT = 8
IN_TRANSIT = 13
libpypsg.docker.set_url_and_run()
Saved settings to /home/runner/.libpypsg/settings.json
Reloading settings...
Create the needed configurations#
Moran et al. [2023] observed super-Earth GJ 486b with JWST NIRSpec/G395H
# Instrument
inst = params.InstrumentParameters(
telescope=params.SingleDishParameters.jwst(),
bandpass=params.BandpassParameters(
wl_blue=2.87*u.um,
wl_red=5.14*u.um,
resolving_power=200,
wavelength_unit=u.um,
flux_unit=u.Unit('W m-2 um-1')
),
detector=params.DetectorParameters(
beam_width=5*u.arcsec,
integration_time=0.5*u.s,
ccd=params.ccdParameters(
pixel_sampling=8,
read_noise=16.8*u.electron,
dark_current=0.005*u.electron/u.s,
throughput=0.3,
emissivity=0.1,
temperature=50*u.K
)
)
)
# Observation
observation = params.ObservationParameters(
observation_time=3.53*u.hr,
integration_time=8*u.min
)
# PSG
psg_kwargs = dict(
gcm_binning=200,
phase_binning=1,
nmax=0,
lmax=0,
use_continuum_stellar=True,
continuum=['Rayleigh', 'Refraction'],
)
psg_params = params.psgParameters(
use_molecular_signatures=True,
**psg_kwargs
)
psg_no_atm = params.psgParameters(
use_molecular_signatures=False,
**psg_kwargs
)
# Star and Planet
star_teff = 3343*u.K
star_rad = 0.339*u.R_sun
a_rstar = 11.229 # Moran+23 Table 1
rp_rstar = 0.03709 # Moran+23 Table 1
planet_rad = star_rad * rp_rstar
orbit_rad = star_rad * a_rstar
orbit_period = 1.467119*u.day
planet_rot_period = orbit_period
star_rot_period = np.inf*u.s # assume the star does not change.
planet_mass = 2.82*u.M_earth
star_mass = 0.323*u.M_sun
inclination = 89.06*u.deg
planet_norm_factor = {
'rock_quiet': 1,
'rock_spotted': 0.97,
'water_quiet': 0.98,
} # We will multiply the radius by these scalars rather than fit
# later (because I have run this code before and know how much
# each model is off by).
print(f'Planet radius: {planet_rad.to_value(u.R_earth)} R_earth')
print(f'Semimajor axis: {orbit_rad.to_value(u.AU)} AU')
observation_angle = (2*np.pi*u.rad * observation.observation_time/orbit_period).to(u.deg)
initial_phase = 180*u.deg - 0.5*observation_angle
planet_kwargs = dict(
name='GJ486b',
gravity=params.GravityParameters('kg',planet_mass),
semimajor_axis=orbit_rad,
orbit_period=orbit_period,
rotation_period=planet_rot_period,
eccentricity=0,
obliquity=0*u.deg,
obliquity_direction=0*u.deg,
init_phase=initial_phase,
init_substellar_lon=0*u.deg
)
quiet_rock_planet = params.PlanetParameters(
radius=planet_rad*planet_norm_factor['rock_quiet'],
**planet_kwargs
)
quiet_water_planet = params.PlanetParameters(
radius=planet_rad*planet_norm_factor['water_quiet'],
**planet_kwargs
)
spotted_rock_planet = params.PlanetParameters(
radius=planet_rad*planet_norm_factor['rock_spotted'],
**planet_kwargs
)
system_params = params.SystemParameters(
distance=8.07*u.pc,
inclination=inclination,
phase_of_periastron=0*u.deg
)
star_dict = {
'teff': star_teff,
'radius': star_rad
}
planet_dict = {'semimajor_axis': orbit_rad}
gcm_dict = {
'nlayer': 40,
'nlon': 60,
'nlat': 30,
'epsilon': 100,
'albedo': 0.3,
'emissivity': 1.0,
'lat_redistribution': 0.0,
'gamma': 1.0,
'psurf': 1*u.bar,
'ptop': 1e-10*u.bar,
'wind': {'U': '0 m/s','V':'0 m/s'},
'molecules':{'H2O':0.99}
}
# Create two sets of GCM Parameters
h2o_atm = {'molecules':{'H2O':0.99}}
gcm_h2o = params.gcmParameters.from_dict({
'star':star_dict,
'planet':planet_dict,
'gcm':{'vspec':dict(gcm_dict,**h2o_atm),'mean_molec_weight':18}
})
star_kwargs = dict(
psg_star_template='M',
teff=star_teff,
mass=star_mass,
radius=star_rad,
period=star_rot_period,
misalignment=0*u.deg,
misalignment_dir=0*u.deg,
ld=params.LimbDarkeningParameters.lambertian(),
faculae=params.FaculaParameters.none(),
flares=params.FlareParameters.none(),
granulation=params.GranulationParameters.none(),
grid_params=100000,
)
quiet_star = params.StarParameters(
spots=params.SpotParameters.none(),
**star_kwargs
)
spotted_star = params.StarParameters(
spots=params.SpotParameters(
distribution='iso',
initial_coverage=0.11,
area_mean=500*MSH,
area_logsigma=0.2,
teff_umbra=2700*u.K,
teff_penumbra=2700*u.K,
equillibrium_coverage=0.0,
burn_in=0*u.s,
growth_rate=0.0/u.day,
decay_rate=0*MSH/u.day,
initial_area=10*MSH
),
**star_kwargs
)
# Set parameters for simulation
header_kwargs = dict(
seed = SEED,
spec_grid = params.VSPECGridParameters(
max_teff=3400*u.K,min_teff=2300*u.K,
impl_bin='rust',impl_interp='scipy',fail_on_missing=False
)
)
internal_params_kwargs = dict(
system=system_params,
obs=observation,
gcm=gcm_h2o,
inst=inst
)
# Make the three cases
params_rock_quiet = params.InternalParameters(
header=params.Header(data_path=Path('.vspec/rock_quiet'),**header_kwargs),
star = quiet_star,
psg=psg_no_atm,
planet=quiet_rock_planet,
**internal_params_kwargs
)
params_h2o_quiet = params.InternalParameters(
header=params.Header(data_path=Path('.vspec/h2o_quiet'),**header_kwargs),
star = quiet_star,
psg=psg_params,
planet=quiet_water_planet,
**internal_params_kwargs
)
params_rock_spotted = params.InternalParameters(
header=params.Header(data_path=Path('.vspec/rock_spotted'),**header_kwargs),
star = spotted_star,
psg=psg_no_atm,
planet=spotted_rock_planet,
**internal_params_kwargs
)
Planet radius: 1.371472837835719 R_earth
Semimajor axis: 0.017702612840063636 AU
Run VSPEC for the simplest case#
We read in the config file and run the model.
model_rock_quiet = ObservationModel(params_rock_quiet)
model_rock_quiet.build_planet()
model_rock_quiet.build_spectra()
Starting at phase 161.9544290544939 deg, observe for 3.53 h in 26 steps
Phases = [161.95 163.32 164.68 166.04 167.41 168.77 170.13 171.5 172.86 174.22
175.59 176.95 178.31 179.68 181.04 182.4 183.77 185.13 186.49 187.86
189.22 190.58 191.95 193.31 194.67 196.03 197.4 ] deg
Build Planet: 0%| | 0/27 [00:00<?, ?it/s]
Build Planet: 4%|▎ | 1/27 [00:00<00:21, 1.20it/s]
Build Planet: 7%|▋ | 2/27 [00:01<00:20, 1.20it/s]
Build Planet: 11%|█ | 3/27 [00:02<00:19, 1.22it/s]
Build Planet: 15%|█▍ | 4/27 [00:03<00:18, 1.22it/s]
Build Planet: 19%|█▊ | 5/27 [00:04<00:18, 1.22it/s]
Build Planet: 22%|██▏ | 6/27 [00:04<00:17, 1.23it/s]
Build Planet: 26%|██▌ | 7/27 [00:05<00:16, 1.22it/s]
Build Planet: 30%|██▉ | 8/27 [00:06<00:15, 1.23it/s]
Build Planet: 33%|███▎ | 9/27 [00:07<00:14, 1.22it/s]
Build Planet: 37%|███▋ | 10/27 [00:08<00:13, 1.22it/s]
Build Planet: 41%|████ | 11/27 [00:09<00:13, 1.22it/s]
Build Planet: 44%|████▍ | 12/27 [00:09<00:12, 1.22it/s]
Build Planet: 48%|████▊ | 13/27 [00:10<00:11, 1.22it/s]
Build Planet: 52%|█████▏ | 14/27 [00:11<00:10, 1.22it/s]
Build Planet: 56%|█████▌ | 15/27 [00:12<00:09, 1.23it/s]
Build Planet: 59%|█████▉ | 16/27 [00:13<00:09, 1.22it/s]
Build Planet: 63%|██████▎ | 17/27 [00:13<00:08, 1.21it/s]
Build Planet: 67%|██████▋ | 18/27 [00:14<00:07, 1.22it/s]
Build Planet: 70%|███████ | 19/27 [00:15<00:06, 1.22it/s]
Build Planet: 74%|███████▍ | 20/27 [00:16<00:05, 1.22it/s]
Build Planet: 78%|███████▊ | 21/27 [00:17<00:04, 1.22it/s]
Build Planet: 81%|████████▏ | 22/27 [00:18<00:04, 1.22it/s]
Build Planet: 85%|████████▌ | 23/27 [00:18<00:03, 1.23it/s]
Build Planet: 89%|████████▉ | 24/27 [00:19<00:02, 1.22it/s]
Build Planet: 93%|█████████▎| 25/27 [00:20<00:01, 1.23it/s]
Build Planet: 96%|█████████▋| 26/27 [00:21<00:00, 1.22it/s]
Build Planet: 100%|██████████| 27/27 [00:22<00:00, 1.22it/s]
Build Planet: 100%|██████████| 27/27 [00:22<00:00, 1.22it/s]
Creating interpolators:
thermal
thermal, combined
thermal, combined, stellar
thermal, combined, stellar, photon noise
thermal, combined, stellar, photon noise, detector noise
thermal, combined, stellar, photon noise, detector noise, telescope noise
thermal, combined, stellar, photon noise, detector noise, telescope noise, background noise
thermal, combined, stellar, photon noise, detector noise, telescope noise, background noise, transit
Finished!
Build Spectra: 0%| | 0/26 [00:00<?, ?it/s]
Loading Spectra: 0%| | 0/12 [00:00<?, ?it/s]
Loading Spectra: 8%|▊ | 1/12 [00:00<00:03, 3.39it/s]
Loading Spectra: 17%|█▋ | 2/12 [00:00<00:02, 3.38it/s]
Loading Spectra: 25%|██▌ | 3/12 [00:00<00:02, 3.38it/s]
Loading Spectra: 33%|███▎ | 4/12 [00:01<00:02, 3.37it/s]
Loading Spectra: 42%|████▏ | 5/12 [00:01<00:02, 3.37it/s]
Loading Spectra: 50%|█████ | 6/12 [00:01<00:01, 3.38it/s]
Loading Spectra: 58%|█████▊ | 7/12 [00:02<00:01, 3.38it/s]
Loading Spectra: 67%|██████▋ | 8/12 [00:02<00:01, 3.39it/s]
Loading Spectra: 75%|███████▌ | 9/12 [00:02<00:00, 3.39it/s]
Loading Spectra: 83%|████████▎ | 10/12 [00:02<00:00, 3.39it/s]
Loading Spectra: 92%|█████████▏| 11/12 [00:03<00:00, 3.38it/s]
Loading Spectra: 100%|██████████| 12/12 [00:03<00:00, 3.38it/s]
Loading Spectra: 100%|██████████| 12/12 [00:03<00:00, 3.38it/s]
Build Spectra: 4%|▍ | 1/26 [00:03<01:33, 3.72s/it]
Build Spectra: 8%|▊ | 2/26 [00:03<00:39, 1.63s/it]
Build Spectra: 12%|█▏ | 3/26 [00:04<00:22, 1.04it/s]
Build Spectra: 15%|█▌ | 4/26 [00:04<00:14, 1.53it/s]
Build Spectra: 19%|█▉ | 5/26 [00:04<00:10, 2.08it/s]
Build Spectra: 23%|██▎ | 6/26 [00:04<00:07, 2.67it/s]
Build Spectra: 27%|██▋ | 7/26 [00:04<00:05, 3.25it/s]
Build Spectra: 31%|███ | 8/26 [00:04<00:04, 3.79it/s]
Build Spectra: 35%|███▍ | 9/26 [00:05<00:03, 4.27it/s]
Build Spectra: 38%|███▊ | 10/26 [00:05<00:03, 4.65it/s]
Build Spectra: 42%|████▏ | 11/26 [00:05<00:03, 4.36it/s]
Build Spectra: 46%|████▌ | 12/26 [00:05<00:03, 4.20it/s]
Build Spectra: 50%|█████ | 13/26 [00:06<00:03, 4.08it/s]
Build Spectra: 54%|█████▍ | 14/26 [00:06<00:03, 4.00it/s]
Build Spectra: 58%|█████▊ | 15/26 [00:06<00:02, 3.95it/s]
Build Spectra: 62%|██████▏ | 16/26 [00:06<00:02, 3.93it/s]
Build Spectra: 65%|██████▌ | 17/26 [00:07<00:02, 3.91it/s]
Build Spectra: 69%|██████▉ | 18/26 [00:07<00:02, 3.96it/s]
Build Spectra: 73%|███████▎ | 19/26 [00:07<00:01, 4.39it/s]
Build Spectra: 77%|███████▋ | 20/26 [00:07<00:01, 4.75it/s]
Build Spectra: 81%|████████ | 21/26 [00:07<00:00, 5.02it/s]
Build Spectra: 85%|████████▍ | 22/26 [00:08<00:00, 5.24it/s]
Build Spectra: 88%|████████▊ | 23/26 [00:08<00:00, 5.41it/s]
Build Spectra: 92%|█████████▏| 24/26 [00:08<00:00, 5.49it/s]
Build Spectra: 96%|█████████▌| 25/26 [00:08<00:00, 5.60it/s]
Build Spectra: 100%|██████████| 26/26 [00:08<00:00, 5.68it/s]
Build Spectra: 100%|██████████| 26/26 [00:08<00:00, 2.99it/s]
Load in the data#
We can use VSPEC to read in the synthetic data we just created.
data_rock_quiet = PhaseAnalyzer(model_rock_quiet.directories['all_model'])
/opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/site-packages/VSPEC/analysis.py:116: RuntimeWarning: No Layer info, maybe globes or molecular signatures are off
warnings.warn(
Plot the transit#
Lets plot the lightcurve of each channel.
def plot_transit(data:PhaseAnalyzer,title:str,color:str):
time_from_mid_transit = (data.time-0.5*(data.time[-1]+data.time[0])).to(u.hour)
fig,axes = plt.subplots(2,1,tight_layout=True)
axes[0].scatter(
time_from_mid_transit,
data.lightcurve('total',(0,-1),normalize=0,noise=False),
label = 'white light curve',c=color
)
axes[0].set_xlabel('Time past mid-transit (hour)')
axes[0].set_ylabel('Transit depth (ppm)')
axes[0].legend()
axes[0].set_title(title)
# standardize the epochs to use for analysis
pre_transit = PRE_TRANSIT
in_transit = IN_TRANSIT
unocculted_spectrum = data.spectrum('total',pre_transit)
occulted_spectrum = data.spectrum('total',in_transit)
lost_to_transit = unocculted_spectrum-occulted_spectrum
transit_depth = (lost_to_transit/unocculted_spectrum).to_value(u.dimensionless_unscaled)
axes[1].plot(data.wavelength, 1e6*(transit_depth),c=color)
axes[1].set_xlabel(f'Wavelength ({data.wavelength.unit})')
axes[1].set_ylabel('Transit depth (ppm)')
ylo,yhi = axes[1].get_ylim()
if yhi-ylo < 0.5:
mean = 0.5*(ylo+yhi)
axes[1].set_ylim(mean-0.25,mean+0.25)
return fig
plot_transit(data_rock_quiet,'Spotless Star and Bare Rock','xkcd:lavender').show()
Run the other models#
Let’s do the same analysis for the other cases.
Spotless Star, H2O Planet#
model_h2o_quiet = ObservationModel(params_h2o_quiet)
model_h2o_quiet.build_planet()
model_h2o_quiet.build_spectra()
data_h2o_quiet = PhaseAnalyzer(model_h2o_quiet.directories['all_model'])
plot_transit(data_h2o_quiet,'Spotless Star and 1 bar H2O Atmosphere','xkcd:azure').show()
Starting at phase 161.9544290544939 deg, observe for 3.53 h in 26 steps
Phases = [161.95 163.32 164.68 166.04 167.41 168.77 170.13 171.5 172.86 174.22
175.59 176.95 178.31 179.68 181.04 182.4 183.77 185.13 186.49 187.86
189.22 190.58 191.95 193.31 194.67 196.03 197.4 ] deg
Build Planet: 0%| | 0/27 [00:00<?, ?it/s]
Build Planet: 4%|▎ | 1/27 [00:01<00:28, 1.11s/it]
Build Planet: 7%|▋ | 2/27 [00:02<00:28, 1.14s/it]
Build Planet: 11%|█ | 3/27 [00:03<00:26, 1.11s/it]
Build Planet: 15%|█▍ | 4/27 [00:04<00:25, 1.11s/it]
Build Planet: 19%|█▊ | 5/27 [00:05<00:24, 1.10s/it]
Build Planet: 22%|██▏ | 6/27 [00:06<00:22, 1.09s/it]
Build Planet: 26%|██▌ | 7/27 [00:07<00:22, 1.12s/it]
Build Planet: 30%|██▉ | 8/27 [00:08<00:21, 1.11s/it]
Build Planet: 33%|███▎ | 9/27 [00:10<00:20, 1.14s/it]
Build Planet: 37%|███▋ | 10/27 [00:11<00:19, 1.12s/it]
Build Planet: 41%|████ | 11/27 [00:12<00:18, 1.15s/it]
Build Planet: 44%|████▍ | 12/27 [00:13<00:17, 1.16s/it]
Build Planet: 48%|████▊ | 13/27 [00:14<00:16, 1.18s/it]
Build Planet: 52%|█████▏ | 14/27 [00:15<00:15, 1.18s/it]
Build Planet: 56%|█████▌ | 15/27 [00:17<00:14, 1.19s/it]
Build Planet: 59%|█████▉ | 16/27 [00:18<00:13, 1.19s/it]
Build Planet: 63%|██████▎ | 17/27 [00:19<00:11, 1.19s/it]
Build Planet: 67%|██████▋ | 18/27 [00:20<00:10, 1.18s/it]
Build Planet: 70%|███████ | 19/27 [00:21<00:09, 1.16s/it]
Build Planet: 74%|███████▍ | 20/27 [00:22<00:07, 1.14s/it]
Build Planet: 78%|███████▊ | 21/27 [00:24<00:06, 1.13s/it]
Build Planet: 81%|████████▏ | 22/27 [00:25<00:05, 1.12s/it]
Build Planet: 85%|████████▌ | 23/27 [00:26<00:04, 1.11s/it]
Build Planet: 89%|████████▉ | 24/27 [00:27<00:03, 1.11s/it]
Build Planet: 93%|█████████▎| 25/27 [00:28<00:02, 1.11s/it]
Build Planet: 96%|█████████▋| 26/27 [00:29<00:01, 1.10s/it]
Build Planet: 100%|██████████| 27/27 [00:30<00:00, 1.11s/it]
Build Planet: 100%|██████████| 27/27 [00:30<00:00, 1.13s/it]
Creating interpolators:
thermal
thermal, combined
thermal, combined, stellar
thermal, combined, stellar, photon noise
thermal, combined, stellar, photon noise, detector noise
thermal, combined, stellar, photon noise, detector noise, telescope noise
thermal, combined, stellar, photon noise, detector noise, telescope noise, background noise
thermal, combined, stellar, photon noise, detector noise, telescope noise, background noise, transit
Finished!
Build Spectra: 0%| | 0/26 [00:00<?, ?it/s]
Loading Spectra: 0%| | 0/12 [00:00<?, ?it/s]
Loading Spectra: 8%|▊ | 1/12 [00:00<00:03, 3.36it/s]
Loading Spectra: 17%|█▋ | 2/12 [00:00<00:02, 3.35it/s]
Loading Spectra: 25%|██▌ | 3/12 [00:00<00:02, 3.36it/s]
Loading Spectra: 33%|███▎ | 4/12 [00:01<00:02, 3.35it/s]
Loading Spectra: 42%|████▏ | 5/12 [00:01<00:02, 3.35it/s]
Loading Spectra: 50%|█████ | 6/12 [00:01<00:01, 3.35it/s]
Loading Spectra: 58%|█████▊ | 7/12 [00:02<00:01, 3.35it/s]
Loading Spectra: 67%|██████▋ | 8/12 [00:02<00:01, 3.35it/s]
Loading Spectra: 75%|███████▌ | 9/12 [00:02<00:00, 3.36it/s]
Loading Spectra: 83%|████████▎ | 10/12 [00:02<00:00, 3.36it/s]
Loading Spectra: 92%|█████████▏| 11/12 [00:03<00:00, 3.36it/s]
Loading Spectra: 100%|██████████| 12/12 [00:03<00:00, 3.35it/s]
Loading Spectra: 100%|██████████| 12/12 [00:03<00:00, 3.35it/s]
Build Spectra: 4%|▍ | 1/26 [00:03<01:35, 3.81s/it]
Build Spectra: 8%|▊ | 2/26 [00:04<00:40, 1.71s/it]
Build Spectra: 12%|█▏ | 3/26 [00:04<00:23, 1.03s/it]
Build Spectra: 15%|█▌ | 4/26 [00:04<00:15, 1.39it/s]
Build Spectra: 19%|█▉ | 5/26 [00:04<00:11, 1.84it/s]
Build Spectra: 23%|██▎ | 6/26 [00:04<00:08, 2.29it/s]
Build Spectra: 27%|██▋ | 7/26 [00:05<00:07, 2.70it/s]
Build Spectra: 31%|███ | 8/26 [00:05<00:05, 3.06it/s]
Build Spectra: 35%|███▍ | 9/26 [00:05<00:05, 3.37it/s]
Build Spectra: 38%|███▊ | 10/26 [00:05<00:04, 3.62it/s]
Build Spectra: 42%|████▏ | 11/26 [00:06<00:04, 3.44it/s]
Build Spectra: 46%|████▌ | 12/26 [00:06<00:04, 3.28it/s]
Build Spectra: 50%|█████ | 13/26 [00:06<00:04, 3.20it/s]
Build Spectra: 54%|█████▍ | 14/26 [00:07<00:03, 3.17it/s]
Build Spectra: 58%|█████▊ | 15/26 [00:07<00:03, 3.16it/s]
Build Spectra: 62%|██████▏ | 16/26 [00:07<00:03, 3.14it/s]
Build Spectra: 65%|██████▌ | 17/26 [00:08<00:02, 3.12it/s]
Build Spectra: 69%|██████▉ | 18/26 [00:08<00:02, 3.16it/s]
Build Spectra: 73%|███████▎ | 19/26 [00:08<00:02, 3.42it/s]
Build Spectra: 77%|███████▋ | 20/26 [00:08<00:01, 3.64it/s]
Build Spectra: 81%|████████ | 21/26 [00:09<00:01, 3.80it/s]
Build Spectra: 85%|████████▍ | 22/26 [00:09<00:01, 3.93it/s]
Build Spectra: 88%|████████▊ | 23/26 [00:09<00:00, 4.02it/s]
Build Spectra: 92%|█████████▏| 24/26 [00:09<00:00, 4.09it/s]
Build Spectra: 96%|█████████▌| 25/26 [00:10<00:00, 4.13it/s]
Build Spectra: 100%|██████████| 26/26 [00:10<00:00, 4.17it/s]
Build Spectra: 100%|██████████| 26/26 [00:10<00:00, 2.51it/s]
Spotted Star, CO2 Planet#
model_rock_spotted = ObservationModel(params_rock_spotted)
model_rock_spotted.build_planet()
model_rock_spotted.build_spectra()
data_rock_spotted = PhaseAnalyzer(model_rock_spotted.directories['all_model'])
plot_transit(data_rock_spotted,'Spotted Star and Bare Rock','xkcd:golden yellow').show()
Starting at phase 161.9544290544939 deg, observe for 3.53 h in 26 steps
Phases = [161.95 163.32 164.68 166.04 167.41 168.77 170.13 171.5 172.86 174.22
175.59 176.95 178.31 179.68 181.04 182.4 183.77 185.13 186.49 187.86
189.22 190.58 191.95 193.31 194.67 196.03 197.4 ] deg
Build Planet: 0%| | 0/27 [00:00<?, ?it/s]
Build Planet: 4%|▎ | 1/27 [00:00<00:21, 1.20it/s]
Build Planet: 7%|▋ | 2/27 [00:01<00:20, 1.20it/s]
Build Planet: 11%|█ | 3/27 [00:02<00:19, 1.22it/s]
Build Planet: 15%|█▍ | 4/27 [00:03<00:18, 1.22it/s]
Build Planet: 19%|█▊ | 5/27 [00:04<00:18, 1.22it/s]
Build Planet: 22%|██▏ | 6/27 [00:04<00:17, 1.22it/s]
Build Planet: 26%|██▌ | 7/27 [00:05<00:16, 1.22it/s]
Build Planet: 30%|██▉ | 8/27 [00:06<00:15, 1.22it/s]
Build Planet: 33%|███▎ | 9/27 [00:07<00:14, 1.22it/s]
Build Planet: 37%|███▋ | 10/27 [00:08<00:13, 1.22it/s]
Build Planet: 41%|████ | 11/27 [00:09<00:13, 1.22it/s]
Build Planet: 44%|████▍ | 12/27 [00:09<00:12, 1.22it/s]
Build Planet: 48%|████▊ | 13/27 [00:10<00:11, 1.21it/s]
Build Planet: 52%|█████▏ | 14/27 [00:11<00:10, 1.21it/s]
Build Planet: 56%|█████▌ | 15/27 [00:12<00:09, 1.22it/s]
Build Planet: 59%|█████▉ | 16/27 [00:13<00:09, 1.22it/s]
Build Planet: 63%|██████▎ | 17/27 [00:13<00:08, 1.22it/s]
Build Planet: 67%|██████▋ | 18/27 [00:14<00:07, 1.22it/s]
Build Planet: 70%|███████ | 19/27 [00:15<00:06, 1.21it/s]
Build Planet: 74%|███████▍ | 20/27 [00:16<00:05, 1.21it/s]
Build Planet: 78%|███████▊ | 21/27 [00:17<00:04, 1.22it/s]
Build Planet: 81%|████████▏ | 22/27 [00:18<00:04, 1.22it/s]
Build Planet: 85%|████████▌ | 23/27 [00:18<00:03, 1.22it/s]
Build Planet: 89%|████████▉ | 24/27 [00:19<00:02, 1.22it/s]
Build Planet: 93%|█████████▎| 25/27 [00:20<00:01, 1.22it/s]
Build Planet: 96%|█████████▋| 26/27 [00:21<00:00, 1.22it/s]
Build Planet: 100%|██████████| 27/27 [00:22<00:00, 1.22it/s]
Build Planet: 100%|██████████| 27/27 [00:22<00:00, 1.22it/s]
Generated 94 mature spots
Creating interpolators:
thermal
thermal, combined
thermal, combined, stellar
thermal, combined, stellar, photon noise
thermal, combined, stellar, photon noise, detector noise
thermal, combined, stellar, photon noise, detector noise, telescope noise
thermal, combined, stellar, photon noise, detector noise, telescope noise, background noise
thermal, combined, stellar, photon noise, detector noise, telescope noise, background noise, transit
Finished!
Build Spectra: 0%| | 0/26 [00:00<?, ?it/s]
Loading Spectra: 0%| | 0/12 [00:00<?, ?it/s]
Loading Spectra: 8%|▊ | 1/12 [00:00<00:03, 3.34it/s]
Loading Spectra: 17%|█▋ | 2/12 [00:00<00:03, 3.32it/s]
Loading Spectra: 25%|██▌ | 3/12 [00:00<00:02, 3.33it/s]
Loading Spectra: 33%|███▎ | 4/12 [00:01<00:02, 3.34it/s]
Loading Spectra: 42%|████▏ | 5/12 [00:01<00:02, 3.33it/s]
Loading Spectra: 50%|█████ | 6/12 [00:01<00:01, 3.34it/s]
Loading Spectra: 58%|█████▊ | 7/12 [00:02<00:01, 3.35it/s]
Loading Spectra: 67%|██████▋ | 8/12 [00:02<00:01, 3.35it/s]
Loading Spectra: 75%|███████▌ | 9/12 [00:02<00:00, 3.35it/s]
Loading Spectra: 83%|████████▎ | 10/12 [00:02<00:00, 3.35it/s]
Loading Spectra: 92%|█████████▏| 11/12 [00:03<00:00, 3.36it/s]
Loading Spectra: 100%|██████████| 12/12 [00:03<00:00, 3.36it/s]
Loading Spectra: 100%|██████████| 12/12 [00:03<00:00, 3.35it/s]
Build Spectra: 4%|▍ | 1/26 [00:03<01:38, 3.93s/it]
Build Spectra: 8%|▊ | 2/26 [00:04<00:43, 1.82s/it]
Build Spectra: 12%|█▏ | 3/26 [00:04<00:26, 1.15s/it]
Build Spectra: 15%|█▌ | 4/26 [00:04<00:18, 1.21it/s]
Build Spectra: 19%|█▉ | 5/26 [00:05<00:13, 1.53it/s]
Build Spectra: 23%|██▎ | 6/26 [00:05<00:10, 1.82it/s]
Build Spectra: 27%|██▋ | 7/26 [00:05<00:09, 2.08it/s]
Build Spectra: 31%|███ | 8/26 [00:06<00:07, 2.29it/s]
Build Spectra: 35%|███▍ | 9/26 [00:06<00:06, 2.45it/s]
Build Spectra: 38%|███▊ | 10/26 [00:07<00:06, 2.58it/s]
Build Spectra: 42%|████▏ | 11/26 [00:07<00:06, 2.48it/s]
Build Spectra: 46%|████▌ | 12/26 [00:07<00:05, 2.41it/s]
Build Spectra: 50%|█████ | 13/26 [00:08<00:05, 2.38it/s]
Build Spectra: 54%|█████▍ | 14/26 [00:08<00:05, 2.37it/s]
Build Spectra: 58%|█████▊ | 15/26 [00:09<00:04, 2.35it/s]
Build Spectra: 62%|██████▏ | 16/26 [00:09<00:04, 2.33it/s]
Build Spectra: 65%|██████▌ | 17/26 [00:10<00:03, 2.33it/s]
Build Spectra: 69%|██████▉ | 18/26 [00:10<00:03, 2.35it/s]
Build Spectra: 73%|███████▎ | 19/26 [00:10<00:02, 2.50it/s]
Build Spectra: 77%|███████▋ | 20/26 [00:11<00:02, 2.61it/s]
Build Spectra: 81%|████████ | 21/26 [00:11<00:01, 2.70it/s]
Build Spectra: 85%|████████▍ | 22/26 [00:11<00:01, 2.75it/s]
Build Spectra: 88%|████████▊ | 23/26 [00:12<00:01, 2.80it/s]
Build Spectra: 92%|█████████▏| 24/26 [00:12<00:00, 2.84it/s]
Build Spectra: 96%|█████████▌| 25/26 [00:12<00:00, 2.86it/s]
Build Spectra: 100%|██████████| 26/26 [00:13<00:00, 2.88it/s]
Build Spectra: 100%|██████████| 26/26 [00:13<00:00, 1.97it/s]
/opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/site-packages/VSPEC/analysis.py:116: RuntimeWarning: No Layer info, maybe globes or molecular signatures are off
warnings.warn(
Plot the observed spectra#
Let’s compare the transits. We also load in the actual JWST data.
Get the data#
Reduced data from Moran et al. [2023] is publicly available. However, you must download it from the figure caption of the online version.
import pandas as pd
try:
filename = Path(__file__).parent / 'moran2023_fig3.txt'
except NameError:
filename = 'moran2023_fig3.txt'
df = pd.read_fwf(filename,colspecs=[(0,8),(9,14),(15,20),(21,25),(26,28)],
header=20,names=['Reduction','Wave','Width','Depth','e_Depth'])
used_eureka = df['Reduction']=='Eureka'
moran_x = df.loc[used_eureka,'Wave']
moran_y = df.loc[used_eureka,'Depth']
moran_yerr = df.loc[used_eureka,'e_Depth']
moran_mean = np.mean(moran_y)
Make the figure#
fig, ax = plt.subplots(1,1)
for data,label,color in zip(
[data_rock_quiet,data_h2o_quiet,data_rock_spotted],
['Rock', 'H2O', 'Rock+Spots'],
['xkcd:lavender','xkcd:azure','xkcd:golden yellow']
):
pre_transit = PRE_TRANSIT
in_transit = IN_TRANSIT
unocculted_spectrum = data.spectrum('total',pre_transit)
occulted_spectrum = data.spectrum('total',in_transit)
lost_to_transit = unocculted_spectrum-occulted_spectrum
transit_depth = (lost_to_transit/unocculted_spectrum).to_value(u.dimensionless_unscaled)*1e6
depth_mean = np.mean(transit_depth)
shift = moran_mean/depth_mean
shift = 1
ax.plot(data.wavelength,transit_depth*shift,label=label,color=color)
ax.errorbar(moran_x,moran_y,yerr=moran_yerr,
fmt='o',color='k',label='Moran+23 Eureka!')
ax.set_xlabel('Wavelength (um)')
ax.set_ylabel('Transit depth (ppm)')
ax.legend()
<matplotlib.legend.Legend object at 0x7fd410499c50>
Make the figure again#
This time without the bare rock
fig, ax = plt.subplots(1,1,figsize=(5,3.5),tight_layout=True)
for data,label,color in zip(
[data_h2o_quiet,data_rock_spotted],
['H2O', 'Rock+Spots'],
['xkcd:azure','xkcd:golden yellow']
):
pre_transit = PRE_TRANSIT
in_transit = IN_TRANSIT
unocculted_spectrum = data.spectrum('total',pre_transit)
occulted_spectrum = data.spectrum('total',in_transit)
lost_to_transit = unocculted_spectrum-occulted_spectrum
transit_depth = (lost_to_transit/unocculted_spectrum).to_value(u.dimensionless_unscaled)
ax.plot(data.wavelength,transit_depth*1e6,label=label,color=color)
ax.errorbar(moran_x,moran_y,yerr=moran_yerr,
fmt='o',color='k',label='Moran+23 Eureka!')
ax.set_xlabel('Wavelength (um)')
ax.set_ylabel('Transit depth (ppm)')
ax.legend()
<matplotlib.legend.Legend object at 0x7fd418530f90>
Total running time of the script: (2 minutes 40.511 seconds)