Note
Go to the end to download the full example code.
Observe the MIRECLE Phase Curve of Proxima Centauri b#
This example observes the closest exoplanet with the Mid-Infrared Exoplanet CLimate Explorer.
Saved settings to /home/runner/.libpypsg/settings.json
Reloading settings...
Create the needed configurations#
MIRECLE is described in Mandell et al. [2022]
# Instrument
inst = params.InstrumentParameters.mirecle()
# Observation
observation = params.ObservationParameters(
observation_time=12*u.day,
integration_time=8*u.hr
)
# PSG
psg_params = params.psgParameters(
use_molecular_signatures=True,
gcm_binning=200,
phase_binning=1,
use_continuum_stellar=True,
nmax=0,
lmax=0,
continuum=['Rayleigh', 'Refraction','CIA_all'],
)
# Star and Planet
star_teff = 2900*u.K
star_rad = 0.141*u.R_sun
inclination = 85*u.deg
planet_mass = 1.07*u.M_earth/np.sin(inclination)
planet_rad = 1*u.R_earth * planet_mass.to_value(u.M_earth)**3
orbit_rad = 0.04856*u.AU
orbit_period = 11.18*u.day
planet_rot_period = orbit_period
star_rot_period = 90*u.day
planet_mass = 1.07*u.M_earth
star_mass = 0.122*u.M_sun
initial_phase = 180*u.deg
planet_params = params.PlanetParameters(
name='proxcenb',
radius=planet_rad,
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
)
system_params = params.SystemParameters(
distance=1.3*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': 30,
'nlon': 30,
'nlat': 15,
'epsilon': 1.5,
'albedo': 0.3,
'emissivity': 1.0,
'lat_redistribution': 0.5,
'gamma': 1.4,
'psurf': 1*u.bar,
'ptop': 1e-8*u.bar,
'wind': {'U': '0 m/s','V':'0 m/s'},
'molecules':{'CO2':1e-4}
}
gcm_params = params.gcmParameters.from_dict({
'star':star_dict,
'planet':planet_dict,
'gcm':{'vspec':gcm_dict,'mean_molec_weight':28}
})
quiet_star = params.StarParameters(
spots=params.SpotParameters.none(),
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.proxima(),
faculae=params.FaculaParameters.none(),
flares=params.FlareParameters.none(),
granulation=params.GranulationParameters.none(),
grid_params=(500, 1000),
)
# Set parameters for simulation
internal_params = params.InternalParameters(
header=params.Header(
data_path=Path('.vspec/proxcenb'),
spec_grid = params.VSPECGridParameters(
max_teff=3400*u.K,min_teff=2300*u.K,
impl_bin='rust',impl_interp='scipy',fail_on_missing=False
),
seed = SEED),
star = quiet_star,
psg=psg_params,
planet=planet_params,
system=system_params,
obs=observation,
gcm=gcm_params,
inst=inst
)
Run VSPEC#
We read in the config file and run the model.
model = ObservationModel(internal_params)
model.build_planet()
model.build_spectra()
Starting at phase 180.0 deg, observe for 12.0 d in 36 steps
Phases = [180. 190.73 201.47 212.2 222.93 233.67 244.4 255.13 265.87 276.6
287.33 298.07 308.8 319.53 330.27 341. 351.74 2.47 13.2 23.94
34.67 45.4 56.14 66.87 77.6 88.34 99.07 109.8 120.54 131.27
142. 152.74 163.47 174.2 184.94 195.67 206.4 ] deg
Build Planet: 0%| | 0/37 [00:00<?, ?it/s]
Build Planet: 3%|▎ | 1/37 [00:01<00:41, 1.15s/it]
Build Planet: 5%|▌ | 2/37 [00:02<00:40, 1.15s/it]
Build Planet: 8%|▊ | 3/37 [00:03<00:38, 1.14s/it]
Build Planet: 11%|█ | 4/37 [00:04<00:38, 1.18s/it]
Build Planet: 14%|█▎ | 5/37 [00:05<00:37, 1.16s/it]
Build Planet: 16%|█▌ | 6/37 [00:06<00:35, 1.16s/it]
Build Planet: 19%|█▉ | 7/37 [00:08<00:34, 1.15s/it]
Build Planet: 22%|██▏ | 8/37 [00:09<00:33, 1.14s/it]
Build Planet: 24%|██▍ | 9/37 [00:10<00:32, 1.15s/it]
Build Planet: 27%|██▋ | 10/37 [00:11<00:30, 1.15s/it]
Build Planet: 30%|██▉ | 11/37 [00:12<00:29, 1.14s/it]
Build Planet: 32%|███▏ | 12/37 [00:13<00:28, 1.15s/it]
Build Planet: 35%|███▌ | 13/37 [00:14<00:27, 1.15s/it]
Build Planet: 38%|███▊ | 14/37 [00:16<00:26, 1.14s/it]
Build Planet: 41%|████ | 15/37 [00:17<00:25, 1.15s/it]
Build Planet: 43%|████▎ | 16/37 [00:18<00:23, 1.14s/it]
Build Planet: 46%|████▌ | 17/37 [00:19<00:22, 1.15s/it]
Build Planet: 49%|████▊ | 18/37 [00:20<00:21, 1.15s/it]
Build Planet: 51%|█████▏ | 19/37 [00:21<00:20, 1.14s/it]
Build Planet: 54%|█████▍ | 20/37 [00:22<00:19, 1.15s/it]
Build Planet: 57%|█████▋ | 21/37 [00:24<00:18, 1.15s/it]
Build Planet: 59%|█████▉ | 22/37 [00:25<00:17, 1.14s/it]
Build Planet: 62%|██████▏ | 23/37 [00:26<00:16, 1.15s/it]
Build Planet: 65%|██████▍ | 24/37 [00:27<00:14, 1.14s/it]
Build Planet: 68%|██████▊ | 25/37 [00:28<00:13, 1.14s/it]
Build Planet: 70%|███████ | 26/37 [00:29<00:12, 1.14s/it]
Build Planet: 73%|███████▎ | 27/37 [00:31<00:11, 1.16s/it]
Build Planet: 76%|███████▌ | 28/37 [00:32<00:10, 1.15s/it]
Build Planet: 78%|███████▊ | 29/37 [00:33<00:09, 1.15s/it]
Build Planet: 81%|████████ | 30/37 [00:34<00:07, 1.14s/it]
Build Planet: 84%|████████▍ | 31/37 [00:35<00:06, 1.14s/it]
Build Planet: 86%|████████▋ | 32/37 [00:36<00:05, 1.14s/it]
Build Planet: 89%|████████▉ | 33/37 [00:37<00:04, 1.13s/it]
Build Planet: 92%|█████████▏| 34/37 [00:38<00:03, 1.13s/it]
Build Planet: 95%|█████████▍| 35/37 [00:40<00:02, 1.13s/it]
Build Planet: 97%|█████████▋| 36/37 [00:41<00:01, 1.13s/it]
Build Planet: 100%|██████████| 37/37 [00:42<00:00, 1.13s/it]
Build Planet: 100%|██████████| 37/37 [00:42<00:00, 1.14s/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/36 [00:00<?, ?it/s]
Loading Spectra: 0%| | 0/12 [00:00<?, ?it/s]
Loading Spectra: 8%|▊ | 1/12 [00:00<00:03, 3.23it/s]
Loading Spectra: 17%|█▋ | 2/12 [00:00<00:03, 3.25it/s]
Loading Spectra: 25%|██▌ | 3/12 [00:00<00:02, 3.26it/s]
Loading Spectra: 33%|███▎ | 4/12 [00:01<00:02, 3.27it/s]
Loading Spectra: 42%|████▏ | 5/12 [00:01<00:02, 3.27it/s]
Loading Spectra: 50%|█████ | 6/12 [00:01<00:01, 3.28it/s]
Loading Spectra: 58%|█████▊ | 7/12 [00:02<00:01, 3.27it/s]
Loading Spectra: 67%|██████▋ | 8/12 [00:02<00:01, 3.27it/s]
Loading Spectra: 75%|███████▌ | 9/12 [00:02<00:00, 3.28it/s]
Loading Spectra: 83%|████████▎ | 10/12 [00:03<00:00, 3.28it/s]
Loading Spectra: 92%|█████████▏| 11/12 [00:03<00:00, 3.28it/s]
Loading Spectra: 100%|██████████| 12/12 [00:03<00:00, 3.28it/s]
Loading Spectra: 100%|██████████| 12/12 [00:03<00:00, 3.27it/s]
Build Spectra: 3%|▎ | 1/36 [00:04<02:21, 4.04s/it]
Build Spectra: 6%|▌ | 2/36 [00:04<01:03, 1.88s/it]
Build Spectra: 8%|▊ | 3/36 [00:04<00:39, 1.19s/it]
Build Spectra: 11%|█ | 4/36 [00:05<00:27, 1.15it/s]
Build Spectra: 14%|█▍ | 5/36 [00:05<00:21, 1.45it/s]
Build Spectra: 17%|█▋ | 6/36 [00:05<00:17, 1.72it/s]
Build Spectra: 19%|█▉ | 7/36 [00:06<00:14, 1.96it/s]
Build Spectra: 22%|██▏ | 8/36 [00:06<00:13, 2.15it/s]
Build Spectra: 25%|██▌ | 9/36 [00:06<00:11, 2.30it/s]
Build Spectra: 28%|██▊ | 10/36 [00:07<00:10, 2.39it/s]
Build Spectra: 31%|███ | 11/36 [00:07<00:10, 2.48it/s]
Build Spectra: 33%|███▎ | 12/36 [00:08<00:09, 2.55it/s]
Build Spectra: 36%|███▌ | 13/36 [00:08<00:08, 2.60it/s]
Build Spectra: 39%|███▉ | 14/36 [00:08<00:08, 2.65it/s]
Build Spectra: 42%|████▏ | 15/36 [00:09<00:07, 2.68it/s]
Build Spectra: 44%|████▍ | 16/36 [00:09<00:07, 2.70it/s]
Build Spectra: 47%|████▋ | 17/36 [00:09<00:06, 2.72it/s]
Build Spectra: 50%|█████ | 18/36 [00:10<00:06, 2.73it/s]
Build Spectra: 53%|█████▎ | 19/36 [00:10<00:06, 2.73it/s]
Build Spectra: 56%|█████▌ | 20/36 [00:11<00:05, 2.73it/s]
Build Spectra: 58%|█████▊ | 21/36 [00:11<00:05, 2.74it/s]
Build Spectra: 61%|██████ | 22/36 [00:11<00:05, 2.73it/s]
Build Spectra: 64%|██████▍ | 23/36 [00:12<00:04, 2.71it/s]
Build Spectra: 67%|██████▋ | 24/36 [00:12<00:04, 2.70it/s]
Build Spectra: 69%|██████▉ | 25/36 [00:12<00:04, 2.71it/s]
Build Spectra: 72%|███████▏ | 26/36 [00:13<00:03, 2.71it/s]
Build Spectra: 75%|███████▌ | 27/36 [00:13<00:03, 2.71it/s]
Build Spectra: 78%|███████▊ | 28/36 [00:13<00:02, 2.71it/s]
Build Spectra: 81%|████████ | 29/36 [00:14<00:02, 2.73it/s]
Build Spectra: 83%|████████▎ | 30/36 [00:14<00:02, 2.73it/s]
Build Spectra: 86%|████████▌ | 31/36 [00:15<00:01, 2.74it/s]
Build Spectra: 89%|████████▉ | 32/36 [00:15<00:01, 2.74it/s]
Build Spectra: 92%|█████████▏| 33/36 [00:15<00:01, 2.73it/s]
Build Spectra: 94%|█████████▍| 34/36 [00:16<00:00, 2.74it/s]
Build Spectra: 97%|█████████▋| 35/36 [00:16<00:00, 2.74it/s]
Build Spectra: 100%|██████████| 36/36 [00:16<00:00, 2.74it/s]
Build Spectra: 100%|██████████| 36/36 [00:16<00:00, 2.13it/s]
Load in the data#
We can use VSPEC to read in the synthetic data we just created.
data = PhaseAnalyzer(model.directories['all_model'])
Plot the phase curve#
fig,ax = plt.subplots(1,1,figsize=(4,4),tight_layout=True)
emission = (data.thermal/data.total).to_value(u.dimensionless_unscaled)*1e6
noise = (data.noise/data.total).to_value(u.dimensionless_unscaled)*1e6
sim_noise = model.rng.normal(loc=0,scale=noise)
sim_data = emission + sim_noise
time = (data.time - data.time[0]).to_value(u.day)
wl = data.wavelength.to_value(u.um)
im = ax.pcolormesh(time,wl,sim_data,cmap='viridis')
fig.colorbar(im,ax=ax,label='Emission (ppm)')
ax.set_xlabel('Time (days)')
ax.set_ylabel('Wavelength ($\\mu m$)')
Text(33.00000000000001, 0.5, 'Wavelength ($\\mu m$)')
Plot the integrated spectrum#
true = np.mean(emission,axis=1)
observed = np.mean(sim_data,axis=1)
err = np.sqrt(np.sum(noise**2,axis=1))/noise.shape[1]
fig,ax = plt.subplots(1,1,figsize=(4,3),tight_layout=True)
ax.plot(wl,true,c='xkcd:azure',label='True')
ax.errorbar(wl,observed,yerr=err,fmt='o',color='xkcd:rose pink',label='Observed',markersize=2)
ax.set_xlabel('Wavelength ($\\mu m$)')
ax.set_ylabel('Planetary Emission (ppm)')
ax.legend()
<matplotlib.legend.Legend object at 0x7fd4185d9c50>
Total running time of the script: (1 minutes 11.684 seconds)