Data files and code to work with data from Koposov et al. (2010)

koposov10.py

#Module to read and process the Koposov et al. (2010) data
import os
import numpy
from galpy.util import bovy_coords
_DATADIR='../data/'
# Koposov et al. (2010) transformation matrix
_TKOP= numpy.zeros((3,3))
_TKOP[0,:]= [-0.4776303088,-0.1738432154,0.8611897727]
_TKOP[1,:]= [0.510844589,-0.8524449229,0.111245042]
_TKOP[2,:]= [0.7147776536,0.4930681392,0.4959603976]
def read_vr():
    return numpy.loadtxt(os.path.join(_DATADIR,'koposov10_rv.txt'),
                         comments='#',delimiter='|')

def read_pos():
    return numpy.loadtxt(os.path.join(_DATADIR,'koposov10_pos.txt'),
                         comments='#',delimiter='|')

def read_dist():
    return numpy.loadtxt(os.path.join(_DATADIR,'koposov10_dist.txt'),
                         comments='#',delimiter='|')

def read_pm():
    return numpy.loadtxt(os.path.join(_DATADIR,'koposov10_pm.txt'),
                         comments='#',delimiter='|')

def read_vxvv():
    #First read all of the data
    rvdata= read_vr()
    posdata= read_pos()
    distdata= read_dist()
    pmdata= read_pm()
    ndata= rvdata.shape[0]+posdata.shape[0]+distdata.shape[0]+pmdata.shape[0]
    #Set up output
    missing= -9999.99
    large_err= 1000.
    vxvv= numpy.zeros((ndata,6))+missing
    vxvv_err= numpy.zeros((ndata,6))+large_err
    #Load all of the data
    #RV
    istart= 0
    vxvv[istart:istart+rvdata.shape[0],0]= rvdata[:,0]
    vxvv[istart:istart+rvdata.shape[0],5]= rvdata[:,2]
    vxvv_err[istart:istart+rvdata.shape[0],0]= 0.01
    vxvv_err[istart:istart+rvdata.shape[0],5]= rvdata[:,3]
    istart+= rvdata.shape[0]
    #Position
    vxvv[istart:istart+posdata.shape[0],0]= posdata[:,0]
    vxvv[istart:istart+posdata.shape[0],1]= posdata[:,1]
    vxvv_err[istart:istart+posdata.shape[0],0]= 0.01
    vxvv_err[istart:istart+posdata.shape[0],1]= posdata[:,2]
    istart+= posdata.shape[0]
    #Distance
    vxvv[istart:istart+distdata.shape[0],0]= distdata[:,0]
    vxvv[istart:istart+distdata.shape[0],2]= distdata[:,1]
    vxvv_err[istart:istart+distdata.shape[0],0]= 0.01
    vxvv_err[istart:istart+distdata.shape[0],2]= distdata[:,2]
    istart+= distdata.shape[0]
    #PM
    vxvv[istart:istart+pmdata.shape[0],0]= pmdata[:,0]
    vxvv[istart:istart+pmdata.shape[0],3]= pmdata[:,1] #cos(phi2) doesn't matter
    vxvv[istart:istart+pmdata.shape[0],4]= pmdata[:,2]
    vxvv_err[istart:istart+pmdata.shape[0],0]= 0.01
    vxvv_err[istart:istart+pmdata.shape[0],3]= pmdata[:,3]
    vxvv_err[istart:istart+pmdata.shape[0],4]= pmdata[:,3]
    return (vxvv,vxvv_err)

@bovy_coords.scalarDecorator
@bovy_coords.degreeDecorator([0,1],[0,1])
def lb_to_phi12(l,b,degree=False):
    #Convert l and b to phi1,phi2 coordinates
    #First convert to ra and dec
    radec= bovy_coords.lb_to_radec(l,b)
    ra= radec[:,0]
    dec= radec[:,1]
    XYZ= numpy.array([numpy.cos(dec)*numpy.cos(ra),
                      numpy.cos(dec)*numpy.sin(ra),
                      numpy.sin(dec)])
    phiXYZ= numpy.dot(_TKOP,XYZ)
    phi2= numpy.arcsin(phiXYZ[2])
    phi1= numpy.arctan2(phiXYZ[1],phiXYZ[0])
    phi1[phi1<0.]+= 2.*numpy.pi
    return numpy.array([phi1,phi2]).T

