This notebook creates the table '../group_prob_dv.ecsv'


In [ ]:
from os import path

# Third-party
from astropy.io import ascii
from astropy.table import Table
import astropy.coordinates as coord
import astropy.units as u
from astropy.constants import G, c
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
import numpy as np
plt.style.use('apw-notebook')
%matplotlib inline
from scipy.stats import scoreatpercentile
import sqlalchemy

import corner
import emcee
from scipy.integrate import quad
from scipy.misc import logsumexp

from gwb.data import TGASData

from comoving_rv.log import logger
from comoving_rv.db import Session, Base, db_connect
from comoving_rv.db.model import (Run, Observation, TGASSource, SimbadInfo, PriorRV,
                                  SpectralLineInfo, SpectralLineMeasurement, RVMeasurement,
                                  GroupToObservations)

APW RV's


In [ ]:
# base_path = '/Volumes/ProjectData/gaia-comoving-followup/'
base_path = '../../data/'
db_path = path.join(base_path, 'db.sqlite')
engine = db_connect(db_path)
session = Session()

In [ ]:
# Produced by/for scripts/mixture_model.py - see Create-tgas-fits.ipynb
prob_gids = Table.read('../../data/tgas_apw1.fits')['group_id']
pair_probs = np.load('../../data/pair_probs_apw.npy')

In [ ]:
group_ids = np.array([x[0] 
                      for x in session.query(Observation.group_id).distinct().all() 
                      if x[0] is not None and x[0] > 0 and x[0] != 10])
len(group_ids)

In [ ]:
base_q = session.query(Observation).join(RVMeasurement).filter(RVMeasurement.rv != None)

names = ['group_id', 'dv', 'prob']
rows = []
for gid in np.unique(group_ids):
    try:
        gto = session.query(GroupToObservations).filter(GroupToObservations.group_id == gid).one()
        obs1 = base_q.filter(Observation.id == gto.observation1_id).one()
        obs2 = base_q.filter(Observation.id == gto.observation2_id).one()
    except sqlalchemy.orm.exc.NoResultFound:
        print('Skipping group {0}'.format(gid))
        continue
    
    raw_rv_diff = (obs1.measurements[0].x0 - obs2.measurements[0].x0) / 6563. * c.to(u.km/u.s)    
    mean_rv = np.mean([obs1.rv_measurement.rv.value, 
                       obs2.rv_measurement.rv.value]) * obs2.rv_measurement.rv.unit
    
    rv1 = mean_rv + raw_rv_diff/2.
    rv_err1 = obs1.measurements[0].x0_error / 6563. * c.to(u.km/u.s)
    rv2 = mean_rv - raw_rv_diff/2.
    rv_err2 = obs2.measurements[0].x0_error / 6563. * c.to(u.km/u.s)
    
    # Compute point-estimate difference in 3D velocity
    icrs1 = obs1.icrs(with_rv=rv1, lutz_kelker=False)
    icrs2 = obs2.icrs(with_rv=rv2, lutz_kelker=False)
    
    icrs1.set_representation_cls(coord.CartesianRepresentation, coord.CartesianDifferential)
    icrs2.set_representation_cls(coord.CartesianRepresentation, coord.CartesianDifferential)
    
    dv = np.sqrt((icrs1.v_x-icrs2.v_x)**2 + 
                 (icrs1.v_y-icrs2.v_y)**2 + 
                 (icrs1.v_z-icrs2.v_z)**2)
    
    prob = pair_probs[prob_gids == gid]
    assert len(prob) == 1
    
    rows.append((gid, dv.value, prob[0]))
    
dtype = dict(names=names, formats=['i4']+['f8']*(len(names)-1))
tbl = np.array(rows, dtype)
tbl = Table(tbl)
tbl['dv'].unit = u.km/u.s

In [ ]:
base_q = session.query(Observation).join(RVMeasurement).filter(RVMeasurement.rv != None)

dv_15s = []
dv_meds = []
dv_85s = []

