"""
DFD-CLASS — Background module.

Computes H(z), distances, the sound horizon r_drag, the CMB acoustic scale
theta*, and BAO/SNe/H0 observables for the DFD background, with all densities
DERIVED (see dfd_params.py). Validates numerically against real data:
  - SH0ES local H0 (Riess+ 2022)
  - BAO consensus (6dFGS, MGS, BOSS DR12, eBOSS QSO, Lya)
  - Planck 2018 LambdaCDM as a reference cosmology

Honest grading: real chi^2 numbers, no fudge. Where DFD's derived background
mismatches a datum, it is reported as a tension for later modules (the optical-psi
distance law and the psi-gravity perturbation closure are the levers).
"""

import numpy as np
from scipy.integrate import quad
import dfd_params as P


# ------------------------------------------------------------------ background
def E(z, par=P):
    """Dimensionless Hubble rate H(z)/H0 for w=-1 (optical-psi DE)."""
    zp1 = 1.0 + z
    return np.sqrt(par.OMEGA_R * zp1**4 + par.OMEGA_M * zp1**3 + par.OMEGA_LAMBDA)


def Hofz(z, par=P):
    return par.H0 * E(z, par)


def _Elcdm(z, pl=P.PLANCK):
    zp1 = 1.0 + z
    Or = (P.OMEGA_R_H2) / pl["h"]**2  # same radiation
    Om = pl["Omega_m"]
    OL = 1.0 - Om - Or
    return np.sqrt(Or * zp1**4 + Om * zp1**3 + OL)


def Hofz_lcdm(z, pl=P.PLANCK):
    return pl["H0"] * _Elcdm(z)


# ------------------------------------------------------------------ distances
def comoving_distance(z, par=P, Efun=None):
    """D_C(z) [Mpc] (flat)."""
    Efun = Efun or (lambda zz: E(zz, par))
    integ, _ = quad(lambda zz: 1.0 / Efun(zz), 0.0, z, limit=200)
    return (P.C_KMS / par.H0) * integ


def D_M(z, par=P):
    return comoving_distance(z, par)


def D_A(z, par=P):
    return D_M(z, par) / (1.0 + z)


def D_L(z, par=P):
    return D_M(z, par) * (1.0 + z)


def D_H(z, par=P):
    return P.C_KMS / Hofz(z, par)


def D_V(z, par=P):
    dm = D_M(z, par)
    return (z * dm**2 * D_H(z, par))**(1.0 / 3.0)


# ------------------------------------------------------------- sound horizon
def _R_baryon(z, par=P):
    """Baryon-photon momentum ratio R = 3 rho_b / 4 rho_gamma."""
    return (3.0 * par.OMEGA_B_H2) / (4.0 * par.OMEGA_G_H2) / (1.0 + z)


def _cs_over_c(z, par=P):
    return 1.0 / np.sqrt(3.0 * (1.0 + _R_baryon(z, par)))


def sound_horizon(z_end, par=P, z_max=1.0e6):
    """r_s(z_end) = int_{z_end}^inf c_s/H dz  [Mpc]."""
    def integrand(z):
        return P.C_KMS * _cs_over_c(z, par) / Hofz(z, par)
    integ, _ = quad(integrand, z_end, z_max, limit=400)
    return integ


def r_drag(par=P):
    return sound_horizon(P.Z_DRAG, par)


def theta_star(par=P):
    """CMB acoustic angular scale theta* = r_s(z*) / D_M(z*)."""
    rs = sound_horizon(P.Z_STAR, par)
    dm = D_M(P.Z_STAR, par)
    return rs / dm


def age_Gyr(par=P):
    """Age of universe today [Gyr]."""
    integ, _ = quad(lambda z: 1.0 / ((1.0 + z) * Hofz(z, par)), 0.0, 1.0e5, limit=400)
    return integ * (977.792)  # (1 Mpc/km/s) -> Gyr conversion: 1/(km/s/Mpc)=977.79 Gyr


