In [1]:
%matplotlib inline
#Lets you play a movie inline in the ipython notebook
%matplotlib nbagg
from astropy import units as u
from astroquery.xmatch import XMatch
from astroquery.vizier import Vizier
import matplotlib.pyplot as plt
import numpy as np
from astropy.table import Table, hstack, join
from astropy.io import ascii
import pyfits
In [2]:
#No limit on rows returned
Vizier.ROW_LIMIT = -1
In [3]:
#main Hipparcos catalog
dat = Vizier.get_catalogs('I/311/hip2')
dat = dat[0]
print dat
In [4]:
#Late G and K stars with parallax errors less than 100%
keep, = np.where((dat['B-V'] > 0.65) & (dat['B-V'] < 0.80 )& ( dat['e_Plx']/dat['Plx'] < 1.0) )# &
print keep.size
In [5]:
#rename columns for querying in Xquery
cut = dat[keep]
cut.rename_column('_RAJ2000','ra')
cut.rename_column('_DEJ2000','dec')
#Now Match to 2Mass Catalog
finder = cut['ra','dec']
#Mask needs to be False or Xquery won't work
finder = Table(finder,masked=False)
In [6]:
#use Xmatch to find the 2MASS ID and IR photometry for these HIP stars
tab = XMatch.query(cat1=finder,cat2='vizier:II/246/out',max_distance=0.5*u.arcsec,colRA1='ra',colDec1='dec')
print 'Total matches = {0:4.0f}, Unique Matches {1:4.0f}'.format(tab['2MASS'].size,np.unique(tab['2MASS']).size)
In [7]:
#print tab
#remove objects without 2MASS phot. (which there are none in this iteration)
keepers, = np.where(tab['2MASS'].mask == False)
tab = tab[keepers]
In [8]:
#Unfortunately Xmatch cuts the rows it doesn't find
from astropy.coordinates import match_coordinates_sky, SkyCoord
c1 = SkyCoord(ra=cut['ra'],dec=cut['dec'],unit='deg')
c2 = SkyCoord(ra=tab['ra'],dec=tab['dec'],unit='deg')
#Get the indexes, which we can match
idx1,idx2, idis1, _ = c1.search_around_sky(c2, 0.5*u.arcsec)
#Create empty string array
cut['2MASS'] = ' '*18
#Fill string with values match in the previous lines using the index match from c2 to c1
cut['2MASS'][idx2] = tab['2MASS']
#Join Tables on 2MASS ID
tab = join(cut,tab,keys=['2MASS'])
#Make sure RA and Dec line up after match (there is an ra and dec component in cut and tab,
# so joining appends _1 and _2, respectively)
ddec = (tab['dec_1']-tab['dec_2'])
dra = tab['ra_1']-tab['ra_2']
print (np.sqrt(ddec**2+(dra*np.cos(np.radians(ddec)))**2)).max()
In [9]:
#Now let us check if the stars have apogee spectra (Apogee targeted HIP stars for calibration early on the 1m)
apdir = 'http://data.sdss3.org/sas/dr12/apogee/spectro/redux/r5/stars/apo1m/hip/'
#This is the location of the Apogee HIP stars list on the Apogee d12 server
apohip = ascii.read(apdir+'r5_stars_apo1m_hip.sha1sum',
names=['num','2MASSfits'],exclude_names=['num'])
In [10]:
#remove all repeat spectra because we are only looking for certain reduced apogee spectra
reducer = np.core.defchararray.find(apohip['2MASSfits'],'apStar-r5')
#Only keep file names which star with 'apStar-r5'
reduced, = np.where(reducer == 0)
apohipfile = apohip[reduced]
In [11]:
#Create 2MASS ID list from Apogee file names
apohipfile['2MASS'] = apohipfile['2MASSfits']
for j,i in enumerate(apohipfile['2MASSfits']):
apohipfile['2MASS'][j] = str(i[12:-11])
#Fix formating so we can match to tab['2MASS'], which is a string column
apohipfile = Table(apohipfile,names=['2MASSfits','2MASS',],dtype=('S48','S16'))
In [12]:
#Combine the Apogee and Hipparcos list
fhipap = join(tab,apohipfile,keys=['2MASS'])
print 'Final number of HIP stars with Apogee Spectra {0:4.0f}'.format(fhipap['2MASS'].size)
#You can now also double check 2MASS and HIP IDs in simbad and they are the same (i.e. the Xmatch query worked)
print fhipap['HIP','2MASS']
In [13]:
#Now let's grab the spectra from the SDSS archive
spec = []
for j,i in enumerate(fhipap['2MASSfits']):
spec.append(pyfits.open(apdir+i)[1])
In [14]:
#add tick marks to plots
def fancy_plot(ax):
#Turn minor ticks on
ax.minorticks_on()
#set the width of the ticks
ax.tick_params(which='both',width=1)
#set the length of the major ticks
ax.tick_params(which='major',length=7)
#set length of the minor ticks
ax.tick_params(which='minor',length=3)
In [15]:
#And for the final scene, a movie
import matplotlib.animation as an
#First create a plot with nothing on it
fig, ax = plt.subplots()
l, = ax.plot(10**(spec[0].header['CRVAL1']+spec[0].header['CDELT1']*np.arange(spec[0].data.size)),
spec[0].data,color=None)
ax.set_xlabel('Wavelength [$\AA$]')
ax.set_ylabel('Flux [10$^{-17}$ erg/s/cm$^2$/$\AA$]')
fancy_plot(ax)
#create a function, which updates the plot information
def update_movie(num):
spec0 = spec[num].data
wave0 = 10.**(spec[num].header['CRVAL1'] + spec[num].header['CDELT1']*np.arange(spec0.size))
l.set_data(wave0,spec0)
l.set_color('black')
ax.set_ylim([0,spec0.max()])
ani = an.FuncAnimation(fig,update_movie,frames=7)
#must have ffmpeg or mencoder installed to save (this step can be nontrivial)
ani.save('test.mp4',writer='ffmpeg',fps=.5)
#press the X at the top right hand corner to stop video playing
In [ ]: