(* | Evolution of the FLRW Universe | Simon Tyran, Vienna | flrw.yukterez.net | *) set = {"GlobalAdaptive", "MaxErrorIncreases"->100, Method->"GaussKronrodRule"}; (* Integration Rule *) n = 100; (* Recursion Depth *) int[f_, {x_, xmin_, xmax_}] := (* Integral *) NIntegrate[f, {x, xmin, xmax}, Method->set, MaxRecursion->n, WorkingPrecision->wp]; wp = 24; (* Working Precision *) im = 640; (* Image Size *) ηmax = 63045845/1000000; ηMax=Round[ηmax]; (* Conformal Plot Range *) c = 299792458 m/sek; (* Lightspeed *) G = 667384*^-16 m^3 kg^-1 sek^-2; (* Newton Constant *) Gyr = 10^7*36525*24*3600 sek; (* Billion Years *) Glyr = Gyr*c; (* Billion Lightyears *) Mpc = 30856775777948584200000 m; (* Megaparsec *) kB = 13806488*^-30 kg m^2/sek^2/K; (* Boltzmann Constant *) h = 662606957*^-42 kg m^2/sek; (* Planck Constant *) ρc[H_] := 3H^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 *) kg = m = sek = 1; (* SI Units *) ΩR = 168132/100000 ρR 8 π G/3/H0^2; (* Radiation Proportion including Neutrinos *) ΩM = 315/1000; (* Matter Proportion including Dark Matter *) ΩΛ = 1-ΩM-ΩR; (* Dark Energy Proportion *) ΩT = ΩR+ΩM+ΩΛ; (* Total Density over Critical Density *) ΩK = 1-ΩT; (* Curvature Density *) H0 = 67150 m/Mpc/sek; (* Hubble Constant *) H[a_] := H0 Sqrt[ΩR/a^4+ΩM/a^3+ΩK/a^2+ΩΛ] (* Hubble Parameter *) sol = Quiet[NDSolve[{A'[t]/A[t] == H[A[t]], A[0] == 1*^-15}, A, {t, 0, 500 Gyr}, Method->{"ImplicitRungeKutta", "DifferenceOrder"->20}, MaxSteps->Infinity, WorkingPrecision->wp]]; a[t_] := Evaluate[(A[t]/.sol)[[1]]]; (* Scale Factor a by Time t *) т[a_] := int[1/A/H[A], {A, 0, a}]; (* Time t by Scale Factor a *) rP[t_] := a[t] int[c/a[т], {т, 0, t}]; (* Particle Horizon by t *) rp[a_] := a int[c/A^2/H[A], {A, 0, a}]; (* Particle Horizon by a *) rE[t_] := a[t] int[c/a[т], {т, t, Infinity}]; (* Event Horizon by t *) re[a_] := a int[c/A^2/H[A], {A, a, Infinity}]; (* Event Horizon by a *) rL[t0_, t_] := a[t] int[c/a[т], {т, t, t0}]; (* Light Cone by t *) rl[a0_, a_] := a int[c/A^2/H[A], {A, a, a0}]; (* Light Cone by a *) rH[t_] := c/H[a[t]]; (* Hubble Radius by t *) rh[a_] := c/H[a]; (* Hubble Radius by a *) ã[η_] := Quiet[FindRoot[rp[Ã]/Ã/Glyr-η, (* Scale Factor a by Conformal Time η *) {Ã, 0.00001}, WorkingPrecision->wp, MaxIterations->1000][[1, 2]]]; t0 = t/.FindRoot[a[t]-1, {t, 10 Gyr}]; ti = t Gyr; τi = τ Gyr; "t0"->t0/Gyr "Gyr" (* Current Time *) "PROPER DISTANCES, f(t)" pt = Quiet[Plot[Evaluate[ {rH[τi]/Glyr, rP[τi]/Glyr, rE[τi]/Glyr}], {τ, 0, 70}, Frame->True, AspectRatio->1, PlotRange->{{0, 70}, {0, 70}}, PlotStyle->{{Thickness[0.005]}, {Darker[Green], Thickness[0.005]}, {Purple, Thickness[0.005]}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; plot1[t_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[Evaluate[ {rL[ti, τi]/Glyr, -rL[ti, τi]/Glyr}], {τ, 0, 70}, Frame->True, AspectRatio->1, PlotRange->{{0, 70}, {0, 70}}, PlotStyle->{{Orange, Thickness[0.005]}, {{Orange, Thickness[0.005]}, Dashed}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}], pt]], 90 Degree]}}]]; Do[Print[plot1[t]], {t, 1t0/Gyr, 5t0/Gyr, 1t0/Gyr}] plot2 = Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[ Join[{0}, Table[a[τ Gyr] n^(7/2)/250, {n, 4, 55}]]], {τ, 0, 70}, Frame->True, AspectRatio->1, PlotRange->{{0, 70}, {0, 70}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Lighter@LightGray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] "COMOVING DISTANCES, f(t)" ct = Quiet[Plot[Evaluate[ {rH[τi]/(a[τi]Glyr), rP[τi]/(a[τi]Glyr), rE[τi]/(a[τi]Glyr)}], {τ, 0, 70}, Frame->True, AspectRatio->1, PlotRange->{{0, 70}, {0, 70}}, PlotStyle->{{Thickness[0.005]}, {Darker[Green], Thickness[0.005]}, {Purple, Thickness[0.005]}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; plot3[t_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[Evaluate[ {rL[ti, τi]/(a[τi]Glyr), -rL[ti, τi]/(a[τi]Glyr)}], {τ, 0, 70}, Frame->True, AspectRatio->1, PlotRange->{{0, 70}, {0, 70}}, PlotStyle->{{Orange, Thickness[0.005]}, {{Orange, Thickness[0.005]}, Dashed}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}], ct]], 90 Degree]}}]]; Do[Print[plot3[t]], {t, 1t0/Gyr, 5t0/Gyr, 1t0/Gyr}] plot4 = Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[ Join[{0}, Table[ n^(7/2)/250, {n, 4, 55}]]], {τ, 0, 70}, Frame->True, AspectRatio->1, PlotRange->{{0, 70}, {0, 70}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Lighter@LightGray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] "PROPER DISTANCES, f(a)" pa = Quiet[Plot[Evaluate[ {rh[α]/Glyr, rp[α]/Glyr, re[α]/Glyr}], {α, 0, 70 Gyr/t0}, Frame->True, AspectRatio->1, PlotRange->{{0, 70 Gyr/t0}, {0, 70}}, PlotStyle->{{Thickness[0.005]}, {Darker[Green], Thickness[0.005]}, {Purple, Thickness[0.005]}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; plot5[å_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[Evaluate[ {rl[å, α]/Glyr, -rl[å, α]/Glyr}], {α, 0, 70 Gyr/t0}, Frame->True, AspectRatio->1, PlotRange->{{0, 70 Gyr/t0}, {0, 70}}, PlotStyle->{{Orange, Thickness[0.005]}, {{Orange, Thickness[0.005]}, Dashed}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}], pa]], 90 Degree]}}]]; Do[Print[plot5[å]], {å, 1, 5, 1}] plot6 = Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[ Join[{0}, Table[α n^(7/2)/250, {n, 4, 55}]]], {α, 0, 70 Gyr/t0}, Frame->True, AspectRatio->1, PlotRange->{{0, 70 Gyr/t0}, {0, 70}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Lighter@LightGray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] "COMOVING DISTANCES, f(a)" ca = Quiet[Plot[Evaluate[ {rh[α]/Glyr/α, rp[α]/Glyr/α, re[α]/Glyr/α}], {α, 0, 70 Gyr/t0}, Frame->True, AspectRatio->1, PlotRange->{{0, 70 Gyr/t0}, {0, 70}}, PlotStyle->{{Thickness[0.005]}, {Darker[Green], Thickness[0.005]}, {Purple, Thickness[0.005]}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; plot7[å_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[Evaluate[ {rl[å, α]/Glyr/α, -rl[å, α]/Glyr/α}], {α, 0, 70 Gyr/t0}, Frame->True, AspectRatio->1, PlotRange->{{0, 70 Gyr/t0}, {0, 70}}, PlotStyle->{{Orange, Thickness[0.005]}, {{Orange, Thickness[0.005]}, Dashed}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}], ca]], 90 Degree]}}]]; Do[Print[plot7[å]], {å, 1, 5, 1}] plot8 = Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[ Join[{0}, Table[n^(7/2)/250, {n, 4, 55}]]], {α, 0, 70 Gyr/t0}, Frame->True, AspectRatio->1, PlotRange->{{0, 70 Gyr/t0}, {0, 70}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Lighter@LightGray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] (* Interpolating Function for faster Output *) ā = Quiet[Interpolation[ParallelTable[{η, ã[η]}, {η, ηmax/im, ηmax, ηmax/im}]]]; rpN = rp[1]/Glyr; "PROPER DISTANCES, f(η)" pη = Quiet[Plot[rh[ā[Ct]]/Glyr, {Ct, 0, ηMax}, Frame->True, AspectRatio->70/ηMax, PlotRange->{{0, ηMax}, {0, 70}}, PlotStyle->{Thickness[0.005]}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; Pη = Quiet[Plot[Evaluate[ ā[Ct] {Ct, (ηmax-Ct)}], {Ct, 0, ηMax}, Frame->True, AspectRatio->70/ηMax, PlotRange->{{0, ηMax}, {0, 70}}, PlotStyle->{{Darker[Green], Thickness[0.005]}, {Purple, Thickness[0.005]}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; plot9[η_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[Evaluate[ ā[Ct] {η-Ct, Ct-η}], {Ct, 0, ηMax}, Frame->True, AspectRatio->70/ηMax, PlotRange->{{0, ηMax}, {0, 70}}, PlotStyle->{{Orange, Thickness[0.005]}, {{Orange, Thickness[0.005]}, Dashed}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}], pη, Pη]], 90 Degree]}}]]; Do[Print[plot9[η]], {η, rpN/5, rpN, rpN/5}] plot10 = Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[ Join[{0}, Table[ā[Ct] n^(7/2)/250, {n, 4, 55}]]], {Ct, 0, ηMax}, Frame->True, AspectRatio->70/ηMax, PlotRange->{{0, ηMax}, {0, 70}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Lighter@LightGray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] "COMOVING DISTANCES, f(η)" cη = Quiet[Plot[rh[ā[Ct]]/Glyr/ā[Ct], {Ct, 0, ηMax}, Frame->True, AspectRatio->70/ηMax, PlotRange->{{0, ηMax}, {0, 70}}, PlotStyle->{Thickness[0.005]}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; Cη = Quiet[Plot[Evaluate[{Ct, ηmax-Ct}], {Ct, 0, ηMax}, Frame->True, AspectRatio->70/ηMax, PlotRange->{{0, ηMax}, {0, 70}}, PlotStyle->{{Darker[Green], Thickness[0.005]}, {Purple, Thickness[0.005]}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; plot11[η_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[Evaluate[{η-Ct, Ct-η}], {Ct, 0, ηMax}, Frame->True, AspectRatio->70/ηMax, PlotRange->{{0, ηMax}, {0, 70}}, PlotStyle->{{Orange, Thickness[0.005]}, {{Orange, Thickness[0.005]}, Dashed}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}], cη, Cη]], 90 Degree]}}]]; Do[Print[plot11[η]], {η, rpN/5, rpN, rpN/5}] plot12 = Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[ Join[{0}, Table[n^(7/2)/250, {n, 4, 55}]]], {Ct, 0, ηMax}, Frame->True, AspectRatio->70/ηMax, PlotRange->{{0, ηMax}, {0, 70}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Lighter@LightGray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] Quit[] (* The method below evaluates everything, even the functions that are straight lines, therefore it is slower than the method above where therefore straight lines are given as simple functions *) "PROPER DISTANCES, f(η)" pη = Quiet[Plot[Evaluate[ {rh[ā[Ct]]/Glyr, rp[ā[Ct]]/Glyr, re[ā[Ct]]/Glyr}], {Ct, 0, ηMax}, Frame->True, AspectRatio->70/ηMax, PlotRange->{{0, ηMax}, {0, 70}}, PlotStyle->{{Thickness[0.005]}, {Darker[Green], Thickness[0.005]}, {Purple, Thickness[0.005]}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; plot9[η_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[Evaluate[ {If[Ct<=η, rl[ā[η], ā[Ct]]/Glyr, -1], If[Ct>=η, -rl[ā[η], ā[Ct]]/Glyr, -1]}], {Ct, 0, ηMax}, Frame->True, AspectRatio->70/ηMax, PlotRange->{{0, ηMax}, {0, 70}}, PlotStyle->{{Orange, Thickness[0.005]}, {{Orange, Thickness[0.005]}, Dashed}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}], pη]], 90 Degree]}}]]; Do[Print[plot9[η]], {η, rpN/5, rpN, rpN/5}] plot10 = Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[ Join[{0}, Table[ā[Ct] n^(7/2)/250, {n, 4, 55}]]], {Ct, 0, ηMax}, Frame->True, AspectRatio->70/ηMax, PlotRange->{{0, ηMax}, {0, 70}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Lighter@LightGray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] "COMOVING DISTANCES, f(η)" cη = Quiet[Plot[Evaluate[ {rh[ā[Ct]]/Glyr/ā[Ct], rp[ā[Ct]]/Glyr/ā[Ct], re[ā[Ct]]/Glyr/ā[Ct]}], {Ct, 0, ηMax}, Frame->True, AspectRatio->70/ηMax, PlotRange->{{0, ηMax}, {0, 70}}, PlotStyle->{{Thickness[0.005]}, {Darker[Green], Thickness[0.005]}, {Purple, Thickness[0.005]}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; plot11[η_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[Evaluate[ {If[Ct<=η, rl[ā[η], ā[Ct]]/Glyr/ā[Ct], -1], If[Ct>=η, -rl[ā[η], ā[Ct]]/Glyr/ā[Ct], -1]}], {Ct, 0, ηMax}, Frame->True, AspectRatio->70/ηMax, PlotRange->{{0, ηMax}, {0, 70}}, PlotStyle->{{Orange, Thickness[0.005]}, {{Orange, Thickness[0.005]}, Dashed}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}], cη]], 90 Degree]}}]]; Do[Print[plot11[η]], {η, rpN/5, rpN, rpN/5}] plot12 = Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[ Join[{0}, Table[n^(7/2)/250, {n, 4, 55}]]], {Ct, 0, ηMax}, Frame->True, AspectRatio->70/ηMax, PlotRange->{{0, ηMax}, {0, 70}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Lighter@LightGray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]]