import camb, numpy as np
def get(omch2):
    p=camb.CAMBparams(); p.set_cosmology(H0=67.4,ombh2=0.0224,omch2=omch2,mnu=0.06,tau=0.054)
    p.InitPower.set_params(As=2.1e-9,ns=0.965); p.set_for_lmax(2600)
    r=camb.get_results(p); cl=r.get_cmb_power_spectra(p,CMB_unit='muK',spectra=['total'])['total'][:,0]
    return cl
def peaks(TT,lo=120,hi=2500):
    P=[]
    for l in range(lo,hi):
        if TT[l]>TT[l-1] and TT[l]>TT[l+1] and TT[l]>TT[l-8] and TT[l]>TT[l+8]:
            if not P or l-P[-1]>80: P.append(l)
    return P
TT=get(1e-5)                       # baryon-only
pk=peaks(TT); H=np.array([TT[l] for l in pk]); 
print("baryon-only peaks l=",pk[:6])
print("  H/H1 =",np.round(H[:6]/H[0],3))
# fit Silk scale l_D from the damping tail: ln(H_n) ~ const - (l_n/l_D)^2 over peaks 3..6
import numpy as np
ls=np.array(pk[2:6]); hs=np.log(np.array([TT[l] for l in ls]))
# assume undamped envelope ~ const*l^p ; fit ln H = a + b*ln l - (l/lD)^2 ; simpler: fit ln H vs l^2 on peaks 3-6
A=np.vstack([np.ones_like(ls,dtype=float), ls.astype(float)**2]).T
coef,*_=np.linalg.lstsq(A, hs, rcond=None); lD=np.sqrt(-1/coef[1])
print(f"  fitted Silk scale l_D = {lD:.0f}")
# UNDAMP: multiply peak heights by exp(+(l/lD)^2), renormalize to peak1
undamp=lambda l,h: h*np.exp((l/lD)**2)
Hu=np.array([undamp(pk[i],H[i]) for i in range(6)]); Hu=Hu/Hu[0]
print("  UNDAMPED (Silk removed) H/H1 =",np.round(Hu,3))
print(f"  => undamped peak3/peak1 = {Hu[2]:.3f}   (Planck target 0.44 ; baryon-only damped 0.205)")
print()
print("CONTROL: at FIXED Omega_b (fixed Silk), does peak3 track z_eq/driving?")
for om in [0.0,0.04,0.08,0.12]:
    T=get(max(om,1e-5)); p=peaks(T); h=np.array([T[l] for l in p]); h=h/h[0]
    print(f"  omch2={om:.2f}: peaks{p[:3]} H/H1={np.round(h[:3],3)}  (Silk ~unchanged; peak3 set by driving)")
