"""
FOUR-CASE DECOMPOSITION (the referee's exact program, "no rhetoric").
The real theorem is NOT "conformal screen can't move theta*"; it is "ONLY relative optical/material scaling can
move theta*." DFD's optical metric screens PHOTON propagation; the acoustic ruler r_s is set by baryon-trace
matter physics. So we must check the four cases explicitly and see whether the consistent TWO-SECTOR model
survives, with actual numbers. Honesty-guard BOTH ways: don't fake a cancellation, don't fake a survival.

theta* = r_s(z*) / D_M(z*).   r_s forms at z>z* (early); the optical-DE screen lives at low z (late).
"""
import numpy as np
from scipy.integrate import quad, cumulative_trapezoid
from scipy.optimize import brentq
import dfd_params as P
import background as B

C, H0 = P.C_KMS, P.H0
ZSTAR, ZDRAG = P.Z_STAR, P.Z_DRAG
ALPHA = getattr(P, 'ALPHA', 1/137.035999084)
THETA_PL, THETA_ERR = 1.04109, 0.00030

# ---- raw DFD background (no screen) ----
rs_star = B.sound_horizon(ZSTAR)         # ruler, formed at z>z*
rd      = B.r_drag()
DMstar  = B.D_M(ZSTAR)
theta_raw = 100*rs_star/DMstar
print("="*86)
print("FOUR-CASE DECOMPOSITION OF theta*  (DFD optical/material split)")
print("="*86)
print(f"\n  raw DFD: r_s(z*)={rs_star:.3f} Mpc,  D_M(z*)={DMstar:.1f} Mpc,  100 theta*={theta_raw:.5f}")
print(f"  Planck 100 theta*={THETA_PL} +/- {THETA_ERR}  ->  raw DFD off by "
      f"{(theta_raw-THETA_PL)/THETA_ERR:+.1f} sigma")

# screen shape: late-time optical-DE energy fraction Omega_Lambda(z), amplitude A
OmL0 = P.OMEGA_LAMBDA
def shape(z):   # ~1 at z=0, ->0 by z~3, EXACTLY 0 at recombination
    return OmL0/(OmL0 + P.OMEGA_M*(1+z)**3 + P.OMEGA_R*(1+z)**4)

# integration grids
zlo = np.linspace(0, ZSTAR, 200000)              # photon path 0..z*
cH_lo = C/(H0*np.array([B.E(z) for z in zlo]))
DM_lo = cumulative_trapezoid(cH_lo, zlo, initial=0.0)   # = D_M(z), ends at DMstar

print("\n  screen profile delta psi(z) = A*Omega_Lambda(z):  (A=0.0667 shown)")
for z in (0.0,0.5,1.0,2.0,5.0,100.0,ZSTAR):
    print(f"     z={z:8.1f}:  Omega_L(z)={shape(z):.3e}   ->  screen at the RULER epoch is ~0")

# ================================================================= CASE 1
# Pure conformal single-sector: SAME e^{psi} scales BOTH the r_s integrand (z>z*) AND the D_M integrand (z<z*).
# A universal screen must give delta theta* = 0. We apply the SAME profile to both sectors over ALL z.
print("\n" + "-"*86)
print("CASE 1 — pure conformal single-sector (screen scales BOTH ruler and distance). Expect d theta*=0.")
def case1(A):
    # D_M with screen on the photon path
    DM_s = np.trapezoid(np.exp(A*shape(zlo))*cH_lo, zlo)
    # r_s with the SAME screen on the ruler integrand (z>z*)
    rs_s = quad(lambda z: np.exp(A*shape(z))*C*B._cs_over_c(z)/B.Hofz(z), ZSTAR, 1e6, limit=400)[0]
    return 100*rs_s/DM_s
for A in (0.0, 0.0667, 0.20):
    th = case1(A)
    print(f"     A={A:6.4f}:  100 theta*={th:.5f}   (shift {(th-theta_raw)/theta_raw*100:+.4f}%  "
          f"= {(th-theta_raw)/ (THETA_ERR):+.2f} sigma)")
print("     -> the screen at z>z* is ~0, so r_s is untouched; the shift comes ONLY from the low-z D_M leg.")
print("        A *universal* zoom (constant in z) would cancel exactly; this z-shaped DE screen does NOT")
print("        touch the ruler, so it is already an optical/material split, NOT a universal conformal zoom.")

# Demonstrate the true universal-zoom cancellation explicitly (constant screen):
def case1_const(c):  # constant conformal factor everywhere
    return 100*(np.exp(c)*rs_star)/(np.exp(c)*DMstar)
print(f"     CONSTANT zoom check: c=0.2 -> 100 theta*={case1_const(0.2):.5f} (identical to raw {theta_raw:.5f}) "
      f"=> universal zoom cancels EXACTLY, as the no-go requires.")

# ================================================================= CASE 2
# Optical-only propagation screen: photons see D_M^opt, r_s physical. (= the kernel split.)
print("\n" + "-"*86)
print("CASE 2 — optical-only propagation screen on D_M (r_s physical). Legit iff DFD screens light vs rulers.")
def theta_case2(A):
    DM_s = np.trapezoid(np.exp(A*shape(zlo))*cH_lo, zlo)
    return 100*rs_star/DM_s
