Use Filon quadrature for line-of-sight integration

[hersle]

using SymBoltz, Plots
M = ΛCDM()
pars = parameters_Planck18(M)
prob = CosmologyProblem(M, pars)

ls = 25:25:2500
jl = SphericalBesselCache(ls)
jlint = SphericalBesselIntegralCache(ls)

Dls1 = spectrum_cmb(:TT, prob, jl; normalization = :Dl, xs = nothing)
Dls2 = spectrum_cmb(:TT, prob, jlint; normalization = :Dl, xs = nothing)
plot(jl.l, Dls1)
plot!(jlint.l, Dls2)

[hersle]

Potential for time cut approximation:

k = 100.0
sol = solve(prob, k)
τs = sol.bg.t
Ss = sol(M.ST, τs, k)
xs = k .* (τs[end] .- τs)
plot(xs, 100Ss ./ maximum(abs.(Ss)); ylims = (-1, +1), label = "100S(x)", xlabel = "x")
il = 5
plot!(jl.x, jl.y[:,il] ./ maximum(abs.(jl.y[:,il])); xlims = extrema(xs), ylims = (-1, +1), label = "jₗ(x)")

[hersle]

Potential for multipole cut approximation (setting $\Delta_l(k)$ = 0 without integrating):

τ0 = τs[end]
ks = range(0.0, 1*jl.l[end]/τ0, length = 100)[2:end] # skip 0
sol = solve(prob, ks)
Ss = sol(M.ST, τs, ks)
Ss = reshape(Ss, (1, size(Ss)...))
Δs = SymBoltz.los_integrate(Ss, jl.l, τs, ks, jl)
Δs = Δs[1, :, :]
lgΔs = max.(log10.(abs.(Δs)), -300)
heatmap(ls, ks.*τ0, lgΔs; xlabel = "l", ylabel = "k τ₀", title = "log10(abs(Δₗ(k)))")

Possible algorithm:

1. Start with highest $k$ and assume $\Delta_l(k) \neq 0$ for all $l$.
2. LOS integrate, but skip those $l$ for which $\Delta_l(k) \approx 0$ at the previous (higher) $k$-value.

How to multithread?