@bovy_coords.scalarDecorator
@bovy_coords.degreeDecorator([0,1],[0,1])
def phi12_to_lb(phi1,phi2,degree=False):
    #Convert phi1,phi2 to l,b coordinates
    phiXYZ= numpy.array([numpy.cos(phi2)*numpy.cos(phi1),
                         numpy.cos(phi2)*numpy.sin(phi1),
                         numpy.sin(phi2)])
    eqXYZ= numpy.dot(_TKOP.T,phiXYZ)
    #Get ra dec
    dec= numpy.arcsin(eqXYZ[2])
    ra= numpy.arctan2(eqXYZ[1],eqXYZ[0])
    ra[ra<0.]+= 2.*numpy.pi
    #Now convert to l,b
    return bovy_coords.radec_to_lb(ra,dec)

def test_lb_to_phi12():
    #Test the coordinate conversion, according to Willett et al. (2009)
    # (l,b) = (216,32) is approximately (phi12) = (-70,-2.)
    phi1, phi2= lb_to_phi12(216.,32.,degree=True)
    assert numpy.fabs(phi1-290.) < 5., 'phi1 at the high-v end of the stream does not agree between Koposov and Willett'
    assert numpy.fabs(phi2+2.) < 2., 'phi2 at the high-v end of the stream does not agree between Koposov and Willett'
    # (l,b) = (160,60) is approximately (phi12) = (-28,0)
    phi1, phi2= lb_to_phi12(160.,60.,degree=True)
    assert numpy.fabs(phi1-332.) < 7., 'phi1 at the lower-v end of the stream does not agree between Koposov and Willett'
    assert numpy.fabs(phi2) < 2., 'phi2 at the lower-v end of the stream does not agree between Koposov and Willett'
    # Also test for arrays
    s= numpy.ones(3)
    phi1phi2= lb_to_phi12(160.*s,60.*s,degree=True)
    phi1= phi1phi2[:,0]
    phi2= phi1phi2[:,1]
    assert numpy.all(numpy.fabs(phi1-332.) < 7.), 'phi1 at the lower-v end of the stream does not agree between Koposov and Willett'
    assert numpy.all(numpy.fabs(phi2) < 2.), 'phi2 at the lower-v end of the stream does not agree between Koposov and Willett'
    return None

def test_phi12_to_lb():
    # Test that lb_to_phi12 and phi12_to_lb are each other's inverse
    phi1, phi2= lb_to_phi12(216.,32.,degree=True)
    l,b= phi12_to_lb(phi1,phi2,degree=True)
    assert numpy.fabs(l-216.) < 0.001, 'phi12_to_lb not the inverse of lb_to_phi12'
    assert numpy.fabs(b-32.) < 0.001, 'phi12_to_lb not the inverse of lb_to_phi12'
    phi1, phi2= lb_to_phi12(160.,60.,degree=True)
    l,b= phi12_to_lb(phi1,phi2,degree=True)
    assert numpy.fabs(l-160.) < 0.001, 'phi12_to_lb not the inverse of lb_to_phi12'
    assert numpy.fabs(b-60.) < 0.001, 'phi12_to_lb not the inverse of lb_to_phi12'
    # Also test for arrays
    s= numpy.ones(3)
    lb= phi12_to_lb(phi1*s,phi2*s,degree=True)
    l= lb[:,0]
    b= lb[:,1]
    assert numpy.all(numpy.fabs(l-160.) < 0.001), 'phi12_to_lb not the inverse of lb_to_phi12'
    assert numpy.all(numpy.fabs(b-60.) < 0.001), 'phi12_to_lb not the inverse of lb_to_phi12'
    return None