A2 = brentq(lambda A: theta_case2(A)-THETA_PL, -0.5, 0.5)
print(f"     amplitude to PIN Planck theta*:  A={A2:.4f}   -> 100 theta*={theta_case2(A2):.5f}")
print(f"     derived gate value 2 sqrt(alpha) R_d/(1+R_d) = {2*np.sqrt(ALPHA)*B._R_baryon(ZDRAG)/(1+B._R_baryon(ZDRAG)):.4f}")

# ================================================================= CASE 3
# Baryon-trace acoustic-ruler correction: which WEIGHTING of f_b? endpoint vs sound-horizon vs visibility.
print("\n" + "-"*86)
print("CASE 3 — baryon-trace ruler correction: the f_b WEIGHTING question (external review).")
Rd = B._R_baryon(ZDRAG); Rstar = B._R_baryon(ZSTAR)
fb_drag = Rd/(1+Rd); fb_star = Rstar/(1+Rstar)
# sound-horizon-integral weighted <Pi_b>_{r_s}: weight by the r_s integrand c_s/H over z>z*
def w_rs(z): return C*B._cs_over_c(z)/B.Hofz(z)
def Pi(z): R=B._R_baryon(z); return R/(1+R)
num = quad(lambda z: Pi(z)*w_rs(z), ZSTAR, 1e6, limit=400)[0]
den = quad(lambda z: w_rs(z),       ZSTAR, 1e6, limit=400)[0]
fb_rs_weighted = num/den
# visibility-ish weighting: recombination is sharply peaked near z*, so ~ endpoint value
print(f"     Pi_b(z_drag)           = {fb_drag:.4f}   (endpoint, drag surface) <- used in A=2 sqrt(a) f_b")
print(f"     Pi_b(z_star)           = {fb_star:.4f}   (endpoint, recombination)")
print(f"     <Pi_b>_{{r_s integral}}  = {fb_rs_weighted:.4f}   (sound-horizon-weighted over z>z*)")
print(f"     -> the weighting MATTERS: {fb_drag:.3f} vs {fb_rs_weighted:.3f} "
      f"({100*(fb_rs_weighted/fb_drag-1):+.0f}% ). The drag-surface value is a CHOICE that must be justified;")
print(f"        the honest gate is bracketed [{min(fb_rs_weighted,fb_drag):.3f}, {max(fb_rs_weighted,fb_drag):.3f}].")

# ================================================================= CASE 4
# Correct two-sector model: theta* = r_s^{trace-corrected} / D_M^{optical}. Frozen amplitude, check theta*+BAO.
print("\n" + "-"*86)
print("CASE 4 — consistent TWO-SECTOR model: D_M optically screened + r_s baryon-trace corrected. FROZEN.")
A_frozen = 2*np.sqrt(ALPHA)*fb_drag
def two_sector(A, fb_weight):
    DM_s = np.trapezoid(np.exp(A*shape(zlo))*cH_lo, zlo)            # optical screen on photon path
    # baryon-trace ruler correction: r_s -> r_s*(1 + fb_weight*delta), delta from the trace coupling at the ruler.
    # screen at z>z* is ~0, so the ruler correction is tiny -> r_s ~ physical (this is the key point).
    rs_corr = quad(lambda z: np.exp(fb_weight*A*shape(z))*C*B._cs_over_c(z)/B.Hofz(z), ZSTAR, 1e6, limit=400)[0]
    return 100*rs_corr/DM_s, rs_corr
th4, rs4 = two_sector(A_frozen, fb_drag)
print(f"     frozen A = 2 sqrt(alpha) f_b = {A_frozen:.4f}")
print(f"     r_s^corrected = {rs4:.3f} Mpc (vs physical {rs_star:.3f}: trace correction = "
      f"{100*(rs4/rs_star-1):+.4f}%, negligible b/c screen~0 at z>z*)")
print(f"     100 theta*(two-sector, frozen) = {th4:.5f}   "
      f"(Planck {THETA_PL}, {(th4/100-THETA_PL/100)/(THETA_ERR/100):+.2f} sigma)")

# BAO under the same frozen two-sector screen (kernel split: transverse integrated, radial local)
def bao_chi2_screen(A):
    DM_opt = cumulative_trapezoid(np.exp(A*shape(zlo))*cH_lo, zlo, initial=0.0)
    chi2 = 0.0
    for b in B.BAO:
        z = b["z"]
        DMo = np.interp(z, zlo, DM_opt)
        DHo = np.exp(A*shape(z))*C/B.Hofz(z)
        pred = DMo/rd if b["kind"]=="DM" else DHo/rd if b["kind"]=="DH" else (z*DMo**2*DHo)**(1/3)/rd
        chi2 += ((pred-b["val"])/b["err"])**2
    return chi2
print(f"     BAO chi^2 (frozen two-sector, 12 pts) = {bao_chi2_screen(A_frozen):.2f}  "
      f"(chi^2/N={bao_chi2_screen(A_frozen)/12:.2f})")
print(f"     theta* chi^2 = {((th4/100-THETA_PL/100)/(THETA_ERR/100))**2:.2f}")

print("\n" + "="*86)
print("READ: Case 1 shows a universal zoom cancels (no-go holds) BUT the DFD DE screen is z-shaped and")
print("vanishes at the ruler epoch -> it is already optical-only. Case 2 pins theta* with the screen on")
print("the photon path. Case 3 shows the f_b weighting is a real open choice (drag vs integral). Case 4 is")
print("the consistent two-sector model; if its FROZEN theta*+BAO chi^2 is good, the cosmology fit SURVIVES")
print("the correct physics and was NOT a bookkeeping artifact.")
