"""Line them up. Eclipses computed from orbital mechanics (astronomy-engine, no calendar assumed) against
eclipse records as an annal states them: reign-year within a sequence of rulers, month, and the
day's name in the sixty-day cycle. Nothing is taken from a concordance. Two free parameters are
searched exhaustively: which astronomical year the annal's first year is, and the phase of the
sixty-day cycle. For each candidate pair, count records whose stated cycle-day coincides with a
computed solar eclipse in the right year. A real alignment is a lone spike; a false one is flat.

Records file (JSON): {"name":..., "site": [lat, lon], "rulers": [{"name":..., "years": N}, ...],
 "records": [{"ruler":..., "year": y, "month": m, "ganzhi": "癸未", "text": "..."}]}
"""
import argparse, json, math, sys
import numpy as np
from paths import ROOT, topic_dir
STEMS = "甲乙丙丁戊己庚辛壬癸"; BRANCHES = "子丑寅卯辰巳午未申酉戌亥"
def ganzhi_index(gz):
    s, b = STEMS.index(gz[0]), BRANCHES.index(gz[1])
    for i in range(60):
        if i % 10 == s and i % 12 == b: return i
    raise ValueError(gz)

def solar_eclipses(y0, y1, site=None):
    """Solar eclipses between astronomical years y0..y1 from orbital mechanics (astronomy-engine), no
    calendar assumed. With a site, only eclipses visible there (Sun up, obscuration > 0), and the
    time returned is the local peak. Returns (JD_UT of peak, kind, magnitude)."""
    import astronomy as A
    t = A.Time.Make(y0, 1, 1, 0, 0, 0); end = A.Time.Make(y1 + 1, 1, 1, 0, 0, 0); out = []
    if site is None:
        while True:
            e = A.SearchGlobalSolarEclipse(t)
            if e.peak.ut >= end.ut: break
            out.append((e.peak.ut + 2451545.0, str(e.kind).split(".")[-1], round(e.obscuration or 0, 2))); t = e.peak.AddDays(10)
    else:
        obs = A.Observer(site[0], site[1], 50)
        while True:
            e = A.SearchLocalSolarEclipse(t, obs)
            if e.peak.time.ut >= end.ut: break
            if e.peak.altitude > -0.5 and e.obscuration > 0.01:
                out.append((e.peak.time.ut + 2451545.0, str(e.kind).split(".")[-1], round(e.obscuration, 2)))
            t = e.peak.time.AddDays(10)
    return out

def main():
    ap = argparse.ArgumentParser(); ap.add_argument("--records", required=True); ap.add_argument("--from", dest="y0", type=int, default=-1000); ap.add_argument("--to", dest="y1", type=int, default=-200)
    ap.add_argument("--out", default=None); ap.add_argument("--local", action="store_true", help="only eclipses visible at the annal's site"); ap.add_argument("--slack", type=int, default=0, help="days of slack on the cycle day (0 = exact)"); ap.add_argument("--check", type=int, default=None, help="also report this first-year's best cell (e.g. the conventional one)"); a = ap.parse_args()
    R = json.loads(open(a.records).read()); rulers = R["rulers"]; start = {}; acc = 0
    for r in rulers: start[r["name"]] = acc; acc += r["years"]
    recs = [dict(x, offset=start[x["ruler"]] + x["year"] - 1, gz=ganzhi_index(x["ganzhi"])) for x in R["records"] if x.get("ganzhi")]
    span = acc; print(f"{len(recs)} dated eclipse records over {span} reign-years; searching first-year = {a.y0}..{a.y1 - span}", file=sys.stderr)
    site = R.get("site") if a.local else None
    ecl = solar_eclipses(a.y0, a.y1, site); print(f"{len(ecl)} computed solar eclipses ({'visible at ' + str(site) if site else 'anywhere on Earth'}) in {a.y0}..{a.y1}", file=sys.stderr)
    lon = (R.get("site") or [0, 0])[1]
    jd = np.array([e[0] + lon / 360.0 for e in ecl]); jd_int = np.floor(jd + 0.5).astype(np.int64)   # LOCAL civil date at the site
    def jd_year(j):   # proleptic Julian calendar year from a Julian day number (Meeus)
        z = int(math.floor(j + 0.5)); a_ = z; b = a_ + 1524; c = int((b - 122.1) / 365.25); d = int(365.25 * c); e_ = int((b - d) / 30.6001)
        m = e_ - 1 if e_ < 14 else e_ - 13; return c - 4716 if m > 2 else c - 4715
    ecl_year = np.array([jd_year(float(j)) for j in jd])
    def report_cell(Y0, p, label):
        h = 0; lines = []
        for r in recs:
            ry = Y0 + r["offset"]; in_year = (ecl_year == ry) | ((ecl_year == ry - 1) & (ecl_month >= 11))
            cand = np.where(in_year & np.array([month_ok(r["month"], m) for m in ecl_month]))[0]
            gzs = [(int(jd_int[c]), int((jd_int[c] + p) % 60), ecl[c][1], ecl[c][2]) for c in cand]
            ok = any(g == r["gz"] or (g + a.slack) % 60 == r["gz"] for _, g, _, _ in gzs); h += ok
            lines.append(f"  {'ok ' if ok else 'NO '} {r['ruler']} {r['year']:>2} month {str(r['month'] or '?'):>2} {r['ganzhi']} → year {ry}: {[(k, STEMS[g % 10] + BRANCHES[g % 12], kind[:3], mag) for k, g, kind, mag in gzs] or 'no eclipse in that window'}")
        print(f"\n{label}: first year {Y0} ({'BCE ' + str(1 - Y0)}), phase {p}: {h}/{len(recs)} record by record"); print("\n".join(lines))
    def jd_month(j):
        z = int(math.floor(j + 0.5)); b = z + 1524; c = int((b - 122.1) / 365.25); d = int(365.25 * c); e_ = int((b - d) / 30.6001)
        return e_ - 1 if e_ < 14 else e_ - 13
    ecl_month = np.array([jd_month(float(j)) for j in jd])
    # The annal's month M is a Zhou lunisolar month; Zhou month 1 holds the winter solstice, so it
    # sits at Julian December/January: Zhou M ≈ Julian M − 1, give or take one month.
    def month_ok(zhou_m, jul_m):
        if not zhou_m: return True
        target = (zhou_m - 2) % 12 + 1; diff = (jul_m - target) % 12   # Zhou month 1 ≈ December/January
        return diff == 0   # exact Zhou month (the cycle-proof test)
    # joint grid over (first year Y0, cycle phase p): matches[Y0, p]
    Y0s = list(range(a.y0, a.y1 - span + 1)); grid = np.zeros((len(Y0s), 60), dtype=np.int16)
    gz_tab = np.stack([(jd_int + p) % 60 for p in range(60)])          # 60 × n_ecl
    gz_tab_next = np.stack([(jd_int + a.slack + p) % 60 for p in range(60)])   # with --slack 0 this is the same table
    for yi, Y0 in enumerate(Y0s):
        for r in recs:
            ry = Y0 + r["offset"]
            # the Lu year runs from about the preceding December: allow the last two Julian months of ry-1
            in_year = (ecl_year == ry) | ((ecl_year == ry - 1) & (ecl_month >= 11))
            cand = np.where(in_year & np.array([month_ok(r["month"], m) for m in ecl_month]))[0]
            if cand.size == 0: continue
            hit = ((gz_tab[:, cand] == r["gz"]) | (gz_tab_next[:, cand] == r["gz"])).any(axis=1)   # per phase
            grid[yi] += hit
    flat = grid.ravel(); mean, sd = flat.mean(), flat.std(); order = np.argsort(-flat)[:8]
    print(f"\nrecords {len(recs)} · grid {len(Y0s)} first-years × 60 phases = {flat.size} cells · null mean {mean:.2f} · sd {sd:.2f} · max {flat.max()}")
    print(f"cells scoring ≥ 10: {(flat >= 10).sum()} · ≥ 20: {(flat >= 20).sum()}")
    for k in order:
        yi, p = divmod(int(k), 60); Y0 = Y0s[yi]; h = int(flat[k])
        print(f"  first year = {Y0:5d} ({'BCE ' + str(1 - Y0) if Y0 <= 0 else 'CE ' + str(Y0)})  phase {p:2d}  matches {h:2d}/{len(recs)}  z={(h - mean) / (sd or 1):.1f}")
    yi, p = divmod(int(order[0]), 60); report_cell(Y0s[yi], p, "best cell")
    if a.check is not None and a.check in Y0s:
        row = grid[Y0s.index(a.check)]; pc = int(row.argmax()); print(f"\ncheck: first year {a.check} scores {int(row[pc])}/{len(recs)} at phase {pc}; rank among all cells: {int((flat > row[pc]).sum()) + 1}"); report_cell(a.check, pc, "conventional")
    results = [(Y0s[int(k) // 60], int(k) % 60, int(flat[k])) for k in np.argsort(-flat)[:200]]
    if a.out: json.dump({"records": len(recs), "span": span, "search": [a.y0, a.y1], "null_mean": float(mean), "null_sd": float(sd), "top": results}, open(a.out, "w"))
if __name__ == "__main__": main()