#!/usr/bin/python3
#
# Christian Dersch <christian.dersch@physik.uni-marburg.de>
# License: MIT
#

from pyvo.dal import TAPService

tap_service = TAPService("http://dc.g-vo.org/tap")

# Get star position and proper motion from catalog using TAP
query = "SELECT * FROM gaia.dr3lite WHERE source_id=4152407474508172160"
star = tap_service.search(query)
star_ra = star["ra"][0]
star_dec = star["dec"][0]
star_pmra = star["pmra"][0]
star_pmdec = star["pmdec"][0]

# Get close NGC objects from OpenNGC
query = """
    SELECT
    TOP 10
    name,DISTANCE(%f, %f,raj2000, dej2000) as d,raj2000,dej2000,pmra,pmdec,min_ax_deg,maj_ax_deg,obj_type, ic_cross, other_ngc
    FROM openngc.data
    WHERE 
    1=CONTAINS(POINT('ICRS', raj2000, dej2000),CIRCLE('ICRS', %f, %f, 1 ))
    AND name LIKE 'NGC____'
    ORDER BY d
    """ % (star_ra,star_dec,star_ra,star_dec)

ngc = tap_service.search(query)
# Assume that closest object is the right one. This is safe, as proper
# motion calculations below will not match in both coordinates in case of
# mismatch
ngc_ra = ngc["raj2000"][0]
ngc_dec = ngc["dej2000"][0]
ngc_pmra = ngc["pmra"][0]
ngc_pmdec = ngc["pmdec"][0]

# Stupid calculation: Get time by simple linear velocity in ra
t = 3.6e6*(star_ra-ngc_ra)/(star_pmra-ngc_pmra)

# Now check for dec, should be within NGC objects size, ignoring angular orientation...
d0 = star_dec - t*star_pmdec/3.6e6
d1 = ngc_dec - t*ngc_pmdec/3.6e6
diff = abs(d0-d1)
match = (diff < max(ngc["min_ax_deg"][0],ngc["maj_ax_deg"][0]))

# If match: print result
if match:
    print("Object: %s" %(ngc["name"][0]))
    print("Ejection about %d years ago" %(t))
else:
    print("Closest NGC object did not match... Don't be lazy and investigate!")