#!/usr/bin/env python3
"""The age of the universe at the redshifts Lumina covers, in the simulation's own cosmology.
Fact-Checked Science, "S0 News, October 28, 2026", section 4 ("Our analysis"). Standard library only.

Run it
  python3 -I cosmic_time.py
  It needs no input file.

How
  The age of a flat expanding universe at redshift z is
      t(z) = (1/H0) * integral from z to infinity of dz' / [(1+z') E(z')],
      E(z) = sqrt(Omega_r (1+z)^4 + Omega_m (1+z)^3 + Omega_L),
  the lookback-time integral of Hogg (1999), "Distance measures in cosmology", arXiv:astro-ph/9905116, section 10,
  equation 30, taken from z to infinity instead of from 0 to z. The parameters are Lumina's, from Table 2 of Zier et
  al. (2026), arXiv:2605.15310 (Planck 2018, TT,TE,EE+lowE+lensing+BAO): Omega_m = 0.3096, Omega_L = 0.6903,
  H0 = 67.66 km/s/Mpc, CMB temperature 2.7255 K, and three massless neutrino species at 1.9524 K, which the paper
  says Lumina includes in its expansion history along with the photons. Radiation:
      Omega_gamma h^2 = 2.4728e-5 (T_CMB / 2.7255 K)^4,  Omega_nu = 3 (7/8) (T_nu / T_CMB)^4 Omega_gamma.
  The integral runs over ln(a), the logarithm of the scale factor, with Simpson's rule.

What it prints
  The age at z = 1100 (the cosmic microwave background), 49 (Lumina's start), 23, 10, 7, 6, 5.2 (hydrogen fully
  ionized in Lumina's median sub-volume), 4.75, 4 and 3 (Lumina's end), and today, each as a share of the present
  age; then four checks: the same cosmology without radiation; Planck 2018's TT,TE,EE+lowE+lensing parameters
  (Omega_m = 0.3153, H0 = 67.36, whose Table 2 gives an age of 13.797 +/- 0.023 Gyr); its TT,TE,EE+lowE+lensing+BAO
  parameters, the column Lumina's come from (Omega_m = 0.3111, H0 = 67.66; Table 2: 13.787 +/- 0.020 Gyr); and a
  matter-only universe, whose age is exactly 2/(3 H0). Planck's Omega_m includes one neutrino of mass 0.06 eV
  (its Table 2 notes); Lumina counts its three neutrinos as massless (Zier et al., Table 2), which is why its
  Omega_m, 0.3096, is lower and its age today, 13.805 Gyr, higher.
"""
import math
import sys

MPC_KM = 3.0856775814913673e19   # kilometres in a megaparsec
GYR_S = 3.15576e16               # seconds in a billion Julian years


def age(z, H0, Om, OL, Or):
    """Age of the universe at redshift z, in billions of years."""
    t_hubble = MPC_KM / H0 / GYR_S
    x1, x0, n = -math.log1p(z), -60.0, 20000      # ln(a) from about 1e-26 to a = 1/(1+z)
    h = (x1 - x0) / n
    s = 0.0
    for i in range(n + 1):
        a = math.exp(x0 + i * h)
        f = 1.0 / math.sqrt(Or / a ** 4 + Om / a ** 3 + OL)
        s += (1 if i in (0, n) else 4 if i % 2 else 2) * f
    return t_hubble * s * h / 3.0


def main():
    sys.stdout.reconfigure(encoding="utf-8")
    H0, Om, OL = 67.66, 0.3096, 0.6903
    h = H0 / 100
    t_cmb, t_nu = 2.7255, 1.9524
    Og = 2.4728e-5 * (t_cmb / 2.7255) ** 4 / h ** 2
    Onu = 3 * 7 / 8 * (t_nu / t_cmb) ** 4 * Og
    Or = Og + Onu
    print(f"Lumina's cosmology: H0 = {H0}, Omega_m = {Om}, Omega_L = {OL}, Omega_r = {Or:.3e} "
          f"(photons {Og:.3e}, neutrinos {Onu:.3e})")
    t0 = age(0, H0, Om, OL, Or)
    print(f"Age today: {t0:.3f} billion years")
    for z in (1100, 49, 23, 10, 7, 6, 5.2, 4.75, 4, 3):
        t = age(z, H0, Om, OL, Or)
        print(f"  z = {z:>6}: {t * 1000:8.1f} million years after the Big Bang, {100 * t / t0:5.2f}% of the present "
              f"age, {t0 - t:6.2f} billion years before today")
    t3, t49 = age(3, H0, Om, OL, Or), age(49, H0, Om, OL, Or)
    print(f"Lumina's full-physics run, z = 49 to 3: {t49 * 1000:.0f} million to {t3:.2f} billion years after the "
          f"Big Bang; it stops {t0 - t3:.2f} billion years before today, at {100 * t3 / t0:.1f}% of cosmic history.")
    print("\nChecks")
    print(f"  without radiation: z = 3 at {age(3, H0, Om, OL, 0.0):.3f}, today {age(0, H0, Om, OL, 0.0):.3f} billion years")
    Om2, H02 = 0.3153, 67.36
    Or2 = 2.4728e-5 / (H02 / 100) ** 2 * (1 + 3 * 7 / 8 * (4 / 11) ** (4 / 3) * 3.046 / 3)
    print(f"  Planck 2018 TT,TE,EE+lowE+lensing: z = 3 at {age(3, H02, Om2, 1 - Om2 - Or2, Or2):.3f}, today "
          f"{age(0, H02, Om2, 1 - Om2 - Or2, Or2):.3f} billion years (Planck's Table 2: 13.797 +/- 0.023)")
    Om3, H03 = 0.3111, 67.66
    Or3 = 2.4728e-5 / (H03 / 100) ** 2 * (1 + 3 * 7 / 8 * (4 / 11) ** (4 / 3) * 3.046 / 3)
    print(f"  Planck 2018 TT,TE,EE+lowE+lensing+BAO: z = 3 at {age(3, H03, Om3, 1 - Om3 - Or3, Or3):.3f}, today "
          f"{age(0, H03, Om3, 1 - Om3 - Or3, Or3):.3f} billion years (Planck's Table 2: 13.787 +/- 0.020)")
    print(f"  matter only, H0 = 67.66: {age(0, 67.66, 1.0, 0.0, 0.0):.3f} billion years; "
          f"2/(3 H0) = {2 / 3 * MPC_KM / 67.66 / GYR_S:.3f}")


if __name__ == "__main__":
    main()