"""
LOS AQUAL SIMULATION (the referee's ask: pin xbar from the real density field, don't infer it).
The screen boost is Q=<1/mu(x)>, x=g_pec/a0, a0=2sqrt(alpha)cH0=1.197e-10 m/s^2. KEY physics: the peculiar
GRAVITATIONAL acceleration g_pec is a LARGE-SCALE quantity (g_k ~ delta_k/k, Poisson), so it is set by the
environment at the smoothing scale R the screen samples, NOT by small-scale clumps. Compute x_rms(R) and the
volume-weighted <1/mu> from the real CAMB P(k); find what R the BAO-need (Q~5, xbar~0.24) corresponds to, and
whether that R is the physical void/large-scale environment galaxy sightlines traverse. honesty-guard both ways:
report x(R) honestly; do not tune R to hit the answer.
"""
import numpy as np
from scipy.integrate import quad
import camb
import dfd_params as P

ALPHA=P.ALPHA; twosqrta=2*np.sqrt(ALPHA)
c_si=2.998e8; Mpc=3.0857e22; G=6.674e-11
H0_si=P.H0*1e3/Mpc
a0=twosqrta*c_si*H0_si
rho_crit=3*H0_si**2/(8*np.pi*G); rho_m=P.OMEGA_M*rho_crit

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, 0.5], kmax=50.0)
res=camb.get_results(pars)
PK=camb.get_matter_power_interpolator(pars, nonlinear=True, kmax=50, zmax=1)
s8=res.get_sigma8_0()

# peculiar accel: |g_k| = (3/2)Om H0^2 delta_k/k (phys, a=1).  Smoothed at scale R (Gaussian W=exp(-k^2R^2/2)).
# sigma_g(R)^2 = [(3/2)Om H0^2]^2 * (1/2pi^2) int P(k) W(kR)^2 dk   (the 1/k^2 and k^2 cancel)
pref=(1.5*P.OMEGA_M*H0_si**2)   # (3/2)Om H0^2  [1/s^2], SI
h=P.H0/100
def sigma_g(R_Mpch):
    # sigma_g^2 = pref^2 * (1/2pi^2) * int P(k) W(kR)^2 dk ; k in 1/Mpc, P in Mpc^3 -> J in Mpc^2
    R_Mpc=R_Mpch/h
    def ip(kM):  # kM in 1/Mpc
        W=np.exp(-(kM*R_Mpc)**2/2)
        return PK.P(0.5,kM)*W**2
    J=quad(ip,1e-3,50,limit=400)[0]               # [Mpc^2]
    J_m2=J*Mpc**2                                  # -> m^2
    sg2=pref**2*(1/(2*np.pi**2))*J_m2              # (1/s^2)^2 * m^2 = (m/s^2)^2
    return np.sqrt(max(sg2,0))

def Qmu_volume(xrms):
    # volume-weighted <1/mu(x)>, x = xrms*|chi| with chi Maxwell (3D gaussian accel magnitude)
    gg=np.linspace(1e-3,8,4000); pdf=np.sqrt(2/np.pi)*gg**2*np.exp(-gg**2/2)  # Maxwell, <chi^2>=3 -> rescale
    pdf/=np.trapezoid(pdf,gg)
    x=xrms*gg/np.sqrt(3)        # so that rms(x)=xrms
    inv_mu=(1+x)/x
    return np.trapezoid(inv_mu*pdf,gg)

print("="*86)
print("LOS AQUAL SIM — pin xbar from the real peculiar-acceleration field (CAMB P(k), z=0.5)")
print("="*86)
print(f"  a0=2sqrt(alpha)cH0={a0:.3e} m/s^2 (MOND a0=1.2e-10)  sigma8={s8:.3f}")
print(f"  BAO-need: Q~5.1, xbar~0.24 (from anisotropic_screen.py)\n")
print(f"  {'R[Mpc/h]':>9s} {'sigma_g[m/s^2]':>14s} {'x_rms=sig_g/a0':>15s} {'<1/mu>_vol':>11s}")
for R in (2.,5.,8.,12.,20.,30.,50.):
    sg=sigma_g(R); xr=sg/a0; Q=Qmu_volume(xr)
    print(f"  {R:9.0f} {sg:14.3e} {xr:15.3f} {Q:11.2f}")
print("\n  READ: find the R where x_rms ~ 0.2-0.3 (Q~5). If that R is the physical large-scale/void")
print("  environment (R~10-30 Mpc/h, where galaxy LOS spend most path), the boost is DERIVED-in-range,")
print("  not tuned. If x_rms only reaches 0.24 at unphysically large R, the boost is not supported.")
print("  HONEST: g_pec is large-scale; whether the SCREEN samples this volume-filling low-g field (boost")
print("  lives) or the mass (clumps, boost dies) is the source-weighting question = the sky-map test.")
