Add nonzero massless neutrino initial conditions

src/components/species.jl

@@ -130,7 +130,7 @@ function photons(g; polarization = true, lmax = 6, name = :γ, kwargs...)
         F[0] ~ -2*g.Ψ # Dodelson (7.89) # TODO: derive automatically
         F[1] ~ 2//3 * k/g.ℰ*g.Ψ # Dodelson (7.95)
         F[2] ~ (polarization ? -8//15 : -20//45) * k/τ̇ * F[1] # depends on whether polarization is included
-        [F[l] ~ -l/(2*l+1) * k/τ̇ * F[l-1] for l in 3:lmax]...
+        [F[l] ~ -l//(2*l+1) * k/τ̇ * F[l-1] for l in 3:lmax]...
     ] .|> O(ϵ^1)
     if polarization
         append!(eqs1, [
@@ -142,7 +142,7 @@ function photons(g; polarization = true, lmax = 6, name = :γ, kwargs...)
             G[0] ~ 5//16 * F[2],
             G[1] ~ -1//16 * k/τ̇ * F[2],
             G[2] ~ 1//16 * F[2],
-            [G[l] ~ -l/(2*l+1) * k/τ̇ * G[l-1] for l in 3:lmax]...
+            [G[l] ~ -l//(2*l+1) * k/τ̇ * G[l-1] for l in 3:lmax]...
         ] .|> O(ϵ^1))
     else
         append!(eqs1, [collect(G .~ 0)...] .|> O(ϵ^1)) # pin to zero
@@ -176,7 +176,7 @@ function massless_neutrinos(g; lmax = 6, name = :ν, kwargs...)
         δ ~ -2 * g.Ψ # adiabatic: δᵢ/(1+wᵢ) == δⱼ/(1+wⱼ) (https://cmb.wintherscoming.no/theory_initial.php#adiabatic)
         θ ~ 1//2 * (k^2/g.ℰ) * g.Ψ
         σ ~ 1//15 * (k/g.ℰ)^2 * g.Ψ
-        [F[l] ~ 0 #=1/(2*l+1) * k/g.ℰ * Θ[l-1]=# for l in 3:lmax]... # TODO: nonzero
+        [F[l] ~ +l//(2*l+1) * k/g.ℰ * F[l-1] for l in 3:lmax]...
     ] .|> O(ϵ^1)
     description = "Massless neutrinos"
     return extend(ν, ODESystem(eqs1, t, vars, pars; initialization_eqs=ics1, name, kwargs...); description)