# ------------------------------------------------------------------ BAO data
# Consensus BAO (dimensionless ratios vs r_drag). type: 'DV', or 'DM'/'DH' pair.
BAO = [
    dict(name="6dFGS",        z=0.106, kind="DV", val=2.976, err=0.133),
    dict(name="SDSS MGS",     z=0.15,  kind="DV", val=4.466, err=0.168),
    dict(name="BOSS z1 DM",   z=0.38,  kind="DM", val=10.23, err=0.17),
    dict(name="BOSS z1 DH",   z=0.38,  kind="DH", val=25.00, err=0.76),
    dict(name="BOSS z2 DM",   z=0.51,  kind="DM", val=13.36, err=0.21),
    dict(name="BOSS z2 DH",   z=0.51,  kind="DH", val=22.33, err=0.58),
    dict(name="BOSS z3 DM",   z=0.61,  kind="DM", val=15.45, err=0.22),
    dict(name="BOSS z3 DH",   z=0.61,  kind="DH", val=20.86, err=0.55),
    dict(name="eBOSS QSO DM", z=1.48,  kind="DM", val=30.21, err=0.79),
    dict(name="eBOSS QSO DH", z=1.48,  kind="DH", val=13.23, err=0.47),
    dict(name="Lya DM",       z=2.33,  kind="DM", val=37.50, err=1.10),
    dict(name="Lya DH",       z=2.33,  kind="DH", val=8.99,  err=0.19),
]


def bao_predict(z, kind, rd, par=P):
    if kind == "DV":
        return D_V(z, par) / rd
    if kind == "DM":
        return D_M(z, par) / rd
    if kind == "DH":
        return D_H(z, par) / rd
    raise ValueError(kind)


def bao_chi2(par=P, rd=None):
    rd = rd if rd is not None else r_drag(par)
    chi2 = 0.0
    rows = []
    for b in BAO:
        pred = bao_predict(b["z"], b["kind"], rd, par)
        r = (pred - b["val"]) / b["err"]
        chi2 += r * r
        rows.append((b["name"], b["z"], b["val"], b["err"], pred, r))
    return chi2, rows


# ------------------------------------------------------------------- report
def main():
    P.summary()
    print("\n" + "=" * 64)
    print("BACKGROUND OBSERVABLES (DFD derived background)")
    print("=" * 64)

    rd = r_drag()
    rs_star = sound_horizon(P.Z_STAR)
    th = theta_star()
    age = age_Gyr()

    print(f"  r_drag                = {rd:8.3f} Mpc   (Planck {P.PLANCK['r_drag']:.2f})")
    print(f"  r_s(z*)               = {rs_star:8.3f} Mpc")
    print(f"  100*theta*            = {100*th:8.5f}     (Planck {100*P.PLANCK['theta_star']:.5f})")
    print(f"  age                   = {age:8.3f} Gyr   (Planck {P.PLANCK['age_Gyr']:.3f})")
    print(f"  H0 (DFD)              = {P.H0:8.3f}       SH0ES 73.04 +/- 1.04")
    h0_pull = (P.H0 - 73.04) / 1.04
    print(f"    -> H0 pull vs SH0ES = {h0_pull:+.2f} sigma  "
          f"(Planck-LCDM pull = {(P.PLANCK['H0']-73.04)/1.04:+.2f} sigma)")

    # theta* tension
    dth = 100*(th - P.PLANCK['theta_star'])/P.PLANCK['theta_star']
    print(f"  theta* vs Planck      = {dth:+.2f}%   "
          f"({'MATCH' if abs(dth)<0.1 else 'TENSION -> optical-psi distance lever'})")

    print("\n" + "-" * 64)
    print("BAO fit (DFD background, derived r_drag):")
    print(f"{'datum':14s} {'z':>5s} {'obs':>8s} {'err':>6s} {'DFD':>8s} {'pull':>7s}")
    chi2, rows = bao_chi2()
    for name, z, val, err, pred, r in rows:
        print(f"{name:14s} {z:5.2f} {val:8.3f} {err:6.3f} {pred:8.3f} {r:+7.2f}")
    print(f"  BAO chi^2 (DFD)   = {chi2:.2f}  over {len(rows)} points "
          f"(chi^2/N = {chi2/len(rows):.2f})")

    # LCDM BAO for contrast (Planck cosmology, its own r_drag)
    class _PL:
        H0 = P.PLANCK["H0"]; OMEGA_M = P.PLANCK["Omega_m"]
        OMEGA_R = P.OMEGA_R_H2 / P.PLANCK["h"]**2
        OMEGA_LAMBDA = 1 - P.PLANCK["Omega_m"] - OMEGA_R
        OMEGA_B_H2 = P.OMEGA_B_H2; OMEGA_G_H2 = P.OMEGA_G_H2
    chi2_l, _ = bao_chi2(par=_PL, rd=P.PLANCK["r_drag"])
    print(f"  BAO chi^2 (LCDM)  = {chi2_l:.2f}  over {len(rows)} points "
          f"(chi^2/N = {chi2_l/len(rows):.2f})")

    return dict(r_drag=rd, theta_star=th, age=age, bao_chi2=chi2, bao_chi2_lcdm=chi2_l)


if __name__ == "__main__":
    main()
