python中arcsec_Python astropy.units 模块，degree() 实例源码 - 编程字典

def xmatch(cat1,cat2,maxdist=2,

colRA1='RA',colDec1='DEC',epoch1=2000.,

colRA2='RA',colDec2='DEC',epoch2=2000.,

colpmRA2='pmra',colpmDec2='pmdec',

swap=False):

"""

NAME:

xmatch

PURPOSE:

cross-match two catalogs (incl. proper motion in cat2 if epochs are different)

INPUT:

cat1 - First catalog

cat2 - Second catalog

maxdist= (2) maximum distance in arcsec

colRA1= ('RA') name of the tag in cat1 with the right ascension in degree in cat1 (assumed to be ICRS)

colDec1= ('DEC') name of the tag in cat1 with the declination in degree in cat1 (assumed to be ICRS)

epoch1= (2000.) epoch of the coordinates in cat1

colRA2= ('RA') name of the tag in cat2 with the right ascension in degree in cat2 (assumed to be ICRS)

colDec2= ('DEC') name of the tag in cat2 with the declination in degree in cat2 (assumed to be ICRS)

epoch2= (2000.) epoch of the coordinates in cat2

colpmRA2= ('pmra') name of the tag in cat2 with the proper motion in right ascension in degree in cat2 (assumed to be ICRS; includes cos(Dec)) [only used when epochs are different]

colpmDec2= ('pmdec') name of the tag in cat2 with the proper motion in declination in degree in cat2 (assumed to be ICRS) [only used when epochs are different]

swap= (False) if False, find closest matches in cat2 for each cat1 source, if False do the opposite (important when one of the catalogs has duplicates)

OUTPUT:

(index into cat1 of matching objects,

index into cat2 of matching objects,

angular separation between matching objects)

HISTORY:

2016-09-12 - Written - Bovy (UofT)

2016-09-21 - Account for Gaia epoch 2015 - Bovy (UofT)

"""

if ('ref_epoch' in cat1.dtype.fields and numpy.fabs(epoch1-2015.) > 0.01)\

or ('ref_epoch' in cat2.dtype.fields and \

numpy.fabs(epoch2-2015.) > 0.01):

warnings.warn("You appear to be using a Gaia catalog, but are not setting the epoch to 2015., which may lead to incorrect matches")

depoch= epoch2-epoch1

if depoch != 0.:

# Use proper motion to get both catalogs at the same time

dra=cat2[colpmRA2]/numpy.cos(cat2[colDec2]/180.*numpy.pi)\

/3600000.*depoch

ddec= cat2[colpmDec2]/3600000.*depoch

else:

dra= 0.

ddec= 0.

mc1= acoords.SkyCoord(cat1[colRA1],cat1[colDec1],

unit=(u.degree, u.degree),frame='icrs')

mc2= acoords.SkyCoord(cat2[colRA2]-dra,cat2[colDec2]-ddec,

unit=(u.degree, u.degree),frame='icrs')

if swap:

idx,d2d,d3d = mc2.match_to_catalog_sky(mc1)

m1= numpy.arange(len(cat2))

else:

idx,d2d,d3d = mc1.match_to_catalog_sky(mc2)

m1= numpy.arange(len(cat1))

mindx= d2d < maxdist*u.arcsec

m1= m1[mindx]

m2= idx[mindx]

if swap:

return (m2,m1,d2d[mindx])

else:

return (m1,m2,d2d[mindx])