#!/usr/bin/env python3 """ LifeAstroAPI astronomy accuracy benchmark. Compares LifeAstroAPI with independent public references: planets Geocentric apparent ecliptic longitude of Sun..Pluto at 24 moments (1950-2050) against NASA JPL Horizons (DE441). riseset Sunrise, sunset, moonrise and moonset for 6 cities on 4 dates in 2027 against the moment JPL Horizons puts the upper limb of the disc on the refracted horizon (the standard civil definition). eclipses Every solar and lunar eclipse 2021-2030: time of greatest eclipse and type against NASA's Five Millennium Canon (Espenak & Meeus). Python 3.8+, standard library only. Credits used: planets 480, riseset 480, eclipses 80. The free tier allows 500 a day, so run one part per day with --only. export LIFEASTRO_API_KEY=your_key python3 accuracy_benchmark.py # all parts python3 accuracy_benchmark.py --only planets Writes planets.csv, riseset.csv, eclipses.csv and summary.json to --out. """ import argparse import csv import datetime as dt import html import json import math import os import re import statistics import sys import time import urllib.error import urllib.parse import urllib.request HORIZONS = "https://ssd.jpl.nasa.gov/api/horizons.api" NASA = { "solar": "https://eclipse.gsfc.nasa.gov/SEcat5/SE2001-2100.html", "lunar": "https://eclipse.gsfc.nasa.gov/LEcat5/LE2001-2100.html", } # Horizons ids. Jupiter..Pluto are system barycentres, which is what the # Swiss Ephemeris returns by default; the offset to the body centre is # below 0.02 arcsec at these distances. JPL_ID = {"Sun": "10", "Moon": "301", "Mercury": "199", "Venus": "299", "Mars": "499", "Jupiter": "5", "Saturn": "6", "Uranus": "7", "Neptune": "8", "Pluto": "9"} VEDIC = ["Sun", "Moon", "Mercury", "Venus", "Mars", "Jupiter", "Saturn"] CITIES = [ # name, lat, lon, fixed UTC offset ("Delhi", 28.6139, 77.2090, "+05:30"), ("Chennai", 13.0827, 80.2707, "+05:30"), ("Singapore", 1.3521, 103.8198, "+08:00"), ("London", 51.5074, -0.1278, "+00:00"), ("Toronto", 43.6532, -79.3832, "-05:00"), ("Sydney", -33.8688, 151.2093, "+10:00"), ] DATES = ["2027-03-20", "2027-06-21", "2027-09-23", "2027-12-21"] EPOCH = dt.datetime(1970, 1, 1, tzinfo=dt.timezone.utc) def jd_of(t): return 2440587.5 + (t - EPOCH).total_seconds() / 86400.0 def moments(): """24 moments from 1950 to 2050, spread over months and hours.""" out = [] for i, y in enumerate(1950 + round(k * 100 / 23) for k in range(24)): out.append(dt.datetime(y, (i * 5) % 12 + 1, (i * 7) % 28 + 1, (i * 11) % 24, (i * 17) % 60, tzinfo=dt.timezone.utc)) return out class Api: def __init__(self, base, key, delay): self.base, self.key, self.delay = base.rstrip("/"), key, delay def get(self, path, **params): url = f"{self.base}{path}?{urllib.parse.urlencode(params)}" req = urllib.request.Request(url, headers={"Authorization": f"Bearer {self.key}"}) for _ in range(6): try: with urllib.request.urlopen(req, timeout=60) as r: time.sleep(self.delay) return json.loads(r.read()) except urllib.error.HTTPError as e: if e.code == 429: time.sleep(int(e.headers.get("Retry-After", "20"))) continue raise SystemExit(f"{path}: HTTP {e.code} {e.read()[:200]!r}") raise SystemExit(f"{path}: still rate limited") def horizons(params): q = {"format": "json", "OBJ_DATA": "'NO'", "MAKE_EPHEM": "'YES'", "EPHEM_TYPE": "'OBSERVER'", "ANG_FORMAT": "'DEG'", "EXTRA_PREC": "'YES'", "CSV_FORMAT": "'YES'"} q.update({k: f"'{v}'" for k, v in params.items()}) for attempt in range(4): try: with urllib.request.urlopen(f"{HORIZONS}?{urllib.parse.urlencode(q)}", timeout=120) as r: res = json.loads(r.read())["result"] break except (urllib.error.URLError, TimeoutError): time.sleep(5 * (attempt + 1)) else: raise SystemExit("JPL Horizons did not answer") body = res[res.index("$$SOE") + 5:res.index("$$EOE")] time.sleep(0.5) return [[c.strip() for c in line.split(",")] for line in body.strip().splitlines()] def wrap_arcsec(a, b): d = (a - b + 180.0) % 360.0 - 180.0 return d * 3600.0 def stats(values): v = [abs(x) for x in values] if not v: return {"n": 0} v.sort() return {"n": len(v), "median": round(statistics.median(v), 3), "p95": round(v[min(len(v) - 1, int(0.95 * len(v)))], 3), "max": round(v[-1], 3)} # ── planets ────────────────────────────────────────────────────────────────── def run_planets(api, out): ms = moments() ours = {} # (i, body) -> (tropical lon, source) for i, t in enumerate(ms): q = dict(date=t.strftime("%Y-%m-%d"), time=t.strftime("%H:%M"), tz="UTC", lat=0, lon=0) v = api.get("/v1/chart/planets", **q) for p in v["planets"]: if p["name"] in VEDIC: ours[(i, p["name"])] = ((p["sidereal_lon"] + v["ayanamsa_value"]) % 360, "chart/planets (sidereal + ayanamsa)") w = api.get("/v1/western/natal/planets", **q) for p in w["planets"]: if p["name"] in ("Uranus", "Neptune", "Pluto"): ours[(i, p["name"])] = (p["longitude"], "western/natal/planets") rows = [] tlist = " ".join(f"{jd_of(t):.6f}" for t in ms) for body, cid in JPL_ID.items(): ref = horizons({"COMMAND": cid, "CENTER": "500@399", "TLIST": tlist, "TLIST_TYPE": "JD", "TIME_TYPE": "UT", "QUANTITIES": "31"}) for i, r in enumerate(ref): lon = float(r[3]) mine, src = ours[(i, body)] rows.append({"moment_utc": ms[i].strftime("%Y-%m-%d %H:%M"), "body": body, "endpoint": src, "lifeastroapi_deg": f"{mine:.7f}", "jpl_horizons_deg": f"{lon:.7f}", "diff_arcsec": f"{wrap_arcsec(mine, lon):.3f}"}) write_csv(out, "planets.csv", rows) by_body = {b: stats([float(r["diff_arcsec"]) for r in rows if r["body"] == b]) for b in JPL_ID} return {"all": stats([float(r["diff_arcsec"]) for r in rows]), "by_body": by_body, "unit": "arcsec"} # ── rise / set ─────────────────────────────────────────────────────────────── def limb_crossing(body, lat, lon, jd0, rising): """JD when the upper limb meets the refracted horizon, from JPL elevations.""" rows = horizons({"COMMAND": JPL_ID[body], "CENTER": "coord@399", "COORD_TYPE": "GEODETIC", "SITE_COORD": f"{lon},{lat},0", "START_TIME": f"JD{jd0 - 0.005:.6f}", "STOP_TIME": f"JD{jd0 + 0.005:.6f}", "STEP_SIZE": "48", "APPARENT": "REFRACTED", "QUANTITIES": "4,13", "CAL_FORMAT": "JD"}) pts = [(float(r[0]), float(r[4]) + float(r[5]) / 7200.0) for r in rows] # elevation + semi-diameter for (t1, f1), (t2, f2) in zip(pts, pts[1:]): if (rising and f1 < 0 <= f2) or (not rising and f1 >= 0 > f2): return t1 + (t2 - t1) * f1 / (f1 - f2) return None def run_riseset(api, out): rows = [] for city, lat, lon, tz in CITIES: for d in DATES: sun = api.get("/v1/panchang/sunrise-sunset", date=d, tz=tz, lat=lat, lon=lon) moon = api.get("/v1/panchang/moonrise-moonset", date=d, tz=tz, lat=lat, lon=lon) for event, body, jd, rising in (("sunrise", "Sun", sun.get("sunrise_jd"), True), ("sunset", "Sun", sun.get("sunset_jd"), False), ("moonrise", "Moon", moon.get("moonrise_jd"), True), ("moonset", "Moon", moon.get("moonset_jd"), False)): if jd is None: continue ref = limb_crossing(body, lat, lon, jd, rising) rows.append({"city": city, "date": d, "event": event, "lifeastroapi_jd": f"{jd:.6f}", "jpl_limb_jd": f"{ref:.6f}" if ref else "", "diff_seconds": f"{(jd - ref) * 86400:.1f}" if ref else "not within ±7 min"}) write_csv(out, "riseset.csv", rows) res = {} for e in ("sunrise", "sunset", "moonrise", "moonset"): sel = [r for r in rows if r["event"] == e] res[e] = stats([float(r["diff_seconds"]) for r in sel if r["jpl_limb_jd"]]) res[e]["outside_window"] = sum(1 for r in sel if not r["jpl_limb_jd"]) res["unit"] = "seconds" return res # ── eclipses ───────────────────────────────────────────────────────────────── MON = {m: i for i, m in enumerate(["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"], 1)} TYPE = {"solar": {"T": "total", "A": "annular", "H": "hybrid", "P": "partial"}, "lunar": {"T": "total", "P": "partial", "N": "penumbral"}} def nasa_catalog(kind): with urllib.request.urlopen(NASA[kind], timeout=120) as r: text = html.unescape(re.sub(r"<[^>]+>", "", r.read().decode("latin-1"))) out = [] for line in text.splitlines(): m = re.match(r"\s*\d+\s+(\d{4}) (\w{3}) (\d\d)\s+(\d\d):(\d\d):(\d\d)\s+(-?\d+)\s+-?\d+\s+\d+\s+(\w+)", line) if not m: continue y, mon, d, hh, mm, ss, dT, typ = m.groups() if not 2021 <= int(y) <= 2030: continue td = dt.datetime(int(y), MON[mon], int(d), int(hh), int(mm), int(ss), tzinfo=dt.timezone.utc) out.append({"td": td, "ut": td - dt.timedelta(seconds=int(dT)), "delta_t": int(dT), "type": TYPE[kind][typ[0]]}) return out def run_eclipses(api, out): rows = [] for kind in ("solar", "lunar"): ours = [] for a, b in (("2021-01-01", "2025-12-31"), ("2026-01-01", "2030-12-31")): ours += api.get(f"/v1/eclipses/{kind}", start_date=a, end_date=b)["eclipses"] for ref in nasa_catalog(kind): jd_ref = jd_of(ref["ut"]) match = min(ours, key=lambda e: abs(e["peak_jd"] - jd_ref), default=None) ok = match is not None and abs(match["peak_jd"] - jd_ref) < 1 rows.append({"kind": kind, "nasa_greatest_ut": ref["ut"].strftime("%Y-%m-%d %H:%M:%S"), "nasa_delta_t_s": ref["delta_t"], "nasa_type": ref["type"], "lifeastroapi_peak_utc": match["peak_utc"] if ok else "", "lifeastroapi_type": match["type"] if ok else "", "diff_seconds": f"{(match['peak_jd'] - jd_ref) * 86400:.1f}" if ok else "missing", "type_matches": str(ok and match["type"] == ref["type"]).lower()}) write_csv(out, "eclipses.csv", rows) found = [r for r in rows if r["diff_seconds"] != "missing"] return {"nasa_eclipses": len(rows), "found": len(found), "type_matches": sum(1 for r in rows if r["type_matches"] == "true"), "time": stats([float(r["diff_seconds"]) for r in found]), "unit": "seconds"} def write_csv(out, name, rows): with open(os.path.join(out, name), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=list(rows[0].keys())) w.writeheader() w.writerows(rows) def main(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--api", default=os.environ.get("LIFEASTRO_API_URL", "https://api.lifeastroapi.com")) ap.add_argument("--key", default=os.environ.get("LIFEASTRO_API_KEY")) ap.add_argument("--only", choices=["planets", "riseset", "eclipses"]) ap.add_argument("--out", default="benchmark-results") ap.add_argument("--delay", type=float, default=0.7, help="seconds between API calls") a = ap.parse_args() if not a.key: sys.exit("Set LIFEASTRO_API_KEY or pass --key (free key: https://lifeastroapi.com/signup/)") os.makedirs(a.out, exist_ok=True) api = Api(a.api, a.key, a.delay) path = os.path.join(a.out, "summary.json") summary = {} if a.only and os.path.exists(path): # parts run on different days add up with open(path) as f: summary = json.load(f) summary.update({"generated_at": dt.datetime.now(dt.timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ"), "api": a.api, "api_version": api.get("/v1/meta/version")}) parts = {"planets": run_planets, "riseset": run_riseset, "eclipses": run_eclipses} for name, fn in parts.items(): if a.only in (None, name): print(f"running {name} …", file=sys.stderr) summary[name] = fn(api, a.out) with open(path, "w") as f: json.dump(summary, f, indent=2) print(json.dumps(summary, indent=2)) if __name__ == "__main__": main()