Skip to content

Commit

Permalink
cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
jmsull committed Aug 8, 2024
1 parent ca20359 commit c37469e
Showing 1 changed file with 2 additions and 60 deletions.
62 changes: 2 additions & 60 deletions examples/ttteee_sfint_conv.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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")
= (15/ π^2 *Bolt.ρ_crit(𝕡) * 𝕡.Ω_r)^(1/4)
norm_fac = ( (Tγ * Bolt.Kelvin_natural_unit_conversion * 1e6)^2 / (2π) * 4^2)
if Dℓ
Expand All @@ -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ᴱᴱ
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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)


0 comments on commit c37469e

Please sign in to comment.