# Licensed under the MIT License - see LICENSE.rst
from __future__ import (absolute_import, division, print_function,
unicode_literals)
import os
import numpy as np
from astropy.io import ascii
from astropy.table import Column, join
from astropy.coordinates import SkyCoord
import astropy.units as u
import matplotlib.pyplot as plt
from astropy.constants import R_sun
from .gaia import Nprime_fov, sigma_fov
__all__ = ['get_table_ms']
environment_variable = 'MRSPOC_DATA_DIR'
data_dir_path = os.getenv(environment_variable, '.')
tgas_path = os.path.join(os.path.abspath(data_dir_path),
'tgas_bright_g_lt_12.tsv')
hipparcos_path = os.path.join(os.path.abspath(data_dir_path),
'hipparcos.tsv')
boyajian_path = os.path.join(os.path.abspath(data_dir_path),
'boyajian2012.csv')
mamajek_path = os.path.join(os.path.abspath(data_dir_path),
'spt_mamajek.txt')
[docs]def get_table_ms(plot=True, ax=None):
"""
Open the TGAS catalog for all stars brighter than G < 12,
make CMD cuts to flag just the main sequence stars.
Parameters
----------
plot : bool (optional)
Make a plot of the CMD and color/mag cuts.
ax : `~matplotlib.pyplot.Axes` (optional)
If ``ax`` is not `None`, make the plot on ``ax``.
Returns
-------
table : `~astropy.table.Table`
Table of TGAS sources, magnitudes, parallaxes, distances, and
anticipated astrometric uncertainties.
ms : `~numpy.ndarray`
Boolean array, ``True`` for rows of ``table`` where the star is a
main sequence star within the color/mag cuts.
"""
table = ascii.read(tgas_path, delimiter=';', data_start=3)
# floatify:
table['BTmag'] = table['BTmag'].astype(float)
table['VTmag'] = table['VTmag'].astype(float)
# Compute the galactic latitude of each star, add to table
coords = SkyCoord(ra=table['RA_ICRS'] * u.deg,
dec=table['DE_ICRS'] * u.deg, frame='icrs')
galactic_coords = coords.transform_to('galactic')
abs_galactic_latitude = abs(galactic_coords.b).degree
table.add_column(Column(data=abs_galactic_latitude, name='b'))
# Compute distance, CMD
def color_cut(b_minus_v):
return -9. + 4.9 * b_minus_v
parallax_mas = table['Plx']
Vmag = table['VTmag']
bt_minus_vt = table['BTmag'] - table['VTmag']
parallax_arcsec = parallax_mas / 1000
dist_pc = 1. / parallax_arcsec
# Add astrometric uncertainty column to table
table.add_column(Column(data=sigma_fov(table['<Gmag>']), name='sigma_fov'))
# Add a distance column to the table:
table.add_column(Column(data=dist_pc * u.pc, name='distance'))
# Add a Nfov column to the table:
table.add_column(Column(data=Nprime_fov(abs_galactic_latitude),
name='N_fov'))
M_V = Vmag - 5 * (np.log10(dist_pc) + 1)
b_minus_v_lower = 0.6 # 0.64 # (B-V)_sun = 0.65
b_minus_v_upper = 2
main_sequence = ((np.abs(M_V - color_cut(bt_minus_vt)) < 1.) &
(bt_minus_vt > b_minus_v_lower) &
(bt_minus_vt < b_minus_v_upper))
main_sequence_table = table[main_sequence]
# Now match the B-V color table from HIPPARCOS to the main sequence TGAS table
hipparcos_table = ascii.read(hipparcos_path, delimiter=';', header_start=0,
data_start=3)
hipparcos_table.add_index("HIP")
main_sequence_table['HIP'][main_sequence_table['HIP'].mask] = 0
main_sequence_color_table = join(main_sequence_table, hipparcos_table,
keys='HIP')
# Cut again by the color cuts, this time with the real Johnson B and V,
# rather than Tycho magnitudes:
main_sequence = ((main_sequence_color_table['B-V'].data.data < b_minus_v_upper) &
(main_sequence_color_table['B-V'].data.data > b_minus_v_lower))
main_sequence_color_table = main_sequence_color_table[main_sequence]
# Add in stellar radii with color-radius relation from Boyajian 2012
R_star = bv_to_radius(main_sequence_color_table['B-V'].data.data)
main_sequence_color_table.add_column(Column(data=R_star, name='R_star'))
# Add in a column of interferometric angular diameters from
# Boyajian 2012 where available:
boyajian = ascii.read(boyajian_path)
ang_diams = np.zeros(len(main_sequence_color_table))
for row in boyajian:
ang_diams[row['HIP'] == main_sequence_color_table['HIP']] = row['D(UD)']
main_sequence_color_table.add_column(Column(data=ang_diams,
name='angular_diameter'))
boyajian_radii = main_sequence_color_table['angular_diameter'] != 0
half_angle = (main_sequence_color_table['angular_diameter'][boyajian_radii]
* u.marcsec/2)
distance_pc = (main_sequence_color_table['Plx_1'][
boyajian_radii].data.data / 1000)**-1 * u.pc
measured_radii = distance_pc * np.tan(half_angle)
R_star[boyajian_radii] = measured_radii
# In radius reference column, `1`==color-radius estimate;
# `2`==interferometric measurement
refs = np.ones(len(R_star))
refs[boyajian_radii] = 2
main_sequence_color_table.add_column(Column(data=refs, name='rstar_ref'))
# Add column containing approximate stellar effective temperatures based
# on B-V -> T_eff table from Eric Mamajek:
# http://www.pas.rochester.edu/~emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt
mamajek = ascii.read(mamajek_path, format='commented_header')
bv_to_teff = lambda bv: np.interp(bv, mamajek['B-V'], mamajek['Teff'])
approx_teffs = bv_to_teff(main_sequence_color_table['B-V'])
main_sequence_color_table.add_column(Column(data=approx_teffs, name='Teff'))
if plot:
if ax is None:
ax = plt.gca()
polygon_x = [0.6, 0.6, 2.0, 2.0, 0.6]
polygon_y = [color_cut(0.6) - 1, color_cut(0.6) + 1,
color_cut(2) + 1, color_cut(2) - 1,
color_cut(0.6) - 1]
H, xedges, yedges = np.histogram2d(bt_minus_vt[abs(bt_minus_vt) > 1e-3],
M_V[abs(bt_minus_vt) > 1e-3],
bins=1000)
extent = [xedges.min(), xedges.max(), yedges.max(), yedges.min()]
ax.imshow(np.log10(H.T), extent=extent, cmap=plt.cm.Greys, aspect=0.2)
ax.plot(polygon_x, polygon_y, lw=2, color='r', ls='--')
ax.set(xlim=[-0.5, 3], ylim=[2, -15],
ylabel='$M_{VT}$', xlabel="BT - VT")
return main_sequence_color_table
def bv_to_radius(b_minus_v):
"""
Estimate radii for stars on the main sequence using their ``B-V`` color,
using a simple relation calibrated on interferometry by Boyajian et al. 2012
Parameters
----------
b_minus_v : float
B-V color.
Returns
-------
radius : `~astropy.units.Quantity`
Stellar radius.
"""
# Boyajian 2012
X = b_minus_v
a0 = 0.3830
a1 = 0.9907
a2 = -0.6038
Y = 0
# Ignore metallicity
a3 = 0
a4 = 0
a5 = 0
return (a0 + a1 * X + a2 * X ** 2 + a3 * X * Y +
a4 * Y + a5 * Y ** 2) * R_sun