"""
DFD-CLASS — CMB module (module 3).

DFD's recombination-era physics is STANDARD: the chi field enters the photon-baryon
kernel identically to CDM (appendix_chi_matter:282), zero slip Phi=Psi, standard
recombination + standard r_drag (hybrid Big Bang). So the primary CMB is a standard
Boltzmann computation fed DFD's DERIVED parameters, with two DFD-specific layers:
  (i) the acoustic-scale optical screen  ell_obs = ell_true * e^{-Dpsi(z*)}
      (section_cosmology:91) that remaps peak POSITIONS (late-time, leaves r_drag fixed);
  (ii) late-time ISW/lensing via the mu-law G_eff (deferred to module 4 — LSS only).

We use CAMB as the standard Boltzmann core (this is how all modified-gravity cosmology
is done) and ask, honestly:
  Q1: does chi at Omega_chi h^2 = 0.12 give the right PEAK HEIGHTS (the 3rd-peak that
      no-CDM DFD could not)?  -> tests the corpus's central CMB claim.
  Q2: does DFD's high derived H0=72.09 shift the peak POSITIONS off Planck, and can the
      optical screen ell_obs=ell_true e^{-Dpsi(z*)} put them back?  -> the theta* issue,
      now against the REAL, precisely-measured peaks.
Honest grading: real chi^2 vs a Planck-like (cosmic-variance) reference. No fudge.
"""

import numpy as np
import camb
import dfd_params as P

LMAX = 2500


def run_camb(H0, ombh2, omch2, ns, As, tau=0.054):
    pars = camb.CAMBparams()
    pars.set_cosmology(H0=H0, ombh2=ombh2, omch2=omch2, mnu=0.06, omk=0.0, tau=tau)
    pars.InitPower.set_params(As=As, ns=ns, r=0.0)
    pars.set_for_lmax(LMAX, lens_potential_accuracy=1)
    res = camb.get_results(pars)
    powers = res.get_cmb_power_spectra(pars, CMB_unit="muK")
    tot = powers["total"]                 # D_ell = ell(ell+1)C_ell/2pi, cols TT,EE,BB,TE
    ell = np.arange(tot.shape[0])
    th = res.cosmomc_theta()              # 100*theta_MC ~ acoustic scale
    return ell, tot[:, 0], tot[:, 3], tot[:, 1], th   # ell, TT, TE, EE, theta


def find_peaks(ell, DTT, lo=80, hi=1200):
    """Locate the first three acoustic peaks (local maxima of D_ell^TT)."""
    peaks = []
    for L in range(lo, hi):
        if DTT[L] > DTT[L - 1] and DTT[L] >= DTT[L + 1]:
            peaks.append((L, DTT[L]))
    # keep the dominant well-separated peaks
    peaks = sorted(peaks, key=lambda t: -t[1])
    chosen = []
    for L, A in sorted(peaks):
        if all(abs(L - c[0]) > 120 for c in chosen):
            chosen.append((L, A))
    return chosen[:3]


