Das Prinzip des Seins

(* Syntax: Mathematica ||| lcdm.yukterez.net *)

c       = 299792458 m/sek; (* Lightspeed *)
ca      = 1; (* Perturbation Velocity *)
G       = 667384*^-16 m^3 kg^-1 sek^-2;  (* Newton Constant *)
Gyr     = 10^7*36525*24*3600*sek; (* Billion Year Scale *)
Glyr    = Gyr*c; (* Billion Lightyear Scale *)
Mpc     = 30856775777948584200000 m; (* Megaparsec *)
kB      = 13806488*^-30 kg m^2/sek^2/K; (* Boltzmann *)
h       = 662606957*^-42 kg m^2/sek; (* Planck *)
ρR      = 8 π^5 kB^4 T^4/15/c^5/h^3; (* Radiation Density *)
ρΛ      = ρc[H0] ΩΛ; (* Dark Energy Density *)
T       = 2725/1000 K; (* CMB Temperature *)

ρc[H_] := 3 H^2/8/π/G; (* Critical Density *)

H0 = 67150 m/Mpc/sek; (* Hubble Constant *)
ΩR = ρR/ρc[H0]; ΩM = 317/1000; ΩΛ = 683/1000 - ΩR; ΩT = ΩR + ΩM + ΩΛ; ΩK = 1 - ΩT; (* Density Parameters *)

kg = 1; m = 1; sek = 1;  K = 1; (* Units *)
set  = {"GlobalAdaptive", "MaxErrorIncreases" -> 100, Method -> "GaussKronrodRule"}; (* Integration Rule *)
n    = 100; (* Recursion Depth *)
tE   = 300 Gyr;  (* Eventhorizon Limit *)

aE[t_]     := Power[(Sqrt[ΩM/ΩΛ] Sinh[(3 H0 Sqrt[ΩΛ])/2 t])^2, (3)^-1]; (* Solving Region *)
w[a_, w0_] := (1 + w0) (Sqrt[1 + (ΩΛ^-1 - 1) a^-3] - (ΩΛ^-1 -1) a^-3 Tanh[1/Sqrt[1 + (ΩΛ^-1 - 1) a^-3]]^-1)^2 (1/Sqrt[ΩΛ] - (ΩΛ^-1 - 1) Tanh[Sqrt[ΩΛ]]^-1)^-2 -1;

F[a_, w0_] := Sqrt[ΩR a^-4 + ΩM a^-3 + ΩK a^-2 + ΩΛ a^(-3 (w[a, w0] + 1))]; (* Density Function by Scalefactor*)
φ[z_, w0_] := Sqrt[ΩR (z + 1)^4 + ΩM (z + 1)^3 + ΩK (z + 1)^2 + ΩΛ ((1 + z)^(3 (w[1/(z + 1), w0] + 1))) ]; (* Density Function by Redshift *)
H[a_, w0_] := H0 F[a, w0]; (* Hubble Parameter by Scalefactor *)
ε[z_, w0_] := H0 φ[z, w0]; (* Hubble Parameter by Redshift *)

int[f_, {x_, xmin_, xmax_}] := Quiet[NIntegrate[f, {x, xmin, xmax}, Method -> set, MaxRecursion -> n]];

ta[A_, w0_] := int[1/a/ H[a, w0], {a, 0, A}]; (* Time by Scalefactor *)
α[τ_, w0_]  := Quiet[A /.FindRoot[ta[A, w0] - τ, {A, 1}]] (* Scalefactor by Time *)

tz[Z_, w0_] := int[1/(1 + z)/ ε[z, w0], {z, Z, \[Infinity]}]; (* Time by Redshift *)
χ[τ_, w0_]  := Z /. Quiet[FindRoot[tz[Z, w0] - τ, {Z, 0}]] (* Redshift by Time *)

rH[τ_, w0_] := c/H[α[τ, w0], w0]; (* Hubble Radius *)

lC[τ_, w0_]     := int[-c α[τ, w0]/a^2/H[a, w0], {a, 1, α[τ, w0]}]; (* Light Cone of t0 *)
Lc[τ_, t_, w0_] := int[-c α[τ, w0]/a^2/H[a, w0], {a, α[t, w0], α[τ, w0]}]; (* Light Cone of t *)
eH[τ_, w0_]     := α[τ, w0] int[c/(α[time, w0]), {time, τ, tE}]; (* Event Horizon *)
pH[τ_, w0_]     := int[-α[τ, w0] c/a^2/H[a, w0], {a, α[τ, w0], 0}]; (* Particle Horizon *)
g[τ_, w0_]      := tc /. Quiet[FindRoot[pH[tc, w0]/c - τ, {tc, τ}]]; (* Conformal Time *)

ωR[τ_, w0_] := ΩR α[τ, w0]^-4/ρc[H[α[τ, w0]]]; (* Radiation Evolution *)
ωM[τ_, w0_] := ΩM α[τ, w0]^-3/ρc[H[α[τ, w0]]]; (* Matter Evolution *)
ωK[τ_, w0_] := ΩK α[τ, w0]^-2/ρc[H[α[τ, w0]]]; (* Curvature Evolution *)
ωΛ[τ_, w0_] := ΩΛ α[τ, w0]^(-3 (w[α[τ, w0], w0] + 1))/ρc[H[α[τ, w]]]; (* Dark Energy Evolution *)

t0[w0_] := ta[1, w0]/Gyr; (* Age of the Universe, now *)
"t0 in Gyr" -> t0[-1]

Plot[(tz[0, -1] - tz[z2, -1])/Gyr, {z2, 0, 4}, Frame -> True,
PlotRange -> All, PlotStyle -> {Red, Thick}]
(* Light Travel Time *)

Plot[lC[tz[z2, -1], -1]/Glyr, {z2, 0, 4}, Frame -> True,
PlotRange -> All, PlotStyle -> {Green, Thick}]
(* Distance at Emission *)

Plot[lC[tz[z2, -1], -1] (z2 + 1)/Glyr, {z2, 0, 4}, Frame -> True,
PlotRange -> All, PlotStyle -> {Blue, Thick}]
(* Distance at  Absorption *)

Plot[ε[z2, -1] lC[tz[z2, -1], -1]/c, {z2, 0, 4}, Frame -> True,
PlotRange -> All, PlotStyle -> {Orange, Thick}]
(* Recessional Velocity at Emission *)

Plot[H0 lC[tz[z2, -1], -1] (z2 + 1)/c, {z2, 0, 4}, Frame -> True,
PlotRange -> All, PlotStyle -> {Brown, Thick}]
(* Recessional Velocity at Absorption *)

v = (2 c z + c z^2)/(2 + 2 z + z^2)/c; Plot[v, {z, 0, 4}, Frame -> True,
PlotRange -> All, PlotStyle -> {Magenta, Thick}]
(* Local Minkowski Redshift *)

z1  = 4; (* Redshift *)
t1  = tz[z1, -1]; "Age at Emission" -> t1/Gyr "Gyr"
tr  = tz[0, -1] - tz[z1, -1]; "Light Travel Time" -> tr/Gyr "Gyr"
lC1 = lC[t1, -1]; "Distance at Emission" -> lC1/Glyr "Glyr"
lC2 = lC[t1, -1] (z1 + 1); "Distance at Absorption" -> lC2/Glyr "Glyr"
v1  = ε[z1, -1] lC1; "Recessional Velocity at  Emission" -> v1/c "c"
v2  = H0 lC2; "Recessional Velocity at  Absorption" -> v2/c "c"

                   Age at Emission -> 1.53583 Gyr
                 Light Travel Time -> 12.2838 Gyr
              Distance at Emission -> 4.79045 Glyr
            Distance at Absorption -> 23.9522 Glyr
  Recessional Velocity at Emission -> 2.08956 c
Recessional Velocity at Absorption -> 1.64492 c