#!/usr/bin/env python

# (c) Ken Harris, GPL, 2022

'''
Basic (command line) celestial almanac app :
Given a time and position, calculate the position of 2 bodies and the angle between them (minus the semi-diameter)

Default time is "now" and can be set on the command line w/ the --time flag
Default position is near San Francisco (w/ height of 6 meters) and can be set w/ the --geo flag
Default bodies are Aldebran and the moon, and can be set w/ the --name0 & --name1 flags

Example command line arguments :

python cel-nav-skyfield.py -t 2022-03-10T00:41:38 --geo 23.1553,-109.725,7 -n Sun

NB :
The "height" in the --geo argument is used when calculating "altaz", so includes the "dip"
The --geo info is used for calculations, so includes parallax
'''

import skyfield.api # pip install skyfield
import numpy
import datetime

def semidiameter(app0): # Apparent
    if app0.target == 10: # sun
        radius = 696340.0 # km
    elif app0.target == 301: # moon
        radius = 1737.1 # km
    else:
        radius = 0.0
    return skyfield.units.Angle(radians=numpy.arcsin(radius / app0.distance().km))
    # XXX use math.arcsin ??

############

sfload = skyfield.api.Loader('~/skyfield-data')

eph0 = sfload('de421.bsp') # 16.7 MB download

#print(f'# eph0.names() : {eph0.names()}')
#print(f'# eph0 : {eph0}')

############

'''
import skyfield.data.stellarium
url2 = ('https://raw.githubusercontent.com/Stellarium/stellarium/master'
       '/skycultures/western_SnT/star_names.fab')

with sfload.open(url2) as f2:
    star_names = skyfield.data.stellarium.parse_star_names(f2)

print(f'# star_names : {star_names}')
'''

# late bind ? hippa = None ; ...
import skyfield.data.hipparcos
# hip_main.dat # 53 MB download
with sfload.open(skyfield.data.hipparcos.URL) as f0:
    hippa = skyfield.data.hipparcos.load_dataframe(f0)

############

def name2body(name2):
    # XXX 'aries' ??
    name2 = name2.lower()
    if name2 in ['sun', 'earth', 'moon']:
        body0 = eph0[name2]
    elif name2 in ['venus', 'mars', 'jupiter', 'saturn']:
        body0 = eph0[name2 + ' barycenter'] # faster
    else:
        if name2 == 'hamal': hipnum0 = 9884
        elif name2 == 'aldebaran': hipnum0 = 21421
        elif name2 == 'betelgeuse': hipnum0 = 27989
        elif name2 == 'sirius': hipnum0 = 32349
        elif name2 == 'pollux': hipnum0 = 37826
        elif name2 == 'regulus': hipnum0 = 49669
        elif name2 == 'spica': hipnum0 = 65474
        elif name2 == 'antares': hipnum0 = 80763
        elif name2 == 'altair': hipnum0 = 97649
        elif name2 == 'fomalhaut': hipnum0 = 113368
        elif name2 == 'markab': hipnum0 = 113963
        # else: print star names
        body0 = skyfield.api.Star.from_dataframe(hippa.loc[hipnum0])

    return body0

def myprint(astrometric, apparent): # print lots of debug info

    #print(f'## astrometric : {astrometric}')
    #print(f'## apparent : {apparent}')

    #print(f'## apparent.center : {apparent.center}')
    #print(f'## apparent.target : {apparent.target}')

    #print(f'## eph0.names()[apparent.center] : {eph0.names()[apparent.center]}') # doesn't work w/ eph0 + latlon
    #print(f'## eph0.names()[apparent.target] : {eph0.names()[apparent.target]}') # doesn't work w/ stars

    #print(f'# astrometric.radec("date") : {astrometric.radec("date")}')

    print(f'# apparent.radec("date") : {apparent.radec("date")}')
    # if altaz-able ...
    print(f'# apparent.altaz("standard") : {apparent.altaz("standard")}')
    print(f'# semidiameter(apparent) : {semidiameter(apparent).arcminutes():.1f}')

def proc0(args1):
    args0 = vars(args1) # usings vars() to create dict

    s0 = args0['geo'].split(',')
    # use "parse" ? : r0 = '{lat:g},{:s}{lng:g}{}' # ignore Altitude
    geo0 = skyfield.api.wgs84.latlon(float(s0[0]), float(s0[1]), float(s0[2])) # ITRS ECEF # XXX height used for "dip" calc (in altaz) ?

    ############ time

    #ts = sfload.timescale(builtin=False) # downloads finals2000A.all : Earth Orientation Parameters data file
    ts = sfload.timescale() # timescale object with built-in UT1-tables

    if args0['time'] is None:
        dt0 = datetime.datetime.utcnow()
    else:
        dt0 = datetime.datetime.strptime(args0['time'], '%Y-%m-%dT%H:%M:%S')

    # import dateparser
    # t = dateparser.parse('12/12/12')

    # XXX ts.utc / ts.ut1 / ts.from_datetime ?? utc_strftime / ut1_strftime ?? TAI / TT / TDB / UT1 ?
    t0 = ts.ut1(dt0.year, dt0.month, dt0.day,
            dt0.hour, dt0.minute, dt0.second + dt0.microsecond / 1e6)

    t0._nutation_angles = skyfield.nutationlib.iau2000b(t0.tt) # XXX ??

    print(f'## t0 : {t0.ut1_strftime()}') # utc_strftime

    ############

    print('================================')

    observer0 = name2body('earth') + geo0 # parallax
    print(f'# observer0 : {observer0}')
    barycentric1 = observer0.at(t0) # Barycentric BCRS position and velocity
    #print(f'## barycentric1 : {barycentric1}')

    print('================================')

    name0 = args0['name0'] ; print(f'# name0 : {name0}')
    body0 = name2body(name0) ; print(f'# body0 : {body0}')

    astrometric0 = barycentric1.observe(body0) # Astrometric ICRS / ICRF position and velocity
    apparent0 = astrometric0.apparent() # deflection & aberration : Apparent GCRS position and velocity
    myprint(astrometric0, apparent0)

    print('================================')

    name1 = args0['name1'] ; print(f'# name1 : {name1}')
    body1 = name2body(name1) ; print(f'# body1 : {body1}')

    astrometric1 = barycentric1.observe(body1) # Astrometric ICRS / ICRF position and velocity
    apparent1 = astrometric1.apparent() # deflection & aberration : Apparent GCRS position and velocity
    myprint(astrometric1, apparent1)

    print('================================')

    #print(f'# astrometric0.separation_from(astrometric1) : {astrometric0.separation_from(astrometric1)}')
    print(f'# apparent0.separation_from(apparent1) : {apparent0.separation_from(apparent1)}')

    hc0 = apparent0.separation_from(apparent1).degrees - semidiameter(apparent0).degrees - semidiameter(apparent1).degrees
    print('# separation_from minus semi-diameters : ', skyfield.units.Angle(degrees=hc0))

def main():
    import argparse
    parser = argparse.ArgumentParser(description="almanac")

    # sun, moon, venus, mars, saturn, jupiter
    # Capella Procyon Sirius Rigel Betelgeuse
    # Elnath Menkar
    # lunar stars : Hamal Aldebaran Pollux Regulus Spica Antares Altair Fomalhaut Markab

    parser.add_argument("-n", "--name0", default = 'Aldebaran', type = str, help="heavenly body")
    parser.add_argument("-N", "--name1", default = 'moon', type = str, help="heavenly body")
    parser.add_argument("-t", "--time", default = None, help="time : %Y-%m-%dT%H:%M:%S")
    parser.add_argument("-g", "--geo", default = '38.0,-122.0,6', help="assumed position : lat,lng,meters") # SF area
    #parser.add_argument("-T", "--temp", default = "standard", help="temperature_C (for refraction)")

    proc0(parser.parse_args())

import sys
if __name__ == '__main__':
    sys.exit(main())