#!/usr/bin/env python3
"""Reproduce asterwise.com/accuracy: Asterwise vs NASA JPL Horizons vs Swiss Ephemeris.

    pip install pyswisseph requests
    ASTERWISE_API_KEY=aw_... python3 compare.py

Without ASTERWISE_API_KEY the script still compares Swiss Ephemeris with Horizons.
Point SE_EPHE_PATH at a directory holding the Swiss Ephemeris data files
(sepl_18.se1, semo_18.se1, sepl_24.se1 from https://github.com/aloistr/swisseph/tree/master/ephe);
pyswisseph ships without them and otherwise falls back to its built-in Moshier
ephemeris, which is less precise than the files Asterwise runs on.
Horizons is queried live (https://ssd.jpl.nasa.gov/api/horizons.api), so results
can move by hundredths of an arcsecond if JPL updates an ephemeris.
"""
import os, sys, requests
import swisseph as swe

INSTANTS = [(1950,1,1,0,0),(1969,7,20,20,17),(1985,11,12,1,15),(2000,1,1,12,0),
            (2012,12,21,11,11),(2024,4,8,18,17),(2026,9,5,0,0),(2050,6,30,12,0)]
BODIES = [("Sun","10",swe.SUN),("Moon","301",swe.MOON),("Mercury","199",swe.MERCURY),("Venus","299",swe.VENUS),
          ("Mars","499",swe.MARS),("Jupiter","599",swe.JUPITER),("Saturn","699",swe.SATURN),
          ("Uranus","799",swe.URANUS),("Neptune","899",swe.NEPTUNE),("Pluto","999",swe.PLUTO)]
API = "https://api.asterwise.com/v1/western/natal"
KEY = os.environ.get("ASTERWISE_API_KEY")
if os.environ.get("SE_EPHE_PATH"):
    swe.set_ephe_path(os.environ["SE_EPHE_PATH"])

def horizons(body_id, jds):
    r = requests.get("https://ssd.jpl.nasa.gov/api/horizons.api", params={
        "format":"text","COMMAND":f"'{body_id}'","OBJ_DATA":"'NO'","MAKE_EPHEM":"'YES'","EPHEM_TYPE":"'OBSERVER'",
        "CENTER":"'500@399'","TLIST":"'"+" ".join(f"{j:.6f}" for j in jds)+"'","QUANTITIES":"'31'",
        "CAL_FORMAT":"'JD'","ANG_FORMAT":"'DEG'","APPARENT":"'AIRLESS'"}, timeout=60)
    body = r.text.split("$$SOE")[1].split("$$EOE")[0].strip().splitlines()
    return {round(float(l.split()[0]),6): float(l.split()[1]) for l in body}

def asterwise(y,m,d,hh,mm):
    r = requests.post(API, headers={"X-API-Key": KEY}, json={"name":"check","date":f"{y:04d}-{m:02d}-{d:02d}",
        "time":f"{hh:02d}:{mm:02d}","latitude":0.0,"longitude":0.0,"timezone":"UTC"}, timeout=60)
    r.raise_for_status()
    return {p["name"]: p["longitude"] for p in r.json()["data"]["planets"]}

def arcsec(a, b): return ((a - b + 180) % 360 - 180) * 3600

jds = [swe.julday(y,m,d,hh+mm/60) for y,m,d,hh,mm in INSTANTS]
hz = {name: horizons(hid, jds) for name,hid,_ in BODIES}
aw = {jd: asterwise(*inst) for jd, inst in zip(jds, INSTANTS)} if KEY else {}
worst = 0.0
print(f"{'instant (JD UT)':>16} {'body':8} {'swisseph':>12} {'horizons':>12} {'d(arcsec)':>10}" + ("  asterwise  d(arcsec)" if KEY else ""))
for jd, inst in zip(jds, INSTANTS):
    for name, hid, pid in BODIES:
        s = swe.calc_ut(jd, pid, swe.FLG_SWIEPH | swe.FLG_SPEED)[0][0]
        h = hz[name][round(jd,6)]
        line = f"{jd:16.5f} {name:8} {s:12.6f} {h:12.6f} {arcsec(s,h):+10.3f}"
        if KEY:
            a = aw[jd][name]; line += f"  {a:10.6f} {arcsec(a,h):+10.3f}"; worst = max(worst, abs(arcsec(a,h)))
        else:
            worst = max(worst, abs(arcsec(s,h)))
        print(line)
print(f"largest difference vs Horizons: {worst:.3f} arcsec")