"""Errors, statistics and boundary counts for the study's three engines. Horospire ephemeris accuracy study, version 1.0. MIT licence (LICENSE-CODE.txt). Usage: python summarise.py WORK_DIR OUT_DIR MEASURED_ON Reads WORK_DIR/moments.json and reference.json (reference.py), house.json (house_engine.php), swiss.json (swiss_ephemeris.py) and astronomy.json (astronomy_engine.mjs). Writes to OUT_DIR: summary.csv Horospire's engine, per body and for all bodies daily-errors.csv.gz Horospire's engine minus the reference, every placement libraries.csv the three engines, per body, the same way study.json the figures the study page prints Error: the shortest angular difference in ecliptic longitude, engine minus reference, in arcseconds. A placement is in a different sign, decan or whole degree when the engine and the reference put it on different sides of a 30, 10 or 1 degree boundary. """ import csv import gzip import io import json import math import os import sys BODIES = ["sun", "moon", "mercury", "venus", "mars", "jupiter", "saturn", "uranus", "neptune", "pluto"] GRIDS = {"sign": 30.0, "decan": 10.0, "degree": 1.0} VERSION = "1.0" REFERENCE = "JPL DE440s through Skyfield 1.55" ENGINES = { "house": ("Horospire EphemerisEngine", "dfb56ab71", "Sun: Meeus ch. 25; Moon: Meeus ch. 47, 13 longitude terms; Mercury to Saturn: Meeus Table 31.A mean elements; Uranus to Pluto: JPL approximate Keplerian elements (Standish); Kepler solve, IAU 1976 precession to the equinox of date, Jupiter-Saturn terms after Schlyter; civil time used as dynamical time, no Delta T"), "astronomy": ("astronomy-engine", "2.1.19", "SunPosition; Ecliptic(GeoVector(body, time, aberration true)); its own Delta T"), "swiss": ("Swiss Ephemeris", "2.10.03 (pyswisseph 2.10.3.2)", "Moshier mode: calc_ut with FLG_MOSEPH | FLG_SPEED, no ephemeris files; default apparent geocentric, true ecliptic and equinox of date; its own Delta T"), } def sep(a, b): return (a - b + 180.0) % 360.0 - 180.0 def pct(sorted_vals, q): k = (len(sorted_vals) - 1) * q f, c = math.floor(k), math.ceil(k) return sorted_vals[f] + (sorted_vals[c] - sorted_vals[f]) * (k - f) def stats(moments, ref, lon): per = {} for bi, body in enumerate(BODIES): signed, edge = [], [] counts = {g: 0 for g in GRIDS} worst, worst_at = -1.0, None for mi, civil in enumerate(moments): t = ref[mi][bi] e = lon[mi][bi] % 360.0 d = sep(e, t) * 3600.0 signed.append(d) if abs(d) > worst: worst, worst_at = abs(d), civil for g, size in GRIDS.items(): if math.floor(e / size) != math.floor(t / size): counts[g] += 1 if g == "sign": edge.append(min(t % 30.0, 30.0 - t % 30.0) * 3600.0) a = sorted(abs(x) for x in signed) per[body] = { "mean_arcsec": round(sum(a) / len(a), 2), "median_arcsec": round(pct(a, 0.5), 2), "p95_arcsec": round(pct(a, 0.95), 2), "max_arcsec": round(worst, 2), "max_at": worst_at, "signed_mean_arcsec": round(sum(signed) / len(signed), 2), "wrong_sign": counts["sign"], "wrong_decan": counts["decan"], "wrong_degree": counts["degree"], "sign_edge_max_arcsec": round(max(edge), 2) if edge else None, } return per def main(): work, out, measured_on = sys.argv[1:4] os.makedirs(out, exist_ok=True) moments = json.load(open(os.path.join(work, "moments.json"))) reference = json.load(open(os.path.join(work, "reference.json"))) ref = reference["lon"] scales = reference["time_scale"] lon = {k: json.load(open(os.path.join(work, f"{k}.json")))["lon"] for k in ENGINES} n = len(moments) placements = n * len(BODIES) results = {k: stats(moments, ref, lon[k]) for k in ENGINES} # summary.csv: Horospire's engine. head = ["version", "measured_on", "engine", "engine_version", "reference", "sample", "moments", "placements", "body", "mean_arcsec", "median_arcsec", "p95_arcsec", "max_arcsec", "max_at", "signed_mean_arcsec", "wrong_sign", "wrong_decan", "wrong_degree"] with open(os.path.join(out, "summary.csv"), "w", newline="") as h: w = csv.writer(h, lineterminator="\r\n") w.writerow(head) base = [VERSION, measured_on, ENGINES["house"][0], ENGINES["house"][1], REFERENCE, "one moment a day, 1900-01-01 to 2050-12-31", n] for body in BODIES: r = results["house"][body] w.writerow(base + [n, body.capitalize()] + [f"{r[k]:.2f}" for k in ("mean_arcsec", "median_arcsec", "p95_arcsec", "max_arcsec")] + [r["max_at"], f"{r['signed_mean_arcsec']:.2f}", r["wrong_sign"], r["wrong_decan"], r["wrong_degree"]]) w.writerow(base + [placements, "All bodies", "", "", "", "", "", ""] + [sum(results["house"][b][g] for b in BODIES) for g in ("wrong_sign", "wrong_decan", "wrong_degree")]) # daily-errors.csv.gz: every placement, deterministic bytes (gzip mtime 0). text = io.StringIO() w = csv.writer(text, lineterminator="\r\n") w.writerow(["civil_time", "time_scale"] + [f"{b}_error_arcsec" for b in BODIES]) for mi, civil in enumerate(moments): w.writerow([civil, scales[mi]] + [f"{sep(lon['house'][mi][bi] % 360.0, ref[mi][bi]) * 3600.0:.3f}" for bi in range(len(BODIES))]) with open(os.path.join(out, "daily-errors.csv.gz"), "wb") as raw: with gzip.GzipFile(filename="daily-errors.csv", mode="wb", fileobj=raw, mtime=0, compresslevel=9) as gz: gz.write(text.getvalue().encode("utf-8")) # libraries.csv: the three engines. with open(os.path.join(out, "libraries.csv"), "w", newline="") as h: w = csv.writer(h, lineterminator="\r\n") w.writerow(["version", "measured_on", "engine", "engine_version", "how_called", "body", "mean_arcsec", "median_arcsec", "p95_arcsec", "max_arcsec", "wrong_sign", "wrong_degree"]) for k, (name, ver, how) in ENGINES.items(): for body in BODIES: r = results[k][body] w.writerow([VERSION, measured_on, name, ver, how, body.capitalize()] + [f"{r[x]:.2f}" for x in ("mean_arcsec", "median_arcsec", "p95_arcsec", "max_arcsec")] + [r["wrong_sign"], r["wrong_degree"]]) totals = {k: {g: sum(results[k][b][f"wrong_{g}"] for b in BODIES) for g in GRIDS} for k in ENGINES} json.dump({ "version": VERSION, "measured_on": measured_on, "reference": REFERENCE, "moments": n, "placements": placements, "from": moments[0][:10], "to": moments[-1][:10], "engines": {k: {"name": v[0], "version": v[1], "how": v[2]} for k, v in ENGINES.items()}, "bodies": BODIES, "results": results, "totals": totals, }, open(os.path.join(out, "study.json"), "w"), indent=1) for k in ENGINES: print(k, "wrong sign", totals[k]["sign"], "decan", totals[k]["decan"], "degree", totals[k]["degree"], "of", placements) if __name__ == "__main__": main()