for gid in np.unique(group_ids):
    try:
        gto = session.query(GroupToObservations).filter(GroupToObservations.group_id == gid).one()
        obs1 = base_q.filter(Observation.id == gto.observation1_id).one()
        obs2 = base_q.filter(Observation.id == gto.observation2_id).one()
    except sqlalchemy.orm.exc.NoResultFound:
        print('Skipping group {0}'.format(gid))
        continue
    
    raw_rv_diff = (obs1.measurements[0].x0 - obs2.measurements[0].x0) / 6563. * c.to(u.km/u.s)    
    mean_rv = np.mean([obs1.rv_measurement.rv.value, 
                       obs2.rv_measurement.rv.value]) * obs2.rv_measurement.rv.unit
    
    rv1 = mean_rv + raw_rv_diff/2.
    rv_err1 = obs1.measurements[0].x0_error / 6563. * c.to(u.km/u.s)
    rv2 = mean_rv - raw_rv_diff/2.
    rv_err2 = obs2.measurements[0].x0_error / 6563. * c.to(u.km/u.s)
    
    # Compute point-estimate difference in 3D velocity
    icrs1 = obs1.icrs_samples(size=2**16, custom_rv=(rv1,rv_err1))
    icrs2 = obs2.icrs_samples(size=2**16, custom_rv=(rv2,rv_err2))
    
    icrs1.set_representation_cls(coord.CartesianRepresentation, coord.CartesianDifferential)
    icrs2.set_representation_cls(coord.CartesianRepresentation, coord.CartesianDifferential)
    
    dv = np.sqrt((icrs1.v_x-icrs2.v_x)**2 + 
                 (icrs1.v_y-icrs2.v_y)**2 + 
                 (icrs1.v_z-icrs2.v_z)**2)
    
    dv_15, dv_med, dv_85 = scoreatpercentile(dv.value, [15, 50, 85])
    dv_15s.append(dv_15)
    dv_meds.append(dv_med)
    dv_85s.append(dv_85)
    
tbl['dv_15'] = dv_15s*u.km/u.s
tbl['dv_50'] = dv_meds*u.km/u.s
tbl['dv_85'] = dv_85s*u.km/u.s

In [ ]:
base_q = session.query(Observation).join(RVMeasurement).filter(RVMeasurement.rv != None)

sep_2d = []
sep_3d = []
sep_tan = []
d_mins = []
for gid in tbl['group_id']:
    try:
        gto = session.query(GroupToObservations).filter(GroupToObservations.group_id == gid).one()
        obs1 = base_q.filter(Observation.id == gto.observation1_id).one()
        obs2 = base_q.filter(Observation.id == gto.observation2_id).one()
    except sqlalchemy.orm.exc.NoResultFound:
        print('Skipping group {0}'.format(gid))
        continue
    
    icrs1 = obs1.icrs(lutz_kelker=False)
    icrs2 = obs2.icrs(lutz_kelker=False)
    
    sep_2d.append(icrs1.separation(icrs2))
    sep_3d.append(icrs1.separation_3d(icrs2))
    
    R = min(icrs1.distance.value, icrs2.distance.value) * u.pc
    sep_tan.append(2*R*np.sin(icrs1.separation(icrs2)/2))
    d_mins.append(min(icrs1.distance.value, icrs2.distance.value))
    
tbl['sep_2d'] = u.Quantity(sep_2d)
tbl['sep_3d'] = u.Quantity(sep_3d)
tbl['chord_length'] = u.Quantity(sep_tan)
tbl['d_min'] = d_mins * icrs1.distance.unit

In [ ]:
tbl.write('group_prob_dv.ecsv', format='ascii.ecsv', overwrite=True)

RAVE


In [ ]:
# Produced by/for scripts/mixture_model.py - see Create-tgas-fits.ipynb
rave_prob_gids = Table.read('../../data/tgas_rave1.fits')['group_id']
rave_pair_probs = np.load('../../data/pair_probs_rave.npy')

In [ ]:
tgas = TGASData('../../../gaia-comoving-stars/data/stacked_tgas.fits')

In [ ]:
star = ascii.read('../../../gaia-comoving-stars/paper/t1-1-star.txt')
rave_stars = star[(star['group_size'] == 2) & (~star['rv'].mask)]
rave_stars = rave_stars.group_by('group_id')

Get only ones where both stars have RV measurements


In [ ]:
group_idx = np.array([i for i,g in enumerate(rave_stars.groups) if len(g) > 1])
rave_stars = rave_stars.groups[group_idx]