@bovy_coords.scalarDecorator
@bovy_coords.degreeDecorator([2,3],[])
def pmllpmbb_to_pmphi12(pmll,pmbb,l,b,degree=False):
    #First go to ra and dec
    radec= bovy_coords.lb_to_radec(l,b)
    ra= radec[:,0]
    dec= radec[:,1]
    pmradec= bovy_coords.pmllpmbb_to_pmrapmdec(pmll,pmbb,l,b,degree=False)
    pmra= pmradec[:,0]
    pmdec= pmradec[:,1]
    #Now transform ra,dec pm to phi1,phi2
    phi12= lb_to_phi12(l,b,degree=False)
    phi1= phi12[:,0]
    phi2= phi12[:,1]
    #Build A and Aphi matrices
    A= numpy.zeros((3,3,len(ra)))
    A[0,0]= numpy.cos(ra)*numpy.cos(dec)
    A[0,1]= -numpy.sin(ra)
    A[0,2]= -numpy.cos(ra)*numpy.sin(dec)
    A[1,0]= numpy.sin(ra)*numpy.cos(dec)
    A[1,1]= numpy.cos(ra)
    A[1,2]= -numpy.sin(ra)*numpy.sin(dec)
    A[2,0]= numpy.sin(dec)
    A[2,1]= 0.
    A[2,2]= numpy.cos(dec)
    AphiInv= numpy.zeros((3,3,len(ra)))
    AphiInv[0,0]= numpy.cos(phi1)*numpy.cos(phi2)
    AphiInv[0,1]= numpy.cos(phi2)*numpy.sin(phi1)
    AphiInv[0,2]= numpy.sin(phi2)
    AphiInv[1,0]= -numpy.sin(phi1)
    AphiInv[1,1]= numpy.cos(phi1)
    AphiInv[1,2]= 0.
    AphiInv[2,0]= -numpy.cos(phi1)*numpy.sin(phi2)
    AphiInv[2,1]= -numpy.sin(phi1)*numpy.sin(phi2)
    AphiInv[2,2]= numpy.cos(phi2)
    TA= numpy.dot(_TKOP,numpy.swapaxes(A,0,1))
    #Got lazy...
    trans= numpy.zeros((2,2,len(ra)))
    for ii in range(len(ra)):
        trans[:,:,ii]= numpy.dot(AphiInv[:,:,ii],TA[:,:,ii])[1:,1:]
    return (trans*numpy.array([[pmra,pmdec],[pmra,pmdec]])).sum(1).T

@bovy_coords.scalarDecorator
@bovy_coords.degreeDecorator([2,3],[])
def pmphi12_to_pmllpmbb(pmphi1,pmphi2,phi1,phi2,degree=False):
    #First go from phi12 to ra and dec
    lb= phi12_to_lb(phi1,phi2)
    radec= bovy_coords.lb_to_radec(lb[:,0],lb[:,1])
    ra= radec[:,0]
    dec= radec[:,1]
    #Build A and Aphi matrices
    AInv= numpy.zeros((3,3,len(ra)))
    AInv[0,0]= numpy.cos(ra)*numpy.cos(dec)
    AInv[0,1]= numpy.sin(ra)*numpy.cos(dec)
    AInv[0,2]= numpy.sin(dec)
    AInv[1,0]= -numpy.sin(ra)
    AInv[1,1]= numpy.cos(ra)
    AInv[1,2]= 0.
    AInv[2,0]= -numpy.cos(ra)*numpy.sin(dec)
    AInv[2,1]= -numpy.sin(ra)*numpy.sin(dec)
    AInv[2,2]= numpy.cos(dec)
    Aphi= numpy.zeros((3,3,len(ra)))
    Aphi[0,0]= numpy.cos(phi1)*numpy.cos(phi2)
    Aphi[0,1]= -numpy.sin(phi1)
    Aphi[0,2]= -numpy.cos(phi1)*numpy.sin(phi2)
    Aphi[1,0]= numpy.sin(phi1)*numpy.cos(phi2)
    Aphi[1,1]= numpy.cos(phi1)
    Aphi[1,2]= -numpy.sin(phi1)*numpy.sin(phi2)
    Aphi[2,0]= numpy.sin(phi2)
    Aphi[2,1]= 0.
    Aphi[2,2]= numpy.cos(phi2)
    TAphi= numpy.dot(_TKOP.T,numpy.swapaxes(Aphi,0,1))
    #Got lazy...
    trans= numpy.zeros((2,2,len(ra)))
    for ii in range(len(ra)):
        trans[:,:,ii]= numpy.dot(AInv[:,:,ii],TAphi[:,:,ii])[1:,1:]
    pmradec= (trans*numpy.array([[pmphi1,pmphi2],[pmphi1,pmphi2]])).sum(1).T
    pmra= pmradec[:,0]
    pmdec= pmradec[:,1]
    #Now convert to pmll
    return bovy_coords.pmrapmdec_to_pmllpmbb(pmra,pmdec,ra,dec)

