"""
DFD-CLASS — module 11: MEASURE x-bar = g_psi/a0 from first principles (the decisive number).

Module 10 reduced BAO closure to one quantity: the LOS acceleration environment x_bar = g/a0.
  x_bar ~ 1   -> screen ~ 0.022 -> CLOSES
  x_bar << 1  -> 1/mu blows up -> screen >> 0.022 -> OVERSHOOTS (breaks ALL distances)
  x_bar >> 1  -> 1/mu -> 1 -> screen ~ 0.007 -> SHORT
So we don't guess x_bar -- we COMPUTE it from the measured matter power spectrum. The peculiar
gravitational acceleration field follows from P(k):
   g(x) = grad(potential),  |g_k| = 4 pi G rho_bar |delta_k| / k
   sigma_g^2(R) = (4 pi G rho_bar)^2 * int dk/(2pi^2) P(k) W^2(kR)
and a0(z) = 2 sqrt(alpha) c H(z). x_bar(R) = sigma_g(R)/a0. Honest: report whatever it gives,
including OVERSHOOT (which would be as fatal as undershoot). Honesty-guard both ways.
"""

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

# SI constants
G = 6.674e-11; c = 2.998e8; Mpc = 3.0857e22
H0_si = P.H0 * 1e3 / Mpc                     # 1/s
rho_crit0 = 3 * H0_si**2 / (8 * np.pi * G)   # kg/m^3
ALPHA = P.ALPHA


def get_pk(z):
    pars = camb.CAMBparams()
    pars.set_cosmology(H0=P.H0, ombh2=P.OMEGA_B_H2, omch2=P.OMEGA_CHI_H2, mnu=0.06, tau=0.054)
    pars.InitPower.set_params(As=P.A_S, ns=P.N_S)
    pars.set_matter_power(redshifts=[z], kmax=50.0)
    PKlin = camb.get_matter_power_interpolator(pars, nonlinear=False, kmax=50, zmax=z+0.1)
    PKnl = camb.get_matter_power_interpolator(pars, nonlinear=True, kmax=50, zmax=z+0.1)
    return PKlin, PKnl


def Hofz(z):
    return H0_si * np.sqrt(P.OMEGA_M*(1+z)**3 + P.OMEGA_R*(1+z)**4 + P.OMEGA_LAMBDA)


def a0(z):
    return 2*np.sqrt(ALPHA) * c * Hofz(z)


def sigma_g(R_Mpc, z, PK):
    """RMS peculiar gravitational acceleration smoothed on scale R [Mpc], SI m/s^2."""
    rho_bar = P.OMEGA_M * rho_crit0 * (1+z)**3            # kg/m^3
    pref = (4*np.pi*G*rho_bar)**2
    def integ(lnk):
        k = np.exp(lnk)            # 1/Mpc
        Pk = PK.P(z, k) * Mpc**3   # m^3
        k_si = k / Mpc             # 1/m
        W = np.exp(-0.5*(k*R_Mpc)**2)
        # int dk_si P_si W^2 = int dlnk_mpc * k_si * P_si * W^2  (k_si converts dlnk_mpc -> dk_si)
        return k_si * (Pk / (2*np.pi**2)) * W**2
    I = quad(integ, np.log(1e-3), np.log(40), limit=120)[0]   # m^2 (SI: P[m^3]*k_si[1/m])
    return np.sqrt(pref * I)       # m/s^2


def main():
    z = 0.5
    PKlin, PKnl = get_pk(z)
    print("=" * 72)
    print(f"MODULE 11 — measuring x_bar = g_psi/a0 from P(k) at z={z}")
    print("=" * 72)
    print(f"  a0(z={z}) = {a0(z):.3e} m/s^2  (a0(0) = {a0(0):.3e}, known 1.2e-10)")
    print(f"\n  x_bar(R) = sigma_g(R)/a0 vs smoothing scale R:")
    print(f"{'R [Mpc]':>8s} {'sigma_g (lin)':>14s} {'x_bar lin':>10s} {'x_bar NL':>10s}")
    for R in [1, 2, 5, 10, 20, 50]:
        sg_l = sigma_g(R, z, PKlin)
        sg_n = sigma_g(R, z, PKnl)
        print(f"{R:8.0f} {sg_l:14.3e} {sg_l/a0(z):10.3f} {sg_n/a0(z):10.3f}")

    # the screen-relevant x_bar: nonlinear, small-scale (galaxy environment ~1-2 Mpc)
    sg_small = sigma_g(1.5, z, PKnl)
    x_small = sg_small/a0(z)
    sg_large = sigma_g(30, z, PKlin)
    x_large = sg_large/a0(z)

    print("\n" + "-" * 72)
    print("HONEST VERDICT (module 11 — the decisive number):")
    print(f"  Large-scale (bulk, ~30 Mpc) x_bar = {x_large:.3f}  (deep-MOND)")
    print(f"  Small-scale (galaxy env, ~1.5 Mpc, nonlinear) x_bar = {x_small:.3f}")
    print(f"  Closure needs x_bar ~ 1 (module 10). The peculiar-acceleration field gives a")
    print(f"  SCALE-DEPENDENT x_bar spanning {x_large:.2f} (large scale) to {x_small:.2f} (galaxy env).")
    if x_small < 0.5:
        print(f"  => x_bar < 1 on most scales -> the mu-boost 1/mu={1/(x_small/(1+x_small)):.1f} is LARGE ->")
        print(f"     the screen would OVERSHOOT 0.022 substantially. That is AS FATAL as undershoot:")
        print(f"     it would break the monopole/SNe/all distances, not just miss BAO. HONEST: this")
        print(f"     calc does NOT support a clean closure; it points to OVER-correction unless the")
        print(f"     screen is sourced by the higher-g (small-scale/halo) regions where x_bar -> 1.")
    elif 0.5 <= x_small <= 2.0:
        print(f"  => x_bar ~ {x_small:.1f} ~ 1 on galaxy-environment scales -> screen ~ 0.022 -> CLOSURE")
        print(f"     is SUPPORTED at the natural scale. Encouraging, consistent with module 10's x~1.")
    else:
        print(f"  => x_bar = {x_small:.1f} >> 1 -> 1/mu -> 1 -> screen ~ 0.007 -> SHORT. BAO does not close.")
    print(f"  KEY UNCERTAINTY now SHARP: WHICH scale sources the optical-psi screen? Large-scale")
    print(f"  (bulk, x~{x_large:.2f}) -> overshoot; galaxy-env (x~{x_small:.2f}) -> {'close' if 0.5<=x_small<=2 else 'over/under'}. The corpus")
    print(f"  monopole Delta-psi~0.18 normalizes this; the differential needs the sourcing-scale calc.")
    return dict(x_large=x_large, x_small=x_small, a0=a0(z))


if __name__ == "__main__":
    main()
