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 pypsg
from VSPEC import ObservationModel,PhaseAnalyzer
from VSPEC import params
from VSPEC.config import MSH
SEED = 10
PRE_TRANSIT = 8
IN_TRANSIT = 13
pypsg.docker.set_url_and_run()
Saved settings to /home/ted/.pypsg/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,
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_periasteron=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,
spectral_grid='default'
)
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(
teff_min=2300*u.K,teff_max=3400*u.K,
seed = SEED
)
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/26 [00:00<?, ?it/s]
Build Planet: 8%|▊ | 2/26 [00:00<00:02, 10.15it/s]
Build Planet: 15%|█▌ | 4/26 [00:00<00:02, 9.83it/s]
Build Planet: 23%|██▎ | 6/26 [00:00<00:02, 9.91it/s]
Build Planet: 27%|██▋ | 7/26 [00:00<00:01, 9.80it/s]
Build Planet: 31%|███ | 8/26 [00:00<00:01, 9.67it/s]
Build Planet: 35%|███▍ | 9/26 [00:00<00:01, 9.58it/s]
Build Planet: 38%|███▊ | 10/26 [00:01<00:01, 9.48it/s]
Build Planet: 42%|████▏ | 11/26 [00:01<00:01, 8.88it/s]
Build Planet: 46%|████▌ | 12/26 [00:01<00:01, 8.84it/s]
Build Planet: 50%|█████ | 13/26 [00:01<00:01, 8.66it/s]
Build Planet: 54%|█████▍ | 14/26 [00:01<00:01, 8.49it/s]
Build Planet: 58%|█████▊ | 15/26 [00:01<00:01, 8.19it/s]
Build Planet: 62%|██████▏ | 16/26 [00:01<00:01, 8.23it/s]
Build Planet: 65%|██████▌ | 17/26 [00:01<00:01, 8.31it/s]
Build Planet: 69%|██████▉ | 18/26 [00:02<00:00, 8.40it/s]
Build Planet: 73%|███████▎ | 19/26 [00:02<00:00, 8.59it/s]
Build Planet: 77%|███████▋ | 20/26 [00:02<00:00, 8.75it/s]
Build Planet: 81%|████████ | 21/26 [00:02<00:00, 8.81it/s]
Build Planet: 85%|████████▍ | 22/26 [00:02<00:00, 8.86it/s]
Build Planet: 88%|████████▊ | 23/26 [00:02<00:00, 8.48it/s]
Build Planet: 92%|█████████▏| 24/26 [00:02<00:00, 8.79it/s]
Build Planet: 96%|█████████▌| 25/26 [00:02<00:00, 8.70it/s]
Build Planet: 100%|██████████| 26/26 [00:02<00:00, 8.68it/s]
Build Planet: 27it [00:03, 8.70it/s]
Build Planet: 27it [00:03, 8.90it/s]
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.28it/s]
Loading Spectra: 17%|█▋ | 2/12 [00:00<00:03, 3.27it/s]
Loading Spectra: 25%|██▌ | 3/12 [00:00<00:02, 3.24it/s]
Loading Spectra: 33%|███▎ | 4/12 [00:01<00:02, 3.16it/s]
Loading Spectra: 42%|████▏ | 5/12 [00:01<00:02, 3.15it/s]
Loading Spectra: 50%|█████ | 6/12 [00:01<00:01, 3.21it/s]
Loading Spectra: 58%|█████▊ | 7/12 [00:02<00:01, 3.28it/s]
Loading Spectra: 67%|██████▋ | 8/12 [00:02<00:01, 3.33it/s]
Loading Spectra: 75%|███████▌ | 9/12 [00:02<00:00, 3.32it/s]
Loading Spectra: 83%|████████▎ | 10/12 [00:03<00:00, 3.37it/s]
Loading Spectra: 92%|█████████▏| 11/12 [00:03<00:00, 3.38it/s]
Loading Spectra: 100%|██████████| 12/12 [00:03<00:00, 3.41it/s]
Loading Spectra: 100%|██████████| 12/12 [00:03<00:00, 3.31it/s]
Build Spectra: 4%|▍ | 1/26 [00:11<04:39, 11.16s/it]
Build Spectra: 8%|▊ | 2/26 [00:11<01:52, 4.69s/it]
Build Spectra: 12%|█▏ | 3/26 [00:11<01:00, 2.62s/it]
Build Spectra: 15%|█▌ | 4/26 [00:11<00:36, 1.65s/it]
Build Spectra: 19%|█▉ | 5/26 [00:11<00:23, 1.11s/it]
Build Spectra: 23%|██▎ | 6/26 [00:11<00:15, 1.28it/s]
Build Spectra: 27%|██▋ | 7/26 [00:12<00:10, 1.73it/s]
Build Spectra: 31%|███ | 8/26 [00:12<00:07, 2.26it/s]
Build Spectra: 35%|███▍ | 9/26 [00:12<00:06, 2.83it/s]
Build Spectra: 38%|███▊ | 10/26 [00:12<00:04, 3.41it/s]
Build Spectra: 42%|████▏ | 11/26 [00:12<00:04, 3.65it/s]
Build Spectra: 46%|████▌ | 12/26 [00:13<00:03, 3.84it/s]
Build Spectra: 50%|█████ | 13/26 [00:13<00:03, 3.99it/s]
Build Spectra: 54%|█████▍ | 14/26 [00:13<00:02, 4.11it/s]
Build Spectra: 58%|█████▊ | 15/26 [00:13<00:02, 4.20it/s]
Build Spectra: 62%|██████▏ | 16/26 [00:13<00:02, 4.24it/s]
Build Spectra: 65%|██████▌ | 17/26 [00:14<00:02, 4.30it/s]
Build Spectra: 69%|██████▉ | 18/26 [00:14<00:01, 4.43it/s]
Build Spectra: 73%|███████▎ | 19/26 [00:14<00:01, 4.89it/s]
Build Spectra: 77%|███████▋ | 20/26 [00:14<00:01, 5.27it/s]
Build Spectra: 81%|████████ | 21/26 [00:14<00:00, 5.51it/s]
Build Spectra: 85%|████████▍ | 22/26 [00:15<00:00, 5.73it/s]
Build Spectra: 88%|████████▊ | 23/26 [00:15<00:00, 5.91it/s]
Build Spectra: 92%|█████████▏| 24/26 [00:15<00:00, 6.06it/s]
Build Spectra: 96%|█████████▌| 25/26 [00:15<00:00, 6.16it/s]
Build Spectra: 100%|██████████| 26/26 [00:15<00:00, 6.25it/s]
Build Spectra: 100%|██████████| 26/26 [00:15<00:00, 1.66it/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'])
/home/ted/github/VSPEC/VSPEC/analysis.py:131: 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/26 [00:00<?, ?it/s]
Build Planet: 4%|▍ | 1/26 [00:00<00:04, 5.04it/s]
Build Planet: 8%|▊ | 2/26 [00:00<00:04, 5.92it/s]
Build Planet: 12%|█▏ | 3/26 [00:00<00:03, 6.07it/s]
Build Planet: 15%|█▌ | 4/26 [00:00<00:03, 6.31it/s]
Build Planet: 19%|█▉ | 5/26 [00:00<00:03, 6.45it/s]
Build Planet: 23%|██▎ | 6/26 [00:00<00:03, 6.51it/s]
Build Planet: 27%|██▋ | 7/26 [00:01<00:02, 6.62it/s]
Build Planet: 31%|███ | 8/26 [00:01<00:02, 6.68it/s]
Build Planet: 35%|███▍ | 9/26 [00:01<00:02, 6.67it/s]
Build Planet: 38%|███▊ | 10/26 [00:01<00:02, 6.54it/s]
Build Planet: 42%|████▏ | 11/26 [00:01<00:02, 5.97it/s]
Build Planet: 46%|████▌ | 12/26 [00:01<00:02, 5.58it/s]
Build Planet: 50%|█████ | 13/26 [00:02<00:02, 5.37it/s]
Build Planet: 54%|█████▍ | 14/26 [00:02<00:02, 5.26it/s]
Build Planet: 58%|█████▊ | 15/26 [00:02<00:02, 5.18it/s]
Build Planet: 62%|██████▏ | 16/26 [00:02<00:01, 5.03it/s]
Build Planet: 65%|██████▌ | 17/26 [00:02<00:01, 5.00it/s]
Build Planet: 69%|██████▉ | 18/26 [00:03<00:01, 5.26it/s]
Build Planet: 73%|███████▎ | 19/26 [00:03<00:01, 5.50it/s]
Build Planet: 77%|███████▋ | 20/26 [00:03<00:01, 5.69it/s]
Build Planet: 81%|████████ | 21/26 [00:03<00:00, 5.80it/s]
Build Planet: 85%|████████▍ | 22/26 [00:03<00:00, 5.84it/s]
Build Planet: 88%|████████▊ | 23/26 [00:03<00:00, 5.89it/s]
Build Planet: 92%|█████████▏| 24/26 [00:04<00:00, 5.87it/s]
Build Planet: 96%|█████████▌| 25/26 [00:04<00:00, 6.03it/s]
Build Planet: 100%|██████████| 26/26 [00:04<00:00, 6.00it/s]
Build Planet: 27it [00:04, 6.13it/s]
Build Planet: 27it [00:04, 5.84it/s]
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.43it/s]
Loading Spectra: 17%|█▋ | 2/12 [00:00<00:02, 3.42it/s]
Loading Spectra: 25%|██▌ | 3/12 [00:00<00:02, 3.43it/s]
Loading Spectra: 33%|███▎ | 4/12 [00:01<00:02, 3.43it/s]
Loading Spectra: 42%|████▏ | 5/12 [00:01<00:02, 3.42it/s]
Loading Spectra: 50%|█████ | 6/12 [00:01<00:01, 3.41it/s]
Loading Spectra: 58%|█████▊ | 7/12 [00:02<00:01, 3.40it/s]
Loading Spectra: 67%|██████▋ | 8/12 [00:02<00:01, 3.40it/s]
Loading Spectra: 75%|███████▌ | 9/12 [00:02<00:00, 3.40it/s]
Loading Spectra: 83%|████████▎ | 10/12 [00:02<00:00, 3.41it/s]
Loading Spectra: 92%|█████████▏| 11/12 [00:03<00:00, 3.41it/s]
Loading Spectra: 100%|██████████| 12/12 [00:03<00:00, 3.41it/s]
Loading Spectra: 100%|██████████| 12/12 [00:03<00:00, 3.41it/s]
Build Spectra: 4%|▍ | 1/26 [00:11<04:43, 11.35s/it]
Build Spectra: 8%|▊ | 2/26 [00:11<01:54, 4.78s/it]
Build Spectra: 12%|█▏ | 3/26 [00:11<01:01, 2.68s/it]
Build Spectra: 15%|█▌ | 4/26 [00:11<00:37, 1.69s/it]
Build Spectra: 19%|█▉ | 5/26 [00:12<00:24, 1.15s/it]
Build Spectra: 23%|██▎ | 6/26 [00:12<00:16, 1.22it/s]
Build Spectra: 27%|██▋ | 7/26 [00:12<00:11, 1.64it/s]
Build Spectra: 31%|███ | 8/26 [00:12<00:08, 2.11it/s]
Build Spectra: 35%|███▍ | 9/26 [00:12<00:06, 2.60it/s]
Build Spectra: 38%|███▊ | 10/26 [00:12<00:05, 3.09it/s]
Build Spectra: 42%|████▏ | 11/26 [00:13<00:04, 3.31it/s]
Build Spectra: 46%|████▌ | 12/26 [00:13<00:03, 3.52it/s]
Build Spectra: 50%|█████ | 13/26 [00:13<00:03, 3.67it/s]
Build Spectra: 54%|█████▍ | 14/26 [00:13<00:03, 3.78it/s]
Build Spectra: 58%|█████▊ | 15/26 [00:14<00:02, 3.87it/s]
Build Spectra: 62%|██████▏ | 16/26 [00:14<00:02, 3.93it/s]
Build Spectra: 65%|██████▌ | 17/26 [00:14<00:02, 3.93it/s]
Build Spectra: 69%|██████▉ | 18/26 [00:14<00:02, 3.97it/s]
Build Spectra: 73%|███████▎ | 19/26 [00:15<00:01, 4.34it/s]
Build Spectra: 77%|███████▋ | 20/26 [00:15<00:01, 4.67it/s]
Build Spectra: 81%|████████ | 21/26 [00:15<00:01, 4.92it/s]
Build Spectra: 85%|████████▍ | 22/26 [00:15<00:00, 5.13it/s]
Build Spectra: 88%|████████▊ | 23/26 [00:15<00:00, 5.28it/s]
Build Spectra: 92%|█████████▏| 24/26 [00:16<00:00, 5.40it/s]
Build Spectra: 96%|█████████▌| 25/26 [00:16<00:00, 5.48it/s]
Build Spectra: 100%|██████████| 26/26 [00:16<00:00, 5.53it/s]
Build Spectra: 100%|██████████| 26/26 [00:16<00:00, 1.59it/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/26 [00:00<?, ?it/s]
Build Planet: 4%|▍ | 1/26 [00:00<00:02, 9.42it/s]
Build Planet: 8%|▊ | 2/26 [00:00<00:02, 9.31it/s]
Build Planet: 15%|█▌ | 4/26 [00:00<00:02, 9.90it/s]
Build Planet: 19%|█▉ | 5/26 [00:00<00:02, 9.87it/s]
Build Planet: 23%|██▎ | 6/26 [00:00<00:02, 9.48it/s]
Build Planet: 27%|██▋ | 7/26 [00:00<00:01, 9.53it/s]
Build Planet: 31%|███ | 8/26 [00:00<00:01, 9.50it/s]
Build Planet: 35%|███▍ | 9/26 [00:00<00:01, 9.54it/s]
Build Planet: 42%|████▏ | 11/26 [00:01<00:01, 9.52it/s]
Build Planet: 46%|████▌ | 12/26 [00:01<00:01, 9.05it/s]
Build Planet: 50%|█████ | 13/26 [00:01<00:01, 9.11it/s]
Build Planet: 54%|█████▍ | 14/26 [00:01<00:01, 8.65it/s]
Build Planet: 58%|█████▊ | 15/26 [00:01<00:01, 8.86it/s]
Build Planet: 62%|██████▏ | 16/26 [00:01<00:01, 8.99it/s]
Build Planet: 65%|██████▌ | 17/26 [00:01<00:01, 8.99it/s]
Build Planet: 69%|██████▉ | 18/26 [00:01<00:00, 9.03it/s]
Build Planet: 73%|███████▎ | 19/26 [00:02<00:00, 9.07it/s]
Build Planet: 77%|███████▋ | 20/26 [00:02<00:00, 9.12it/s]
Build Planet: 81%|████████ | 21/26 [00:02<00:00, 9.14it/s]
Build Planet: 85%|████████▍ | 22/26 [00:02<00:00, 9.09it/s]
Build Planet: 88%|████████▊ | 23/26 [00:02<00:00, 9.06it/s]
Build Planet: 92%|█████████▏| 24/26 [00:02<00:00, 9.06it/s]
Build Planet: 96%|█████████▌| 25/26 [00:02<00:00, 9.11it/s]
Build Planet: 100%|██████████| 26/26 [00:02<00:00, 8.80it/s]
Build Planet: 27it [00:02, 8.93it/s]
Build Planet: 27it [00:02, 9.16it/s]
Generated 94 mature spots
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.31it/s]
Loading Spectra: 17%|█▋ | 2/12 [00:00<00:02, 3.34it/s]
Loading Spectra: 25%|██▌ | 3/12 [00:00<00:02, 3.35it/s]
Loading Spectra: 33%|███▎ | 4/12 [00:01<00:02, 3.36it/s]
Loading Spectra: 42%|████▏ | 5/12 [00:01<00:02, 3.38it/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.40it/s]
Loading Spectra: 75%|███████▌ | 9/12 [00:02<00:00, 3.41it/s]
Loading Spectra: 83%|████████▎ | 10/12 [00:02<00:00, 3.41it/s]
Loading Spectra: 92%|█████████▏| 11/12 [00:03<00:00, 3.40it/s]
Loading Spectra: 100%|██████████| 12/12 [00:03<00:00, 3.40it/s]
Loading Spectra: 100%|██████████| 12/12 [00:03<00:00, 3.39it/s]
Build Spectra: 4%|▍ | 1/26 [00:11<04:38, 11.12s/it]
Build Spectra: 8%|▊ | 2/26 [00:11<01:53, 4.73s/it]
Build Spectra: 12%|█▏ | 3/26 [00:11<01:01, 2.69s/it]
Build Spectra: 15%|█▌ | 4/26 [00:11<00:38, 1.73s/it]
Build Spectra: 19%|█▉ | 5/26 [00:12<00:25, 1.20s/it]
Build Spectra: 23%|██▎ | 6/26 [00:12<00:17, 1.14it/s]
Build Spectra: 27%|██▋ | 7/26 [00:12<00:12, 1.48it/s]
Build Spectra: 31%|███ | 8/26 [00:12<00:09, 1.84it/s]
Build Spectra: 35%|███▍ | 9/26 [00:13<00:07, 2.20it/s]
Build Spectra: 38%|███▊ | 10/26 [00:13<00:06, 2.54it/s]
Build Spectra: 42%|████▏ | 11/26 [00:13<00:05, 2.68it/s]
Build Spectra: 46%|████▌ | 12/26 [00:14<00:05, 2.78it/s]
Build Spectra: 50%|█████ | 13/26 [00:14<00:04, 2.87it/s]
Build Spectra: 54%|█████▍ | 14/26 [00:14<00:04, 2.93it/s]
Build Spectra: 58%|█████▊ | 15/26 [00:15<00:03, 2.98it/s]
Build Spectra: 62%|██████▏ | 16/26 [00:15<00:03, 3.01it/s]
Build Spectra: 65%|██████▌ | 17/26 [00:15<00:02, 3.03it/s]
Build Spectra: 69%|██████▉ | 18/26 [00:16<00:02, 3.09it/s]
Build Spectra: 73%|███████▎ | 19/26 [00:16<00:02, 3.29it/s]
Build Spectra: 77%|███████▋ | 20/26 [00:16<00:01, 3.45it/s]
Build Spectra: 81%|████████ | 21/26 [00:16<00:01, 3.56it/s]
Build Spectra: 85%|████████▍ | 22/26 [00:17<00:01, 3.65it/s]
Build Spectra: 88%|████████▊ | 23/26 [00:17<00:00, 3.72it/s]
Build Spectra: 92%|█████████▏| 24/26 [00:17<00:00, 3.77it/s]
Build Spectra: 96%|█████████▌| 25/26 [00:17<00:00, 3.80it/s]
Build Spectra: 100%|██████████| 26/26 [00:18<00:00, 3.82it/s]
Build Spectra: 100%|██████████| 26/26 [00:18<00:00, 1.44it/s]
/home/ted/github/VSPEC/VSPEC/analysis.py:131: 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)
/home/ted/github/VSPEC/examples/end_to_end/plot_spotted_transit.py:352: DeprecationWarning:
Pyarrow will become a required dependency of pandas in the next major release of pandas (pandas 3.0),
(to allow more performant data types, such as the Arrow string type, and better interoperability with other libraries)
but was not found to be installed on your system.
If this would cause problems for you,
please provide us feedback at https://github.com/pandas-dev/pandas/issues/54466
import pandas as pd
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 0x7f0fc01c60d0>
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 0x7f0f9817f310>
Total running time of the script: (1 minutes 10.749 seconds)