#!/usr/bin/env python3
# Angkor / Draco precession check used on truthinstone.com
# Published by Truth in Stone (https://truthinstone.com/claims/angkor-draco-10500-bc/).
# License: CC0 1.0 Universal (public domain dedication). https://creativecommons.org/publicdomain/zero/1.0/
# Requires: Python 3, numpy, pyerfa, matplotlib. Star data (from SIMBAD, queried 2026-10-23) are included below.
"""
Method (same as our Sphinx/Leo and Goebekli Tepe Pillar 43 scripts)
* Precession: IAU SOFA/ERFA long-term precession of Vondrak, Capitaine & Wallace (2011, A&A 534, A22),
  via pyerfa (erfa.ltp). Stars propagated with SIMBAD J2000 proper motions, parallaxes and radial
  velocities using erfa.pmsafe.
* Site: Angkor (Bayon), lat 13.441 N. Refraction ignored except a standard -0.833 deg for sunrise.
* March equinox: Sun at RA 0h, Dec 0 of date (mean equinox). Calendar date and Delta-T do not matter
  because only the Sun's position relative to the stars is used.
* "Dawn": Sun altitude -18 (astronomical), -12 (nautical), -6 (civil), -0.833 (sunrise).
* Draco's orientation: we project the stars onto a ground-plan-like chart (azimuth as bearing, zenith
  distance as radius, north up, east right) and fit one epoch's chart to another with a similarity transform
  (shift, scale, rotation; no reflection). The best-fit rotation tells how far the pattern has turned.
Years are astronomical (0 = 1 BC, -10499 = 10,500 BC).
"""
import math, json
import numpy as np
import erfa, warnings
warnings.filterwarnings("ignore", category=erfa.ErfaWarning)

LAT = math.radians(13.441)
D2R = math.pi/180; MAS = D2R/3600e3; J2000JD = 2451545.0
# name: ra, dec (deg, J2000), pmra*cosdec, pmdec (mas/yr), plx (mas), rv (km/s), V
STARS = {
 "Thuban (alf Dra)":   (211.09732332021125, 64.3758696175675, -56.351, 12.329, 12.5156, -6.1, 3.68),
 "Rastaban (bet Dra)": (262.60817373291667, 52.30138870861111, -15.89, 12.28, 8.58, -20.825, 2.81),
 "Eltanin (gam Dra)":  (269.1515411786243, 51.48889561763423, -8.48, -22.79, 21.14, -27.91, 2.23),
 "Altais (del Dra)":   (288.1387557204383, 67.66153865705445, 94.599, 92.374, 33.3499, 24.818, 3.07),
 "eps Dra":            (297.04312691249, 70.26793029824, 79.486, 40.041, 21.3219, 3.12, 3.91),
 "zet Dra":            (257.1966155398917, 65.71464532697362, -16.55, 32.757, 7.8274, -19.2, 3.17),
 "eta Dra":            (245.99785970199997, 61.51420976085001, -16.394, 57.094, 35.7163, -14.734, 2.74),
 "Edasich (iot Dra)":  (231.23239085120377, 58.966065864404726, -8.183, 17.296, 32.5239, -10.912, 3.29),
 "kap Dra":            (188.37059715355997, 69.78823885895, -58.571, 10.11, 7.0093, -12, 3.89),
 "Giausar (lam Dra)":  (172.85092364299916, 69.33107573011195, -41.24, -18.851, 8.6458, 6.6, 3.85),
 "Grumium (ksi Dra)":  (268.38220642831584, 56.87264192915222, 93.698, 78.914, 29.0603, -26.09, 3.75),
 "chi Dra":            (275.2640939542736, 72.73284772361531, 531.21, -349.71, 124.11, 31.9, 3.58),
 "Vega (alf Lyr)":     (279.234734787025, 38.783688956244, 200.94, 286.23, 130.23, -13.5, 0.03),
}
DRACO = [k for k in STARS if "Vega" not in k]

def cvec(ra, dec): return np.array([math.cos(dec)*math.cos(ra), math.cos(dec)*math.sin(ra), math.sin(dec)])
def v2rd(v):
    v = v/np.linalg.norm(v); return math.atan2(v[1], v[0]) % (2*math.pi), math.asin(v[2])

def star_date(name, y):
    ra, dec, pmra, pmdec, plx, rv, _ = STARS[name]; ra, dec = ra*D2R, dec*D2R
    jd2 = J2000JD + (y-2000.0)*365.25
    r2, d2, *_ = erfa.pmsafe(ra, dec, pmra*MAS/math.cos(dec), pmdec*MAS, plx/1000.0, rv, J2000JD, 0.0, jd2, 0.0)
    return v2rd(erfa.ltp(float(y)) @ cvec(float(r2), float(d2)))

def altaz(ha, dec):
    alt = math.asin(math.sin(LAT)*math.sin(dec) + math.cos(LAT)*math.cos(dec)*math.cos(ha))
    az = math.atan2(-math.cos(dec)*math.sin(ha), math.sin(dec)*math.cos(LAT) - math.cos(dec)*math.sin(LAT)*math.cos(ha))
    return math.degrees(alt), math.degrees(az) % 360

def lst_at_sun(h):  # morning, Sun RA 0, Dec 0
    return -math.acos(math.sin(h*D2R)/math.cos(LAT)) % (2*math.pi)

def sky(y, lst):
    out = {}
    for n in STARS:
        ra, dec = star_date(n, y); ha = (lst - ra) % (2*math.pi)
        a, z = altaz(ha, dec)
        out[n] = dict(dec=math.degrees(dec), ha_h=(math.degrees(ha)/15) % 24, alt=a, az=z)
    return out

res = {}
YEARS = {"10,500 BC": -10499, "10,450 BC": -10449, "AD 1150 (Angkor Wat)": 1150, "AD 1200 (Bayon)": 1200}
lowcrit = 90 - math.degrees(LAT)   # dec needed to stay above horizon at lower culmination
print(f"Angkor lat 13.441 N: a star stays above the horizon at lower culmination only if Dec > {lowcrit:.2f} deg")
for lab, y in YEARS.items():
    print(f"\n=== {lab} (year {y}) ===")
    res[lab] = {}
    for h in (-18, -12, -6, -0.833):
        lst = lst_at_sun(h); s = sky(y, lst)
        res[lab][str(h)] = s
        up = [n for n in DRACO if s[n]["alt"] > 0]
        alts = [s[n]["alt"] for n in DRACO]
        print(f" Sun {h:>6}: LST {math.degrees(lst)/15:5.2f} h | Draco stars above horizon {len(up)}/{len(DRACO)} | alt range {min(alts):.1f}..{max(alts):.1f}")
    s = res[lab]["-12"]
    for n in STARS:
        print(f"   {n:20s} dec {s[n]['dec']:6.2f}  HA {s[n]['ha_h']:5.2f} h  alt {s[n]['alt']:6.1f}  az {s[n]['az']:6.1f}")

# Thuban->Eltanin orientation at nautical dawn as a function of epoch: does the head sit below or above the tail?
print("\nHead (Eltanin) minus tail-end (Thuban) altitude at equinox nautical dawn, and Draco 'lower culmination' test")
rows = []
for y in range(-13000, 2001, 250):
    s = sky(y, lst_at_sun(-12))
    # mean hour angle of the Draco stars (circular mean)
    has = [s[n]["ha_h"]*15*D2R for n in DRACO]
    mha = (math.degrees(math.atan2(np.mean(np.sin(has)), np.mean(np.cos(has))))/15) % 24
    nup = sum(1 for n in DRACO if s[n]["alt"] > 0)
    rows.append((y, s["Eltanin (gam Dra)"]["alt"] - s["Thuban (alf Dra)"]["alt"], mha, nup,
                 s["Eltanin (gam Dra)"]["alt"], s["Thuban (alf Dra)"]["alt"]))