def test_pmllpmbb_to_pmphi12():
    #Again test based on comparing Willett et al. to Koposov et al.
    # At (l,b) = (216,32) Willett et al. claim (pmll,pmbb) =~ (8.6,-7.4)
    # From Koposov this should give at approximately (phi12) = (-70,-2.) a
    # proper motion of about (pmphi12) = (-11.5,-3.)
    pmphi1,pmphi2= pmllpmbb_to_pmphi12(8.6,-7.4,216.,32.,degree=True)
    assert numpy.fabs(pmphi1+11.5) < 2., 'pmphi1 at the high-v end of the stream does not agree between Koposov and Willett'
    assert numpy.fabs(pmphi2+3.) < 2., 'pmphi2 at the high-v end of the stream does not agree between Koposov and Willett'
    # At (l,b) = (160,60) Willett et al. claim (pmll,pmbb) =~ (11.5,-0.6)
    # From Koposov this should give at approximately (phi12) = (-28,0) a 
    # proper motion of about (pmphi12) = (-12.5,-3.)
    pmphi1,pmphi2= pmllpmbb_to_pmphi12(11.5,-0.6,160.,60.,degree=True)
    assert numpy.fabs(pmphi1+12.5) < 2., 'pmphi1 at the lower-v end of the stream does not agree between Koposov and Willett'
    assert numpy.fabs(pmphi2+3.) < 2., 'pmphi2 at the lower-v end of the stream does not agree between Koposov and Willett'
    #Also test for arrays
    s= numpy.ones(3)
    pmphi1pmphi2= pmllpmbb_to_pmphi12(11.5*s,-0.6*s,160.*s,60.*s,degree=True)
    pmphi1= pmphi1pmphi2[:,0]
    pmphi2= pmphi1pmphi2[:,1]
    assert numpy.all(numpy.fabs(pmphi1+12.5) < 2.), 'pmphi1 at the lower-v end of the stream does not agree between Koposov and Willett'
    assert numpy.all(numpy.fabs(pmphi2+3.) < 2.), 'pmphi2 at the lower-v end of the stream does not agree between Koposov and Willett'  
    return None

def test_pmphi12_to_pmllpmbb():
    # Test that pmphi12_to_pmllpmbb is the inverse of pmllpmbb_to_pmphi12
    # First test
    phi1, phi2= lb_to_phi12(216.,32.,degree=True)
    pmphi1,pmphi2= pmllpmbb_to_pmphi12(8.6,-7.4,216.,32.,degree=True)
    pmll,pmbb= pmphi12_to_pmllpmbb(pmphi1,pmphi2,phi1,phi2,degree=True)
    assert numpy.fabs(pmll-8.6) < 0.001, 'pmphi12_to_pmllpmbb is not the inverse of pmllpmbb_to_pmphi12'
    assert numpy.fabs(pmbb+7.4) < 0.001, 'pmphi12_to_pmllpmbb is not the inverse of pmllpmbb_to_pmphi12'
    # Second test
    phi1, phi2= lb_to_phi12(160.,60.,degree=True)
    pmphi1,pmphi2= pmllpmbb_to_pmphi12(11.5,-0.6,160.,60.,degree=True)
    pmll,pmbb= pmphi12_to_pmllpmbb(pmphi1,pmphi2,phi1,phi2,degree=True)
    assert numpy.fabs(pmll-11.5) < 0.001, 'pmphi12_to_pmllpmbb is not the inverse of pmllpmbb_to_pmphi12'
    assert numpy.fabs(pmbb+0.6) < 0.001, 'pmphi12_to_pmllpmbb is not the inverse of pmllpmbb_to_pmphi12'
    #Also test for arrays
    s= numpy.ones(3)
    pmllpmbb= pmphi12_to_pmllpmbb(pmphi1*s,pmphi2*s,phi1*s,phi2*s,degree=True)
    pmll= pmllpmbb[:,0]
    pmbb= pmllpmbb[:,1]
    assert numpy.all(numpy.fabs(pmll-11.5) < 0.001), 'pmphi12_to_pmllpmbb is not the inverse of pmllpmbb_to_pmphi12'
    assert numpy.all(numpy.fabs(pmbb+0.6) < 0.001), 'pmphi12_to_pmllpmbb is not the inverse of pmllpmbb_to_pmphi12'
    return None

koposov10_dist.txt

# phi1 | dist | dist_err
-55.00|7.20|0.30
-45.00|7.59|0.40
-35.00|7.83|0.30
-25.00|8.69|0.40
-15.00|8.91|0.40
0.00|9.86|0.50

koposov10_pm.txt

# phi1 | pm_phi1 | pm_phi2 | pm_err
-55.00|-13.60 |-5.70|1.30
-45.00|-13.10 |-3.30|0.70
-35.00|-12.20 |-3.10|1.00
-25.00|-12.60 |-2.70|1.40
-15.00|-10.80 |-2.80|1.00

koposov10_pos.txt

# phi1 | phi2 | phi2_err
-60.00|-0.64|0.15
-56.00|-0.89|0.27
-54.00|-0.45|0.15
-48.00|-0.08|0.13
-44.00|0.01|0.14
-40.00|-0.00|0.09
-36.00|0.04|0.10
-34.00|0.06|0.13
-32.00|0.04|0.06
-30.00|0.08|0.10
-28.00|0.03|0.12
-24.00|0.06|0.05
-22.00|0.06|0.13
-18.00|-0.05|0.11
-12.00|-0.29|0.16
-2.00 |-0.87|0.07

koposov10_rv.txt

# phi1 | phi2  | vrad | vrad err
-45.23|-0.04|28.8|6.9
-43.17|-0.09|29.3|10.2
-39.54|-0.07|2.9|8.7
-39.25|-0.22|-5.2|6.5
-37.95|0.00 |1.1|5.6
-37.96|-0.00|-11.7|11.2
-35.49|-0.05|-50.4|5.2
-35.27|-0.02|-30.9|12.8
-34.92|-0.15|-35.3|7.5
-34.74|-0.08|-30.9|9.2
-33.74|-0.18|-74.3|9.8
-32.90|-0.15|-71.5|9.6
-32.25|-0.17|-71.5|9.2
-29.95|0.00|-92.7|8.7
-26.61|-0.11|-114.2|7.3
-25.45|-0.14|-67.8|7.1
-24.86|0.01 |-111.2|17.8
-21.21|-0.02|-144.4|10.5
-14.47|-0.15|-179.0|10.0
-13.73|-0.28|-191.4|7.5
-13.02|-0.21|-162.9|9.6
-12.68|-0.26|-217.2|10.7
-12.55|-0.23|-172.2|6.6
#-56.|-9999.99|39.|14.
#-47.|-9999.99|-7.|10.
#-28.|-9999.99|-61.|6.
#-24.|-9999.99|-83.|9.