In [ ]:
import os
from os import path
# Third-party
from astropy.io import fits
from astropy.stats import median_absolute_deviation
from astropy.table import Table, QTable, join
from astropy.time import Time
import astropy.units as u
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import numpy as np
%matplotlib inline
import h5py
import pandas as pd
from sqlalchemy import func
import tqdm
from thejoker import JokerSamples
from twoface.config import TWOFACE_CACHE_PATH
from twoface.samples_analysis import unimodal_P, MAP_sample
from twoface.db import (db_connect, AllStar, AllVisit, AllVisitToAllStar, NessRG,
StarResult, Status, JokerRun)
from twoface.plot import plot_two_panel, plot_phase_fold, plot_data_orbits, _RV_LBL
from twoface.mass import get_m2_min, mf, asini, a2sini, stellar_radius
In [ ]:
plot_path = '../../paper/1-catalog/figures/'
table_path = '../../paper/1-catalog/tables/'
In [ ]:
Session, _ = db_connect(path.join(TWOFACE_CACHE_PATH, 'apogee.sqlite'))
session = Session()
samples_file = path.join(TWOFACE_CACHE_PATH, 'apogee-jitter.hdf5')
mcmc_samples_file = path.join(TWOFACE_CACHE_PATH, 'apogee-jitter-mcmc.hdf5')
In [ ]:
run = session.query(JokerRun).limit(1).one()
joker_pars = run.get_joker_params()
In [ ]:
high_K_stars = session.query(AllStar).join(StarResult).filter(StarResult.status_id>0).filter(StarResult.high_K).all()
len(high_K_stars)
In [ ]:
mean_R = []
max_R = []
with h5py.File(mcmc_samples_file, 'r') as f:
for k in f.keys():
R = f[k]['chain-stats/gelman_rubin'][:]
mean_R.append(np.mean(R))
max_R.append(np.max(R))
In [ ]:
# mean_R_old = []
# max_R_old = []
# with h5py.File('../../cache/apogee-jitter-mcmc-old.hdf5', 'r') as f:
# for k in f.keys():
# R = f[k]['chain-stats/gelman_rubin'][:]
# mean_R_old.append(np.mean(R))
# max_R_old.append(np.max(R))
In [ ]:
# bins = np.linspace(0.9, 2.2, 32)
# plt.hist(mean_R, bins=bins, alpha=0.3)
# plt.hist(mean_R_old, bins=bins, alpha=0.3);
# plt.axvline(1.1)
# print((np.array(mean_R) < 1.1).sum(),
# (np.array(mean_R_old) < 1.1).sum())
# plt.figure()
# plt.hist(max_R, bins=bins, alpha=0.3)
# plt.hist(max_R_old, bins=bins, alpha=0.3);
# plt.axvline(1.1)
# print((np.array(max_R) < 1.1).sum(),
# (np.array(max_R_old) < 1.1).sum())
For all high-K stars, classify as unimodal or not based on TheJoker samples. Then do same for MCMC samples, AND the selections:
In [ ]:
unimodal_thejoker = []
with h5py.File(samples_file, 'r') as f:
for star in tqdm.tqdm(high_K_stars):
samples = JokerSamples.from_hdf5(f[star.apogee_id])
data = star.apogeervdata()
unimodal_thejoker.append(unimodal_P(samples, data))
unimodal_thejoker = np.array(unimodal_thejoker)
unimodal_thejoker.sum()
In [ ]:
unimodal_mcmc = []
converged_mcmc = []
with h5py.File(mcmc_samples_file, 'r') as f:
for star in tqdm.tqdm(high_K_stars):
if star.apogee_id not in f:
unimodal_mcmc.append(False)
converged_mcmc.append(True)
continue
R = f[star.apogee_id]['chain-stats/gelman_rubin'][:]
converged_mcmc.append(np.mean(R) <= 1.1)
samples = JokerSamples.from_hdf5(f[star.apogee_id])
data = star.apogeervdata()
unimodal_mcmc.append(unimodal_P(samples, data))
unimodal_mcmc = np.array(unimodal_mcmc)
converged_mcmc = np.array(converged_mcmc)
unimodal_mcmc.sum(), converged_mcmc.sum()
In [ ]:
unimodal_mask = unimodal_thejoker | unimodal_mcmc
unimodal_converged_mask = unimodal_thejoker & (unimodal_mcmc & converged_mcmc)
unimodal_converged_idx = np.where(unimodal_converged_mask)[0]
unimodal_mask.sum(), unimodal_converged_mask.sum()
In [ ]:
unimodal_stars = np.array(high_K_stars)[unimodal_mask]
unimodal_converged = converged_mcmc[unimodal_mask]
In [ ]:
rows = dict()
rows['APOGEE_ID'] = []
for k in JokerSamples._valid_keys:
rows[k] = []
rows[k + '_err'] = []
rows['t0'] = []
rows['converged'] = []
rows['Gelman-Rubin'] = []
with h5py.File(mcmc_samples_file, 'r') as mcmc_f, h5py.File(samples_file, 'r') as joker_f:
for i, star in tqdm.tqdm(enumerate(unimodal_stars)):
data = star.apogeervdata()
if star.apogee_id in mcmc_f: # and unimodal_converged[i]:
samples = JokerSamples.from_hdf5(mcmc_f[star.apogee_id])
R = mcmc_f[star.apogee_id]['chain-stats/gelman_rubin'][:]
else:
samples = JokerSamples.from_hdf5(joker_f[star.apogee_id])
R = np.full(7, np.nan)
rows['APOGEE_ID'].append(star.apogee_id)
MAP = MAP_sample(data, samples, joker_pars)
for k in samples.keys():
rows[k].append(MAP[k])
# if unimodal_converged[i]:
# rows[k+'_err'].append(1.5 * median_absolute_deviation(samples[k]))
# else:
# rows[k+'_err'].append(np.nan * samples[k].unit)
rows[k+'_err'].append(1.5 * median_absolute_deviation(samples[k]))
rows['t0'].append(data.t0.tcb.mjd)
rows['converged'].append(unimodal_converged[i])
rows['Gelman-Rubin'].append(R)
for k in rows:
if hasattr(rows[k][0], 'unit'):
rows[k] = u.Quantity(rows[k])
rows['t0'] = Time(rows['t0'], format='mjd', scale='tcb')
In [ ]:
tbl = Table(rows, masked=True)
In [ ]:
ness_tbl = Table.read('../../data/NessRG.fits')
ness_tbl.rename_column('2MASS', 'APOGEE_ID')
ness_tbl = ness_tbl[np.isin(ness_tbl['APOGEE_ID'], tbl['APOGEE_ID'])]
# trim the duplicates...
_, unq_idx = np.unique(ness_tbl['APOGEE_ID'], return_index=True)
ness_tbl = ness_tbl[unq_idx]
In [ ]:
def stddev(vals):
return 1.5 * median_absolute_deviation(vals, ignore_nan=True)
In [ ]:
rnd = np.random.RandomState(seed=42)
N = rnd.normal
tbl['M1'] = np.full(len(tbl), np.nan) * u.Msun
tbl['M1_err'] = np.full(len(tbl), np.nan) * u.Msun
tbl['M2_min'] = np.full(len(tbl), np.nan) * u.Msun
tbl['M2_min_err'] = np.full(len(tbl), np.nan) * u.Msun
tbl['q_min'] = np.full(len(tbl), np.nan)
tbl['q_min_err'] = np.full(len(tbl), np.nan)
tbl['R1'] = np.full(len(tbl), np.nan) * u.Rsun
tbl['R1_err'] = np.full(len(tbl), np.nan) * u.Rsun
tbl['a_sini'] = np.full(len(tbl), np.nan) * u.au
tbl['a_sini_err'] = np.full(len(tbl), np.nan) * u.au
tbl['a2_sini'] = np.full(len(tbl), np.nan) * u.au
tbl['a2_sini_err'] = np.full(len(tbl), np.nan) * u.au
n_samples = 8192
for i, row in tqdm.tqdm(enumerate(tbl)):
ness_row = ness_tbl[ness_tbl['APOGEE_ID'] == row['APOGEE_ID']]
if len(ness_row) == 0:
continue
star = AllStar.get_apogee_id(session, row['APOGEE_ID'])
m1_samples = np.exp(N(ness_row['lnM'], ness_row['e_logM'], size=n_samples)) * u.Msun
loggs = N(star.logg, star.logg_err, n_samples)
Ps = N(row['P'], row['P_err'], n_samples) * tbl['P'].unit
Ks = N(row['K'], row['K_err'], n_samples) * tbl['K'].unit
es = N(row['e'], row['e_err'], n_samples)
# else:
# Ps = ([row['P']] * n_samples) * tbl['P'].unit
# Ks = ([row['K']] * n_samples) * tbl['K'].unit
# es = np.array([row['e']] * n_samples)
mass_func = mf(P=Ps, K=Ks, e=es)
m2_mins = get_m2_min(m1_samples, mass_func)
asinis = asini(Ps, es, Ks, m1_samples, m2_mins)
a2sinis = a2sini(Ps, es, Ks, m1_samples, m2_mins)
R1s = stellar_radius(loggs, m1_samples).to(u.Rsun)
tbl['M1'][i] = np.median(m1_samples).to(u.Msun).value
tbl['M2_min'][i] = np.nanmedian(m2_mins).to(u.Msun).value
tbl['a_sini'][i] = np.nanmedian(asinis).to(u.au).value
tbl['a2_sini'][i] = np.nanmedian(a2sinis).to(u.au).value
tbl['R1'][i] = np.nanmedian(R1s).to(u.Rsun).value
tbl['M1_err'][i] = stddev(m1_samples).to(u.Msun).value
tbl['M2_min_err'][i] = stddev(m2_mins).to(u.Msun).value
tbl['a_sini_err'][i] = stddev(asinis).to(u.au).value
tbl['a2_sini_err'][i] = stddev(a2sinis).to(u.au).value
tbl['R1_err'][i] = stddev(R1s).to(u.Rsun).value
tbl['q_min'] = (u.Quantity(tbl['M2_min']) / u.Quantity(tbl['M1'])).decompose()
tbl['q_min_err'] = tbl['q_min'] * \
np.sqrt((tbl['M2_min_err']/tbl['M2_min'])**2 +
(tbl['M1_err']/tbl['M1'])**2)
mask_ = np.isnan(tbl['M1']) | np.isnan(tbl['M2_min'])
tbl['M1'].mask = mask_
tbl['M1_err'].mask = mask_
tbl['M2_min'].mask = mask_
tbl['M2_min_err'].mask = mask_
Add Ness columns following our columns:
In [ ]:
tbl_with_ness = join(tbl, ness_tbl, keys='APOGEE_ID', join_type='outer')
assert len(tbl_with_ness) == len(tbl)
In [ ]:
allstar_tbl = fits.getdata('/Users/adrian/data/APOGEE_DR14/allStar-l31c.2.fits')
allstar_tbl = allstar_tbl[np.isin(allstar_tbl['APOGEE_ID'], tbl['APOGEE_ID'])]
# trim the duplicates...
_, unq_idx = np.unique(allstar_tbl['APOGEE_ID'], return_index=True)
allstar_tbl = allstar_tbl[unq_idx]
assert len(allstar_tbl) == len(tbl)
allstar_tbl = Table(allstar_tbl)
allstar_tbl.rename_column('K', 'KS')
allstar_tbl.rename_column('K_ERR', 'KS_ERR')
In [ ]:
full_catalog = join(tbl_with_ness, allstar_tbl, keys='APOGEE_ID')
full_catalog[:1]
In [ ]:
from astropy.io import ascii
In [ ]:
rcdr14 = Table.read('/Users/adrian/data/APOGEE_DR14/apogee-rc-DR14.fits')
rcting = ascii.read('../../data/ting-2018.txt')
In [ ]:
(rcting['Classification'] == 'RC_Pristine').sum()
In [ ]:
full_catalog['DR14RC'] = np.isin(full_catalog['APOGEE_ID'], rcdr14['APOGEE_ID'])
full_catalog['TINGRC'] = np.isin(full_catalog['APOGEE_ID'], rcting[rcting['Classification'] == 'RC_Pristine']['Designation'])
# full_catalog['TINGRC'] = np.isin(full_catalog['APOGEE_ID'], rcting['Designation'])
In [ ]:
len(full_catalog), full_catalog['DR14RC'].sum(), full_catalog['TINGRC'].sum()
In [ ]:
full_catalog['M1'][full_catalog['M1'].mask] = np.nan
full_catalog['M2_min'][full_catalog['M2_min'].mask] = np.nan
In [ ]:
for name in full_catalog.colnames[:30]:
c1 = '\\texttt{{{0}}}'.format(name.replace('_', '\\_'))
try:
c2 = '{0:latex_inline}'.format(full_catalog[name].unit)
except TypeError:
c2 = ''
except AttributeError:
c2 = ''
if len(c1) < 26:
c1 = c1 + ' '*(26 - len(c1))
if len(c2) < 24:
c2 = c2 + ' '*(24 - len(c2))
print('{0} & {1} & <description> \\\\'.format(c1, c2))
In [ ]:
# _path = '../../plots/unimodal/'
# os.makedirs(_path, exist_ok=True)
# units = dict()
# for c in full_catalog.colnames:
# if full_catalog[c].unit is not None:
# units[c] = full_catalog[c].unit
# else:
# units[c] = 1.
# for row in full_catalog:
# apogee_id = row['APOGEE_ID']
# star = AllStar.get_apogee_id(session, apogee_id)
# data = star.apogeervdata()
# row = row[JokerSamples._valid_keys]
# sample = JokerSamples(**{c: row[c]*units[c] for c in row.colnames})
# sample.t0 = data.t0
# fig, axes = plt.subplots(1, 2, figsize=(12, 5), sharey=True)
# plot_data_orbits(data, sample[None], highlight_P_extrema=False,
# ax=axes[0], plot_kwargs=dict(alpha=1., linewidth=1.))
# plot_phase_fold(data, sample, ax=axes[1], label=False)
# axes[1].set_xlabel('phase')
# axes[0].set_title(apogee_id)
# fig.tight_layout()
# fig.savefig(path.join(_path, '{0}.png'.format(apogee_id)), dpi=200)
# plt.close(fig)
In [ ]:
# unimodal:
check = np.array([
'2M05224382+4300425',
'2M08505498+1156503',
'2M08510723+1153019',
'2M08512530+1202563',
'2M09522871+3811487',
'2M10264342+1340172',
'2M10513288-0250550',
'2M13011859+2844170',
'2M13162279+1739074',
'2M13175687+7151180',
'2M13484871+1913474',
'2M14574438+2106271',
'2M15054553+2220325',
'2M15101168+6708289',
'2M16342938-1248117',
'2M18012240-0920302',
'2M18343302+1949166',
'2M18481414-0251133',
'2M17223366+4850318',
'2M15184139+0206004',
'2M21260907+1100178',
'2M17105698+4301117'
])
# Suspect:
# SUSPECT_BROAD_LINES, or SUSPECT_RV_COMBINATIONS
suspect = full_catalog['APOGEE_ID'][(full_catalog['STARFLAG'] & np.sum(2**np.array([16]))) != 0]
check = check[~np.isin(check, suspect)]
print(len(suspect), len(check))
clean_flag = np.zeros(len(full_catalog), dtype=int)
clean_flag[np.isin(full_catalog['APOGEE_ID'], check)] = 1
clean_flag[np.isin(full_catalog['APOGEE_ID'], suspect)] = 2
full_catalog['clean_flag'] = clean_flag
In [ ]:
(full_catalog['clean_flag'] == 0).sum()
In [ ]:
full_catalog.write(path.join(table_path, 'highK-unimodal.fits'), overwrite=True)
In [ ]:
test = QTable.read(path.join(table_path, 'highK-unimodal.fits'),
astropy_native=True, character_as_bytes=False)
In [ ]:
full_catalog = Table.read(path.join(table_path, 'highK-unimodal.fits'))
In [ ]:
arr = np.array(full_catalog[full_catalog['converged'] & np.isfinite(full_catalog['Gelman-Rubin'][:, 0])]['APOGEE_ID'],
dtype='U20')
In [ ]:
np.random.seed(42)
rc = {
'axes.labelsize': 18,
'xtick.labelsize': 14,
'ytick.labelsize': 14
}
subset = full_catalog[full_catalog['converged'] & np.isfinite(full_catalog['Gelman-Rubin'][:, 0])]
rand_subset = np.random.choice(len(subset), size=8, replace=False)
rand_subset = rand_subset[np.argsort(subset['e'][rand_subset])]
with h5py.File(samples_file, 'r') as jok_f, h5py.File(mcmc_samples_file, 'r') as mcmc_f:
with mpl.rc_context(rc):
fig, axes = plt.subplots(4, 2, figsize=(8, 10), sharex=True)
for i, idx in enumerate(rand_subset):
ax = axes.flat[i]
apogee_id = subset[idx]['APOGEE_ID']
star = AllStar.get_apogee_id(session, apogee_id)
data = star.apogeervdata()
if apogee_id in mcmc_f:
f = mcmc_f
print('mcmc')
else:
f = jok_f
print('thejoker')
samples = JokerSamples.from_hdf5(f[star.apogee_id])
samples.t0 = data.t0
if len(samples) > 1:
sample = MAP_sample(data, samples, joker_pars)
else:
sample = samples[0]
fig = plot_phase_fold(data, sample, ax=ax,
jitter_errorbar=True, label=False)
xlim = ax.get_xlim()
ylim = (data.rv.value.min(), data.rv.value.max())
yspan = ylim[1]-ylim[0]
ylim = ax.set_ylim(ylim[0]-0.35*yspan, ylim[1]+0.35*yspan)
text = ('{0}, '.format(star.apogee_id) +
'$P = {0.value:.2f}$ {0.unit:latex}, '.format(sample['P']) +
'$e = {0:.2f}$'.format(sample['e']))
ax.text(xlim[0] + (xlim[1]-xlim[0])/15,
ylim[1] - (ylim[1]-ylim[0])/20,
text, fontsize=10, va='top', ha='left')
# _ = plot_two_panel(data, samples)
ax.set_xlim(-0.02, 1.02)
for i in [0,1]:
axes[-1, i].set_xlabel(r'phase, $\frac{M-M_0}{2\pi}$')
for i in range(4):
axes[i, 0].set_ylabel(_RV_LBL.format(u.km/u.s))
fig.suptitle('High-$K$, unimodal',
x=0.55, y=0.96, fontsize=18)
fig.tight_layout()
fig.subplots_adjust(top=0.92)
fig.savefig(path.join(plot_path, 'highK-unimodal.pdf'))
For my own sake, make the same for unconverged stars:
In [ ]:
np.random.seed(123)
rc = {
'axes.labelsize': 18,
'xtick.labelsize': 14,
'ytick.labelsize': 14
}
subset = full_catalog[np.logical_not(full_catalog['converged'])]
rand_subset = np.random.choice(len(subset), size=8, replace=False)
rand_subset = rand_subset[np.argsort(subset['e'][rand_subset])]
with h5py.File(samples_file, 'r') as jok_f, h5py.File(mcmc_samples_file, 'r') as mcmc_f:
with mpl.rc_context(rc):
fig, axes = plt.subplots(4, 2, figsize=(8, 10), sharex=True)
for i, idx in enumerate(rand_subset):
ax = axes.flat[i]
star = AllStar.get_apogee_id(session, subset[idx]['APOGEE_ID'])
data = star.apogeervdata()
if apogee_id in mcmc_f:
f = mcmc_f
print('mcmc')
else:
f = jok_f
print('thejoker')
samples = JokerSamples.from_hdf5(jok_f[star.apogee_id])
samples.t0 = data.t0
if len(samples) > 1:
sample = MAP_sample(data, samples, joker_pars)
else:
sample = samples[0]
fig = plot_phase_fold(data, sample, ax=ax,
jitter_errorbar=True, label=False)
xlim = ax.get_xlim()
ylim = (data.rv.value.min(), data.rv.value.max())
yspan = ylim[1]-ylim[0]
ylim = ax.set_ylim(ylim[0]-0.35*yspan, ylim[1]+0.35*yspan)
text = ('{0}, '.format(star.apogee_id) +
'$P = {0.value:.2f}$ {0.unit:latex}, '.format(sample['P']) +
'$e = {0:.2f}$'.format(sample['e']))
ax.text(xlim[0] + (xlim[1]-xlim[0])/15,
ylim[1] - (ylim[1]-ylim[0])/20,
text, fontsize=10, va='top', ha='left')
# _ = plot_two_panel(data, samples)
ax.set_xlim(-0.02, 1.02)
for i in [0,1]:
axes[-1, i].set_xlabel(r'phase, $\frac{M-M_0}{2\pi}$')
for i in range(4):
axes[i, 0].set_ylabel(_RV_LBL.format(u.km/u.s))
fig.suptitle('Example stars from the high-$K$, unimodal sample',
x=0.55, y=0.96, fontsize=18)
fig.tight_layout()
fig.subplots_adjust(top=0.92)
In [ ]:
full_catalog['converged'].sum(), len(full_catalog)-full_catalog['converged'].sum()
In [ ]:
# plt.hist(full_catalog['e'][~full_catalog['converged']], bins='auto');
plt.hist(full_catalog['e'], bins='auto');
In [ ]:
emcee_converged = full_catalog[full_catalog['emcee_converged']]
In [ ]:
_path = '../../plots/emcee_converged'
os.makedirs(_path, exist_ok=True)
In [ ]:
with h5py.File(mcmc_samples_file, 'r') as mcmc_f, h5py.File(samples_file, 'r') as f:
for row in emcee_converged:
star = AllStar.get_apogee_id(session, row['APOGEE_ID'])
data = star.apogeervdata()
if star.apogee_id in mcmc_f:
samples = JokerSamples.from_hdf5(mcmc_f[star.apogee_id])
print('mcmc')
else:
samples = JokerSamples.from_hdf5(f[star.apogee_id])
print('thejoker')
samples.t0 = data.t0
fig = plot_two_panel(data, samples,
plot_data_orbits_kw=dict(n_times=16384,
highlight_P_extrema=False))
fig.axes[0].set_title(star.apogee_id)
fig.tight_layout()
fig.savefig(path.join(_path, '{0}.png'.format(star.apogee_id)), dpi=200)
plt.close(fig)
By-eye vetting: these ones are suspicious
In [ ]:
suspicious_ids = ['2M05224382+4300425',
'2M08505498+1156503',
'2M10264342+1340172',
'2M10513288-0250550',
'2M14574438+2106271',
'2M16131259+5043080',
'2M17121495+3211467',
'2M17212080+6003296',
'2M18571262-0328064',
'2M21260907+1100178',
'2M21374395+4304268']
In [ ]:
derp = emcee_converged[~np.isin(emcee_converged['APOGEE_ID'], suspicious_ids)]
In [ ]:
derp = full_catalog
fig, ax = plt.subplots(1, 1, figsize=(6,6))
ax.errorbar(derp['P'], derp['LOGG'],
xerr=derp['P_err'], yerr=derp['LOGG_ERR'],
marker='o', linestyle='none', alpha=0.8)
ax.set_xscale('log')
ax.set_xlim(0.8, 2000)
ax.set_ylim(4., 0)
ax.set_xlabel('P')
ax.set_ylabel('logg')
# -----
fig, ax = plt.subplots(1, 1, figsize=(6,6))
ax.errorbar(derp['P'], derp['e'],
xerr=derp['P_err'], yerr=derp['e_err'],
marker='o', linestyle='none', alpha=0.8)
ax.set_xscale('log')
ax.set_xlim(0.8, 2000)
ax.set_ylim(0, 1)
ax.set_xlabel('P')
ax.set_ylabel('e')
# -----
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
ax = axes[0]
ax.errorbar(derp['M1'], derp['M2_min']/derp['M1'],
xerr=derp['M1_err'], yerr=np.sqrt(derp['M1_err']**2+derp['M2_min_err']**2),
marker='o', linestyle='none', alpha=0.8)
ax.set_xlabel('M1')
ax.set_ylabel('M2/M1')
ax = axes[1]
mass_ratio = derp['M2_min']/derp['M1']
ax.hist(mass_ratio[np.isfinite(mass_ratio)], bins='auto')
ax.set_xlabel('M2/M1')
In [ ]:
with h5py.File(mcmc_samples_file, 'r') as mcmc_f, h5py.File(samples_file, 'r') as f:
for row in derp[rc_mask & (derp['P'] < 20)]:
star = AllStar.get_apogee_id(session, row['APOGEE_ID'])
data = star.apogeervdata()
if star.apogee_id in mcmc_f:
samples = JokerSamples.from_hdf5(mcmc_f[star.apogee_id])
print('mcmc')
else:
samples = JokerSamples.from_hdf5(f[star.apogee_id])
print('thejoker')
samples.t0 = data.t0
fig = plot_two_panel(data, samples,
plot_data_orbits_kw=dict(n_times=16384,
highlight_P_extrema=False))
fig.axes[0].set_title('P = {0:.2f}'.format(samples['P'][0]))
fig.tight_layout()
In [ ]:
derp[rc_mask & (derp['P'] < 20)]
In [ ]: