import camb, numpy as np
def peaks(omch2,label):
    pars=camb.CAMBparams()
    pars.set_cosmology(H0=67.4,ombh2=0.0224,omch2=omch2,mnu=0.06,tau=0.054)
    pars.InitPower.set_params(As=2.1e-9,ns=0.965)
    pars.set_for_lmax(2500,lens_potential_accuracy=0)
    r=camb.get_results(pars)
    cl=r.get_cmb_power_spectra(pars,CMB_unit='muK',spectra=['total'])['total']
    ell=np.arange(cl.shape[0]); TT=cl[:,0]
    # find peaks in TT (Dl) for l>100
    pk=[]
    for l in range(120,2200):
        if TT[l]>TT[l-1] and TT[l]>TT[l+1] and TT[l]>TT[l-8] and TT[l]>TT[l+8]:
            pk.append((l,TT[l]))
    # dedupe close peaks
    P=[]
    for l,h in pk:
        if not P or l-P[-1][0]>80: P.append((l,h))
    h=np.array([p[1] for p in P[:5]]); h=h/h[0]
    locs=[p[0] for p in P[:5]]
    zeq=r.get_derived_params()['zstar']  # not zeq, but report zstar
    print(f"{label:18s} peaks l={locs}  H/H1={np.round(h,3)}  H3/H2={h[2]/h[1]:.3f}")
    return locs,h
print("CAMB 1.6.6 — TT peak heights (Dl=l(l+1)Cl/2pi)\n")
peaks(0.120,"LCDM (omch2=.12)")
peaks(0.055,"half-CDM")
peaks(0.020,"quarter-CDM")
peaks(1e-5,"BARYON-ONLY")
print("\nPlanck observed: H/H1~[1,0.45,0.44], H3/H2~0.97")
