"""
FIGURE cmb_tt: CMB TT power spectrum D_l = l(l+1)C_l/2pi for DFD-derived params vs LCDM.

DFD params: H0=72.09, ombh2=0.02237, omch2=0.12, ns=0.96667, As=2.08e-9
LCDM (Planck 2018): from P.PLANCK.

Marks first 3 acoustic peaks (cmb.find_peaks), annotates 3rd/1st peak height ratio for each
and the 100*theta_star value. CONFIRM: theta* near Planck 1.04109e-2*100, 3rd/1st ratio
DFD ~0.444 vs LCDM ~0.443.
"""
import sys
sys.path.insert(0, ".")

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

import dfd_params as P
import cmb

OUT = "./figures/cmb_tt.png"

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

pkD = cmb.find_peaks(ellD, TTd)
pkL = cmb.find_peaks(ellL, TTl)

r3_D = pkD[2][1] / pkD[0][1]
r3_L = pkL[2][1] / pkL[0][1]

th_D = 100 * thD
th_L = 100 * thL
TH_PLANCK = 1.04109e-2 * 100  # 1.04109

print("=== cmb_tt verification ===")
print(f"100*theta*  DFD  = {th_D:.5f}")
print(f"100*theta*  LCDM = {th_L:.5f}")
print(f"100*theta*  Planck ref = {TH_PLANCK:.5f}")
for i in range(3):
    print(f"P{i+1}: DFD ell={pkD[i][0]:4d} amp={pkD[i][1]:8.1f} | "
          f"LCDM ell={pkL[i][0]:4d} amp={pkL[i][1]:8.1f}")
print(f"3rd/1st peak ratio: DFD = {r3_D:.4f}   LCDM = {r3_L:.4f}")

# ---- plot ----
fig, ax = plt.subplots(figsize=(10, 6.2))
lo, hi = 2, 2500
L = np.arange(lo, hi + 1)

c_dfd = "#c0392b"
c_lcdm = "#2c3e50"

ax.plot(L, TTl[L], color=c_lcdm, lw=1.8, label=r"$\Lambda$CDM (Planck 2018)")
ax.plot(L, TTd[L], color=c_dfd, lw=1.8, ls="-", label="DFD (derived params)")

# mark peaks
for i, (pk, col, mk) in enumerate([(pkD, c_dfd, "o"), (pkL, c_lcdm, "s")]):
    for j in range(3):
        l, a = pk[j]
        ax.plot(l, a, mk, color=col, ms=7, mec="white", mew=0.8, zorder=5)

# annotate peak numbers near DFD peaks
for j in range(3):
    l, a = pkD[j]
    ax.annotate(f"P{j+1}\n$\\ell$={l}", xy=(l, a), xytext=(l, a + 600),
                ha="center", va="bottom", fontsize=9, color=c_dfd)

ax.set_xlabel(r"Multipole moment $\ell$", fontsize=12)
ax.set_ylabel(r"$\mathcal{D}_\ell^{TT}=\ell(\ell+1)C_\ell^{TT}/2\pi\ \ [\mu\mathrm{K}^2]$", fontsize=12)
ax.set_title("DFD vs $\\Lambda$CDM CMB temperature power spectrum (CAMB Boltzmann core)", fontsize=13)
ax.set_xlim(0, 2500)
ax.set_ylim(0, max(TTl[L].max(), TTd[L].max()) * 1.15)
ax.grid(alpha=0.25)

# annotation box with key numbers
txt = (
    r"$100\,\theta_*$:" + "\n"
    f"   DFD   = {th_D:.5f}\n"
    f"   $\\Lambda$CDM = {th_L:.5f}\n"
    f"   Planck = {TH_PLANCK:.5f}\n"
    "\n"
    "3rd/1st peak height ratio:\n"
    f"   DFD   = {r3_D:.3f}\n"
    f"   $\\Lambda$CDM = {r3_L:.3f}"
)
ax.text(0.975, 0.95, txt, transform=ax.transAxes, ha="right", va="top",
        fontsize=10, family="monospace",
        bbox=dict(boxstyle="round,pad=0.5", fc="#fdf6e3", ec="#888", alpha=0.95))

ax.legend(loc="upper center", fontsize=11, frameon=True)

fig.tight_layout()
fig.savefig(OUT, dpi=150)
print(f"\nsaved {OUT}")
