#!/usr/bin/env python3
"""AP48 -- "The Assembly" -- independent verification script.

Reproduces, from CODATA 2022 alone, the numbers behind AP48's closing of the
H0 seam AP46 opened -- including the part that FIRES KS-45.1, the corpus's
own registered Hubble-constant prediction (AP18, already public). This
script does not hide that: it recomputes the firing on purpose, every run.

    python3 verify_ap48.py

No dependencies beyond the Python 3 standard library. No network access.

Source: AP48_The_Assembly_FINAL_v1_1.md sections 8-10, and the AP18 erratum
issued beside it. Full text of the state of record: lucid.rodeo/prereg/ap48/
"""
from __future__ import annotations

import math

# --------------------------------------------------------------------------- #
# CODATA 2022 -- the only measured inputs
# --------------------------------------------------------------------------- #
ALPHA = 7.2973525643e-3      # fine-structure constant
HBAR = 1.054571817e-34       # reduced Planck constant, J*s
C = 299792458.0              # speed of light, m/s
M_E = 9.1093837139e-31       # electron mass, kg
MPC_KM = 3.0856775814913673e19   # 1 Mpc in km
GYR_S = 1e9 * 365.25 * 86400.0   # 1 Julian Gyr in seconds
PI = math.pi

# The AP46 cycle's own inputs (already independently verified in verify_ap46.py)
OMEGA_M = 0.3115
OMEGA_L = 1 - OMEGA_M


def cycle_gyr() -> float:
    """(21/18) * alpha^-18 * tau_C, the AP46 cycle, in Gyr."""
    tau_c = HBAR / (M_E * C ** 2)
    cycle_s = (21 / 18) * ALPHA ** -18 * tau_c
    return cycle_s / GYR_S


def h0t0(omega_m: float = OMEGA_M) -> float:
    """Closed-form H0*t0 for flat LambdaCDM, matter + cosmological constant."""
    omega_l = 1 - omega_m
    return (2 / (3 * math.sqrt(omega_l))) * math.asinh(math.sqrt(omega_l / omega_m))


def age_gyr_for_h0(h0_kms_mpc: float, omega_m: float = OMEGA_M) -> float:
    """Age of the universe (Gyr) at a given H0 (km/s/Mpc), flat LambdaCDM."""
    h0_per_s = (h0_kms_mpc / MPC_KM)  # km/s/Mpc -> km/s/km = 1/s (Mpc in km cancels)
    t0_s = h0t0(omega_m) / h0_per_s
    return t0_s / GYR_S


def h0_for_age_gyr(target_gyr: float, omega_m: float = OMEGA_M, tol: float = 1e-10) -> float:
    """Root-find H0 (km/s/Mpc) giving a target age, via bisection -- no scipy needed."""
    lo, hi = 1.0, 200.0
    for _ in range(200):
        mid = (lo + hi) / 2
        if age_gyr_for_h0(mid, omega_m) > target_gyr:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2


def sigma(value: float, measured: float, u: float) -> float:
    return (value - measured) / u


def main() -> int:
    bar = "=" * 76
    print(bar)
    print("AP48 -- The Assembly -- independent verification (CODATA 2022, stdlib only)")
    print(bar)

    cycle = cycle_gyr()
    h0t0_val = h0t0()
    cycle_inverse_kms_mpc = MPC_KM / (cycle * GYR_S)  # 1/cycle, converted to km/s/Mpc

    derived_h0 = h0t0_val * cycle_inverse_kms_mpc

    print("\nThe seam, closed -- AP46's cycle read as an age, inverted to a rate:")
    print(f"  AP46 cycle (21/18)*alpha^-18*tau_C = {cycle:.4f} Gyr  (claimed 13.830)")
    print(f"  H0*t0 (flat LambdaCDM, Om={OMEGA_M})   = {h0t0_val:.4f}  (claimed 0.954)")
    print(f"  1/cycle as a rate                   = {cycle_inverse_kms_mpc:.3f} km/s/Mpc  (claimed 70.70)")
    print(f"  derived H0 = {h0t0_val:.4f} * {cycle_inverse_kms_mpc:.2f}      = {derived_h0:.3f} km/s/Mpc  (claimed 67.45)")

    lo_age, hi_age = cycle - 0.66, cycle + 0.66
    h0_hi = h0_for_age_gyr(lo_age)   # younger age -> faster/higher H0
    h0_lo = h0_for_age_gyr(hi_age)   # older age -> slower/lower H0
    print(f"\nThe window -- one lane-time (0.66 Gyr) either side of the cycle:")
    print(f"  age {lo_age:.2f} Gyr -> H0 = {h0_hi:.2f} km/s/Mpc  (claimed 70.8)")
    print(f"  age {hi_age:.2f} Gyr -> H0 = {h0_lo:.2f} km/s/Mpc  (claimed 64.4)")

    # KS-45.1 -- AP18's registered claim, fired against AP48's derived rate.
    ap18_h0, ap18_u = 74.3, 1.2
    fired_sigma = sigma(ap18_h0, derived_h0, ap18_u)
    print(f"\nKS-45.1 -- AP18's Hubble constant. FIRED 2026-09-03.")
    print(f"  AP18 registered: H0 = {ap18_h0} +/- {ap18_u} km/s/Mpc")
    print(f"  AP48 derived   : H0 = {derived_h0:.2f} km/s/Mpc")
    print(f"  sigma = ({ap18_h0} - {derived_h0:.2f}) / {ap18_u} = {fired_sigma:+.2f}  (claimed 5.7)")
    print(f"  -> {'FIRED' if abs(fired_sigma) > 1 else 'would survive'} on its registered 1-sigma bar.")
    print( "  An erratum widens AP18's bar to +/-14.9 (it omitted a systematic uncertainty) --")
    print( "  that does NOT un-fire this switch. Widening a bar after contradiction and not")
    print( "  counting the fire would itself be the repair this corpus refuses. Shown, kept.")

    # The corrected AP18 floor check at the derived rate.
    C_S2 = 1.0445
    h0_si = derived_h0 * 1000 / (MPC_KM * 1000)  # km/s/Mpc -> 1/s
    a0 = C_S2 * C * h0_si / (2 * PI)
    scale, scale_u = 1.20e-10, 0.24e-10
    floor_pct = (a0 - scale) / scale * 100
    floor_sigma = (scale - a0) / scale_u
    print(f"\nThe floor, re-read (AP18 section 4, at the derived rate):")
    print(f"  a0 = 1.0445 * c * H0 / (2*pi) = {a0:.4e} m/s^2  (claimed 1.089e-10)")
    print(f"  vs empirical scale 1.20e-10 +/- 0.24e-10 (systematic-honest): {floor_pct:+.1f}%, {floor_sigma:.2f} sigma  (claimed 9.2%, 0.46 sigma)")
    print(f"  coefficient a0/(c*H0) = 1.0445/(2*pi) = {C_S2/(2*PI):.4f}  (claimed 0.1662)")

    print(f"\n{bar}")
    print(f"KS-45.1: FIRED at {fired_sigma:+.2f} sigma. Stays fired -- not touched by this script.")
    print(f"KS-ASM.1: H0 = {derived_h0:.2f} km/s/Mpc, window [{h0_lo:.1f}, {h0_hi:.1f}]. LOCKED 2026-09-04.")
    print(bar)

    fired_as_expected = abs(fired_sigma) > 1
    return 0 if fired_as_expected else 1


if __name__ == "__main__":
    raise SystemExit(main())