PRO READ_CAMB_CL, file, cl, lmax=lmax, polar=polar, cmb_outputscale=cmb_outputscale, cmb_units=cmb_units ;; Check number of columns GET_LUN, unit OPENR, unit, file var='' READF, unit, var CLOSE, unit FREE_LUN, unit var=STRTRIM(var,2) var=STRCOMPRESS(var) res=STRSPLIT(var) ncol=N_ELEMENTS(res) ;; Problem: the CMB_outputscale used in the call to CAMB is not ;; written in the output file ;; We now assume it is 1 by default, and that CAMB spectra are written ;; in K/KCMB. We check that C_T_T = cl.c_1_1 has the correct magnitude at low ell. ;; CMB_outputscale = t0^2 ;; t0=2.726*1.0d6 IF NOT KEYWORD_SET(cmb_outputscale) THEN BEGIN STOP scaling = 1 ENDIF ELSE BEGIN scaling = cmb_outputscale ENDELSE IF NOT KEYWORD_SET(polar) THEN polar=0 IF polar EQ 1 THEN BEGIN ;; If 6 columns IF ncol EQ 6 THEN BEGIN READCOL, file, ell, ctt, cee, cte, cpp, ctp IF KEYWORD_SET(lmax) THEN BEGIN newlmax=MIN([MAX(ell),lmax]) ENDIF ELSE BEGIN newlmax=MAX(ell) ENDELSE newell=ell[where(ell le newlmax)] newind=LINDGEN(N_ELEMENTS(newell)) cl = CL_STRUCTURE(fieldnames=['T','E','B','P'], lmin=0, lmax=newlmax, /double, units=['K/KCMB','K/KCMB','K/KCMB','unitless']) cl.l[newell] = newell cl.c_1_1[newell] = 2.d0*!pi*ctt[newind]/newell/(newell+1L)/scaling cl.c_1_2[newell] = 2.d0*!pi*cte[newind]/newell/(newell+1L)/scaling cl.c_2_2[newell] = 2.d0*!pi*cee[newind]/newell/(newell+1L)/scaling cl.c_1_4[newell] = ctp[newind]/newell/newell/newell/scaling cl.c_4_4[newell] = cpp[newind]/newell/newell/newell/newell/scaling ENDIF ;; If 5 columns IF ncol EQ 5 THEN BEGIN READCOL, file, ell, ctt, cee, cbb, cte IF KEYWORD_SET(lmax) THEN BEGIN newlmax=MIN([MAX(ell),lmax]) ENDIF ELSE BEGIN newlmax=MAX(ell) ENDELSE newell=ell[where(ell le newlmax)] newind=LINDGEN(N_ELEMENTS(newell)) cl = CL_STRUCTURE(fieldnames=['T','E','B'], lmin=0, lmax=newlmax, /double, units=['K/KCMB','K/KCMB','K/KCMB']) cl.l[newell] = newell cl.c_1_1[newell] = 2.d0*!pi*ctt[newind]/newell/(newell+1L)/scaling cl.c_1_2[newell] = 2.d0*!pi*cte[newind]/newell/(newell+1L)/scaling cl.c_2_2[newell] = 2.d0*!pi*cee[newind]/newell/(newell+1L)/scaling cl.c_3_3[newell] = 2.d0*!pi*cbb[newind]/newell/(newell+1L)/scaling ENDIF ENDIF ELSE BEGIN ;; If 6 columns IF ncol EQ 6 THEN BEGIN READCOL, file, ell, ctt, cee, cte, cpp, ctp IF KEYWORD_SET(lmax) THEN BEGIN newlmax=MIN([MAX(ell),lmax]) ENDIF ELSE BEGIN newlmax=MAX(ell) ENDELSE newell=ell[where(ell le newlmax)] newind=LINDGEN(N_ELEMENTS(newell)) cl = CL_STRUCTURE(fieldnames=['T','P'], lmin=0, lmax=newlmax, /double, units=['K/KCMB','unitless']) cl.l[newell] = newell cl.c_1_1[newell] = 2.d0*!pi*ctt[newind]/newell/(newell+1L)/scaling cl.c_1_2[newell] = ctp[newind]/newell/newell/newell/scaling cl.c_2_2[newell] = cpp[newind]/newell/newell/newell/newell/scaling ENDIF ;; If 5 columns IF ncol EQ 5 THEN BEGIN READCOL, file, ell, ctt, cee, cbb, cte IF KEYWORD_SET(lmax) THEN BEGIN newlmax=MIN([MAX(ell),lmax]) ENDIF ELSE BEGIN newlmax=MAX(ell) ENDELSE newell=ell[where(ell le newlmax)] newind=LINDGEN(N_ELEMENTS(newell)) cl = CL_STRUCTURE(fieldnames=['T'], lmin=0, lmax=newlmax, /double, units=['K/KCMB']) cl.l[newell] = newell cl.c_1_1[newell] = 2.d0*!pi*ctt[newind]/newell/(newell+1L)/scaling ENDIF ENDELSE IF KEYWORD_SET(cmb_units) NE 0 THEN BEGIN CL_UPDATE_UNITS, cl, 'T', cmb_units CL_UPDATE_UNITS, cl, 'E', cmb_units CL_UPDATE_UNITS, cl, 'B', cmb_units ENDIF END