In [ ]:
import os
from os import path
# Third-party
import astropy.coordinates as coord
from astropy.constants import G
from astropy.io import fits, ascii
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 scipy.stats import beta
from thejoker import JokerSamples
from twoface.config import TWOFACE_CACHE_PATH
from twoface.samples_analysis import 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
from twoface.mass import get_m2_min, mf, period_at_surface, asini
from twoface.log import log as logger
logger.setLevel(100)
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 [ ]:
(unimodal['LOGG'] <= 0).sum()
In [ ]:
unimodal = QTable.read(path.join(table_path, 'highK-unimodal.fits'),
astropy_native=True, character_as_bytes=False)
bimodal = QTable.read(path.join(table_path, 'highK-bimodal.fits'),
astropy_native=True, character_as_bytes=False)
clean_uni = unimodal[(unimodal['clean_flag'] == 0)]
clean_conv_uni = clean_uni[clean_uni['converged']]
clean_bi = bimodal[(bimodal['clean_flag'] == 0)]# & (bimodal['LOGG'] > 0) &
# (bimodal['TEFF'] > 0)]
print('{0} unimodal, {1} clean'.format(len(unimodal), len(clean_uni)))
print('{0} bimodal, {1} clean'.format(len(bimodal)//2, len(clean_bi)//2))
In [ ]:
troup = ascii.read('../../../papers/thejoker-paper/data/troup16-dr12.csv', format='commented_header')
In [ ]:
n_uni = len(clean_uni['APOGEE_ID'])
n_uni_troup = np.isin(clean_uni['APOGEE_ID'], troup['APOGEE_ID']).sum()
In [ ]:
n_bi = len(np.unique(clean_bi['APOGEE_ID']))
n_bi_troup = np.isin(np.unique(clean_bi['APOGEE_ID']), troup['APOGEE_ID']).sum()
In [ ]:
assert np.isin(unimodal['APOGEE_ID'], bimodal['APOGEE_ID']).sum() == 0
In [ ]:
print('{0} new unimodal companions'.format(n_uni - n_uni_troup))
print('{0} new bimodal companions'.format(n_bi - n_bi_troup))
In [ ]:
sb9 = fits.getdata('../../data/sb9.fit')
In [ ]:
c_sb9 = coord.SkyCoord(ra=sb9['RAJ2000']*u.deg, dec=sb9['DEJ2000']*u.deg)
c_uni = coord.SkyCoord(ra=clean_uni['RA']*u.deg, dec=clean_uni['DEC']*u.deg)
c_bi = coord.SkyCoord(ra=clean_bi['RA']*u.deg, dec=clean_bi['DEC']*u.deg)
In [ ]:
idx, sep_uni, _ = c_uni.match_to_catalog_sky(c_sb9)
idx, sep_bi, _ = c_bi.match_to_catalog_sky(c_sb9)
(sep_uni < 2*u.arcsec).sum(), (sep_bi < 2*u.arcsec).sum()
In [ ]:
((sep_uni < 2*u.arcsec) & np.logical_not(np.isin(np.unique(clean_uni['APOGEE_ID']), troup['APOGEE_ID']))).sum()
In [ ]:
uni_colors = ['k', '#777777']
uni_sizes = [14, 8]
uni_zorders = [10, -10]
uni_labels = ['unimodal, MCMC', r'unimodal, The Joker']
uni_style = dict(alpha=0.7, marker='o', color='k', linewidth=0, s=14)
bi_style = dict(alpha=0.5, marker='s', color='tab:blue', s=12, zorder=-10, linewidth=0)
bi_line_style = dict(marker='', linestyle='-', alpha=bi_style['alpha']/3,
color=bi_style['color'])
err_alpha = 0.2
In [ ]:
fig, axes = plt.subplots(2, 2, figsize=(9, 9), sharey='row', sharex='col')
ax = axes[0, 0]
uni_cats = [clean_conv_uni, clean_uni[~clean_uni['converged']]]
for cat, color, zorder, label, s in zip(uni_cats, uni_colors, uni_zorders, uni_labels, uni_sizes):
style = uni_style.copy()
style['color'] = color
style['zorder'] = zorder
style['s'] = s
ax.scatter(cat['P'], cat['e'],
label=label, **style)
ax.errorbar(clean_uni['P'].value, clean_uni['e'],
xerr=clean_uni['P_err'].value, yerr=clean_uni['e_err'],
marker='', linestyle='none', alpha=err_alpha, zorder=-10,
color=uni_style['color'])
# ax.scatter(clean_bi['P'], clean_bi['e'], **bi_style)
# ax.plot(np.vstack((clean_bi['P'][::2].value, clean_bi['P'][1::2].value)),
# np.vstack((np.array(clean_bi['e'])[::2], np.array(clean_bi['e'])[1::2])),
# marker='', linestyle='-', alpha=err_alpha,
# color=bi_style['color'])
ax.set_ylim(-0.02, 1)
ax.set_ylabel(r'eccentricity, $e$')
# ---
ax = axes[1, 1]
for cat, color, zorder, label, s in zip(uni_cats, uni_colors, uni_zorders, uni_labels, uni_sizes):
style = uni_style.copy()
style['color'] = color
style['zorder'] = zorder
style['s'] = s
ax.scatter(cat['TEFF'], cat['LOGG'],
label=label, **style)
# ax.scatter(clean_uni[clean_uni['TINGRC']]['TEFF'], clean_uni[clean_uni['TINGRC']]['LOGG'],
# label='unimodal, red clump', **uni_style_r)
ax.scatter(clean_bi['TEFF'][::2], clean_bi['LOGG'][::2],
label='bimodal', **bi_style)
ax.legend(loc='upper left', fontsize=14, bbox_to_anchor=(-0.02, 1.34))
ax.errorbar(clean_uni['TEFF'], clean_uni['LOGG'],
xerr=clean_uni['TEFF_ERR'], yerr=clean_uni['LOGG_ERR'],
marker='', linestyle='none', alpha=err_alpha, zorder=-10,
color=uni_style['color'])
ax.errorbar(clean_bi['TEFF'][::2], clean_bi['LOGG'][::2],
xerr=clean_bi['TEFF_ERR'][::2], yerr=clean_bi['LOGG_ERR'][::2],
marker='', linestyle='none', alpha=err_alpha, zorder=-100,
color=bi_style['color'])
ax.set_xlim(5900, 3600)
ax.set_ylim(4, 0)
ax.set_xlabel(r'temperature, $T_{\rm eff}$ [K]')
# ---
ax = axes[1, 0]
for cat, color, zorder, label, s in zip(uni_cats, uni_colors, uni_zorders, uni_labels, uni_sizes):
style = uni_style.copy()
style['color'] = color
style['zorder'] = zorder
style['s'] = s
ax.scatter(cat['P'], cat['LOGG'],
label=label, **style)
# ax.scatter(clean_uni[clean_uni['TINGRC']]['P'], clean_uni[clean_uni['TINGRC']]['LOGG'],
# **uni_style_r)
ax.errorbar(uni_converged['P'].value, uni_converged['LOGG'],
xerr=uni_converged['P_err'].value, yerr=uni_converged['LOGG_ERR'],
marker='', linestyle='none', alpha=err_alpha, zorder=-10,
color=uni_style['color'])
ax.scatter(clean_bi['P'], clean_bi['LOGG'], **bi_style)
ax.plot(np.vstack((clean_bi['P'][::2].value, clean_bi['P'][1::2].value)),
np.vstack((np.array(clean_bi['LOGG'])[::2], np.array(clean_bi['LOGG'])[1::2])),
**bi_line_style)
loggs = np.linspace(0, 4, 128)
Psurfs = period_at_surface(1.35*u.Msun, logg=loggs, e=0, M2=0*u.Msun)
ax.plot(Psurfs.to(u.day).value, loggs, marker='', zorder=-100,
color='tab:red', linewidth=2, alpha=0.5)
ax.axvline(period_at_surface(1.35*u.Msun, logg=0, e=0, M2=0*u.Msun).value,
linestyle='-', zorder=-100, color='tab:red', alpha=0.3)
ax.axhline(2.3, linestyle='-', zorder=-100, color='tab:red', alpha=0.3)
ax.set_xscale('log')
ax.set_xlabel('period, $P$ [day]')
ax.set_ylabel(r'surface gravity, $\log g$')
ax.set_xlim(0.8, 2000)
fig.tight_layout()
axes[0, 1].set_visible(False)
fig.savefig(path.join(plot_path, 'P-Teff-logg-e.pdf'))
In [ ]:
uni_cat = clean_uni[(clean_uni['LOGG'] > 2) & np.isfinite(clean_uni['M1'])]
uni_converged = uni_cat[uni_cat['converged']]
bi_logg = np.array([AllStar.get_apogee_id(session, aid).logg
for aid in clean_bi['APOGEE_ID']])
bi_cat = clean_bi[(bi_logg > 2) & np.isfinite(clean_bi['M1']) & np.isfinite(clean_bi['M2_min'])]
print(len(uni_cat), len(bi_cat)//2)
fig, axes = plt.subplots(1, 2, figsize=(11, 4.8))
ax = axes[0]
uni_cats = [uni_cat[uni_cat['converged']], uni_cat[~uni_cat['converged']]]
for cat, color, zorder, label, s in zip(uni_cats, uni_colors, uni_zorders, uni_labels, uni_sizes):
style = uni_style.copy()
style['color'] = color
style['zorder'] = zorder
style['s'] = s
ax.scatter(cat['M1'], cat['M2_min'],
label=label, **style)
ax.errorbar(uni_converged['M1'].value, uni_converged['M2_min'].value,
xerr=uni_converged['M1_err'].value, yerr=uni_converged['M2_min_err'].value,
marker='', linestyle='none', alpha=err_alpha, zorder=-10,
color=uni_style['color'])
ax.scatter(bi_cat['M1'].value, bi_cat['M2_min'].value, **bi_style)
ax.plot(np.vstack((bi_cat['M1'].value[::2], bi_cat['M1'].value[1::2])),
np.vstack((bi_cat['M2_min'].value[::2], bi_cat['M2_min'].value[1::2])),
**bi_line_style)
ax.plot(np.logspace(-3, 1, 1024),
np.logspace(-3, 1, 1024),
marker='', color='#aaaaaa', zorder=-100, linestyle='--')
ax.axhline(0.08, marker='', color='#aaaaaa', zorder=-100, linestyle='--')
# ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim(0.5, 3)
ax.xaxis.set_ticks(np.arange(0.5, 3+0.1, 0.5))
ax.set_ylim(0.01, 10)
ax.set_ylabel(r'$M_{2, {\rm min}}$ ' + '[{0:latex_inline}]'.format(u.Msun))
axr = ax.twinx()
ylim = ax.get_ylim()
axr.set_yscale('log')
axr.set_ylim((ylim[0]*u.Msun).to(u.Mjup).value,
(ylim[1]*u.Msun).to(u.Mjup).value)
axr.set_ylabel('[{0:latex_inline}]'.format(u.Mjup))
ax.set_xlabel('$M_1$ ' + '[{0:latex_inline}]'.format(u.Msun))
ax.text(2.95, 3, r'$M_{2, {\rm min}} = M_1$',
ha='right', va='bottom', fontsize=15, color='#555555')
ax.text(2.95, 0.07, r'$M_{2, {\rm min}} = 0.08\,{\rm M}_\odot$',
ha='right', va='top', fontsize=15, color='#555555')
# -- next panel
ax = axes[1]
for cat, color, zorder, label, s in zip(uni_cats, uni_colors, uni_zorders, uni_labels, uni_sizes):
style = uni_style.copy()
style['color'] = color
style['zorder'] = zorder
style['s'] = s
y = (cat['R1']/cat['a_sini']).decompose()
ax.scatter(cat['q_min'], y,
label=label, **style)
cat = uni_converged
y = (cat['R1']/cat['a_sini']).decompose()
yerr = y * np.sqrt((cat['R1_err']/cat['R1'])**2 + (cat['a_sini_err']/cat['a_sini'])).decompose()
ax.errorbar(cat['q_min'], y,
xerr=cat['q_min_err'], yerr=yerr,
marker='', linestyle='none', alpha=err_alpha, zorder=-10,
color=uni_style['color'])
bi_logg = np.array([AllStar.get_apogee_id(session, aid).logg
for aid in clean_bi['APOGEE_ID']])
bi_cat = clean_bi[bi_logg > 2]
bi_logg = bi_logg[bi_logg > 2]
q = bi_cat['M2_min'] / bi_cat['M1']
asini_ = asini(bi_cat['P'], bi_cat['e'], bi_cat['K'], bi_cat['M1'], bi_cat['M2_min'])
R1 = np.sqrt(G * bi_cat['M1'] / (10**bi_logg*u.cm/u.s**2)).to(u.au)
_y = R1 / asini_
ax.scatter(q, _y, label='bimodal', **bi_style)
ax.plot(np.vstack((q[::2], q[1::2])),
np.vstack((_y[::2], _y[1::2])),
**bi_line_style)
ax.legend(loc='upper right', fontsize=10, ncol=2)
qgrid = np.linspace(1E-2, 1E1)
r1 = 0.49*qgrid**(-2/3.) / (0.6*qgrid**(-2/3.) + np.log(1+qgrid**(-1/3)))
ax.plot(qgrid, r1, marker='', color='#aaaaaa', zorder=-100, linestyle='--')
ax.set_xlim(1e-2, 1E1)
ax.set_ylim(6e-3, 2E0)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel(r'$q_{\rm min}$')
ax.set_ylabel(r'$R_1 / (a\,\sin i)$')
fig.tight_layout()
fig.subplots_adjust(wspace=0.45)
fig.savefig(path.join(plot_path, 'mass.pdf'))
In [ ]:
# Both samples at q>1
bi_cat[::2][((bi_cat['M2_min'].value[::2] / bi_cat['M1'].value[::2]) > 1) &
((bi_cat['M2_min'].value[1::2] / bi_cat['M1'].value[1::2]) > 1)]
In [ ]:
# Both samples at M2<0.08 Msun
bi_cat[::2][((bi_cat['M2_min'][::2]) < 0.08*u.Msun) &
((bi_cat['M2_min'][1::2]) < 0.08*u.Msun)]
In [ ]:
bi_cat[1::2][((bi_cat['M2_min'][::2]) < 0.08*u.Msun) &
((bi_cat['M2_min'][1::2]) < 0.08*u.Msun)]
In [ ]:
uni_cat[(uni_cat['M2_min'] / uni_cat['M1']) > 0.95]
In [ ]:
from thejoker.data import RVData
In [ ]:
# apogee_id = '2M04015112+5316431' # unimodal, above m2=m1 line
apogee_id = '2M01415794+8520066' # unimodal, just under line
with h5py.File(mcmc_samples_file) as f:
samples = JokerSamples.from_hdf5(f[apogee_id])
chain_stats = f[apogee_id]['chain-stats']
print(list(chain_stats.keys()))
print(chain_stats['gelman_rubin'][:])
In [ ]:
m2s = u.Quantity([get_m2_min(rows['M1'][0], samples.get_orbit(i).m_f) for i in range(len(samples))])
In [ ]:
plt.hist(m2s.value)
In [ ]:
_plot_path = '../../plots/three-interesting'
os.makedirs(_plot_path, exist_ok=True)
with h5py.File(samples_file) as f:
for apogee_id in ['2M04015112+5316431', # BH/NS uni
'2M01415794+8520066', # NS/WD??
'2M01231070+1801407', # BH/NS bi
'2M04582049+2232000', # BD bi
'2M20130687+2044173']: # BD bi
samples = JokerSamples.from_hdf5(f[apogee_id])
star = AllStar.get_apogee_id(session, apogee_id)
data = star.apogeervdata()
rows = bi_cat[bi_cat['APOGEE_ID'] == apogee_id]
if len(rows) == 0:
rows = uni_cat[uni_cat['APOGEE_ID'] == apogee_id]
a2s = u.Quantity([get_a2sini(r['P'], r['e'], r['K'], r['M1'], r['M2_min'])
for r in rows])
es = np.array([r['e'] for r in rows])
print('in Troup: {0}\n'.format(apogee_id in troup['APOGEE_ID']) +
'logg = {0:.2f}, Teff = {1:.1f}\n'.format(star.logg, star.teff) +
'a2sini = {0}\n'.format(a2s) +
'peri = {0}'.format(a2s * (1 - es)))
# with LAMOST RV: http://dr3.lamost.org/spectrum/view?obsid=268505101
# data = RVData(t=Time(np.append(data.t.mjd, [56980.]), format='mjd'),
# rv=np.append(data.rv.value, 11.1)*u.km/u.s,
# stddev=np.append(data.stddev.value, 1.)*u.km/u.s)
fig = plot_two_panel(data, samples)
fig.axes[0].set_title(star.apogee_id)
fig.tight_layout()
# fig.savefig(path.join(_plot_path, '{0}.png'.format(apogee_id)), dpi=250)
break
In [ ]:
star.k
In [ ]:
samples['P'].min(), samples['P'].max()
In [ ]:
with h5py.File('../../cache/2M04015112+5316431.hdf5') as f:
new_samples = JokerSamples.from_hdf5(f)
# new_samples = samples
In [ ]:
fig = plot_two_panel(data, new_samples)
fig.axes[0].set_title('{0} - {1} samples'.format(star.apogee_id, len(new_samples)))
fig.tight_layout()
fig = plot_two_panel(data, samples)
fig.axes[0].set_title('{0} - {1} samples'.format(star.apogee_id, len(samples)))
fig.tight_layout()
In [ ]:
m1 = np.exp(star.ness_rg.lnM) * u.Msun
m2s = np.zeros(len(new_samples)) * u.Msun
for i in range(len(new_samples)):
mf = new_samples.get_orbit(i).m_f
m2s[i] = get_m2_min(m1, mf)
m2s_2 = np.zeros(len(samples)) * u.Msun
for i in range(len(samples)):
mf = samples.get_orbit(i).m_f
m2s_2[i] = get_m2_min(m1, mf)
In [ ]:
m2s
In [ ]:
_, bins, _ = plt.hist(m2s.value, bins='auto', normed=True);
plt.hist(m2s_2.value, bins=bins, normed=True);
plt.axvline(m1.value)
plt.gca().xaxis.set_ticks(np.arange(1, 6, 0.5))
In [ ]:
# ptf = ascii.read('/Users/adrian/Downloads/irsa_catalog_search_results.tbl')
# ptf = ptf[(ptf['oid'] == 35322100000378) & (ptf['fid'] == 2)]
# plt.errorbar(ptf['obsmjd'], ptf['mag_autocorr'],
# yerr=ptf['magerr_auto'], linestyle='none', marker='o')
# plt.ylim(14.0, 13.8)
In [ ]:
wise = ascii.read('/Users/adrian/Downloads/irsa_catalog_search_results.tbl')
# wise = wise[(wise['source_id_mf'] == '0206p181_ac51-027396')]
In [ ]:
plt.errorbar(((wise['mjd']) / new_samples['P'].value[0]) % 1, wise['w1mpro_ep'], wise['w1sigmpro_ep'],
linestyle='none', marker='o')
plt.ylim(12.4, 12)
plt.xlim(0, 1)
In [ ]:
In [ ]: