In [ ]:
# Third-party
from astropy.io import fits
from astropy.table import Table, vstack
import astropy.coordinates as coord
import astropy.units as u
from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt
import numpy as np
plt.style.use('apw-notebook')
%matplotlib inline
In [ ]:
tgas = Table.read('../../gaia-comoving-stars/data/stacked_tgas.fits')
In [ ]:
lamost = fits.getdata('../data/lamost_dr3_stellar.fits')
lamost = lamost[np.isfinite(lamost['rv'])]
In [ ]:
print(lamost.dtype.names)
In [ ]:
lamost_c = coord.ICRS(ra=lamost['ra']*u.degree, dec=lamost['dec']*u.degree)
In [ ]:
tgas_c = coord.ICRS(ra=tgas['ra']*u.deg - tgas['pmra']*u.mas/u.yr*15*u.yr,
dec=tgas['dec']*u.deg - tgas['pmdec']*u.mas/u.yr*15*u.yr)
In [ ]:
idx, sep2d, _ = coord.match_coordinates_sky(lamost_c, tgas_c)
In [ ]:
plt.hist(sep2d.arcsec, bins=np.logspace(-2, 3, 32))
plt.xscale('log')
plt.yscale('log')
In [ ]:
sep_cut = sep2d < 2.5*u.arcsec
tgas_idx = idx[sep_cut]
In [ ]:
idx, = np.where(sep_cut)
In [ ]:
idx.shape, tgas_idx.shape
In [ ]:
star_tbl = Table.read('../../gaia-comoving-stars/paper/t1-1-star.txt', format='ascii.csv')
pair_tbl = Table.read('../../gaia-comoving-stars/paper/t1-2-pair.txt', format='ascii.csv')
In [ ]:
star_tbl_ix = np.in1d(star_tbl['tgas_source_id'], tgas[tgas_idx]['source_id'])
group_ids = []
for g in star_tbl[star_tbl_ix].group_by('group_id').groups:
if len(g) < 2:
continue
group_ids.append(g['group_id'][0])
group_ids = np.array(group_ids)
In [ ]:
len(group_ids), len(dfehs)
In [ ]:
chord_lens = []
dfehs = []
dfeh_errs = []
drvs = []
drv_errs = []
all_tgas_rows = []
for gid in group_ids[group_ids > 319]: # group index at which the 2-member groups start
members = star_tbl[star_tbl['group_id'] == gid]
lamost_rows = []
tgas_rows = []
for sid in members['tgas_source_id']:
hi = idx[tgas[tgas_idx]['source_id'] == sid]
if len(hi) < 1:
continue
lamost_rows.append(Table(lamost[hi]))
tgas_rows.append(tgas[tgas_idx][tgas[tgas_idx]['source_id'] == sid])
tgas_group = vstack(tgas_rows)
lamost_group = vstack(lamost_rows)
pair_group = pair_tbl[pair_tbl['group_id'] == gid]
if len(pair_group) > 1:
pair_group = pair_group[1:2]
chord_len = (pair_group['angsep']*u.arcmin * np.mean(1000./tgas_group['parallax'])*u.pc).to(u.pc,
u.dimensionless_angles())
sky_sep = (pair_group['angsep']*u.arcmin).to(u.degree)
print("{0:.2f} {1:.2f}".format(sky_sep[0], chord_len[0]))
chord_lens.append(chord_len)
d_feh = lamost_group['feh'][0] - lamost_group['feh'][1]
d_feh_err = np.sqrt(lamost_group['feh_err'][0]**2 + lamost_group['feh_err'][1]**2)
dfehs.append(d_feh)
dfeh_errs.append(d_feh_err)
drvs.append(lamost_group['rv'][0] - lamost_group['rv'][1])
drv_errs.append(np.sqrt(lamost_group['rv_err'][0]**2 + lamost_group['rv_err'][1]**2))
# if abs(rv[0]-rv[1]) < 8:
# all_tgas_rows += tgas_rows
print('∆[Fe/H] = {0:.2f} +/- {1:.2f}'.format(d_feh, d_feh_err))
print('∆rv = {0:.2f} +/- {1:.2f}'.format(drvs[-1], drv_errs[-1]))
print()
chord_lens = u.Quantity(chord_lens)
In [ ]:
fig, axes = plt.subplots(1, 2, figsize=(10,5), sharex=True)
# ax.axvline(0, linewidth=1, alpha=0.2, zorder=-10)
axes[0].axhline(0, linewidth=1, alpha=0.2, zorder=-10)
axes[1].axhline(0, linewidth=1, alpha=0.2, zorder=-10)
axes[0].scatter(chord_lens.to(u.pc).value, drvs, marker='o', s=8, zorder=100)
axes[0].errorbar(chord_lens.to(u.pc).value, drvs, yerr=drv_errs,
marker='None', linestyle='none', alpha=1, c='k', ecolor='#aaaaaa')
axes[0].set_xlim(1e-3, 10)
axes[0].set_ylim(-75, 75)
axes[0].set_xscale('log')
axes[0].set_ylabel(r'$\Delta {\rm RV}\,[{\rm km}\,{\rm s}^{-1}]$')
axes[1].scatter(chord_lens.to(u.pc).value, dfehs, marker='o', s=8, zorder=100)
axes[1].errorbar(chord_lens.to(u.pc).value, dfehs, yerr=dfeh_errs,
marker='None', linestyle='none', alpha=1, c='k', ecolor='#aaaaaa')
axes[1].set_xlim(1e-3, 10)
axes[1].set_ylim(-0.75, 0.75)
axes[1].set_xscale('log')
axes[1].set_ylabel(r'$\Delta [{\rm Fe}/{\rm H}]$')
fig.suptitle('Comoving stars in LAMOST DR3', fontsize=20, y=0.95, x=0.55)
fig.tight_layout()
fig.subplots_adjust(top=0.9)
In [ ]:
from matplotlib import colors
In [ ]:
fig, ax = plt.subplots(1, 1, figsize=(7, 6))
ax.axvline(0, linewidth=1, alpha=0.2, zorder=-10, linestyle='--')
ax.axhline(0, linewidth=1, alpha=0.2, zorder=-10, linestyle='--')
cc = ax.scatter(drvs, dfehs, c=chord_lens, s=10, marker='o', zorder=100, alpha=0.75,
norm=colors.LogNorm(vmin=1E-3, vmax=1E1), cmap='magma_r')
ax.errorbar(drvs, dfehs, xerr=drv_errs, yerr=dfeh_errs,
marker='None', linestyle='none', alpha=1, c='k', ecolor='#777777')
ax.set_xlim(-50, 50)
ax.set_ylim(-0.6, 0.6)
ax.set_xlabel(r'$\Delta {\rm RV}\,[{\rm km}\,{\rm s}^{-1}]$')
ax.set_ylabel(r'$\Delta [{\rm Fe}/{\rm H}]$')
cb = fig.colorbar(cc)
cb.set_label('chord length [pc]')
ax.set_title('Comoving stars in LAMOST DR3', fontsize=20)
fig.tight_layout()
fig.subplots_adjust(top=0.9)
fig.savefig('lamost.pdf')
In [ ]: