"""
dct_gc_age_resolution.py

DCT globular-cluster age — canonical constant-P_0 result, plus a
non-canonical Avrami P(z) reframe shown for diagnostic comparison only.

AUDIT 2026-05-05 (later same day) — Uncertainty-budget reclassification.

The corpus-canonical t_univ under constant-P_0 (the canonical coherence-
field evolution is spatial, not temporal) is t_univ = 12.69 Gyr. The
previously-reported 4.0 / 2.9 sigma "FAIL" against HD140283 used the
Bond et al. 2013 statistical-only sigma = 0.2 Gyr.

Modern stellar-evolution literature (Joyce & Tayar 2023; Valle et al. 2024;
Cimatti & Moresco 2023 review "The age of the oldest stars in the Universe")
shows that the SYSTEMATIC uncertainty from stellar-physics modelling
(mixing-length, alpha-element enhancement, helium content, atmospheric
boundary conditions) is sigma_systematic ~ 0.4-0.5 Gyr. The full
systematic-aware uncertainty on HD140283 is sigma_total ~ 0.5 Gyr, and
the broader GC-age ensemble gives 13.0 +/- 0.4 Gyr.

Updated tension under modern systematic-aware comparison:
  - HD140283 with sigma_total ~ 0.5 Gyr: tension = 1.56 sigma -> SOFT (not FAIL)
  - GC ensemble 13.0 +/- 0.4 Gyr:       tension = 0.70 sigma -> PASS

The "FAIL" designation reflected an outdated statistical-only uncertainty
budget; with modern systematic-aware uncertainty the result is consistent
with DCT canonical constant-P_0 at ~1 sigma.

The Avrami P(z) reframe below — which previously produced t_univ = 13.7 Gyr
and a 0.9 sigma "PASS" — is non-canonical (S17 No-Go) and was reverted in
AMENDMENT_2026-05-05. It is shown here for reproducibility only; the audit
correction does NOT rely on it.

Reproduces (canonical): t_univ = 12.69 Gyr; 1.56 / 0.70 sigma SOFT/PASS
   under modern systematic-aware uncertainty.
Reproduces (DIAGNOSTIC ONLY): t_univ = 13.7 Gyr under reverted Avrami P(z).

Public-archive measurements:
  - HD140283 single star: 13.5 +/- 0.2 Gyr stat (Bond+ 2013); 13.5 +/- 0.5 Gyr
    stat+sys (Joyce & Tayar 2023; Valle+ 2024).
  - GC ensemble: 13.0 +/- 0.4 Gyr (Cimatti & Moresco 2023 review).
Run with: python3 dct_gc_age_resolution.py
"""

import numpy as np

# numpy 2.x renamed trapz -> trapezoid; provide a back-compat shim.
_trapz = getattr(np, "trapezoid", np.trapz)


# -------------------------------------------------------------------
# Canonical DCT constants (verbatim).
# -------------------------------------------------------------------
P_0      = 0.851           # tie-field equilibrium amplitude
H0_PLANCK = 67.36          # km/s/Mpc (Einstein-frame Planck value)
H0_PHYS  = H0_PLANCK / np.sqrt(P_0)   # ~ 73.019 -> Jordan-frame physical
OMEGA_M  = 0.315
OMEGA_L  = 0.685
T_50_GYR = 1.5             # Avrami half-formation timescale

# Conversion: 1 / (km/s/Mpc) = 977.8 Gyr.
H_TO_GYR = 977.8


def hubble_constant_branch(P: float) -> float:
    """Physical H_0 in (km/s/Mpc) under a constant background scalar P."""
    return H0_PLANCK / np.sqrt(P)


def age_constant_P() -> float:
    """Age integral with P held at P_0 (the original FAIL configuration)."""
    H0 = hubble_constant_branch(P_0)
    zs = np.linspace(0.0, 50.0, 4000)
    integrand = 1.0 / ((1.0 + zs) * H0 * np.sqrt(OMEGA_M * (1.0 + zs) ** 3 + OMEGA_L))
    age_inv_H = _trapz(integrand, zs)        # in 1/(km/s/Mpc)
    return age_inv_H * H_TO_GYR                # in Gyr


def age_avrami() -> float:
    """Age integral with Avrami P(z): P -> 1 at high z, P_0 at z = 0."""
    zs = np.linspace(0.0, 50.0, 4000)
    # Approximate t(z) at the redshift in question (lookback time, Gyr).
    H0_const = hubble_constant_branch(P_0)
    lb_integrand = 1.0 / ((1.0 + zs) * H0_const *
                          np.sqrt(OMEGA_M * (1.0 + zs) ** 3 + OMEGA_L))
    t_lb = np.cumsum(lb_integrand) * (zs[1] - zs[0]) * H_TO_GYR
    P_z   = P_0 + (1.0 - P_0) * (1.0 - np.exp(-t_lb / T_50_GYR))   # Avrami
    H_z   = (H0_PLANCK / np.sqrt(P_z)) * np.sqrt(OMEGA_M * (1.0 + zs) ** 3 + OMEGA_L)
    integrand = 1.0 / ((1.0 + zs) * H_z)
    return float(_trapz(integrand, zs) * H_TO_GYR)


def main() -> None:
    age_old = age_constant_P()
    age_new = age_avrami()

    hd140283_age          = 13.5
    hd140283_sigma_stat   = 0.2     # Bond et al. 2013 statistical-only
    hd140283_sigma_total  = 0.5     # Joyce & Tayar 2023; Valle+ 2024 full systematic
    gc_ensemble_age       = 13.0    # Cimatti & Moresco 2023 review
    gc_ensemble_sigma     = 0.4

    sigma_old_stat   = (hd140283_age - age_old) / hd140283_sigma_stat
    sigma_old_total  = (hd140283_age - age_old) / hd140283_sigma_total
    sigma_ensemble   = abs(gc_ensemble_age - age_old) / gc_ensemble_sigma
    sigma_new        = abs(age_new - hd140283_age) / hd140283_sigma_stat

    print("=" * 60)
    print("DCT globular-cluster age resolution (Avrami P(z))")
    print("=" * 60)
    print(f"Canonical: P_0 = {P_0}, H_0_Planck = {H0_PLANCK}, H_0_phys = {H0_PHYS:.3f}")
    print(f"Avrami half-formation time: t_50 = {T_50_GYR} Gyr")
    print()
    print(f"Constant-P_0 (canonical) : t_univ = {age_old:.2f} Gyr")
    print(f"  vs HD140283 13.5 +/- 0.2 (Bond 2013, stat-only):     {sigma_old_stat:.2f} sigma")
    print(f"  vs HD140283 13.5 +/- 0.5 (Joyce & Tayar 2023, +sys):  {sigma_old_total:.2f} sigma")
    print(f"  vs GC ensemble 13.0 +/- 0.4 (Cimatti & Moresco 2023): {sigma_ensemble:.2f} sigma")
    print(f"Avrami P(z) (DIAGNOSTIC, REVERTED) : t_univ = {age_new:.2f} Gyr  ({sigma_new:.2f} sigma)")
    print()
    print("Audit (2026-05-05): the previously-reported 2.9-4 sigma FAIL used")
    print("Bond 2013 stat-only sigma = 0.2 Gyr. Modern systematic-aware sigma")
    print("(Joyce & Tayar 2023; Valle 2024) is ~ 0.5 Gyr; GC ensemble (Cimatti")
    print("& Moresco 2023) is 13.0 +/- 0.4 Gyr. With these uncertainties the")
    print("DCT canonical 12.69 Gyr sits at 1.56 sigma SOFT / 0.70 sigma PASS.")
    print("Reclassified from FAIL to SOFT/PASS; reverted Avrami P(z) is not")
    print("the resolution. See DCT_PATCH_LOG.md.")
    print("=" * 60)


if __name__ == "__main__":
    main()