Godstrument/workers/world_sats.py

#!/usr/bin/env python3
"""world_sats.py — the sky's own schedule, as control signals (SGP4, no polling APIs).

Real satellites, propagated locally. TLEs (two-line element sets) are fetched once
from Celestrak and cached; from them SGP4 gives *any* satellite's position at *any*
rate with no per-fix HTTP — the opposite of world_iss.py, which begs a web API for
every point. We track the space stations plus a few named recon birds (GAOFEN,
COSMOS), turn each into a lat/lon/altitude, and derive its **elevation above the
venue city's horizon** — so "something is passing overhead" becomes a continuous,
playable signal, and a recon satellite cresting the sky rings a bell.

Emits (raw floats; the hub normalizes):
  /gs/sats/best_elev        highest elevation among tracked sats, degrees (can be < 0)
  /gs/sats/overhead_count   how many are above 10 degrees right now
  /gs/sats/iss_lat          ISS geodetic latitude  (the sweep, even without world_iss)
  /gs/sats/iss_lon          ISS geodetic longitude
  /gs/sats/pass.event       impulse (magnitude = elevation) when any sat crosses
                            rising through 30 degrees — a pass begins
  /gs/sats/recon_elev       best elevation among the GAOFEN/COSMOS set only

Depends on `sgp4` (pure-python, `pip install sgp4`). Import-safe: if it is missing,
or there is no cached TLE and no network, the worker prints why and exits cleanly.
"""
from __future__ import annotations
import argparse
import datetime
import math
import os
import sys
import time
import urllib.error
import urllib.request

from pythonosc.udp_client import SimpleUDPClient

try:
    from sgp4.api import Satrec, jday
    HAVE_SGP4 = True
except Exception:                       # pure-python but still optional (house rule: import-safe)
    HAVE_SGP4 = False

NAME = "world_sats"
UA = "Godstrument/1.0 (live instrument; contact monsterrobotparty@gmail.com)"

CACHE_DIR = os.path.expanduser("~/.cache/godverse")
CACHE_FILE = os.path.join(CACHE_DIR, "tles.txt")
CACHE_MAX_AGE = 24 * 3600.0             # Celestrak politeness — TLEs age slowly
REFRESH_EVERY = 6 * 3600.0             # re-check the cache age this often while running

# (url, is_recon) — the stations give us the ISS; the named queries are the watchers.
# A NAME query with no match returns the plain text "No GP data found"; parse_tles
# rejects any group whose first data line doesn't start with "1".
TLE_SOURCES = [
    ("https://celestrak.org/NORAD/elements/gp.php?GROUP=stations&FORMAT=tle", False),
    ("https://celestrak.org/NORAD/elements/gp.php?NAME=GAOFEN&FORMAT=tle", True),
    ("https://celestrak.org/NORAD/elements/gp.php?NAME=COSMOS%202486&FORMAT=tle", True),
    ("https://celestrak.org/NORAD/elements/gp.php?NAME=COSMOS%202506&FORMAT=tle", True),
]

PASS_ELEV = 30.0        # a sat rising through this elevation rings the bell
OVERHEAD_ELEV = 10.0    # counted as "overhead" above this
ISS_CATNR = "25544"     # NORAD catalog number of the ISS (ZARYA)
WGS_A = 6378.137        # WGS-84 semi-major axis, km
WGS_B = 6356.7523142    # semi-minor axis, km


def fetch_url(url, timeout=20):
    req = urllib.request.Request(url, headers={"User-Agent": UA})
    with urllib.request.urlopen(req, timeout=timeout) as resp:
        return resp.read().decode("utf-8", "replace")


def parse_tles(text):
    """Parse TLE text into [(name, line1, line2)]. Tolerant of 2-line groups and of
    Celestrak's plain-text 'No GP data found' (any group whose line1 isn't a '1')."""
    lines = [ln.rstrip() for ln in text.splitlines()
             if ln.strip() and not ln.startswith("#")]
    out = []
    i = 0
    while i < len(lines):
        ln = lines[i]
        if ln.startswith("1 ") and i + 1 < len(lines) and lines[i + 1].startswith("2 "):
            out.append(("", ln, lines[i + 1]))                    # 2-line group, no name
            i += 2
        elif (not ln.startswith(("1 ", "2 "))) and i + 2 < len(lines) \
                and lines[i + 1].startswith("1 ") and lines[i + 2].startswith("2 "):
            out.append((ln.strip(), lines[i + 1], lines[i + 2]))  # named 3-line group
            i += 3
        else:
            i += 1                                                # junk / 'No GP data found' — skip
    return out


def _write_cache(groups, recon_catnrs):
    """Cache the merged 3LE with a recon-catnr header so the tags survive a reload."""
    try:
        os.makedirs(CACHE_DIR, exist_ok=True)
        tmp = CACHE_FILE + ".tmp"
        with open(tmp, "w", encoding="utf-8") as f:
            f.write("# recon " + ",".join(sorted(recon_catnrs)) + "\n")
            for name, l1, l2, _is_recon in groups:
                f.write((name + "\n" if name else "") + l1 + "\n" + l2 + "\n")
        os.replace(tmp, CACHE_FILE)                               # atomic
        print(f"[{NAME}] refreshed TLE cache -> {CACHE_FILE}")
    except OSError as e:
        print(f"[{NAME}] warning: could not write cache: {e}")


def _read_cache():
    """Return (groups, recon_catnrs) from the tagged cache, or ([], set())."""
    if not os.path.isfile(CACHE_FILE):
        return [], set()
    try:
        with open(CACHE_FILE, "r", encoding="utf-8") as f:
            text = f.read()
    except OSError:
        return [], set()
    recon = set()
    for ln in text.splitlines():
        if ln.startswith("# recon "):
            recon = {x for x in ln[len("# recon "):].strip().split(",") if x}
            break
    groups = [(name, l1, l2, l1[2:7].strip() in recon) for name, l1, l2 in parse_tles(text)]
    return groups, recon


def load_grouped_tles(force_fetch=False):
    """Return [(name, l1, l2, is_recon)] — a fresh cache if possible, else a fetch.

    Fetches each source separately so recon birds keep their tag, writes the tagged
    cache, and falls back to a stale cache on network failure. Returns [] only when
    there is no cache and no network."""
    fresh = False
    if not force_fetch:
        try:
            fresh = os.path.isfile(CACHE_FILE) and (time.time() - os.path.getmtime(CACHE_FILE)) < CACHE_MAX_AGE
        except OSError:
            fresh = False
    if fresh:
        groups, _ = _read_cache()
        if groups:
            print(f"[{NAME}] using cached TLEs (< 24 h old)")
            return groups
    merged, recon_catnrs = [], set()
    for url, is_recon in TLE_SOURCES:
        try:
            txt = fetch_url(url)
        except (urllib.error.URLError, OSError, ValueError) as e:
            print(f"[{NAME}] TLE fetch failed for {url.split('?')[-1]}: {e}")
            continue
        for name, l1, l2 in parse_tles(txt):
            catnr = l1[2:7].strip()
            if is_recon:
                recon_catnrs.add(catnr)
            merged.append((name, l1, l2, is_recon))
    if merged:
        _write_cache(merged, recon_catnrs)
        return merged
    groups, _ = _read_cache()                                    # network down — serve stale
    if groups:
        print(f"[{NAME}] network down — using stale TLE cache")
    return groups


def gstime(jdut1):
    """Greenwich Mean Sidereal Time (radians) — the IAU-82 series (matches satellite.js)."""
    t = (jdut1 - 2451545.0) / 36525.0
    g = 67310.54841 + (876600.0 * 3600 + 8640184.812866) * t + 0.093104 * t * t - 6.2e-6 * t * t * t
    g = math.radians(g / 240.0) % (2 * math.pi)                  # seconds -> degrees (/240) -> radians
    return g + 2 * math.pi if g < 0 else g


def eci_to_geodetic(r, gmst):
    """TEME position (km) + GMST -> (lat_deg, lon_deg, alt_km). WGS-84, iterative."""
    x, y, z = r
    a, b = WGS_A, WGS_B
    f = (a - b) / a
    e2 = 2 * f - f * f
    R = math.hypot(x, y)
    lon = math.atan2(y, x) - gmst
    lon = (lon + math.pi) % (2 * math.pi) - math.pi              # wrap to [-pi, pi]
    lat = math.atan2(z, R)
    C = 1.0
    for _ in range(20):                                          # converges in a handful of steps
        C = 1.0 / math.sqrt(1 - e2 * math.sin(lat) ** 2)
        lat = math.atan2(z + a * C * e2 * math.sin(lat), R)
    alt = R / math.cos(lat) - a * C
    return math.degrees(lat), math.degrees(lon), alt


def _ecef(lat_deg, lon_deg, alt_km):
    a, b = WGS_A, WGS_B
    la, lo = math.radians(lat_deg), math.radians(lon_deg)
    f = (a - b) / a
    e2 = 2 * f - f * f
    N = a / math.sqrt(1 - e2 * math.sin(la) ** 2)
    return ((N + alt_km) * math.cos(la) * math.cos(lo),
            (N + alt_km) * math.cos(la) * math.sin(lo),
            (N * (1 - e2) + alt_km) * math.sin(la))


def elevation(sat_lat, sat_lon, sat_alt, obs_lat, obs_lon):
    """Elevation angle (degrees) of a sat above an observer's local horizon."""
    sx, sy, sz = _ecef(sat_lat, sat_lon, sat_alt)
    ox, oy, oz = _ecef(obs_lat, obs_lon, 0.0)
    rx, ry, rz = sx - ox, sy - oy, sz - oz
    rng = math.sqrt(rx * rx + ry * ry + rz * rz)
    if rng <= 0:
        return 0.0
    la, lo = math.radians(obs_lat), math.radians(obs_lon)
    ux, uy, uz = math.cos(la) * math.cos(lo), math.cos(la) * math.sin(lo), math.sin(la)
    dot = (rx * ux + ry * uy + rz * uz) / rng
    return math.degrees(math.asin(max(-1.0, min(1.0, dot))))


class Sat:
    __slots__ = ("name", "catnr", "rec", "is_recon", "prev_elev")

    def __init__(self, name, line1, line2, is_recon):
        self.catnr = line1[2:7].strip()
        self.name = name or f"SAT {self.catnr}"
        self.rec = Satrec.twoline2rv(line1, line2)               # handles no_kozai internally
        self.is_recon = is_recon
        self.prev_elev = None

    def geodetic_now(self, dt):
        jd, fr = jday(dt.year, dt.month, dt.day, dt.hour, dt.minute,
                      dt.second + dt.microsecond * 1e-6)
        e, r, _v = self.rec.sgp4(jd, fr)
        if e != 0:                                               # decayed / propagation error
            return None
        return eci_to_geodetic(r, gstime(jd + fr))


def build_sats(groups):
    sats = []
    for name, l1, l2, is_recon in groups:
        try:
            sats.append(Sat(name, l1, l2, is_recon))
        except Exception as e:                                   # a malformed TLE shouldn't kill the set
            print(f"[{NAME}] skip TLE '{name or l1[2:7]}': {e}")
    return sats


def main():
    ap = argparse.ArgumentParser(description="SGP4 satellites -> OSC (the sky's schedule)")
    ap.add_argument("--host", default="127.0.0.1")
    ap.add_argument("--port", type=int, default=9000)
    ap.add_argument("--interval", type=float, default=2.0)
    ap.add_argument("--lat", type=float, default=0.0, help="venue latitude (run.py retargets this)")
    ap.add_argument("--lon", type=float, default=0.0, help="venue longitude")
    args = ap.parse_args()

    if not HAVE_SGP4:
        print(f"[{NAME}] sgp4 not installed — `pip install sgp4` to bring the sky in. Exiting.")
        sys.exit(0)

    sats = build_sats(load_grouped_tles())
    if not sats:
        print(f"[{NAME}] no TLEs (no cache and no network) — exiting politely.")
        sys.exit(0)

    client = SimpleUDPClient(args.host, args.port)
    n_recon = sum(1 for s in sats if s.is_recon)
    have_recon = n_recon > 0
    print(f"[{NAME}] tracking {len(sats)} satellites ({n_recon} recon) from "
          f"({args.lat:.2f}, {args.lon:.2f}); emitting /gs/sats/* every {args.interval:g}s")

    last_refresh = time.time()
    while True:
        if time.time() - last_refresh > REFRESH_EVERY:           # age-check + refresh TLEs in place
            last_refresh = time.time()
            fresh = build_sats(load_grouped_tles())              # refetches only if the cache is stale
            if fresh:
                sats = fresh
                have_recon = any(s.is_recon for s in sats)

        now = datetime.datetime.now(datetime.timezone.utc)
        best_elev, recon_best, overhead = -90.0, -90.0, 0
        for s in sats:
            geo = s.geodetic_now(now)
            if geo is None:
                continue
            lat, lon, alt = geo
            el = elevation(lat, lon, alt, args.lat, args.lon)
            if el > best_elev:
                best_elev = el
            if s.is_recon and el > recon_best:
                recon_best = el
            if el > OVERHEAD_ELEV:
                overhead += 1
            if s.prev_elev is not None and s.prev_elev < PASS_ELEV <= el:   # rising through 30deg
                client.send_message("/gs/sats/pass.event", float(el))       # -> sats.pass.event (impulse)
                print(f"[{NAME}] pass: {s.name} rising through {PASS_ELEV:.0f}deg (now {el:.1f})")
            s.prev_elev = el
            if s.catnr == ISS_CATNR:
                client.send_message("/gs/sats/iss_lat", float(lat))
                client.send_message("/gs/sats/iss_lon", float(lon))

        client.send_message("/gs/sats/best_elev", float(best_elev))
        client.send_message("/gs/sats/overhead_count", float(overhead))
        if have_recon:
            client.send_message("/gs/sats/recon_elev", float(recon_best))

        time.sleep(args.interval)


if __name__ == "__main__":
    try:
        main()
    except KeyboardInterrupt:
        sys.exit(0)