Interpolate in k-space using Chebyshev polynomials

[hersle]

using SymBoltz, Plots
M = ΛCDM()
pars = parameters_Planck18(M)
prob = CosmologyProblem(M, pars)
ls = 25:25:2500
ls_all = ls[begin]:ls[end]
jl = SphericalBesselCache(ls)

plot()
for korder in 20:10:100
    Dls = spectrum_cmb(:TT, prob, jl, ls_all; normalization = :Dl, kτ0s = 0.1*jl.l[begin]:2π/2:2*jl.l[end], korder)
    plot!(ls_all, Dls; label = "order = $korder")
end
plot!()

[hersle]

using SymBoltz, Plots, StaticArrays
Plots.default(dpi = 300)
M = ΛCDM()
pars = parameters_Planck18(M)
prob = CosmologyProblem(M, pars)
sol = solve(prob)

S = M.ST

τs = sol.bg.t[begin:end-1]
ks_raw = 1.0:1.0:1000.0
Ss_raw = source_grid(prob, SVector(S), τs, ks_raw, sol.bg)
Ss_raw = stack(Ss_raw)[1, :, :]

korder = 30
Ss_cheb, _, _ = source_grid_chebyshev(prob, SVector(S), τs, ks_raw, sol.bg; order = korder)
Ss_cheb = Ss_cheb[1, :, :]

abserr = abs.(Ss_cheb .- Ss_raw)
relerr = abserr ./ abs.(Ss_raw)

xlims = (0.00, 0.20)
p1 = surface(τs, ks_raw, transpose(Ss_raw); xlabel = "τ", ylabel = "k", title = "S(τ, k), direct integration of $(length(ks_raw))×k", xlims , size = (1000, 600))
p2 = surface(τs, ks_raw, transpose(Ss_cheb); xlabel = "τ", ylabel = "k", title = "S(τ, k), interpolation from $(korder)-order Chebyshev", xlims, size = (1000, 600))
p3 = surface(τs, ks_raw, transpose(abserr); xlabel = "τ", ylabel = "k", title = "absolute error", xlims, size = (1000, 600))

plot(p1, p2, p3, layout = (1, 3), size = (2000, 400))

[hersle]

Convergence of maximum absolute error of $S^T(\tau, k)$ for Chebyshev interpolation vs. Cubic splines with uniform $k$-points: