Das Prinzip des Seins

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

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 *)

c    = 299792458 m/sek; (* Lightspeed *)
ca   = 1; (* Perturbation Velocity *)
G    = 667384*^-16 m^3 kg^-1 sek^-2;  (* Newton's 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's Constant *)
h    = 662606957*^-42 kg m^2/sek; (* Planck's Constant *)

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

ρ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 *)

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

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; (* Dark Energy Function *)

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]

w0 = -1;
kr1[f_, cl_] := Quiet[{cl, Thickness[0.006], Circle[{0, 0}, f]}];
kr2[f_, cl_] := Quiet[{cl, Circle[{0, 0}, f]}];

Do[Print[Grid[{
    {Quiet[Graphics[{
         {Table[kr2[α[τ Gyr, -1] (3/2)^32/(3/2)^n, Lighter[Lighter[Gray]]], {n, 1, 42}]},
        {kr1[rH[τ Gyr, w0]/Glyr, Cyan]}, {kr1[pH[τ Gyr, w0]/Glyr, Green]},
        {kr1[eH[τ Gyr, w0]/Glyr, Magenta]},
        {kr1[(Quiet[FindMaximum[{Lc[Tm Gyr, τ Gyr, -1]/Glyr}, {Tm, 0}][[1]]]), Orange]}
        }, PlotRange -> {{-60, 60}, {-60, 60}},
             Frame -> True, PlotRangeClipping -> True, ImageSize -> 440]]},
    {Evaluate[N[τ, 8]]}}]], {τ, 1/10000, 25.6, 25.6/500.0}]