horizons.py

#!/usr/bin/env python3
# SPDX-License-Identifier: AGPL-3.0-or-later
# Copyright (C) 2026 PUDDY Inc. <legal@puddystudios.com>
"""NASA JPL HORIZONS client - the authoritative engine for The Star Map.

Queries apparent ecliptic longitude of date (quantity 31, ObsEcLon) and refines
an event to the second:
  - refine_ingress  : longitude crosses a sign cusp (geocentric)
  - refine_lunation : Moon-Sun elongation crosses 0/90/180/270 (geo OR topocentric)
  - refine_station  : longitude turning point, via parabolic vertex (geocentric)

Topocentric queries use CENTER='coord@399' + SITE_COORD='lon,lat,alt_km'.

Notes learned the hard way:
  - Host is ssd.jpl.nasa.gov (ssd.api.jpl.nasa.gov does not resolve here).
  - Param values with spaces MUST be single-quoted ("Too many constants" else).
  - "SECONDS" is not a valid STEP_SIZE unit. A bare integer = number of equal
    intervals, so an N-second window with STEP_SIZE=N yields 1-second spacing.
"""
import re
import time
import urllib.parse
import urllib.request
from datetime import datetime, timezone, timedelta

HOST = "https://ssd.jpl.nasa.gov/api/horizons.api"
GEO = "500@399"

HB = {
    "Sun": "10", "Moon": "301", "Mercury": "199", "Venus": "299", "Mars": "499",
    "Jupiter": "599", "Saturn": "699", "Uranus": "799", "Neptune": "899",
    "Pluto": "999", "Ceres": "1;", "Pallas": "2;", "Juno": "3;", "Vesta": "4;",
    "Chiron": "2060;",
}
DATE = re.compile(r"(\d{4}-\w{3}-\d{2} \d{2}:\d{2}:\d{2}\.\d+)")
FLOAT = re.compile(r"-?\d+\.\d+")


def _get(params, retries=3):
    url = HOST + "?" + urllib.parse.urlencode(params)
    last = None
    for i in range(retries):
        try:
            with urllib.request.urlopen(url, timeout=90) as r:
                return r.read().decode("utf-8", "replace")
        except Exception as e:  # noqa
            last = e
            time.sleep(1.5 * (i + 1))
    raise last


def fetch_eclon(body, start_utc, stop_utc, nsteps, center=GEO, site_coord=None):
    """Return [(datetime_utc, ecl_lon_deg), ...] for one body.

    site_coord: 'lon,lat,alt_km' string for topocentric (center must be coord@399).
    """
    params = {
        "format": "text",
        "COMMAND": f"'{HB[body]}'",
        "EPHEM_TYPE": "OBSERVER",
        "CENTER": f"'{center}'",
        "START_TIME": f"'{start_utc:%Y-%m-%d %H:%M:%S}'",
        "STOP_TIME": f"'{stop_utc:%Y-%m-%d %H:%M:%S}'",
        "STEP_SIZE": f"'{nsteps}'",
        "QUANTITIES": "'31'",
        "ANG_FORMAT": "'DEG'",
        "CAL_FORMAT": "'CAL'",
        "TIME_DIGITS": "'SECONDS'",
    }
    if site_coord:
        params["SITE_COORD"] = f"'{site_coord}'"
    txt = _get(params)
    if "$$SOE" not in txt:
        raise RuntimeError(f"HORIZONS no ephemeris for {body}:\n{txt[-260:]}")
    block = txt.split("$$SOE")[1].split("$$EOE")[0]
    out = []
    for line in block.splitlines():
        dm = DATE.search(line)
        if not dm:
            continue
        # ObsEcLon is the first float after the timestamp (tolerate flag columns)
        vals = FLOAT.findall(line[dm.end():])
        if not vals:
            continue
        t = datetime.strptime(dm.group(1), "%Y-%b-%d %H:%M:%S.%f").replace(tzinfo=timezone.utc)
        out.append((t, float(vals[0])))
    if not out:
        raise RuntimeError(f"parsed 0 rows for {body}")
    return out


def _signed(a, target):
    return (a - target + 180.0) % 360.0 - 180.0


def _interp_crossing(value_fn, pts):
    for i in range(len(pts) - 1):
        v0 = value_fn(pts[i][1])
        v1 = value_fn(pts[i + 1][1])
        if v0 == 0:
            return pts[i][0]
        if v0 * v1 < 0:
            frac = -v0 / (v1 - v0)
            return pts[i][0] + (pts[i + 1][0] - pts[i][0]) * frac
    raise RuntimeError("no crossing in window")


def refine_ingress(body, target_lon, guess_utc, windows=(45, 600, 7200)):
    """Slow bodies sit seconds-to-minutes off the guess; widen until bracketed."""
    last = None
    for hw in windows:
        start = guess_utc - timedelta(seconds=hw)
        stop = guess_utc + timedelta(seconds=hw)
        step = 2 * hw if hw <= 600 else 3600
        pts = fetch_eclon(body, start, stop, step)
        try:
            return _interp_crossing(lambda l: _signed(l, target_lon), pts)
        except RuntimeError as e:
            last = e
            time.sleep(0.3)
    raise last


def refine_lunation(target_elong, guess_utc, hw=60, center=GEO, site_coord=None):
    start = guess_utc - timedelta(seconds=hw)
    stop = guess_utc + timedelta(seconds=hw)
    sun = fetch_eclon("Sun", start, stop, 2 * hw, center, site_coord)
    time.sleep(0.35)
    moon = fetch_eclon("Moon", start, stop, 2 * hw, center, site_coord)
    n = min(len(sun), len(moon))
    elong = [(sun[i][0], moon[i][1] - sun[i][1]) for i in range(n)]
    return _interp_crossing(lambda e: _signed(e, target_elong), elong)


def refine_station(body, guess_utc, hw_hours=24, npts=240):
    """Station = longitude turning point. Fit a parabola to lon(t) over a window
    and take its vertex; robust to the near-motionless extremum."""
    import numpy as np
    start = guess_utc - timedelta(hours=hw_hours)
    stop = guess_utc + timedelta(hours=hw_hours)
    pts = fetch_eclon(body, start, stop, npts)
    t0 = pts[0][0]
    xs = np.array([(t - t0).total_seconds() for t, _ in pts])
    lons = np.degrees(np.unwrap(np.radians([l for _, l in pts])))
    a, b, c = np.polyfit(xs, lons, 2)
    if a == 0:
        raise RuntimeError("degenerate station fit")
    vertex = -b / (2 * a)
    if vertex < xs[0] or vertex > xs[-1]:
        raise RuntimeError("station vertex outside window")
    return t0 + timedelta(seconds=float(vertex))