instruments / chronology / cyclic/align_exact.py · raw
"""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()