"""
dct_edges_21cm.py

DCT EDGES 21-cm absorption test.
At z = 17 the BEC has not yet condensed: the Avrami profile gives
P(z=17) ~ 0.9999, i.e. P -> 1 at high z. The DCT-induced modification of the
21-cm brightness temperature is therefore negligible and the prediction
reduces to the standard 21-cm picture, so DCT is NEUTRAL on the EDGES vs
SARAS-3 controversy.

Reproduces: P(z=17) ~ 0.9999, T_21 modification ~ 0; status = NEUTRAL.
Public archive: EDGES (Bowman et al. 2018) absorption depth at z ~ 17.
Run with: python3 dct_edges_21cm.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
H0_PLANCK = 67.36
OMEGA_M   = 0.315
OMEGA_L   = 0.685
T_50_GYR  = 1.5     # Avrami half-formation timescale in Gyr
H_TO_GYR  = 977.8   # 1/(km/s/Mpc) -> Gyr


def hubble(z: float) -> float:
    H0 = H0_PLANCK / np.sqrt(P_0)
    return H0 * np.sqrt(OMEGA_M * (1.0 + z) ** 3 + OMEGA_L)


def lookback_gyr(z_lim: float, n: int = 4000) -> float:
    """Lookback time from z=0 to z=z_lim in Gyr."""
    zs = np.linspace(0.0, z_lim, n)
    integrand = 1.0 / ((1.0 + zs) * hubble(zs))
    return float(_trapz(integrand, zs) * H_TO_GYR)


def avrami_P(t_lookback_gyr: float) -> float:
    """Avrami profile: P -> 1 at high z (large lookback), P_0 at z = 0."""
    return P_0 + (1.0 - P_0) * (1.0 - np.exp(-t_lookback_gyr / T_50_GYR))


def main() -> None:
    z_edges = 17.0
    t_lb    = lookback_gyr(z_edges)
    P_at_z  = avrami_P(t_lb)

    # Standard 21-cm absorption depth at z ~ 17 (EDGES central value):
    T_21_standard = -0.500  # K (ish; standard CMB-temp baseline)
    # DCT modification scales as (P - 1) since P -> 1 means no scalar effect.
    deviation     = (P_at_z - 1.0)            # ~ -1e-4 at z = 17
    T_21_dct      = T_21_standard * (1.0 + deviation)
    delta_T       = T_21_dct - T_21_standard  # ~ a few x 10^-5 K

    print("=" * 60)
    print("DCT EDGES 21-cm absorption test")
    print("=" * 60)
    print(f"At z = {z_edges}  (Avrami profile)")
    print(f"  lookback time t = {t_lb:.2f} Gyr")
    print(f"  P(z=17)         = {P_at_z:.6f}")
    print(f"  (P-1)           = {deviation:+.2e}")
    print()
    print(f"Standard 21-cm absorption depth: {T_21_standard:.3f} K")
    print(f"DCT-modified absorption depth  : {T_21_dct:.6f} K")
    print(f"Delta T                        : {delta_T:+.2e} K  (negligible)")
    print()
    print("Status: NEUTRAL.")
    print("DCT inactive at z=17, neutral pending SARAS-3 vs EDGES resolution.")
    print("=" * 60)


if __name__ == "__main__":
    main()