"""
DFD-CLASS — module 8: the ANISOTROPIC / density-dependent screen (the density-anisotropy lever + the referee bridge).

Module 6 only proved a SKY-AVERAGED MONOTONE screen can't fit BAO+CMB. But the corpus screen is
Delta-psi(z, n_hat) = LOS integral of the density field, so each probe X samples it with a DIFFERENT
window S_X:  D_X = exp(<Delta-psi>_X) D_bare,  <Delta-psi>_X = (int dOmega S_X Delta-psi)/(int dOmega S_X).
CMB = near-full-sky (monopole). BAO/SNe = galaxy-selected, structure-correlated -> different <Delta-psi>.

Bridge (the referee, correct): a LambdaCDM fitter of an optically-dressed distance infers
   w_obs(z) = -1 + (1+z)/3 * d<Delta-psi>/dz.

THE DECISIVE QUESTION (honest, honesty-guard both ways): is the DFD screen big enough?
  - Standard (linear, no mu-law) galaxy-weighted excess ~ 0.007 (module-7 first cut) -> 3x short of the
    ~0.022 BAO differential.
  - DFD AMPLIFIER: the optical field obeys div[mu(g/a0) grad psi] = (8piG/c^2)rho, so the response is
    1/mu(xbar) enhanced. In low-acceleration (void/large-scale) regions xbar is small -> 1/mu large.
  We compute what enhancement Q_mu = 1/mu(xbar) closes BAO, the xbar it implies, and whether THAT xbar is
  CONSISTENT with the S8 win (which needs only MODEST growth enhancement at 8 Mpc/h).
"""

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

ALPHA = P.ALPHA
A0_OVER_cH = 2.0 * np.sqrt(ALPHA)     # a0(z) = 2 sqrt(alpha) c H(z)


def camb_pk():
    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=[0., 0.3, 0.5, 0.7, 1.0], kmax=30.0)
    res = camb.get_results(pars)
    PK = camb.get_matter_power_interpolator(pars, nonlinear=True, kmax=30, zmax=2)
    s8 = res.get_sigma8_0()
    return PK, s8


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


def chi(z):
    return quad(lambda zz: P.C_KMS / Hofz(zz), 0, z, limit=100)[0]


def sigma_kappa(zs, PK):
    """RMS lensing convergence to source at zs (linear/nonlinear P(k)); the screen fluctuation scale."""
    chis = chi(zs); c = P.C_KMS; Om = P.OMEGA_M; H0 = P.H0
    def W(chil):
        zl = np.interp(chil, [chi(z) for z in np.linspace(0, zs, 20)], np.linspace(0, zs, 20))
        a = 1/(1+zl)
        return 1.5*Om*(H0/c)**2 * chil*(chis-chil)/chis / a, zl
    def Cl(l):
        def integ(chil):
            if chil <= 0 or chil >= chis: return 0.0
            w, zl = W(chil)
            k = l/chil
            if k < 1e-4 or k > 30: return 0.0
            return w**2/chil**2 * PK.P(zl, k)
        return quad(integ, 1.0, chis-1.0, limit=60)[0]
    ls = np.logspace(1, 3.3, 40)
    cl = np.array([Cl(l) for l in ls])
    var = np.trapezoid(ls*cl/(2*np.pi), ls)
    return np.sqrt(max(var, 0))


def mu(x):
    return x / (1 + x)


def main():
    print("=" * 74)
    print("MODULE 8 — anisotropic density-dependent screen: does it close BAO? (honest)")
    print("=" * 74)
    PK, s8 = camb_pk()

    # 1) standard galaxy-weighted screen excess (no mu enhancement) at BAO redshifts
    bg = 1.8        # BAO tracer bias
    r = 0.7         # psi-galaxy correlation coefficient
    print("\n[1] Standard (no mu-law) galaxy-weighted screen excess <Dpsi>_gal = b_g r sigma_kappa:")
    print(f"{'z':>6s} {'sigma_kappa':>12s} {'<Dpsi>_gal':>12s} {'BAO need':>10s}")
    needed = {0.38: 0.036, 0.51: 0.041, 0.61: 0.033}   # from module 6 (ln[data/DFD] at DM points)
    std = {}
    for z in [0.38, 0.51, 0.61]:
        sk = sigma_kappa(z, PK)
        dpsi_std = bg * r * sk
        std[z] = dpsi_std
        print(f"{z:6.2f} {sk:12.4f} {dpsi_std:12.4f} {needed[z]:10.3f}")
    z_mid = 0.51
    diff_need = needed[z_mid] - 0.012     # BAO need minus the CMB monopole requirement
    print(f"  BAO-vs-CMB differential needed at z~0.5: ~{diff_need:.3f}")
    print(f"  standard excess ~{std[z_mid]:.3f} -> shortfall factor Q_need = {diff_need/std[z_mid]:.1f}")

    # 2) the DFD mu-law enhancement Q_mu = 1/mu(xbar); what xbar gives Q_need?
    Q_need = diff_need / std[z_mid]
    # 1/mu(x) = Q -> x/(1+x)=1/Q -> x = 1/(Q-1)
    x_need = 1.0/(Q_need - 1) if Q_need > 1 else np.inf
    print(f"\n[2] DFD mu-law amplifier needed: Q_mu = 1/mu(xbar) = {Q_need:.1f} -> xbar_screen = {x_need:.3f}")
    print(f"   (xbar = g_psi/a0; deep-MOND xbar->0 gives huge Q, Hubble-EFE xbar=5.85 gives Q=1.17)")

    # 3) consistency with S8: the growth enhancement at 8 Mpc/h must stay MODEST (module 4: ~1.17)
    #    8 Mpc/h corresponds to higher density / higher g than the void-dominated screen LOS.
    x_sigma8 = 1.0/(2*np.sqrt(ALPHA))     # Hubble-EFE xbar = 5.85 (the sigma8-scale regime, module 4)
    Q_sigma8 = 1.0/mu(x_sigma8)
    print(f"\n[3] S8 consistency: at the sigma8 scale (8 Mpc/h, denser/higher-g), xbar~{x_sigma8:.2f} "
          f"-> Q_growth={Q_sigma8:.2f} (modest, the module-4 value that WON S8).")
    consistent = x_need < x_sigma8   # screen samples LOWER-g (voids) than sigma8 -> smaller xbar -> larger Q: OK
    print(f"   The screen LOS is void/large-scale dominated (LOWER g => smaller xbar => larger Q) than the")
    print(f"   8 Mpc/h sigma8 scale. Required xbar_screen={x_need:.2f} {'<' if consistent else '>'} "
          f"xbar_sigma8={x_sigma8:.2f}: {'SCALE-CONSISTENT (screen boosted, growth modest)' if consistent else 'INCONSISTENT'}.")

    # 4) implied w_obs(z) drift from the bridge, vs the BAO-needed w
    print(f"\n[4] Bridge w_obs(z) = -1 + (1+z)/3 d<Dpsi>/dz. With the enhanced screen reaching the BAO")
    print(f"   need, the IMPLIED apparent w at z~0.5 drifts from -1 by ~(1+z)/3 * d(0.02-0.04)/dz ~ few %.")
    print(f"   This is a PERCENT-LEVEL apparent-w drift from a FIXED topological vacuum -- exactly the")
    print(f"   corpus regime (section_cosmology:991) and the DESI-DR2 mild-dynamical-w window.")

    print("\n" + "-" * 74)
    print("HONEST VERDICT (module 8):")
    if Q_need <= 1:
        print("  Standard screen ALREADY exceeds the need -> closes BAO. (unexpected; recheck)")
    elif consistent and x_need > 0.05:
        print(f"  VIABLE: closing BAO needs a mu-law boost Q~{Q_need:.1f} (xbar_screen~{x_need:.2f}), which is")
        print(f"  a MODERATE deep-MOND enhancement on void/large-scale LOS -- and it is SCALE-CONSISTENT with")
        print(f"  the modest growth (Q~{Q_sigma8:.2f}) that WON S8 at 8 Mpc/h. So one mu-field can plausibly")
        print(f"  give a BIG screen (closes BAO) AND modest growth (keeps S8). NOT proven -- needs the full")
        print(f"  scale-dependent LOS integral with the real void/g_psi distribution -- but NOT excluded.")
    else:
        print(f"  STRAINED: closing BAO needs xbar_screen~{x_need:.2f}, a strong deep-MOND boost that risks")
        print(f"  over-enhancing growth -> tension with the S8 win. Honest: the escape hatch is tight.")
    print(f"  CLEAN FALSIFIABLE DISCRIMINATOR (ΛCDM cannot make it): BAO/SNe distance residuals must")
    print(f"  CORRELATE with foreground structure (lensing kappa / galaxy density), amplitude ~{std[z_mid]*Q_need:.2f}")
    print(f"  at z~0.5; a pure cosmological constant predicts ZERO. Testable now vs DES+Planck-lensing.")
    return dict(std=std, Q_need=Q_need, x_need=x_need, x_sigma8=x_sigma8, consistent=consistent)


if __name__ == "__main__":
    main()
