There are two ways we will compute CMB spectra in this class: pypico and Class. Class directly computes the spectral densities directly (and therefore takes some time to compute). Pypico, in contrast, fits a functional interpolation scheme to precomputed spectra for fast computation (useful when running MCMC on the LCDM parameters).
This notebook will show how to use pypico and use it to explore how the tempurature spectral density changes with respect to LCDM parameters.
The first thing you need to do to get pypico running is to follow the instructions for installation in python (here is a link).
Part of this installation requires you to download the data interpolation files for pypico. You can get the files by following the link here. Once this data file is on your computer change the following variable to point to your downloaded file.
In [1]:
picodata_path = "/Users/ethananderes/Dropbox/Courses/STA250CMB/data/pypico/pico3_tailmonty_v34.dat"
Out[1]:
Now we can load the data file into pypico using PyCall
In [2]:
using PyCall
@pyimport pypico
const picoload = pypico.load_pico(picodata_path)
# -------- wrap pico
function pico(;
omega_b = 0.0224567,
omega_cdm = 0.118489,
tau_reio = 0.128312,
theta_s = 0.0104098,
logA_s_1010 = 3.29056,
n_s = 0.968602
)
plout::Dict{ASCIIString, Array{Float64,1}} = picoload[:get](;
:re_optical_depth => tau_reio,
symbol("scalar_amp(1)") => 1e-10*exp(logA_s_1010),
:theta => theta_s,
:ombh2 => omega_b,
:omch2 => omega_cdm,
symbol("scalar_spectral_index(1)") => n_s,
:massive_neutrinos => 3.04,
:helium_fraction => 0.25,
:omnuh2 => 0.0,
symbol("scalar_nrun(1)") => 0.0,
:force => true
)
return plout
end
Out[2]:
In [3]:
cls = pico()
Out[3]:
The output of pypico is a dictionary with a bunch of spectral densities. Here are the corresponding ell values for the TT spectral density.
In [5]:
ells = collect(0:length(cls["cl_TT"])-1)
Out[5]:
In [6]:
ellsc = 2π ./ ells ./ (ells + 1) |> x -> (x[1] = 0; x);
rawClTT = cls["cl_TT"] .* ellsc
Out[6]:
In [11]:
using PyPlot
plot(ells, rawClTT, linewidth = 4) #<---plotting the raw cls looks like a spike
Out[11]:
In [14]:
semilogy(ells, rawClTT, linewidth = 2)
Out[14]:
In [15]:
plot(ells, ells.*(ells+1).* rawClTT ./ (2π), linewidth = 2)
Out[15]:
In [16]:
semilogx(ells, ells.*(ells+1).* rawClTT ./ (2π), linewidth = 2)
Out[16]:
In [33]:
# Omega matter: omega_b ≈ 0.0223
μ = 0.0223
σ = 0.00014
l_max = 3000
for θ in linspace(μ-10σ, μ+10σ, 4)
cl = pico(omega_b = θ)["cl_TT"]
clμ = pico(omega_b = μ)["cl_TT"]
subplot(2,1,1)
plot(ells[1:l_max], cl[1:l_max], label = "theta = $θ")
subplot(2,1,2)
plot(ells[1:l_max], (cl-clμ)[1:l_max])
end
subplot(2,1,1)
title("omega_b cl response")
legend()
subplot(2,1,2)
title("delta cl = cl - clmu")
Out[33]:
In [34]:
# Omega dark matter: omega_cdm ≈ 0.118489
μ = 0.118489
σ = 0.00014
l_max = 3000
for θ in linspace(μ-20σ, μ+20σ, 4)
cl = pico(omega_cdm = θ)["cl_TT"]
clμ = pico(omega_cdm = μ)["cl_TT"]
subplot(2,1,1)
plot(ells[1:l_max], cl[1:l_max], label = "theta = $θ")
subplot(2,1,2)
plot(ells[1:l_max], (cl-clμ)[1:l_max])
end
subplot(2,1,1)
title("omega_cdm cl response")
legend()
subplot(2,1,2)
title("delta cl = cl - clmu")
Out[34]:
In [35]:
# tau_reio ≈ 0.128312
μ = 0.118489
σ = 0.00014
l_max = 3000
for θ in linspace(μ-20σ, μ+20σ, 4)
cl = pico(tau_reio = θ)["cl_TT"]
clμ = pico(tau_reio = μ)["cl_TT"]
subplot(2,1,1)
plot(ells[1:l_max], cl[1:l_max], label = "theta = $θ")
subplot(2,1,2)
plot(ells[1:l_max], (cl-clμ)[1:l_max])
end
subplot(2,1,1)
title("tau_reio cl response")
legend()
subplot(2,1,2)
title("delta cl = cl - clmu")
Out[35]:
In [39]:
# theta_s ≈ 0.0104098,
μ = 0.0104098
σ = 0.0000014
l_max = 3000
for θ in linspace(μ-15σ, μ+15σ, 4)
cl = pico(theta_s = θ)["cl_TT"]
clμ = pico(theta_s = μ)["cl_TT"]
subplot(2,1,1)
plot(ells[1:l_max], cl[1:l_max], label = "theta = $θ")
subplot(2,1,2)
plot(ells[1:l_max], (cl-clμ)[1:l_max])
end
subplot(2,1,1)
title("theta_s cl response")
legend()
subplot(2,1,2)
title("delta cl = cl - clmu")
Out[39]:
In [46]:
# logA_s_1010 ≈ 3.29056
μ = 3.29056
σ = 0.0014
l_max = 3000
for θ in linspace(μ-10σ, μ+10σ, 4)
cl = pico(logA_s_1010 = θ)["cl_TT"]
clμ = pico(logA_s_1010 = μ)["cl_TT"]
subplot(2,1,1)
plot(ells[1:l_max], cl[1:l_max], label = "theta = $θ")
subplot(2,1,2)
plot(ells[1:l_max], (cl-clμ)[1:l_max])
end
subplot(2,1,1)
title("logA_s_1010 cl response")
legend()
subplot(2,1,2)
title("delta cl = cl - clmu")
Out[46]:
In [45]:
# n_s ≈ 0.968602
μ = 0.968602
σ = 0.0062
l_max = 3000
for θ in linspace(μ-2σ, μ+2σ, 4)
cl = pico(n_s = θ)["cl_TT"]
clμ = pico(n_s = μ)["cl_TT"]
subplot(2,1,1)
plot(ells[1:l_max], cl[1:l_max], label = "theta = $θ")
subplot(2,1,2)
plot(ells[1:l_max], (cl-clμ)[1:l_max])
end
subplot(2,1,1)
title("n_s cl response")
legend()
subplot(2,1,2)
title("delta cl = cl - clmu")
Out[45]:
In [8]:
using PyPlot
plot((1:5626).^2.*cls["scalar_TT"]) # unlensed
plot((1:5626).^2.*cls["cl_TT"]) # lensed.
Out[8]:
In [ ]: