diff --git a/examples/ttteee_sfint_conv.jl b/examples/ttteee_sfint_conv.jl index b03ac7b..59700b9 100644 --- a/examples/ttteee_sfint_conv.jl +++ b/examples/ttteee_sfint_conv.jl @@ -10,7 +10,7 @@ camb_Cᵀᵀ = readdlm(prefix*"unlens_cellTT.dat")[3:ℓ_stride:end,2] #start at camb_Cᵀᴱ = readdlm(prefix*"unlens_cellTE.dat")[3:ℓ_stride:end,2] camb_Cᴱᴱ = readdlm(prefix*"unlens_cellEE.dat")[3:ℓ_stride:end,2] camb_ℓs = readdlm(prefix*"unlens_cellTT.dat")[3:ℓ_stride:end,1] -ℓs = Int.(camb_ℓs); +ℓs = Int.(camb_ℓs); #FIXME silent error issues if Float64 # global ℓ range # ℓmin,ℓmax,nℓ = 2,20,1200 @@ -19,14 +19,10 @@ camb_ℓs = readdlm(prefix*"unlens_cellTT.dat")[3:ℓ_stride:end,1] # subsample the ℓs because 2401 is too many # (TODO check how camb gets so many ℓs??) - function clttteee(ℓs, 𝕡, bg, ih, sf,sf_P; Dℓ=true) Cᵀᵀ = cltt(ℓs, 𝕡, bg, ih, sf) - println("done TT") Cᵀᴱ = clte(ℓs, 𝕡, bg, ih, sf,sf_P) - println("done TE") Cᴱᴱ = clee(ℓs, 𝕡, bg, ih, sf_P) - println("done EE") Tγ = (15/ π^2 *Bolt.ρ_crit(𝕡) * 𝕡.Ω_r)^(1/4) norm_fac = ( (Tγ * Bolt.Kelvin_natural_unit_conversion * 1e6)^2 / (2π) * 4^2) if Dℓ @@ -45,47 +41,21 @@ function setup_sf(nx,nk) xmin,xmax = -20.0,0.0 dx = (xmax-xmin)/nx bg = Background(𝕡; x_grid=xmin:dx:xmax) - println("done bg") 𝕣 = Bolt.RECFAST(bg=bg, Yp=𝕡.Y_p, OmegaB=𝕡.Ω_b) ih = IonizationHistory(𝕣, 𝕡, bg) - println("done ih") kmin,kmax = 0.1bg.H₀, 1000bg.H₀ #WARNING fd issues with these 𝕡-dep. values ks = quadratic_k(kmin,kmax,nk) sf = source_grid(𝕡, bg, ih, ks, BasicNewtonian()) - println("done sf") sf_P = source_grid_P(𝕡, bg, ih, ks, BasicNewtonian()) - println("done sf_P") return 𝕡, bg, ih, sf, sf_P end -testsf_𝕡,testsf_bg,testsf_ih,testsf_sf,testsf_sf_P = setup_sf(200,100); - -tttest = cltt(ℓs, testsf_𝕡,testsf_bg,testsf_ih,testsf_sf) -#^This is taking at least several minutes...why so much slower than basic_usage? - -# Let's instead try the cltt call with the ℓs from "basic_usage" -# those should run in a minute or less if nothing is very wrong -bu_ℓs = 2:20:1200 -bu_tttest = cltt(bu_ℓs, testsf_𝕡,testsf_bg,testsf_ih,testsf_sf) -ℓs -bu_ℓs - - -tetest = clte(ℓs, testsf_𝕡,testsf_bg,testsf_ih,testsf_sf,testsf_sf_P) - -eetest = clee(ℓs, testsf_𝕡,testsf_bg,testsf_ih,testsf_sf_P) - -test_Cᵀᵀ, test_Cᵀᴱ, test_Cᴱᴱ = clttteee(ℓs, testsf_𝕡,testsf_bg,testsf_ih,testsf_sf,testsf_sf_P) - - - function single_experiment(nx,nk) # setup 𝕡, bg, ih, sf, sf_P = setup_sf(nx,nk) # compute Cℓs Cᵀᵀ, Cᵀᴱ, Cᴱᴱ = clttteee(ℓs, 𝕡, bg, ih, sf,sf_P) - println("done clttteee") # compare to CAMB rᵀᵀ,rᵀᴱ,rᴱᴱ = Cᵀᵀ.-camb_Cᵀᵀ, Cᵀᴱ.-camb_Cᵀᴱ, Cᴱᴱ.-camb_Cᴱᴱ @@ -94,11 +64,8 @@ function single_experiment(nx,nk) end # CMB Cᵀᵀ(ℓ) - function plot_experiments(rrs,nxs,nks,percent=false,sv=["TT","TE","EE"]) - # p = plot(ℓs, @.(ℓs^2*camb_Cᵀᵀ), label="CAMB",color=:black,ls=:solid) p = hline([0.0],color=:black,ls=:solid,label=false) - for (i,rr) in enumerate(rrs) for s in sv if percent @@ -112,11 +79,9 @@ function plot_experiments(rrs,nxs,nks,percent=false,sv=["TT","TE","EE"]) error("Invalid spectrum") end label_use = s*" - nx = $(nxs[i]), nk = $(nks[i])" - # label_use = i == 1 ? s*" - nx = $(nxs[i]), nk = $(nks[i])" : "" plot!(ℓs, rr[1]./camb_C, label=label_use,ls=:solid,color=i) else label_use = s*" - nx = $(nxs[i]), nk = $(nks[i])" - # label_use = i == 1 ? s*" - nx = $(nxs[i]), nk = $(nks[i])" : "" plot!(ℓs, rr[1], label=label_use,ls=:solid,color=i) end end @@ -135,9 +100,6 @@ end nxx = [100,200,300,400] nkk = [50,100,150,200] -# test_exp = single_experiment(200,10) -# plot_experiments([test_exp],[200],[10]) - # Structured runs # the scheme is a list with [(nx_1,nk_1), (nx_1,nk_2),...,(nx_2,nk_1),...] rrs = [single_experiment(nx,nk) @@ -159,25 +121,5 @@ savefig(p,save_path*"rel_clee_diffs.png") save_path = "../../plots/" writedlm(save_path*"clttteee_diffs.dat",rrs) -#---------------------- -# Noodling from before: -# Tγtest = (15/ π^2 *Bolt.ρ_crit(testsf_𝕡) *testsf_𝕡.Ω_r)^(1/4) -# MpctoμK = Tγtest * Bolt.Kelvin_natural_unit_conversion * 1e6 -# plot(ℓs, (ℓs.*(ℓs.+1) .* tttest .* MpctoμK^2 ) ./ camb_Cᵀᵀ) -# plot(ℓs, (ℓs.*(ℓs.+1) .* tttest .* MpctoμK^2 ./ (2*π) .* 4^2) ./ ( camb_Cᵀᵀ)) -# # Here the factors are : -# # 1. The conversion from Mpc to μK in the temperature normalization -# # 2. The ℓ factor and factor of 2π from the definition of Dℓ -# # 3. The factor of 4 for each transfer function since we don't use the MB normalization of CLASS (which agrees with CAMB) -# hline!([1.0],color=:black,ls=:solid) - -# plot(ℓs, (ℓs.*(ℓs.+1) .* tetest .* MpctoμK^2 ) ./ camb_Cᵀᴱ) -# plot!(ℓs, (ℓs.*(ℓs.+1) .* tetest .* MpctoμK^2 ./ (2*π) .* 4^2) ./ ( camb_Cᵀᴱ)) -# hline!([1.0],color=:black,ls=:solid) -# ylims!(1-0.5,1+0.5) - -# plot(ℓs, (ℓs.*(ℓs.+1) .* eetest .* MpctoμK^2 ) ./ camb_Cᴱᴱ) -# plot!(ℓs, (ℓs.*(ℓs.+1) .* eetest .* MpctoμK^2 ./ (2*π) .* 4^2) ./ ( camb_Cᴱᴱ)) -# hline!([1.0],color=:black,ls=:solid) -# ylims!(1-0.5,1+0.5) +