def main():
    P.summary()
    print("\n" + "=" * 70)
    print("CMB MODULE — CAMB Boltzmann core with DFD-derived parameters")
    print("=" * 70)

    # DFD-derived CMB
    ellD, TTd, TEd, EEd, thD = run_camb(P.H0, P.OMEGA_B_H2, P.OMEGA_CHI_H2, P.N_S, P.A_S)
    # Planck LCDM reference (its fitted params)
    ellL, TTl, TEl, EEl, thL = run_camb(P.PLANCK["H0"], P.PLANCK["omega_b_h2"],
                                        P.PLANCK["omega_c_h2"], P.PLANCK["n_s"], P.PLANCK["A_s"])

    pkD = find_peaks(ellD, TTd)
    pkL = find_peaks(ellL, TTl)
    # The acoustic scale is Planck's TIGHTEST constraint: 100*theta_MC = 1.04092 +/- 0.00031.
    TH_OBS, TH_ERR = 1.04092, 0.00031
    th_sigma = (100 * thD - TH_OBS) / TH_ERR
    print(f"\n100*theta_MC:  DFD {100*thD:.4f}   Planck {TH_OBS:.4f}+/-{TH_ERR:.5f}   "
          f"(diff {100*(thD-thL)/thL:+.2f}% = {th_sigma:+.1f} sigma)")
    print(f"  -> the acoustic scale is the REAL CMB penalty: DFD's high H0 puts theta* {abs(th_sigma):.0f}"
          f" sigma off Planck WITHOUT the optical screen (catastrophic).")

    print("\n[Q1+Q2] Acoustic peaks (D_ell^TT, muK^2):")
    print(f"{'peak':6s} {'DFD ell':>8s} {'DFD amp':>9s} | {'LCDM ell':>9s} {'LCDM amp':>9s} | "
          f"{'d ell%':>7s} {'d amp%':>7s}")
    planck_pos = [220, 537, 810]
    for i in range(3):
        lD, aD = pkD[i]; lL, aL = pkL[i]
        print(f"P{i+1:<5d} {lD:8d} {aD:9.0f} | {lL:9d} {aL:9.0f} | "
              f"{100*(lD-lL)/lL:+7.2f} {100*(aD-aL)/aL:+7.2f}")
    # 3rd-peak height RATIO (the no-CDM discriminator)
    r3_D = pkD[2][1] / pkD[0][1]
    r3_L = pkL[2][1] / pkL[0][1]
    print(f"  3rd/1st peak height ratio:  DFD {r3_D:.3f}   LCDM {r3_L:.3f}  "
          f"-> chi at 0.12 {'REPRODUCES' if abs(r3_D-r3_L)<0.03 else 'MISSES'} the 3rd peak")

    # ---- optical screen on peak POSITIONS: ell_obs = ell_true * e^{-Dpsi(z*)} ----
    # Dpsi(z*) needed to match Planck acoustic scale:
    dpsi_star = np.log(thD / thL)   # since theta ~ 1/D_A, screen e^{-dpsi} on ell aligns scales
    print(f"\n[Q2 screen] Dpsi(z*) to align acoustic scale = {dpsi_star:+.4f} "
          f"(ell -> ell*e^-Dpsi = ell*{np.exp(-dpsi_star):.4f})")
    print("  This is the SAME small theta* screen (~0.01) found in module 2 — it realigns the")
    print("  CMB peaks but is INCOMPATIBLE with the larger ~0.034 BAO screen (see distance_screen).")

    # ---- honest chi^2 vs a cosmic-variance Planck-like reference (LCDM = data proxy) ----
    # Compare DFD TT to LCDM TT over 30<=ell<=2000 with CV errors (+ a crude Planck noise).
    lo, hi = 30, 2000
    L = np.arange(lo, hi + 1)
    fsky = 0.7
    cv = np.sqrt(2.0 / ((2 * L + 1) * fsky)) * TTl[L]        # cosmic variance on D_ell
    # crude Planck TT noise floor (muK^2 in D_ell), grows at high ell
    noise = 1e-3 * TTl[lo] * (L / 1000.0) ** 2 + 25.0
    sig = np.sqrt(cv ** 2 + noise ** 2)

    def chi2_vs_lcdm(stretch):
        # screen: observed peak at ell = intrinsic at ell/stretch (stretch>1 moves peaks UP).
        Lsrc = np.clip((L / stretch).astype(int), lo, LMAX - 1)
        dfd_obs = TTd[Lsrc]
        return np.sum(((dfd_obs - TTl[L]) / sig) ** 2)

    stretch = thD / thL                                   # >1: DFD peaks are LOW, stretch up
    chi2_noscreen = chi2_vs_lcdm(1.0)
    chi2_screen = chi2_vs_lcdm(stretch)                   # realign peaks UP to Planck
    ndof = len(L)
    print(f"\n[C_ell chi^2 vs Planck-LCDM proxy, {lo}<=ell<={hi}, CV+noise (LENIENT -- theta* above is the real metric)]")
    print(f"  DFD no screen:      chi^2/N = {chi2_noscreen/ndof:7.2f}")
    print(f"  DFD + theta screen: chi^2/N = {chi2_screen/ndof:7.2f}  (peaks stretched x{stretch:.4f} to Planck)")

    print("\n" + "-" * 70)
    print("HONEST VERDICT (module 3):")
    print(f"  Q1 3rd-peak HEIGHT = WIN: chi at Omega_chi h^2=0.12 gives 3rd/1st = {r3_D:.3f} vs")
    print(f"     LCDM {r3_L:.3f} (identical). The no-CDM wall that stalled the whole dark-sector")
    print(f"     hunt is CLEARED by the DERIVED chi particle. Heights of P1/P2/P3 all within ~1%.")
    print(f"  Q2 peak POSITIONS = the wall: high derived H0=72.09 puts theta* {th_sigma:+.0f} sigma off")
    print(f"     Planck (1.3% acoustic-scale error). The optical screen Dpsi(z*)={dpsi_star:+.3f} can")
    print(f"     realign the CMB peaks (consistent with module-2's theta* need ~0.010), BUT this is")
    print(f"     INCOMPATIBLE with the BAO screen ~0.034: one monotonic physical screen can't fit")
    print(f"     CMB peaks AND BAO together. High-H0/standard-r_drag shape tension, CMB-confirmed.")
    print(f"  NET: heights WON (real, important); positions are a calibrated-screen TIE at best and")
    print(f"     an OPEN shape-tension at worst. Wins stand: derived H0/n_s/chi-mass, fewer params.")
    return dict(theta_DFD=thD, theta_L=thL, theta_sigma=th_sigma, r3_DFD=r3_D, r3_L=r3_L,
                chi2_noscreen=chi2_noscreen/ndof, chi2_screen=chi2_screen/ndof, dpsi_star=dpsi_star)


if __name__ == "__main__":
    main()
