'''
Extract a GDR3 reference catalogue for the field specified.

Created on 2 Jan 2022

@author: nch
'''

from astropy.time import Time
radius = 3.1 * 1.414

#gdr3_catalogue_file_name = "/Users/nch/RECONS/SCR1845/gdr3cat.csv"
gdr3_catalogue_file_name = "/Users/nch/BDs/HIP94931/gdr3cat-E0340-Glt18p0.csv"
# Note: Gaia DR3 source ID is 6439125097427143808 for SCR1845, G = 14.021
gdr3_epoch = Time('2016-01-01T12:00:00.0', format='isot', scale='tcb') # from Lindegren et al. 2021

if __name__ == '__main__':

    # define field centre and radius [degrees]
    #from CatExtract import celestial_tangent_point: E0340 field centre:
    celestial_tangent_point = '19h04m35.561 +42d13m35.41'
    # E0281 field centre (J2000): 19h34m44.080 +42d19m39.30 
    #'19h19m00.5489000285 +41d38m04.582441681' HIP 94931 position
    racen = celestial_tangent_point.split()[0]
    deccen = celestial_tangent_point.split()[1]
    
    # Use ESA's astroquery on the Gaia archive to get a reference star catalogue:
    import warnings
    warnings.filterwarnings("ignore")
    from astroquery.utils.tap.core import TapPlus
    gaia = TapPlus(url="https://gea.esac.esa.int/tap-server/tap")
    gaia.login(credentials_file='/Users/nch/JupyterNoteBooks/Gaia/C9VALBT-10X/my_credentials_file')
    
    # create Astropy coords for the Gaia archive query:
    from astropy.coordinates import SkyCoord
    coords_icrs = SkyCoord(racen, deccen, frame = 'fk5')
    alpha = coords_icrs.ra.degree
    delta = coords_icrs.dec.degree
    
    # query around the given coordinates; filter on RUWE and optionally exclude SCR1845:
    adql_query = "SELECT source_id, ra, dec, pmra, pmdec, phot_g_mean_mag, parallax, bp_rp " + \
                "FROM gaiaedr3.gaia_source " + \
                "WHERE 1=CONTAINS(POINT('ICRS', ra, dec), CIRCLE('ICRS', %f, %f, %f))"%\
                (alpha, delta, radius) + \
                " AND pmra IS NOT NULL AND pmdec IS NOT NULL and parallax IS NOT NULL" + \
                " AND ruwe < 1.4" + \
                " AND source_id NOT IN (6439125097427143808, 2101486923382009472, 2101486923385239808)" + \
                " AND phot_g_mean_mag BETWEEN 0.0 AND 18.0" # just in case this is re-run !!!
    # ... note we don't want to include SCR1845 / HIP 94931 in the local astrometric solution to avoid
    # any possibility of biasing the results. Reference star magnitude range is limited to
    # avoid magnitude terms biasing a fit dominated by more numerous, fainter (than SCR1845)
    # stars
    
    # execute the query
    job = gaia.launch_job_async(adql_query)
    results = job.get_results()
    
    results.write(gdr3_catalogue_file_name, format='csv', overwrite=True)
    
    print("Written GDR3 catalogue with %d records."%(len(results)))