In [ ]:
names = ['group_id', 'prob']
all_rows = []
for gid in np.unique(rave_stars['group_id']):
    rows = rave_stars[rave_stars['group_id'] == gid]
    
    if len(rows) != 2:
        print("skipping group {0} ({1})".format(gid, len(rows)))
        continue
        
    i1 = np.where(tgas._data['source_id'] == rows[0]['tgas_source_id'])[0][0]
    i2 = np.where(tgas._data['source_id'] == rows[1]['tgas_source_id'])[0][0]
    
    star1 = tgas[i1]
    star2 = tgas[i2]

    star1._rv = rows[0]['rv']
    star2._rv = rows[1]['rv']
    star1._rv_err = rows[0]['erv']
    star2._rv_err = rows[1]['erv']
    
    prob = rave_pair_probs[rave_prob_gids == gid]
    assert len(prob) == 1
    
    all_rows.append((gid, prob))
    
dtype = dict(names=names, formats=['i4']+['f8']*(len(names)-1))
tbl = np.array(all_rows, dtype)
tbl = Table(tbl)

In [ ]:
dv_pts = []
dv_15s = []
dv_meds = []
dv_85s = []
sep_2d = []
sep_3d = []
sep_tan = []
d_mins = []

for gid in np.unique(rave_stars['group_id']):
    rows = rave_stars[rave_stars['group_id'] == gid]
    
    if len(rows) != 2:
        print("skipping group {0} ({1})".format(gid, len(rows)))
        continue
        
    i1 = np.where(tgas._data['source_id'] == rows[0]['tgas_source_id'])[0][0]
    i2 = np.where(tgas._data['source_id'] == rows[1]['tgas_source_id'])[0][0]
    
    star1 = tgas[i1]
    star2 = tgas[i2]
    
    icrs1 = star1.get_icrs(rows[0]['rv']*u.km/u.s)
    icrs2 = star2.get_icrs(rows[1]['rv']*u.km/u.s)
    d_mins.append(min(icrs1.distance.value, icrs2.distance.value))
    
    R = min(icrs1.distance.value, icrs2.distance.value) * u.pc
    sep_tan.append(2*R*np.sin(icrs1.separation(icrs2)/2))
    
    # separations    
    sep_2d.append(icrs1.separation(icrs2))
    sep_3d.append(icrs1.separation_3d(icrs2))
    
    icrs1.set_representation_cls(coord.CartesianRepresentation, coord.CartesianDifferential)
    icrs2.set_representation_cls(coord.CartesianRepresentation, coord.CartesianDifferential)
    
    dv_pt = np.sqrt((icrs1.v_x-icrs2.v_x)**2 + 
                    (icrs1.v_y-icrs2.v_y)**2 + 
                    (icrs1.v_z-icrs2.v_z)**2)
    dv_pts.append(dv_pt)
    
    # Compute difference in 3D velocity for samples
    icrs1 = star1.get_icrs_samples(size=2**16, rv=rows[0]['rv']*u.km/u.s, 
                                   rv_err=rows[0]['erv']*u.km/u.s)
    icrs2 = star2.get_icrs_samples(size=2**16, rv=rows[1]['rv']*u.km/u.s, 
                                   rv_err=rows[1]['erv']*u.km/u.s)
    
    icrs1.set_representation_cls(coord.CartesianRepresentation, coord.CartesianDifferential)
    icrs2.set_representation_cls(coord.CartesianRepresentation, coord.CartesianDifferential)
    
    dv = np.sqrt((icrs1.v_x-icrs2.v_x)**2 + 
                 (icrs1.v_y-icrs2.v_y)**2 + 
                 (icrs1.v_z-icrs2.v_z)**2)
    
    dv_15, dv_med, dv_85 = scoreatpercentile(dv.value, [15, 50, 85])
    dv_15s.append(dv_15)
    dv_meds.append(dv_med)
    dv_85s.append(dv_85)

In [ ]:
tbl['dv'] = u.Quantity(dv_pts)
tbl['dv_15'] = dv_15s*u.km/u.s
tbl['dv_50'] = dv_meds*u.km/u.s
tbl['dv_85'] = dv_85s*u.km/u.s
tbl['sep_2d'] = u.Quantity(sep_2d)
tbl['sep_3d'] = u.Quantity(sep_3d)
tbl['chord_length'] = u.Quantity(sep_tan)
tbl['d_min'] = d_mins * u.pc

In [ ]:
tbl.write('group_prob_dv_rave.ecsv', format='ascii.ecsv', overwrite=True)

In [ ]: