Rename variables; improve Latexified equations

docs/src/automatic_differentiation.md

@@ -17,7 +17,7 @@ using SymBoltz
 M = ΛCDM()
 
 function P(k, θ)
-   pars = Dict([M.γ.T0, M.c.Ω0, M.b.Ω0, M.ν.Neff, M.g.h, M.b.rec.Yp] .=> θ)
+   pars = Dict([M.γ.T₀, M.c.Ω₀, M.b.Ω₀, M.ν.Neff, M.g.h, M.b.rec.Yp] .=> θ)
    sol = solve(M, pars, k)
    return power_spectrum(sol, k)
 end

docs/src/comparison.md

@@ -58,20 +58,20 @@ function solve_class(pars, k; exec="class", inpath="/tmp/symboltz_class/input.in
         "h" => pars[M.g.h],
 
         # photons
-        "T_cmb" => pars[M.γ.T0],
+        "T_cmb" => pars[M.γ.T₀],
         "l_max_g" => lmax,
         "l_max_pol_g" => lmax,
 
         # baryons
-        "Omega_b" => pars[M.b.Ω0],
+        "Omega_b" => pars[M.b.Ω₀],
         "YHe" => pars[M.b.rec.Yp],
         "recombination" => "recfast", # TODO: HyREC
         "recfast_Hswitch" => 1,
         "recfast_Heswitch" => 6,
         "reio_parametrization" => "reio_none", # TODO: enable
 
         # cold dark matter
-        "Omega_cdm" => pars[M.c.Ω0],
+        "Omega_cdm" => pars[M.c.Ω₀],
 
         # neutrinos # TODO: set neutrino stuff to 0 unless otherwise specified
         "N_ur" => pars[M.ν.Neff],
@@ -137,7 +137,7 @@ sol = Dict(
     # thermodynamics
     "a_th" => (reverse(sol1["th"]["scalefactora"]), sol2[M.g.a]),
     "τ̇" => (reverse(sol1["th"]["kappa'[Mpc^-1]"]), .- sol2[M.b.rec.τ̇] * (SymBoltz.k0 * h)),
-    "csb²" => (reverse(sol1["th"]["c_b^2"]), sol2[M.b.rec.cs²]),
+    "csb²" => (reverse(sol1["th"]["c_b^2"]), sol2[M.b.rec.cₛ²]),
     "Xe" => (reverse(sol1["th"]["x_e"]), sol2[M.b.rec.Xe]),
     "Tb" => (reverse(sol1["th"]["Tb[K]"]), sol2[M.b.rec.Tb]),
     #"Tb′" => (reverse(sol1["th"]["dTb[K]"]), sol2[M.b.rec.DTb] ./ -sol2[M.g.E]), # convert my dT/dt̂ to CLASS' dT/dz = -1/H * dT/dt 

docs/src/extended_models.md

@@ -53,7 +53,7 @@ using SymBoltz: t, D, k, ϵ # load conformal time, derivative, perturbation wave
 g = M1.g # reuse metric component from the original model
 
 # 1. Create parameters (will resurface as tunable numbers in the full model)
-pars = @parameters w0 wa ρ0 Ω0 cs²
+pars = @parameters w0 wa ρ₀ Ω₀ cₛ²
 
 # 2. Create variables
 vars = @variables ρ(t) P(t) w(t) δ(t) θ(t) σ(t)
@@ -62,12 +62,12 @@ vars = @variables ρ(t) P(t) w(t) δ(t) θ(t) σ(t)
 eqs = [
     # Background equations (of order O(ϵ⁰)
     w ~ w0 + wa * (1 - g.a) # equation of state
-    ρ ~ ρ0 * abs(g.a)^(-3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - g.a)) # energy density # TODO: get rid of abs
+    ρ ~ ρ₀ * abs(g.a)^(-3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - g.a)) # energy density # TODO: get rid of abs
     P ~ w * ρ # pressure
 
     # Perturbation equations (mulitiplied by ϵ to mark them as order O(ϵ¹))
-    D(δ) * ϵ ~ (-(1 + w) * (θ - 3*g.Φ) - 3 * g.ℰ * (cs² - w) * δ) * ϵ # energy overdensity
-    D(θ) * ϵ ~ (-g.ℰ * (1 - 3*w) - D(w) / (1 + w) * θ + cs² / (1 + w) * k^2 * δ - k^2 * σ + k^2 * g.Ψ) * ϵ # momentum
+    D(δ) * ϵ ~ (-(1 + w) * (θ - 3*g.Φ) - 3 * g.ℰ * (cₛ² - w) * δ) * ϵ # energy overdensity
+    D(θ) * ϵ ~ (-g.ℰ * (1 - 3*w) - D(w) / (1 + w) * θ + cₛ² / (1 + w) * k^2 * δ - k^2 * σ + k^2 * g.Ψ) * ϵ # momentum
     σ * ϵ ~ 0 # shear stress
 ]
 
@@ -79,7 +79,7 @@ initialization_eqs = [
 
 # 5. Specify parameter relationships
 defaults = [
-    ρ0 => 3/8π * Ω0 # set ρ0 from Ω0
+    ρ₀ => 3/8π * Ω₀ # set ρ₀ from Ω₀
 ]
 
 # 6. Pack into an ODE system called "X"
@@ -108,9 +108,9 @@ To test, let us set some parameters and solve both models with one perturbation
 For the ΛCDM model:
 ```@example ext
 θ1 = Dict(
-    M1.γ.T0 => 2.7,
-    M1.c.Ω0 => 0.27,
-    M1.b.Ω0 => 0.05,
+    M1.γ.T₀ => 2.7,
+    M1.c.Ω₀ => 0.27,
+    M1.b.Ω₀ => 0.05,
     M1.ν.Neff => 3.0,
     M1.g.h => 0.7,
     M1.b.rec.Yp => 0.25
@@ -123,7 +123,7 @@ And for the w₀wₐCDM model:
 θ2 = merge(θ1, Dict(
     M2.X.w0 => -0.9,
     M2.X.wa => 0.2,
-    M2.X.cs² => 1.0
+    M2.X.cₛ² => 1.0
 ))
 sol2 = solve(M2, θ2, ks)
 ```

docs/src/getting_started.md

@@ -42,9 +42,9 @@ Next, we set our desired cosmological parameters and wavenumbers, and solve the
 ```@example getting_started
 using Unitful, UnitfulAstro # for interfacing without internal code units
 pars = Dict(
-    M.γ.T0 => 2.7,
-    M.c.Ω0 => 0.27,
-    M.b.Ω0 => 0.05,
+    M.γ.T₀ => 2.7,
+    M.c.Ω₀ => 0.27,
+    M.b.Ω₀ => 0.05,
     M.ν.Neff => 3.0,
     M.g.h => 0.7,
     M.b.rec.Yp => 0.25

docs/src/models.md

@@ -9,7 +9,7 @@ SymBoltz.RMΛ
 ```@example RMΛ
 using SymBoltz, Unitful, UnitfulAstro, Plots
 M = RMΛ()
-pars = Dict(M.r.Ω0 => 5e-5, M.m.Ω0 => 0.3, M.g.h => 1.0, M.r.T0 => 0.0) # TODO: don't pass h and T0 to avoid infinite loop
+pars = Dict(M.r.Ω₀ => 5e-5, M.m.Ω₀ => 0.3, M.g.h => 1.0, M.r.T₀ => 0.0) # TODO: don't pass h and T₀ to avoid infinite loop
 ks = [1e-3, 1e-2, 1e-1, 1e-0] / u"Mpc"
 sol = solve(M, pars, ks)
 p1 = plot(sol, log10(M.g.a), [M.r.ρ, M.m.ρ, M.Λ.ρ, M.G.ρ] ./ M.G.ρ; N = 10000)
@@ -46,7 +46,7 @@ M = w0waCDM()
 pars = merge(parameters_Planck18(M), Dict(
     M.X.w0 => -0.9,
     M.X.wa => 0.2,
-    M.X.cs² => 1.0
+    M.X.cₛ² => 1.0
 ))
 ks = [1e-3, 1e-2, 1e-1, 1e-0] / u"Mpc"
 sol = solve(M, pars, ks)
@@ -67,8 +67,8 @@ using SymBoltz, ModelingToolkit
 M = BDΛCDM()
 D = Differential(M.t)
 
-pars_fixed = merge(parameters_Planck18(M), Dict(M.G.ω => 100.0, D(M.G.ϕ) => 0.0)) # unspecified: M.Λ.Ω0, M.G.ϕ
-pars_guess = Dict(M.G.ϕ => 0.95, M.Λ.Ω0 => 0.7) # initial guesses for shooting method
+pars_fixed = merge(parameters_Planck18(M), Dict(M.G.ω => 100.0, D(M.G.ϕ) => 0.0)) # unspecified: M.Λ.Ω₀, M.G.ϕ
+pars_guess = Dict(M.G.ϕ => 0.95, M.Λ.Ω₀ => 0.7) # initial guesses for shooting method
 pars_shoot = shoot(M, pars_fixed, pars_guess, [M.g.ℰ ~ 1, M.G.G ~ 1]; thermo = false, backwards = false) # exact solutions
 pars = merge(pars_fixed, pars_shoot) # merge fixed and shooting parameters
 ```
@@ -89,8 +89,8 @@ using SymBoltz, ModelingToolkit, Unitful, UnitfulAstro, Plots
 M = SymBoltz.BDRMΛ()
 D = Differential(M.t)
 
-pars_fixed = Dict(M.r.Ω0 => 5e-5, M.m.Ω0 => 0.3, M.g.h => 1.0, M.r.T0 => 0.0, M.G.ω => 10.0, D(M.G.ϕ) => 0.0) # unspecified: M.Λ.Ω0, M.G.ϕ
-pars_guess = Dict(M.G.ϕ => 0.95, M.Λ.Ω0 => 0.7) # initial guesses for shooting method
+pars_fixed = Dict(M.r.Ω₀ => 5e-5, M.m.Ω₀ => 0.3, M.g.h => 1.0, M.r.T₀ => 0.0, M.G.ω => 10.0, D(M.G.ϕ) => 0.0) # unspecified: M.Λ.Ω₀, M.G.ϕ
+pars_guess = Dict(M.G.ϕ => 0.95, M.Λ.Ω₀ => 0.7) # initial guesses for shooting method
 pars_shoot = shoot(M, pars_fixed, pars_guess, [M.g.ℰ ~ 1, M.G.G ~ 1]; thermo = false, backwards = false) # exact solutions
 pars = merge(pars_fixed, pars_shoot) # merge fixed and shooting parameters
 

docs/src/parameter_fitting.md

@@ -68,17 +68,17 @@ M = RMΛ()
 
 function dL(z, sol::CosmologySolution)
     a = @. 1 / (z + 1)
-    t = sol(M.g.z, z, M.t / M.g.H0)
-    t0 = sol(M.g.z, 0, M.t / M.g.H0)
+    t = sol(M.g.z, z, M.t / M.g.H₀)
+    t0 = sol(M.g.z, 0, M.t / M.g.H₀)
     r = @. SymBoltz.c * (t0 .- t)
     return @. r / a / SymBoltz.Gpc
 end
 
 function dL(z, M::CosmologyModel, Ωm0, h)
     pars = Dict(
-        M.r.Ω0 => 5e-5,
-        M.m.Ω0 => Ωm0,
-        M.g.h => h, M.r.T0 => 0.0 # TODO: don't set
+        M.r.Ω₀ => 5e-5,
+        M.m.Ω₀ => Ωm0,
+        M.g.h => h, M.r.T₀ => 0.0 # TODO: don't set
     )
     sol = solve(M, pars; thermo = false)
     return dL(z, sol)

src/SymBoltz.jl

@@ -18,7 +18,7 @@ using DataInterpolations
 # TODO: try different AD sensitivity algorithms: https://docs.sciml.ai/SciMLSensitivity/stable/getting_started/
 # TODO: connector systems for Compton scattering / recombination etc.
 
-using ModelingToolkit: t_nounits as t, D_nounits as D # t is conformal time in units of 1/H0
+using ModelingToolkit: t_nounits as t, D_nounits as D # t is conformal time in units of 1/H₀
 k = only(GlobalScope.(@parameters k)) # perturbation wavenumber
 ϵ = only(GlobalScope.(@parameters ϵ)) # perturbative expansion parameter
 

src/components/metric.jl

@@ -7,18 +7,18 @@ Create a symbolic component for the perturbed FLRW spacetime metric in the confo
 """
 function metric(; name = :g, kwargs...)
     vars = a, ℰ, E, H, ℋ, Φ, Ψ, Φ′, Ψ′, g11, g22, h11, h22, z = GlobalScope.(@variables a(t) ℰ(t) E(t) H(t) ℋ(t) Φ(t) Ψ(t) Φ′(t) Ψ′(t) g11(t) g22(t) h11(t) h22(t) z(t))
-    pars = H0, h = GlobalScope.(@parameters H0 h)
+    pars = H₀, h = GlobalScope.(@parameters H₀ h)
     defs = [
-        H0 => H100 * h
-        h => H0 / H100
+        H₀ => H100 * h
+        h => H₀ / H100
     ]
     description = "Spacetime FLRW metric in Newtonian gauge"
     return ODESystem([
         z ~ 1/a - 1
-        ℰ ~ D(a) / a # ℰ = ℋ/ℋ0 = ℋ/H0
-        E ~ ℰ / a # E = H/H0
-        ℋ ~ ℰ * H0
-        H ~ E * H0
+        ℰ ~ D(a) / a # ℰ = ℋ/ℋ0 = ℋ/H₀
+        E ~ ℰ / a # E = H/H₀
+        ℋ ~ ℰ * H₀
+        H ~ E * H₀
         g11 ~ -a^2
         g22 ~ +a^2
         h11 * ϵ ~ -2 * a^2 * Ψ * ϵ