for r in rows[::4]: print(f"  {r[0]:6d}  head-tail dAlt {r[1]:6.1f}  mean HA {r[2]:5.2f} h  stars up {r[3]:2d}  Eltanin alt {r[4]:5.1f}  Thuban alt {r[5]:5.1f}")
# (mean hour angle of a circumpolar figure is only indicative)
# epoch where Draco's mean hour angle at equinox nautical dawn is exactly 12h (lower culmination)
fine = []
for y in range(-14000, 2001, 10):
    s = sky(y, lst_at_sun(-12)); has = [s[n]["ha_h"]*15*D2R for n in DRACO]
    mha = (math.degrees(math.atan2(np.mean(np.sin(has)), np.mean(np.cos(has))))/15) % 24
    fine.append((y, mha))
best = min(fine, key=lambda t: abs(t[1]-12))
best0 = min(fine, key=lambda t: min(abs(t[1]), abs(t[1]-24)))
print(f"\nClosest the mean hour angle of Draco's {len(DRACO)} stars gets to lower culmination (12h) at equinox nautical dawn, 14,000 BC-AD 2000: {best[1]:.2f} h in year {best[0]} -> Draco is never at lower culmination at that moment")
print(f"Draco at UPPER culmination (HA 0h) at equinox nautical dawn: year {best0[0]} (HA {best0[1]:.2f})")
rate = (fine[-1][1]-fine[-101][1]) / (fine[-1][0]-fine[-101][0]) * 1000
print(f"Rate of change of Draco's hour angle at fixed equinox dawn near AD 1000: {rate*15:.1f} deg per 1000 yr")
# Declination range over a full cycle
print("\nDec of date range, 25,800-yr cycle (-24000..+1800):")
for n in ("Thuban (alf Dra)", "Eltanin (gam Dra)", "Rastaban (bet Dra)", "Vega (alf Lyr)"):
    ds = [math.degrees(star_date(n, y)[1]) for y in range(-24000, 1801, 100)]
    print(f"  {n:20s} {min(ds):6.2f} .. {max(ds):6.2f}")
# ---- shape comparison (rotation of the whole figure between epochs) ----
def layout(y, h=-12):
    s = sky(y, lst_at_sun(h)); P = []
    for n in DRACO:
        r = 90 - s[n]["alt"]; a = math.radians(s[n]["az"]); P.append((r*math.sin(a), r*math.cos(a)))
    return np.array(P)
def fit(A, B):
    a = A - A.mean(0); b = B - B.mean(0); za = a[:,0]+1j*a[:,1]; zb = b[:,0]+1j*b[:,1]
    k = (np.conj(zb)@za)/(np.conj(zb)@zb); res = za - k*zb
    return math.sqrt(np.mean(abs(res)**2))/math.sqrt(np.mean(abs(za)**2)), math.degrees(np.angle(k))
A = layout(-10499)
print("\nEquinox nautical-dawn chart of Draco compared with 10,500 BC (rotation needed, rms misfit after fit):")
for y in (-8000, -5000, -2500, 0, 800, 1150, 1200, 2000):
    r, ang = fit(A, layout(y)); print(f"  {y:6d}: rotation {ang:7.1f} deg, misfit {r:.3f}")
iE, iT = DRACO.index("Eltanin (gam Dra)"), DRACO.index("Thuban (alf Dra)")
v1, v2 = layout(-10999), layout(-9999)
a1 = math.degrees(math.atan2(*(v1[iT]-v1[iE]))); a2 = math.degrees(math.atan2(*(v2[iT]-v2[iE])))
print(f"Near 10,500 BC the head-to-tail line turns {abs(a2-a1):.1f} deg per 1000 years at the equinox dawn")
ys = list(range(-14000, -7000, 10)); ds = [star_date("Eltanin (gam Dra)", y)[1] for y in ys]
i = int(np.argmax(ds)); print(f"Eltanin (Draco's head) closest to the celestial pole: dec {math.degrees(ds[i]):.2f} in year {ys[i]} ({1-ys[i]:,} BC)")
json.dump({"years": YEARS, "result": {k: v for k, v in res.items()}, "scan": rows, "lower_culm_year": best, "upper_culm_year": best0},
          open("angkor-draco-results.json", "w"), indent=1, default=float)