In [ ]:

    

import os
from os import path
from collections import OrderedDict

import astropy.units as u
from astropy.constants import c
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('apw-notebook')
%matplotlib inline
from scipy.stats import scoreatpercentile

import pandas as pd

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


    

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 [ ]:

    

interesting_group_ids = [1500, 1229,  1515]


    

In [ ]:

    

tbls = OrderedDict()
for gid in interesting_group_ids:
    for f in os.listdir('isochrone_chains/'):
        if f.startswith(str(gid)):
            if gid not in tbls:
                tbls[gid] = OrderedDict()
        
            tbl = pd.read_hdf('isochrone_chains/{0}'.format(f), 
                              key='samples')
            tbls[gid][f.split('.')[0]] = tbl[::64]


    

In [ ]:

    

tbls.keys()


    

In [ ]:

    

def MAD(x):
    return np.median(np.abs(x - np.median(x)))


    

In [ ]:

    

for gid in tbls:
    obs1, obs2 = session.query(Observation).filter(Observation.group_id == gid).all()
    
    print(obs1.simbad_info.preferred_name, obs2.simbad_info.preferred_name)
    print('TGAS ids: {0} {1}'.format(obs1.tgas_source.source_id,
                                     obs2.tgas_source.source_id))
    
    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**15, custom_rv=(rv1, rv_err1))
    icrs2 = obs2.icrs_samples(size=2**15, custom_rv=(rv2, rv_err2))
    
    sep2d = icrs1.separation(icrs2)
    sep3d = icrs1.separation_3d(icrs2)
    R = np.mean([icrs1.distance.value, icrs2.distance.value], axis=0) * u.pc
    chord = 2*R*np.sin(sep2d/2)
    
    print("Angular sep.: {:.2f}".format(sep2d[0]))
    print("Chord length: {:.2f} +/- {:.2f}".format(np.median(chord), 1.5*MAD(chord)))
    print("3D sep.: {:.2f} +/- {:.2f}".format(np.median(sep3d), 1.5*MAD(sep3d)))
    print("parallaxes: {:.2f} +/- {:.2f}, {:.2f} +/- {:.2f}"
          .format(obs1.tgas_source.parallax, obs1.tgas_source.parallax_error,
                  obs2.tgas_source.parallax, obs2.tgas_source.parallax_error))
    
    rep1 = icrs1.represent_as('cartesian')
    rep2 = icrs2.represent_as('cartesian')
    
    dv = (rep1.differentials['s'] - rep2.differentials['s']).norm().to(u.km/u.s)
    print('|∆v| = {:.2f} + {:.2f} - {:.2f}'.format(np.median(dv), 
                                                   np.median(dv.value)-scoreatpercentile(dv.value, 15),
                                                   scoreatpercentile(dv.value, 85)-np.median(dv.value)))
    
    print()


    

In [ ]:

    

d2 = 1000 / np.random.normal(15.68, 0.32, size=100000) 
d1 = 1000 / np.random.normal(17.98, 0.44, size=100000) 
# plt.hist(d1, bins='auto')
# plt.hist(d2, bins='auto', alpha=0.5);
plt.hist(np.abs(d1 - d2), bins='auto');
print(np.median(d1), np.std(d1))
print(np.median(d2), np.std(d2))
print(np.sqrt(np.median(np.abs(d1 - d2))**2 + 3.3**2))


    

In [ ]:

    

for gid in tbls:
    obses = []
    for name in tbls[gid]: 
        obses.append(session.query(Observation).filter(Observation.object == name).one())
    
    fig, axes = plt.subplots(2, 2, figsize=(8,8))
    
    axes[0,0].set_ylabel('$V$ [mag]')
    axes[1,0].set_ylabel('log$g$')
    
    for i,name in enumerate(tbls[gid]):
        tbl = tbls[gid][name]
        common_name = obses[i].simbad_info.preferred_name
        
        axes[0,i].plot(tbl['B_mag'], tbl['V_mag'], marker='.', linestyle='none', alpha=0.5)
        axes[0,i].set_title(common_name)
        axes[0,i].set_xlabel('$B$ [mag]')
        
        axes[1,i].plot(tbl['Teff_0_0'], tbl['logg_0_0'], marker='.', linestyle='none', alpha=0.5)
        axes[1,i].set_xlabel(r'$T_{\rm eff}$ [K]')
        
        print("{0}\n-----------------".format(common_name))
        for par in ['logg', 'Teff', 'mass', 'radius']:
            med = np.median(tbl[par+'_0_0'])
            mad = np.median(np.abs(tbl[par+'_0_0'] - med))
            print('{0} = {1:.2f} +/- {2:.3f}'.format(par, med, 1.5 * mad))
        print("-----------------\n\n")
    fig.tight_layout()


    

In [ ]: