"""
GATE DERIVATION (rigorous, from the optical kernel) + FULL CMB TT/TE/EE confirmation.
Honesty-guard BOTH ways: derive what the physics says, even if it deflates the earlier 2sqrt(alpha)*f_b match.

PART A — which object does f_b gate?
  The optical metric is n=e^psi. Two distinct couplings:
   * PHOTONS follow null geodesics of g_opt = e^{2psi} eta. Distances (D_L, D_A, D_M) are photon paths
     -> they feel the FULL screen Delta-psi. NO baryon gate on the distance screen.
   * The BARYON ACOUSTIC RULER (r_d) is a property of the photon-BARYON fluid at decoupling; baryons
     carry rest mass and load the oscillator (R=3rho_b/4rho_gamma). Only the RULER feels f_b=R/(1+R).
  => f_b gates the RULER r_d, NOT the distance screen. (This is the referee's 'photons conformal, baryons
     feel psi' stated precisely: it makes r_d^b != r_s^gamma, it does NOT gate D_A.)
  Consequence to test honestly: the earlier A = 2sqrt(alpha)*f_b for the DISTANCE-screen amplitude has
  NO mechanism -- f_b does not act there. So is the distance-screen amplitude derived or fitted?

PART B — run CAMB TT/TE/EE for DFD's DERIVED params; check the intrinsic spectrum (heights, damping,
  polarization), which the screen does NOT touch (screen only sets the overall angle theta*).
"""
import numpy as np
from scipy.integrate import quad
import dfd_params as P

C, H0 = P.C_KMS, P.H0
ALPHA = P.ALPHA
twosqrta = 2*np.sqrt(ALPHA)
R_d = 0.640; f_b = R_d/(1+R_d)

def Dc(z, Om, h):                              # comoving distance, flat, matter+Lambda
    Or = P.OMEGA_R_H2/h**2; OL = 1-Om-Or
    return (C/(100*h))*quad(lambda zz: 1/np.sqrt(Om*(1+zz)**3+Or*(1+zz)**4+OL), 0, z, limit=90)[0]
def Dc_EdS(z, h):
    return (C/(100*h))*quad(lambda zz: 1/np.sqrt(1.0*(1+zz)**3), 0, z, limit=90)[0]

print("="*80)
print("PART A — GATE DERIVATION: what does f_b gate?")
print("="*80)
# corpus screen = ln(D_LCDM/D_EdS); H0-independent ratio
def dpsi_corpus(z, Om=0.30):
    h=0.70
    return np.log(Dc(z,Om,h)/Dc_EdS(z,h))
print("\n  Corpus screen Delta-psi(z)=ln(D_LCDM/D_EdS) [DERIVED, H0-independent, w_eff=-1 'by construction']:")
for z in (0.1,0.3,0.5,1.0,2.0):
    print(f"     z={z}: {dpsi_corpus(z):.3f}   (corpus table: 0.049,0.126,0.184,0.275,0.358)")
print(f"     z->inf saturates ~0.50 (ln(3.305/2.000)=0.502)")

print("\n  KERNEL split target = the SMALL RESIDUAL screen (kernel_split_test): ~0.030 at z=0.5.")
print(f"  Corpus full screen at z=0.5 = {dpsi_corpus(0.5):.3f}  >> 0.030. The corpus screen is DUAL to the")
print(f"  real rho_Lambda already in H(z) (NOT additive) -> applying it = LCDM distances, which at")
print(f"  H0=72.09 still leave theta* off (it does NOT resolve H0). The wall-closing residual (0.030)")
print(f"  is a SEPARATE object NOT in the corpus.")

print(f"\n  f_b RULER gate (the ONE derived lever): r_d^b = r_s * exp(-f_b*Delta-psi(z_d)/3)")
dpsi_rec = 0.499
print(f"     f_b={f_b:.3f}, Delta-psi(z_d)~{dpsi_rec}: r_d^b/r_s = exp(-{f_b*dpsi_rec/3:.4f}) = {np.exp(-f_b*dpsi_rec/3):.4f}")
print(f"     -> r_d 147 -> {147*np.exp(-f_b*dpsi_rec/3):.1f} Mpc (the T47 shift). This gates the RULER, derived.")
print(f"\n  HONEST GATE VERDICT: f_b gates the baryon RULER r_d (derived, baryon inertia). The DISTANCE")
print(f"  screen on photon paths is UNGATED. So the earlier A=2sqrt(alpha)*f_b={twosqrta*f_b:.4f} for the")
print(f"  DISTANCE amplitude has NO mechanism -> it was a NUMERICAL COINCIDENCE. The distance-screen")
print(f"  amplitude that pins theta*+BAO is NOT derived; it is a fitted H0=72.09 reconciliation.")

print("\n" + "="*80)
print("PART B — FULL CMB TT/TE/EE (CAMB): intrinsic spectrum on DFD's DERIVED params")
print("="*80)
import camb
def spec(H0_, ombh2, omch2, ns, As, tau=0.054):
    pars = camb.set_params(H0=H0_, ombh2=ombh2, omch2=omch2, ns=ns, As=As, tau=tau,
                           lmax=2500, WantTensors=False)
    r = camb.get_results(pars)
    cl = r.get_cmb_power_spectra(pars, CMB_unit='muK', raw_cl=False)['total']
    th = r.get_derived_params()['thetastar']
    return cl, th     # cl[:,0]=TT, [:,1]=EE, [:,2]=BB, [:,3]=TE  (Dl)
print("\n  running CAMB (DFD vs LCDM-Planck)...")
dfd, th_dfd = spec(P.H0, P.OMEGA_B_H2, P.OMEGA_CHI_H2, P.N_S, P.A_S)
lcdm, th_lcdm = spec(67.36, 0.02237, 0.12, 0.9649, 2.1e-9)
ell = np.arange(dfd.shape[0])
def peaks(TT, lo, hi):
    seg = TT[lo:hi]; i = np.argmax(seg)+lo; return i, TT[i]
print(f"  100*theta*: DFD={100*th_dfd:.4f}  LCDM={100*th_lcdm:.4f}  (Planck 1.04109; screen fixes the SCALE)")
print(f"\n  Intrinsic TT spectrum (screen-independent: heights, ratios, damping):")
print(f"  {'feature':22s} {'DFD':>12s} {'LCDM':>12s} {'DFD/LCDM':>10s}")
for name,(lo,hi) in [('peak1 (l~220)',(180,260)),('peak2 (l~530)',(480,600)),
                     ('peak3 (l~810)',(740,880)),('damping (l=1500)',(1499,1501))]:
    l1,h1 = peaks(dfd[:,0],lo,hi) if 'damping' not in name else (1500,dfd[1500,0])
    l2,h2 = peaks(lcdm[:,0],lo,hi) if 'damping' not in name else (1500,lcdm[1500,0])
    print(f"  {name:22s} {h1:12.1f} {h2:12.1f} {h1/h2:10.4f}")
p1=peaks(dfd[:,0],180,260)[1]; p2=peaks(dfd[:,0],480,600)[1]; p3=peaks(dfd[:,0],740,880)[1]
q1=peaks(lcdm[:,0],180,260)[1]; q2=peaks(lcdm[:,0],480,600)[1]; q3=peaks(lcdm[:,0],740,880)[1]
print(f"\n  Peak-height RATIOS (the no-CDM-wall test, screen-independent):")
print(f"     2nd/1st: DFD={p2/p1:.4f} LCDM={q2/q1:.4f}   3rd/1st: DFD={p3/p1:.4f} LCDM={q3/q1:.4f}")
print(f"     3rd/2nd: DFD={p3/p2:.4f} LCDM={q3/q2:.4f}")
# polarization sanity
te_dfd = dfd[100:1500,3]; te_lcdm = lcdm[100:1500,3]
ee_dfd = dfd[100:1500,1]; ee_lcdm = lcdm[100:1500,1]
print(f"\n  Polarization (l=100-1500, screen-independent shape):")
print(f"     TE rms DFD/LCDM = {np.sqrt(np.mean(te_dfd**2))/np.sqrt(np.mean(te_lcdm**2)):.4f}")
print(f"     EE rms DFD/LCDM = {np.sqrt(np.mean(ee_dfd**2))/np.sqrt(np.mean(ee_lcdm**2)):.4f}")
print("\n  READ: peak-height ratios + damping + TE/EE are set by Omega_b h^2, Omega_chi h^2, n_s, A_s")
print("  (all DERIVED) and are INDEPENDENT of the distance screen. If DFD/LCDM ~ 1 here, the intrinsic")
print("  spectrum is reproduced; the only screen-dependent quantity is the overall angle theta*.")
