In [ ]:
import os
from os import path
from astropy.time import Time
from astropy.io import fits, ascii
import astropy.units as u
from astropy.table import Table
from astropy.constants import G
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
%matplotlib inline
import h5py
import pandas as pd
from sqlalchemy import func
from scipy.optimize import root
from scipy.stats import scoreatpercentile
import tqdm
from thejoker import JokerSamples
from thejoker.sampler import JokerParams, TheJoker
from thejoker.plot import plot_rv_curves
from twoface.config import TWOFACE_CACHE_PATH
from twoface.db import (db_connect, AllStar, AllVisit, AllVisitToAllStar,
StarResult, Status, JokerRun, initialize_db)
from twoface.data import APOGEERVData
from twoface.plot import plot_data_orbits
In [ ]:
TWOFACE_CACHE_PATH = path.abspath('../cache/')
samples_file = path.join(TWOFACE_CACHE_PATH, 'apogee-jitter.hdf5')
In [ ]:
Session, _ = db_connect(path.join(TWOFACE_CACHE_PATH, 'apogee.sqlite'))
session = Session()
In [ ]:
# stars = session.query(AllStar).join(StarResult, Status)\
# .filter(Status.id.in_([1,2]))\
# .all()
# len(stars)
stars = session.query(AllStar).filter(AllStar.logg < 1.).all()
len(stars)
In [ ]:
%%time
rows = None
with h5py.File(samples_file, 'r') as f:
for star in tqdm.tqdm(stars):
# samples = JokerSamples.from_hdf5(f[star.apogee_id])
K = f[star.apogee_id]['K'][:]
P = f[star.apogee_id]['P'][:]
ecc = f[star.apogee_id]['e'][:]
vsini = [star.vsini] * len(P)
loggs = [star.logg] * len(P)
status = [star.results[0].status_id] * len(P)
apogee_ids = [star.apogee_id] * len(P)
this_rows = list(zip(apogee_ids, loggs, vsini, P, ecc, K, status))
if rows is None:
rows = this_rows
else:
rows = rows + this_rows
tbl_ = np.array(rows, dtype=[('apogee_id', 'U25'), ('logg', float), ('vsini', float),
('P', float), ('e', float), ('K', float),
('status', int)])
df = pd.DataFrame(tbl_)
print(len(df))
In [ ]:
huh = df.groupby('apogee_id')
tbl = huh.filter(lambda x: np.percentile(x['K'], 15) > 0.5)
len(tbl)
In [ ]:
def get_P(m1, m2, logg):
logg = np.asarray(logg)
return (2*np.pi * G**(1/4.) * (m1+m2) / (m1**(3/4.)) * (10**logg*u.cm/u.s**2)**(-3/4)).to(u.day)
In [ ]:
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.plot(tbl['P'], tbl['logg'], marker='.',
linestyle='none', alpha=0.2, ms=2, color='k')
# -- P curves --
logg = np.linspace(0, 4, 1024)
for M1 in [0.5, 1., 2] * u.Msun:
P1 = get_P(M1, 1E-5*u.Msun, logg).value
P2 = get_P(M1, M1, logg).value
ax.fill_betweenx(logg, P1, P2, alpha=0.5, linewidth=0,
label='{0:.1f} {1:latex_inline}'.format(M1.value, M1.unit))
ax.legend(loc='upper right')
ax.set_xscale('log')
ax.set_xlim(1, 32768)
ax.set_ylim(3.55, 0)
ax.set_xlabel('$P$ [day]')
ax.set_ylabel('$\log g$')
fig.tight_layout()
In [ ]:
sub_tbl = tbl[tbl['status'] == 2]
In [ ]:
def make_plot(tbl, c=None, clabel='', **kw):
if c is None:
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.plot(tbl['P'], tbl['logg'], linestyle='none',
marker='.', ms=4, alpha=0.3)
else:
fig, ax = plt.subplots(1, 1, figsize=(7.2, 6))
cs = ax.scatter(tbl['P'], tbl['logg'], c=c,
marker='o', alpha=0.65, s=15, linewidth=0,
cmap='magma_r', **kw)
cb = fig.colorbar(cs)
cb.set_label(clabel)
# -- P curves --
logg = np.linspace(0, 4, 1024)
for M1 in [1.] * u.Msun:
P1 = get_P(M1, 1E-5*u.Msun, logg).value
P2 = get_P(M1, M1, logg).value
ax.fill_betweenx(logg, P1, P2, alpha=0.5, linewidth=0, color='#aaaaaa', zorder=-100,
label='$M_1 = ${0:.1f} {1:latex_inline}'.format(M1.value, M1.unit))
ax.legend(loc='upper right')
ax.set_xscale('log')
ax.set_xlim(1, 32768)
ax.set_ylim(3.55, 0)
ax.set_xlabel('$P$ [day]')
ax.set_ylabel('$\log g$')
fig.tight_layout()
In [ ]:
fig = make_plot(sub_tbl, sub_tbl['K'],
'$K$ [{0:latex_inline}]'.format(u.km/u.s),
norm=LogNorm(1., 50.))
fig
In [ ]:
fig = make_plot(sub_tbl, sub_tbl['e'],
'$e$', vmin=0, vmax=1)
fig
In [ ]:
fig = make_plot(sub_tbl[sub_tbl['vsini']>0], sub_tbl['vsini'][sub_tbl['vsini']>0],
r'$v \sin i$', norm=LogNorm(0.1, 100.))
fig
In [ ]:
corner_tmp = tbl[(tbl['logg'] > 0.) & (tbl['logg'] < 2.)].groupby('apogee_id')
corner = corner_tmp.filter(lambda x: np.percentile(np.log(x['P']), 85) < np.log(get_P(1.5*u.Msun, 0.01*u.Msun, x['logg']).value[0]))
In [ ]:
cmap = plt.get_cmap('magma_r')
norm = LogNorm(0.5, 10)
colors = cmap(norm(corner['K']))
for apid in np.unique(corner['apogee_id']):
mask = corner['apogee_id'] == apid
colors[mask, 3] = min(1, np.sqrt(1/mask.sum()))
In [ ]:
fig, ax = plt.subplots(1, 1, figsize=(7.2, 6))
ax.scatter(corner['P'], corner['logg'], color=colors,
marker='o', s=15, linewidth=0, cmap=cmap, norm=norm)
cs = ax.scatter([-9999, -9999], [-9999, -9999], c=[norm.vmin, norm.vmax],
cmap=cmap, norm=norm) # major hack to put in colorbar!
cb = fig.colorbar(cs)
cb.set_label('$K$ [{0:latex_inline}]'.format(u.km/u.s))
# -- P curves --
logg = np.linspace(0, 4, 1024)
M1 = 1.5*u.Msun
P1 = get_P(M1, 1E-5*u.Msun, logg).value
P2 = get_P(M1, M1, logg).value
ax.fill_betweenx(logg, P1, P2, alpha=0.5, linewidth=0, color='#aaaaaa', zorder=-100,
label='$M_1 = ${0:.1f} {1:latex_inline}'.format(M1.value, M1.unit))
ax.legend(loc='lower right')
ax.set_xscale('log')
ax.set_xlim(1, 1000)
ax.set_ylim(2., -0.02)
ax.set_xlabel('$P$ [day]')
ax.set_ylabel('$\log g$')
fig.tight_layout()
In [ ]:
weirdos = corner.groupby('apogee_id').filter(lambda x: np.percentile(np.log(x['P']), 85) < np.log(10.))
weirdos = np.unique(weirdos['apogee_id']).astype('U25')
len(weirdos)
In [ ]:
with h5py.File(samples_file) as f:
for apid in weirdos:
star = AllStar.get_apogee_id(session, apid)
data = star.apogeervdata()
samples = JokerSamples.from_hdf5(f[apid])
# fig = plot_data_orbits(data, samples, xlim_choice='tight')
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
_ = data.plot(ax=ax)
ax.set_xlabel('BMJD')
ax.set_ylabel('RV [{0:latex_inline}]'.format(u.km/u.s))
fig.tight_layout()
fig.savefig('../plots/astero/{0}.png'.format(apid), dpi=250)
plt.close(fig)
In [ ]:
apokasc = fits.getdata('../data/APOKASC_cat_v4.3.4.fits')
In [ ]:
unq_corner = np.unique(corner['apogee_id']).astype('U25')
in_apokasc = np.isin(unq_corner, np.asarray(apokasc['2MASS_ID']))
In [ ]:
in_apokasc.sum()
In [ ]:
with h5py.File(samples_file) as f:
for apid in unq_corner[in_apokasc]:
star = AllStar.get_apogee_id(session, apid)
data = star.apogeervdata()
samples = JokerSamples.from_hdf5(f[apid])
# fig = plot_data_orbits(data, samples, xlim_choice='tight')
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
_ = data.plot(ax=ax)
ax.set_xlabel('BMJD')
ax.set_ylabel('RV [{0:latex_inline}]'.format(u.km/u.s))
fig.tight_layout()
# fig.savefig('../plots/astero/{0}.png'.format(apid), dpi=250)
# plt.close(fig)
In [ ]:
by_eye = ['2M13424150+2819081', '2M06390580+3444271', '2M13423922+2827574', '2M14035445+2914424',
'2M15184139+0206004', '2M16333828+0010147', '2M18544319+0025050', '2M19403294+2438335',
'2M23435214+6122267', '2M23500237+6053208']
In [ ]:
def plot_phase_fold(data, sample, ax=None):
orbit = sample.get_orbit(0)
P = sample['P']
M0 = sample['M0']
# HACK: hard-set getting the median
t0 = data.t0 + (P/(2*np.pi)*M0).to(u.day, u.dimensionless_angles())
phase = data.phase(P=P, t0=t0)
if ax is None:
ax = plt.gca()
# plot the phase-folded data and orbit
rv_unit = u.km/u.s
ax.errorbar(phase, data.rv.to(rv_unit).value,
data.stddev.to(rv_unit).value,
linestyle='none', marker='o', color='k', markersize=5)
phase_grid = np.linspace(0, 1, 1024)
ax.plot(phase_grid, orbit.radial_velocity(t0 + phase_grid*P),
marker='', zorder=-1, color='#aaaaaa')
ax.set_xlabel('phase')
# ax.set_ylabel('radial velocity [{0:latex_inline}]'.format(rv_unit))
# ax.set_title(r'$\chi^2 = {0:.2f}$'.format(chisq))
return ax.figure
In [ ]:
with h5py.File(samples_file) as f:
for apid in by_eye:
star = AllStar.get_apogee_id(session, apid)
data = star.apogeervdata()
samples = JokerSamples.from_hdf5(f[apid])
fig, axes = plt.subplots(1, 2, figsize=(15, 5), sharey=True)
fig = plot_data_orbits(data, samples, xlim_choice='tight',
highlight_P_extrema=False, n_times=16384, ax=axes[0])
j = samples['P'].argmin()
plot_phase_fold(data, samples[j], ax=axes[1])
fig.tight_layout()
# fig.savefig('../plots/astero/{0}.png'.format(apid), dpi=250)
# plt.close(fig)
# break
In [ ]:
In [ ]: