"""The reference checked against NASA JPL's Horizons system directly.

Horospire ephemeris accuracy study, version 1.0. MIT licence (LICENSE-CODE.txt).

Usage:
    python horizons_check.py DE440S_BSP TT_JDS_JSON OUT_CSV

TT_JDS_JSON is a list of Julian dates in Terrestrial Time. For each of the ten
bodies, one read-only request to https://ssd.jpl.nasa.gov/api/horizons.api,
two seconds apart: QUANTITIES 31 (observer ecliptic longitude of date,
apparent), CENTER 500@399 (the geocentre), APPARENT AIRLESS, TIME_TYPE TT,
TLIST the dates. Bodies 10, 301, 199, 299 and the barycentres 4 to 9, as the
reference uses. The same instants are computed with Skyfield from DE440s, so
the two columns differ only by the two programs, not by any time convention.
Horizons uses DE441; over 1900 to 2050 the two ephemerides agree far below
anything this study reports.
"""
import csv
import json
import sys
import time
import urllib.parse
import urllib.request

from skyfield.api import load, load_file
from skyfield.framelib import ecliptic_frame

BODIES = [("sun", "10", "sun"), ("moon", "301", "moon"), ("mercury", "199", "mercury"), ("venus", "299", "venus"),
          ("mars", "4", "mars barycenter"), ("jupiter", "5", "jupiter barycenter"), ("saturn", "6", "saturn barycenter"),
          ("uranus", "7", "uranus barycenter"), ("neptune", "8", "neptune barycenter"), ("pluto", "9", "pluto barycenter")]
AGENT = "HorospireAudit/1.0 (+https://horospire.com)"


def horizons(command, jds):
    params = {
        "format": "json", "COMMAND": f"'{command}'", "OBJ_DATA": "'NO'", "MAKE_EPHEM": "'YES'",
        "EPHEM_TYPE": "'OBSERVER'", "CENTER": "'500@399'", "TIME_TYPE": "'TT'", "TLIST_TYPE": "'JD'",
        "TLIST": "'" + " ".join(f"{jd:.9f}" for jd in jds) + "'", "QUANTITIES": "'31'",
        "CSV_FORMAT": "'YES'", "EXTRA_PREC": "'YES'", "APPARENT": "'AIRLESS'",
    }
    request = urllib.request.Request("https://ssd.jpl.nasa.gov/api/horizons.api?" + urllib.parse.urlencode(params), headers={"User-Agent": AGENT})
    with urllib.request.urlopen(request, timeout=60) as response:
        text = json.load(response)["result"]
    rows = text.split("$$SOE")[1].split("$$EOE")[0].strip().splitlines()
    return [float([c.strip() for c in row.split(",")][3]) for row in rows], text.split("{source: ")[1].split("}")[0]


def main():
    bsp, jds_path, out_path = sys.argv[1:4]
    jds = sorted(json.load(open(jds_path)))
    eph = load_file(bsp)
    ts = load.timescale(builtin=True)
    t = ts.tt_jd(jds)
    observer = eph["earth"].at(t)
    with open(out_path, "w", newline="") as handle:
        w = csv.writer(handle, lineterminator="\r\n")
        w.writerow(["jd_tt", "body", "horizons_lon_deg", "reference_lon_deg", "horizons_minus_reference_arcsec", "horizons_source"])
        for name, command, target in BODIES:
            _, lon, _ = observer.observe(eph[target]).apparent().frame_latlon(ecliptic_frame)
            got, source = horizons(command, jds)
            for jd, h, r in zip(jds, got, lon.degrees):
                w.writerow([f"{jd:.9f}", name.capitalize(), f"{h:.7f}", f"{float(r) % 360.0:.7f}", f"{((h - float(r) + 180.0) % 360.0 - 180.0) * 3600.0:.4f}", source])
            time.sleep(2.0)


if __name__ == "__main__":
    main()