"""Exploratory checks (not pre-registered) on the planetary configuration dates.

1. Lunar phase at each candidate minute: the Sun-Earth-Moon elongation (0 deg new, 180 full).
   Share of dates within 24.4 deg (about two days) of new or full moon; 27.1% if unrelated.
2. Same-side alignment: mean rank of R1 = |mean exp(iL)| of the four inner planets' longitudes.
3. Earthquake days near syzygy, for comparison: share of M>6.5 days within two days of new or full moon.

Usage: python -I explore.py <horizons dir> <results dir> <usgs csv> [<usgs csv> ...] [--out file.json]
"""
import csv, datetime, json, math, pathlib, sys
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))
import numpy as np
import config_dates as cd

WINDOW = 24.4   # deg; the Moon gains about 12.2 deg a day on the Sun


def elong(hdir, year, minutes):
    e = cd.minute_positions(hdir, "399", year); m = cd.minute_positions(hdir, "301", year)
    rows = np.asarray(minutes) + 1
    a, b = m[rows] - e[rows], -e[rows]
    c = (a * b).sum(1) / np.linalg.norm(a, axis=1) / np.linalg.norm(b, axis=1)
    return np.degrees(np.arccos(np.clip(c, -1, 1)))


def near_syzygy(x):
    return (x < WINDOW) | (x > 180 - WINDOW)


def main():
    args = sys.argv[1:]; out_path = None
    if "--out" in args:
        i = args.index("--out"); out_path = args[i + 1]; del args[i:i + 2]
    hdir, rdir = args[0], pathlib.Path(args[1]); csvs = args[2:]
    out = {"window_deg": WINDOW, "share_if_unrelated": round(2 * WINDOW / 180, 3), "blocks": []}
    quake_days = set()
    for p in csvs:
        for r in csv.DictReader(open(p, encoding="utf-8")):
            if r["type"] == "earthquake" and r["mag"] and float(r["mag"]) > 6.5:
                quake_days.add(r["time"][:10])
    expected = 2 * WINDOW / 180
    for tag in ("2010-2012", "2013-2025", "2000-2009"):
        sel = json.load(open(rdir / f"selection-{tag}.json", encoding="utf-8"))
        dg = json.load(open(rdir / f"diagnose-{tag}.json", encoding="utf-8"))
        hits = n = 0; qh = qn = 0
        for yr in sel:
            y = yr["year"]
            first = {}
            for s in yr["stamps"]:
                first.setdefault(s[:10], s)
            mins = [int((datetime.datetime.strptime(s, "%Y-%m-%d %H:%M") - datetime.datetime(y, 1, 1)).total_seconds() // 60)
                    for s in first.values()]
            if mins:
                ns = near_syzygy(elong(hdir, y, mins)); hits += int(ns.sum()); n += len(ns)
            qd = sorted(d for d in quake_days if d.startswith(str(y)))
            qm = [((datetime.date.fromisoformat(d) - datetime.date(y, 1, 1)).days) * 1440 + 720 for d in qd]
            if qm:
                nq = near_syzygy(elong(hdir, y, qm)); qh += int(nq.sum()); qn += len(nq)
        r1 = {}
        for o in dg:
            r1.setdefault(o["stamp"][:10], o["R1_pct_of_year"])
        v = list(r1.values())
        print(f"{tag}: dates near new/full moon {hits}/{n} = {hits/n:.3f} (unrelated: {expected:.3f}); "
              f"quake days near new/full moon {qh}/{qn} = {qh/qn:.3f}; "
              f"mean rank of R1 (same side of the Sun) {sum(v)/len(v):.3f} (se {math.sqrt(1/12/len(v)):.3f})")
        out["blocks"].append({"block": tag, "dates_near_syzygy": hits, "dates": n, "share_dates": round(hits / n, 3),
                              "quake_days_near_syzygy": qh, "quake_days": qn, "share_quake_days": round(qh / qn, 3),
                              "mean_rank_R1": round(sum(v) / len(v), 3), "se_R1": round(math.sqrt(1 / 12 / len(v)), 3)})
    if out_path:
        with open(out_path, "w", encoding="utf-8", newline="\n") as f:
            json.dump(out, f, indent=1); f.write("\n")


if __name__ == "__main__":
    main()