#!/usr/bin/env python3
# Sphinx / Leo precession computation used on truthinstone.com
# Published by Truth in Stone (https://truthinstone.com/claims/sphinx-leo-equinox-10500-bc/).
# License: CC0 1.0 Universal (public domain dedication). https://creativecommons.org/publicdomain/zero/1.0/
# Requires: Python 3, numpy, pyerfa, astropy, matplotlib. Star data (from SIMBAD, queried 2026-10-22) are included below.
# Outputs are written to the current directory.
"""
Sphinx / Leo precession check for truthinstone.com  (2026-10-22)

Method
------
* Precession: IAU SOFA/ERFA long-term precession of Vondrak, Capitaine & Wallace (2011, A&A 534, A22;
  2012 erratum), via pyerfa (erfa.ltp, erfa.ltecm). Valid to a few arcseconds within +/-1,000 yr of J2000
  and to a fraction of a degree (usually much better) over +/-200,000 yr.
* Constellations: official IAU 1930 boundaries (Delporte), as digitised by Roman (1987, PASP 99, 695),
  via astropy.coordinates.get_constellation.
* Stars: SIMBAD J2000 positions, proper motions, parallaxes, radial velocities (queried 2026-10-22),
  propagated rigorously with erfa.pmsafe (space motion).
* Equinox geometry: at the March equinox the Sun is at ecliptic longitude 0 of date = RA 0h, Dec 0 of date
  (mean equinox; nutation ~17" and aberration ~20" are ignored). The date in the (Julian/proleptic) calendar is
  irrelevant for this geometry, and no Delta-T is needed because nothing depends on clock time.
* Giza: lat 29.9792 N. With Sun Dec = 0, the hour angle at true altitude h is H = -arccos(sin h / cos lat).
  The due-east horizon point (az 90, alt 0) has HA = -6h, Dec 0, so its RA of date = LST + 90 deg.
Years: astronomical numbering (year 0 = 1 BC, -10499 = 10,500 BC). Julian epoch ~ year.
"""
import json, math
import numpy as np
import erfa, warnings
warnings.filterwarnings("ignore", category=erfa.ErfaWarning)
from astropy.coordinates import SkyCoord, get_constellation
import astropy.units as u

LAT = math.radians(29.9792)
D2R = math.pi / 180.0
MAS = D2R / 3600e3

def bc(y):  # astronomical year -> label
    return f"{1-y:,} BC" if y <= 0 else f"AD {y}"

def cvec(ra, dec):
    return np.array([math.cos(dec)*math.cos(ra), math.cos(dec)*math.sin(ra), math.sin(dec)])

def vec2radec(v):
    v = v/np.linalg.norm(v)
    return math.atan2(v[1], v[0]) % (2*math.pi), math.asin(v[2])

_cache = {}
def constel(v_icrs):
    ra, dec = vec2radec(v_icrs)
    key = (round(ra, 7), round(dec, 7))
    if key not in _cache:
        _cache[key] = get_constellation(SkyCoord(ra*u.rad, dec*u.rad, frame='icrs'))
    return _cache[key]

def constel_many(vs):
    radec = np.array([vec2radec(v) for v in vs])
    return list(get_constellation(SkyCoord(radec[:,0]*u.rad, radec[:,1]*u.rad, frame='icrs')))

def vernal_point_icrs(epj):
    """x-axis of the ecliptic/equinox of date expressed in ICRS (= equinox direction)."""
    return erfa.ltecm(epj)[0].copy()

def date_eq_to_icrs(epj, ra_d, dec_d):
    """mean equator & equinox of date -> J2000/ICRS vector (ltp excludes 23 mas frame bias: negligible)."""
    return erfa.ltp(epj).T @ cvec(ra_d, dec_d)

# ---- stars (SIMBAD, J2000/ICRS, pmra = mu_alpha*cos(dec) mas/yr) ----
STARS = {
 "Regulus (alf Leo)":   (152.09296243828146, 11.967208776100023, -248.73,   5.59, 41.13,  0.72),
 "Denebola (bet Leo)":  (177.26490975591017, 14.572058064829658, -497.68, -114.67, 90.91, -0.2),
 "Algieba (gam1 Leo)":  (154.9931435916667, 19.841488666666663,  310.8,  -152.9,  0.0,   0.0),
 "Zosma (del Leo)":     (168.52708926588556, 20.523718139047475,  143.42, -129.88, 55.82, -20.9),
 "Chertan (tet Leo)":   (168.56002362722833, 15.429571086424442,  -60.538, -78.675, 20.316, 7.536),
 "Ras Elased (eps Leo)":(146.46280673687792, 23.77425377711222,  -46.154, -10.799, 14.3394, 4.322),
 "Spica (alf Vir)":     (201.2982473615632, -11.161319485111932,  -42.35,  -30.67, 13.06, -3.31),
 "Acubens (alf Cnc)":   (134.62168421127,    11.8576804028,        42.181,  -31.16, 18.3304, -14.79),
}
J2000JD = 2451545.0

def star_icrs(name, epj, use_pm=True):
    ra, dec, pmra, pmdec, plx, rv = STARS[name]
    ra, dec = ra*D2R, dec*D2R
    if not use_pm:
        return cvec(ra, dec)
    jd2 = J2000JD + (epj-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 cvec(float(r2), float(d2))

def star_of_date(name, epj, use_pm=True):
    return vec2radec(erfa.ltp(epj) @ star_icrs(name, epj, use_pm))

def ecl_lon_of_date(v_icrs, epj):
    w = erfa.ltecm(epj) @ v_icrs
    return math.degrees(math.atan2(w[1], w[0])) % 360, math.degrees(math.asin(w[2]/np.linalg.norm(w)))

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 sun_ha(h_deg):
    return -math.acos(math.sin(h_deg*D2R)/math.cos(LAT))   # Sun at Dec 0 (equinox), morning

# ===================== PART 1: vernal-equinox Sun vs constellations =====================
def transitions(y0, y1, step=1):
    ys = np.arange(y0, y1+1, step)
    cs = constel_many([vernal_point_icrs(float(y)) for y in ys])
    out, start = [], ys[0]
    for i in range(1, len(ys)):
        if cs[i] != cs[i-1]:
            out.append((cs[i-1], int(start), int(ys[i-1]))); start = ys[i]
    out.append((cs[-1], int(start), int(ys[-1])))
    return out

# Sidereal zodiac ("Age of Leo"): vernal point in sidereal sign Leo = sidereal longitude 120-150 deg.
# ayanamsa(t) = lambda_of_date(Spica,t) - lambda_sid(Spica). Lahiri: Spica at 180 deg sidereal.
# Fagan-Bradley: Spica at 29 Vir 06' 05" = 179.1014 deg sidereal (standard definition).
def ayanamsa(epj, spica_sid):
    lam, _ = ecl_lon_of_date(star_icrs("Spica (alf Vir)", epj), epj)
    return ((lam - spica_sid + 180) % 360) - 180

def age_bounds(spica_sid, sign_lo, y0=-16000, y1=2500):
    """years when vernal point sidereal longitude crosses sign_lo (deg)."""
    ys = np.arange(y0, y1+1, 1)
    # sidereal longitude of vernal point = -ayanamsa (mod 360); track continuously (unwrap)
    sid = np.array([(-ayanamsa(float(y), spica_sid)) % 360 for y in ys])
    sid_unw = np.degrees(np.unwrap(np.radians(sid)))
    # vernal point moves backward in sidereal longitude as time increases
    hits = []
    for target in [sign_lo, sign_lo+30]:
        for k in range(-2, 3):
            t = target + 360*k
            idx = np.where(np.diff(np.sign(sid_unw - t)) != 0)[0]
            for i in idx: hits.append((target, int(ys[i+1])))
    return hits

# ===================== PART 2: Giza, equinox dawn, due east =====================
def dawn_table(y, sun_alts=(-18, -12, -6, -0.833)):
    epj = float(y)
    rows = []
    for h in sun_alts:
        H = sun_ha(h)
        lst = 0.0 + H                       # LST = RA_sun + HA_sun, RA_sun = 0 at equinox
        east_alt0 = date_eq_to_icrs(epj, (lst + math.pi/2) % (2*math.pi), 0.0)
        # point due east at 5 deg altitude (visible above the haze): convert alt/az -> HA/Dec
        alt, az = 5*D2R, 90*D2R
        dec5 = math.asin(math.sin(LAT)*math.sin(alt) + math.cos(LAT)*math.cos(alt)*math.cos(az))
        ha5 = math.atan2(-math.sin(az)*math.cos(alt),
                         math.cos(LAT)*math.sin(alt) - math.sin(LAT)*math.cos(alt)*math.cos(az))
        east_alt5 = date_eq_to_icrs(epj, (lst - ha5) % (2*math.pi), dec5)
        stars = {}
        for n in ["Regulus (alf Leo)", "Denebola (bet Leo)", "Algieba (gam1 Leo)", "Spica (alf Vir)", "Acubens (alf Cnc)"]:
            ra_d, dec_d = star_of_date(n, epj)
            stars[n] = altaz(lst - ra_d, dec_d)
        rows.append(dict(sun_alt=h, east_alt0=constel(east_alt0), east_alt5=constel(east_alt5),
                         east_lon=ecl_lon_of_date(east_alt0, epj)[0], stars=stars))
    return rows

def regulus_rising_az(y, use_pm=True):
    _, dec = star_of_date("Regulus (alf Leo)", float(y), use_pm)
    return math.degrees(math.acos(math.sin(dec)/math.cos(LAT))), math.degrees(dec)

if __name__ == "__main__":
    res = {}
    print("=== PART 1: constellation containing the March-equinox Sun (IAU boundaries) ===")
    tr = transitions(-16000, 2100)
    for c, a, b in tr:
        print(f"  {c:12s} {bc(a):>12s} -> {bc(b)}")
    res["vernal_point_constellations"] = [(c, a, b) for c, a, b in tr]

    print("\n=== Sidereal 'Age of Leo' (vernal point in sidereal Leo 120-150 deg) ===")
    for lab, s in [("Lahiri (Spica=180.000)", 180.0), ("Fagan-Bradley (Spica=179.1014)", 179.1014)]:
        z = [y for y in range(-500, 1000) if abs(ayanamsa(float(y), s)) < 0.02]
        hits = age_bounds(s, 120)
        print(f"  {lab}: ayanamsa=0 near {bc(int(np.mean(z))) if z else '?'};", sorted(set(hits)))
        res[f"age_leo_{lab}"] = sorted(set(hits))
    # equal-2160-year convention check
    for y in [-10969, -8809, -10499]:
        print(f"  check {bc(y)}: Lahiri ayanamsa {ayanamsa(float(y),180.0):.2f} deg")

    print("\n=== PART 2: Giza, March equinox, what is on the due-east horizon ===")
    tab = {}
    for y in list(range(-12999, -2498, 500)) + [-2569]:
        rows = dawn_table(y)
        tab[y] = rows
        print(f"\n{bc(y)}:  Regulus rising az {regulus_rising_az(y)[0]:.1f} deg (dec {regulus_rising_az(y)[1]:.1f})")
        for r in rows:
            reg = r['stars']['Regulus (alf Leo)']
            print(f"   Sun {r['sun_alt']:6.1f}: E@0deg={r['east_alt0']:9s} E@5deg={r['east_alt5']:9s} "
                  f"Regulus alt {reg[0]:5.1f} az {reg[1]:5.1f} | Denebola alt {r['stars']['Denebola (bet Leo)'][0]:5.1f} "
                  f"| Spica alt {r['stars']['Spica (alf Vir)'][0]:5.1f}")
    res["dawn_table"] = {bc(y): v for y, v in tab.items()}

    # fine scan, 10-yr steps, for year ranges when Leo occupies the due-east horizon point
    print("\n=== Fine scan (10-yr steps): years when the due-east horizon point (alt 0) is in Leo ===")
    ranges = {}
    ys = np.arange(-15000, -1999, 10)
    for h in (-18, -12, -6, -0.833):
        H = sun_ha(h)
        vs = [date_eq_to_icrs(float(y), (H + math.pi/2) % (2*math.pi), 0.0) for y in ys]
        cs = constel_many(vs)
        segs, start = [], None
        for y, c in zip(ys, cs):
            if c == 'Leo' and start is None: start = y
            if c != 'Leo' and start is not None: segs.append((int(start), int(y-10))); start = None
        if start is not None: segs.append((int(start), int(ys[-1])))
        ranges[h] = segs
        print(f"   Sun at {h:6.1f} deg: Leo due east during " + "; ".join(f"{bc(a)} - {bc(b)}" for a, b in segs))
    res["leo_due_east_ranges"] = {str(k): v for k, v in ranges.items()}

    # Regulus: when was it on the eastern sky at equinox dawn with Sun at -12 (visible, alt>=5)?
    print("\n=== Regulus at equinox dawn (Sun -12 deg), 100-yr steps ===")
    reg_scan = []
    for y in range(-15000, -1999, 100):
        ra_d, dec_d = star_of_date("Regulus (alf Leo)", float(y))
        a, z = altaz(sun_ha(-12) - ra_d, dec_d)
        reg_scan.append((y, a, z))
    vis = [y for y, a, z in reg_scan if a >= 5 and 60 <= z <= 120]
    print("   Regulus >=5 deg alt, az 60-120 at Sun -12:", bc(min(vis)) if vis else None, "->", bc(max(vis)) if vis else None)
    res["regulus_visible_east_sun-12"] = (min(vis), max(vis)) if vis else None
    # proper-motion effect
    for y in (-10499,):
        r1 = star_of_date("Regulus (alf Leo)", float(y), True); r0 = star_of_date("Regulus (alf Leo)", float(y), False)
        sep = math.degrees(math.acos(min(1, np.dot(cvec(*r1), cvec(*r0)))))
        print(f"\nRegulus proper-motion displacement at {bc(y)}: {sep:.2f} deg")
        res["regulus_pm_shift_deg_10500BC"] = sep
    json.dump(res, open("leo-sphinx-precession-results.json", "w"), indent=1, default=str)


# ===================== CHART (run: python this.py chart) =====================
def make_chart(path="leo-equinox-chart.png"):
    import matplotlib; matplotlib.use("Agg"); import matplotlib.pyplot as plt
    E0 = erfa.ltecm(2000.0)
    def lonJ2000(v):
        w = E0 @ v; return math.degrees(math.atan2(w[1], w[0])) % 360
    ys = np.arange(-15000, 2101, 20)
    vp = [vernal_point_icrs(float(y)) for y in ys]
    lon = np.degrees(np.unwrap(np.radians([lonJ2000(v) for v in vp])))   # continuous track
    cons = constel_many(vp)
    H = sun_ha(-12)
    ep = [date_eq_to_icrs(float(y), (H + math.pi/2) % (2*math.pi), 0.0) for y in ys]
    elon = np.degrees(np.unwrap(np.radians([lonJ2000(v) for v in ep])))
    econs = constel_many(ep)
    colors = {"Libra":"#8c564b","Virgo":"#2ca02c","Leo":"#d62728","Cancer":"#9467bd","Gemini":"#1f77b4",
              "Taurus":"#ff7f0e","Aries":"#17becf","Pisces":"#7f7f7f","Scorpius":"#bcbd22","Auriga":"#000000","Perseus":"#e377c2","Aquarius":"#aec7e8","Coma Berenices":"#c49c94","Lynx":"#dbdb8d"}
    fig, ax = plt.subplots(figsize=(11, 6.2), dpi=130)
    yrs_bc = 1 - ys
    for c in sorted(set(cons)):
        mk = np.array([x == c for x in cons])
        ax.scatter(np.where(mk, yrs_bc, np.nan), lon, s=6, color=colors.get(c, "k"), label=f"Equinox Sun in {c}" if c in set(cons) else None)
    for c in sorted(set(econs)):
        mk = np.array([x == c for x in econs])
        ax.scatter(np.where(mk, yrs_bc, np.nan), elon, s=2, marker="s", color=colors.get(c, "k"), alpha=0.45)
    ax.plot([], [], "s", color="0.5", alpha=0.5, ms=3, label="Due-east horizon, Sun 12° below (dawn)")
    # sidereal Age of Leo (Lahiri, real LTP rate) and equal-2160-yr convention
    ax.axvspan(10702, 8495, color="#d62728", alpha=0.07, label="'Age of Leo', Lahiri sidereal zodiac (10,702–8,495 BC)")
    ax.axvline(10500, color="k", ls="--", lw=1); ax.text(10450, 215, "10,500 BC", fontsize=8)
    ax.axvline(2570, color="k", ls=":", lw=1); ax.text(2520, 215, "Khafre era\n~2,570 BC", fontsize=8)
    ax.invert_xaxis()
    ax.set_xlabel("Year BC (astronomical calculation; negative = AD)")
    ax.set_ylabel("Ecliptic longitude, fixed J2000 frame (deg)")
    ax.set_title("Which constellation held the March-equinox Sun? (Vondrák 2011 precession, IAU boundaries)\n"
                 "Large dots: equinox Sun. Small squares: point due east on the horizon at dawn (Sun -12°), Giza")
    ax.legend(fontsize=7, loc="lower left", ncol=2); ax.grid(alpha=0.3)
    fig.tight_layout(); fig.savefig(path); print("chart ->", path)

import sys as _sys
if __name__ == "__main__" and len(_sys.argv) > 1 and _sys.argv[1] == "chart":
    make_chart()