Script to crossmatch literature binary IDs with Simbad and Gaia

get_binary_positions.py

"""
Script to crossmatch literature binary IDs with Simbad and Gaia
"""

import os, sys

import matplotlib.pyplot as plt
import numpy as np
import astropy.io.ascii as at
from astropy.coordinates import SkyCoord
from astroquery.simbad import Simbad
from astropy.table import join, Table
from astropy import units as u

def match_to_simbad(bintable):
    """
    match an input table consisting of disparate identifiers to Simbad
    """

    nbin = len(bintable)
    print(nbin,"stars")
    print(bintable.dtype)
    simbad_base = "cl* NGC 2632 "
    bintable["TIC"] = np.zeros(len(bintable),"U20")
    bintable["GaiaDR2"] = np.zeros(len(bintable),"U50")
    bintable["GaiaDR3"] = np.zeros(len(bintable),"U50")
    bintable["SimbadName"] = np.zeros(len(bintable),"U50")
    bintable["Found"] = np.zeros(len(bintable),"int")

    for i in range(nbin):
        if bintable["KW"].mask[i]==False:
            name = bintable["KW"][i]
            simbad_name = f"{simbad_base} KW {name}"
        elif bintable["HSHJ"].mask[i]==False:
            name = bintable["HSHJ"][i]
            simbad_name = f"{simbad_base} HSHJ {name}"
        elif bintable["JC"].mask[i]==False:
            name = bintable["JC"][i]
            simbad_name = f"{simbad_base} JC {name}"
        elif bintable["JS"].mask[i]==False:
            name = bintable["JS"][i]
            simbad_name = f"{simbad_base} JS {name}"
        elif bintable["VL"].mask[i]==False:
            name = bintable["VL"][i]
            simbad_name = f"{simbad_base} VL {name}"
        # elif bintable["2MASS"].mask[i]==False:
        #     name = bintable["2MASS"][i]
        #     simbad_name = f"2MASS {name}"
        elif bintable["EPIC"].mask[i]==False:
            name = bintable["EPIC"][i]
            simbad_name = f"EPIC {name}"
        elif ((bintable["Other name"].mask[i]==False) and
              ("I-" in bintable["Other name"][i])):
            name = bintable["Other name"][i]
            simbad_name = f"[A66] {name}".replace("-"," ")
        elif ((bintable["Other name"].mask[i]==False) and
              ("Gem" in bintable["Other name"][i])):
            simbad_name = bintable["Other name"][i]
        elif ((bintable["Other name"].mask[i]==False) and
              (bintable["Other name"][i]=="A 1565")):
            simbad_name = "HD 74720"
        elif ((bintable["Other name"].mask[i]==False) and
              (bintable["Other name"][i]=="AD 1336")):
            bintable["2MASS"][i] = "J08301213+2313370"
            simbad_name = None
        elif bintable["Art"].mask[i]==False:
            name = bintable["Art"][i]
            simbad_name = f"{simbad_base} ANM {name}"
        elif ((bintable["WEBDA"].mask[i]==False) and
              (bintable["WEBDA"][i]<=577)):
              name = bintable["WEBDA"][i]
              simbad_name = f"{simbad_base} KW {name}"
        else:
            simbad_name = None
            # print("uh oh\n",bintable["Art","WEBDA","Other name",
            #                          "RA(deg)","Dec(deg)","RA(hms)","Dec(dms)"
            #                          ][i])

        found_position = False
        simbad_found = False
        if simbad_name is not None:
            result_table = Simbad.query_objectids(simbad_name)
            # print(result_table)

            bintable["SimbadName"][i] = simbad_name

            if result_table is None:
                print(simbad_name,"Not Found")
                simbad_found = False
                continue

            for id in result_table["ID"]:
                if "TIC" in id:
                    bintable["TIC"][i] = id[4:]
                if "EPIC" in id:
                    bintable["EPIC"][i] = id[4:]
                    bintable["EPIC"].mask[i] = False
                elif "DR2" in id:
                    bintable["GaiaDR2"][i] = id[9:]
                elif "DR3" in id:
                    bintable["GaiaDR3"][i] = id[9:]

            simbad_found = True
            bintable["Found"][i] = 1


        if bintable["RA(hms)"].mask[i]==False:
            # Search by position
            if bintable["Epoch"][i]==1950:
                coord0 = SkyCoord(bintable["RA(hms)"][i],bintable["Dec(dms)"][i],
                                 frame="fk4",#equinox=1950,
                                 unit=(u.hourangle,u.degree))
                coord = coord0.transform_to("icrs")
            else:
                coord = SkyCoord(bintable["RA(hms)"][i],bintable["Dec(dms)"][i],
                                 frame="icrs",
                                 unit=(u.hourangle,u.degree))
            # print(coord)
            found_position = True
        elif bintable["2MASS"].mask[i]==False:
            tm = bintable["2MASS"][i]
            ra_str = tm[1:3]+":"+tm[3:5]+":"+tm[5:7]+"."+tm[7:9]
            de_str = tm[9:12]+":"+tm[12:14]+":"+tm[14:16]+"."+tm[16:]
            coord = SkyCoord(ra_str,de_str,frame="icrs",
                             unit=(u.hourangle,u.degree))
            # print(coord)
            found_position = True
        # else:
        #     print("Unable to get a position!")
        #     print(bintable["Art","WEBDA","Other name","Source",
        #                  "RA(deg)","Dec(deg)","RA(hms)","Dec(dms)"
        #                  ][i])

        if (found_position==True) and (simbad_found==False):
            # print("\n",coord)
            position_results = Simbad.query_region(coord, radius=15 * u.arcsec)
            if position_results is None:
                print("Not Found!")
                print(bintable["Art","WEBDA","Other name",
                             "RA(deg)","Dec(deg)","RA(hms)","Dec(dms)"
                             ][i])
                continue
            # print(position_results.dtype)

            if len(position_results)>1:
                print("multiple results found! Assuming the first is the closest match:")
                print(position_results["MAIN_ID"])
            simbad_name = position_results["MAIN_ID"][0]
            # print(simbad_name)
            bintable["SimbadName"][i] = simbad_name
            bintable["Found"][i] = 1

            result_table = Simbad.query_objectids(simbad_name)

            for id in result_table["ID"]:
                if "TIC" in id:
                    bintable["TIC"][i] = id[4:]
                if "EPIC" in id:
                    bintable["EPIC"][i] = id[4:]
                    bintable["EPIC"].mask[i] = False
                elif "DR2" in id:
                    bintable["GaiaDR2"][i] = id[9:]
                elif "DR3" in id:
                    bintable["GaiaDR3"][i] = id[9:]
        elif (found_position==False) and (simbad_found==False):
            print("Unable to identify!",simbad_name)
            print(bintable["Art","WEBDA","Other name","Source",
                         "RA(deg)","Dec(deg)","RA(hms)","Dec(dms)"
                         ][i])
        # else:
        #     print("IDK what happened here")

    at.write(bintable,"praesepe_consolidated_simbad.csv",delimiter=",",
             overwrite=True)



if __name__=="__main__":

    bintable = at.read("praesepe_consolidated.csv",delimiter=",")

    match_to_simbad(bintable)