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/ted/.pypsg/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,
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_periasteron=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),
spectral_grid='default'
)
# Set parameters for simulation
internal_params = params.InternalParameters(
header=params.Header(
data_path=Path('.vspec/proxcenb'),
teff_min=2300*u.K,teff_max=3400*u.K,
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()
/home/ted/github/VSPEC/VSPEC/gcm/heat_transfer.py:538: RuntimeWarning: Energy balance off by -6.7 %. eps = 1.5, T0 = 0.5
warnings.warn(msg, RuntimeWarning)
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/36 [00:00<?, ?it/s]
Build Planet: 3%|▎ | 1/36 [00:00<00:08, 4.24it/s]
Build Planet: 6%|▌ | 2/36 [00:00<00:07, 4.68it/s]
Build Planet: 8%|▊ | 3/36 [00:00<00:06, 4.87it/s]
Build Planet: 11%|█ | 4/36 [00:00<00:06, 4.93it/s]
Build Planet: 14%|█▍ | 5/36 [00:01<00:06, 4.87it/s]
Build Planet: 17%|█▋ | 6/36 [00:01<00:06, 4.87it/s]
Build Planet: 19%|█▉ | 7/36 [00:01<00:06, 4.78it/s]
Build Planet: 22%|██▏ | 8/36 [00:01<00:05, 4.85it/s]
Build Planet: 25%|██▌ | 9/36 [00:01<00:05, 4.72it/s]
Build Planet: 28%|██▊ | 10/36 [00:02<00:05, 4.76it/s]
Build Planet: 31%|███ | 11/36 [00:02<00:05, 4.68it/s]
Build Planet: 33%|███▎ | 12/36 [00:02<00:05, 4.71it/s]
Build Planet: 36%|███▌ | 13/36 [00:02<00:04, 4.83it/s]
Build Planet: 39%|███▉ | 14/36 [00:02<00:04, 4.79it/s]
Build Planet: 42%|████▏ | 15/36 [00:03<00:04, 4.79it/s]
Build Planet: 44%|████▍ | 16/36 [00:03<00:04, 4.77it/s]
Build Planet: 47%|████▋ | 17/36 [00:03<00:03, 4.82it/s]
Build Planet: 50%|█████ | 18/36 [00:03<00:03, 4.84it/s]
Build Planet: 53%|█████▎ | 19/36 [00:03<00:03, 4.78it/s]
Build Planet: 56%|█████▌ | 20/36 [00:04<00:03, 4.77it/s]
Build Planet: 58%|█████▊ | 21/36 [00:04<00:03, 4.80it/s]
Build Planet: 61%|██████ | 22/36 [00:04<00:02, 4.79it/s]
Build Planet: 64%|██████▍ | 23/36 [00:04<00:02, 4.74it/s]
Build Planet: 67%|██████▋ | 24/36 [00:05<00:02, 4.79it/s]
Build Planet: 69%|██████▉ | 25/36 [00:05<00:02, 4.78it/s]
Build Planet: 72%|███████▏ | 26/36 [00:05<00:02, 4.74it/s]
Build Planet: 75%|███████▌ | 27/36 [00:05<00:01, 4.73it/s]
Build Planet: 78%|███████▊ | 28/36 [00:05<00:01, 4.59it/s]
Build Planet: 81%|████████ | 29/36 [00:06<00:01, 4.66it/s]
Build Planet: 83%|████████▎ | 30/36 [00:06<00:01, 4.68it/s]
Build Planet: 86%|████████▌ | 31/36 [00:06<00:01, 4.71it/s]
Build Planet: 89%|████████▉ | 32/36 [00:06<00:00, 4.73it/s]
Build Planet: 92%|█████████▏| 33/36 [00:06<00:00, 4.69it/s]
Build Planet: 94%|█████████▍| 34/36 [00:07<00:00, 4.64it/s]
Build Planet: 97%|█████████▋| 35/36 [00:07<00:00, 4.66it/s]
Build Planet: 100%|██████████| 36/36 [00:07<00:00, 4.60it/s]
Build Planet: 37it [00:07, 4.55it/s]
Build Planet: 37it [00:07, 4.73it/s]
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, 2.98it/s]
Loading Spectra: 17%|█▋ | 2/12 [00:00<00:03, 2.91it/s]
Loading Spectra: 25%|██▌ | 3/12 [00:01<00:03, 2.89it/s]
Loading Spectra: 33%|███▎ | 4/12 [00:01<00:02, 2.94it/s]
Loading Spectra: 42%|████▏ | 5/12 [00:01<00:02, 2.94it/s]
Loading Spectra: 50%|█████ | 6/12 [00:02<00:02, 2.95it/s]
Loading Spectra: 58%|█████▊ | 7/12 [00:02<00:01, 2.90it/s]
Loading Spectra: 67%|██████▋ | 8/12 [00:02<00:01, 2.90it/s]
Loading Spectra: 75%|███████▌ | 9/12 [00:03<00:01, 2.92it/s]
Loading Spectra: 83%|████████▎ | 10/12 [00:03<00:00, 3.02it/s]
Loading Spectra: 92%|█████████▏| 11/12 [00:03<00:00, 3.08it/s]
Loading Spectra: 100%|██████████| 12/12 [00:04<00:00, 3.13it/s]
Loading Spectra: 100%|██████████| 12/12 [00:04<00:00, 2.99it/s]
Build Spectra: 3%|▎ | 1/36 [00:14<08:33, 14.68s/it]
Build Spectra: 6%|▌ | 2/36 [00:14<03:30, 6.18s/it]
Build Spectra: 8%|▊ | 3/36 [00:15<01:54, 3.47s/it]
Build Spectra: 11%|█ | 4/36 [00:15<01:10, 2.19s/it]
Build Spectra: 14%|█▍ | 5/36 [00:15<00:46, 1.49s/it]
Build Spectra: 17%|█▋ | 6/36 [00:15<00:31, 1.07s/it]
Build Spectra: 19%|█▉ | 7/36 [00:16<00:23, 1.25it/s]
Build Spectra: 22%|██▏ | 8/36 [00:16<00:17, 1.61it/s]
Build Spectra: 25%|██▌ | 9/36 [00:16<00:13, 1.99it/s]
Build Spectra: 28%|██▊ | 10/36 [00:16<00:10, 2.37it/s]
Build Spectra: 31%|███ | 11/36 [00:17<00:09, 2.71it/s]
Build Spectra: 33%|███▎ | 12/36 [00:17<00:07, 3.03it/s]
Build Spectra: 36%|███▌ | 13/36 [00:17<00:06, 3.31it/s]
Build Spectra: 39%|███▉ | 14/36 [00:17<00:06, 3.53it/s]
Build Spectra: 42%|████▏ | 15/36 [00:18<00:05, 3.69it/s]
Build Spectra: 44%|████▍ | 16/36 [00:18<00:05, 3.78it/s]
Build Spectra: 47%|████▋ | 17/36 [00:18<00:04, 3.83it/s]
Build Spectra: 50%|█████ | 18/36 [00:18<00:04, 3.78it/s]
Build Spectra: 53%|█████▎ | 19/36 [00:19<00:04, 3.43it/s]
Build Spectra: 56%|█████▌ | 20/36 [00:19<00:04, 3.39it/s]
Build Spectra: 58%|█████▊ | 21/36 [00:19<00:04, 3.57it/s]
Build Spectra: 61%|██████ | 22/36 [00:19<00:03, 3.73it/s]
Build Spectra: 64%|██████▍ | 23/36 [00:20<00:03, 3.87it/s]
Build Spectra: 67%|██████▋ | 24/36 [00:20<00:03, 3.96it/s]
Build Spectra: 69%|██████▉ | 25/36 [00:20<00:02, 4.04it/s]
Build Spectra: 72%|███████▏ | 26/36 [00:20<00:02, 4.04it/s]
Build Spectra: 75%|███████▌ | 27/36 [00:21<00:02, 3.92it/s]
Build Spectra: 78%|███████▊ | 28/36 [00:21<00:02, 3.99it/s]
Build Spectra: 81%|████████ | 29/36 [00:21<00:01, 4.05it/s]
Build Spectra: 83%|████████▎ | 30/36 [00:21<00:01, 3.97it/s]
Build Spectra: 86%|████████▌ | 31/36 [00:22<00:01, 3.96it/s]
Build Spectra: 89%|████████▉ | 32/36 [00:22<00:01, 3.96it/s]
Build Spectra: 92%|█████████▏| 33/36 [00:22<00:00, 4.02it/s]
Build Spectra: 94%|█████████▍| 34/36 [00:22<00:00, 4.04it/s]
Build Spectra: 97%|█████████▋| 35/36 [00:23<00:00, 3.48it/s]
Build Spectra: 100%|██████████| 36/36 [00:23<00:00, 3.62it/s]
Build Spectra: 100%|██████████| 36/36 [00:23<00:00, 1.53it/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 0x7f1a4ec65550>
Total running time of the script: (0 minutes 34.251 seconds)