🗎 Code for the
Big Crunch space/time-diagrams
(* Radiation Dominated Big Crunch Universe | Simon Tyran | https://yukterez.net *) set = {"GlobalAdaptive", "MaxErrorIncreases"->100}; (* Integration Rule *) n = 100; (* Recursion Depth *) int[f_, {x_, xmin_, xmax_}] := (* Integral *) Quiet[NIntegrate[f, {x, xmin, xmax}, Method->set, MaxRecursion->n, WorkingPrecision->wp]]; wp = MachinePrecision; (* Working Precision *) im = 180; (* Image Size *) prmax = 70; ptmax = tmax; (* Regular 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 *) ρc[H_] := 3H^2/8/π/G; (* Critical Density *) kg = m = sek = 1; (* SI Units *) ΩR = 2; (* Radiation Proportion including Neutrinos *) ΩM = 0; (* Matter Proportion including Dark Matter *) ΩΛ = 0; (* Dark Energy Proportion *) ΩT = ΩR+ΩM+ΩΛ; (* Total Density over Critical Density *) ΩK = 1-ΩT; (* Curvature Density *) rK = c/H0/Sqrt[-ΩK]; (* Curvature Radius *) H0 = 67150 m/Mpc/sek; (* Hubble Constant *) H[a_] := H0 Sqrt[ΩR/a^4+ΩM/a^3+ΩK/a^2+ΩΛ] (* Hubble Parameter *) Ж[t_] := (Sqrt[ΩR]+H0(t-t ΩR))/(t(2Sqrt[ΩR]+H0 (t-t ΩR))); a[t_] := Re[Sqrt[H0 t (2Sqrt[ΩR]+H0 (t-t ΩR))]]; (* Scale Factor a by Time t *) т[a_] := Re[(2H0 Sqrt[ΩR]-Sqrt[-4a^2 H0^2 (-1+ΩR)+4H0^2ΩR])/(2H0^2(ΩR-1))]; rP[t_] := Re[a[t] int[c/a[т], {т, 0, t}]]; (* Proper Particle Horizon by t *) rp[a_] := Re[a int[c/A^2/H[A], {A, 0, a}]]; (* Proper Particle Horizon by a *) RP[t_] := Re[int[c/a[т], {т, 0, t}]]; (* Comoving Particle Horizon by t *) Rp[a_] := Re[int[c/A^2/H[A], {A, 0, a}]]; (* Comoving Particle Horizon by a *) rE[t_] := Nothing; (* Proper Event Horizon by t *) re[a_] := Nothing; (* Proper Event Horizon by a *) RE[t_] := Nothing; (* Comoving Event Horizon by t *) Rε[a_] := Nothing; (* Comoving Event Horizon by a *) rL[t0_, t_] := Re[a[t] int[c/a[т], {т, t, t0}]]; (* Proper Light Cone by t *) rl[a0_, a_] := Re[a int[c/A^2/H[A], {A, a, a0}]]; (* Proper Light Cone by a *) RL[t0_, t_] := Re[int[c/a[т], {т, t, t0}]]; (* Comoving Light Cone by t *) Rl[a0_, a_] := Re[int[c/A^2/H[A], {A, a, a0}]]; (* Comoving Light Cone by a *) rH[t_] := c/Abs[Ж[t]]; (* Proper Hubble Radius by t *) rh[a_] := c/Abs[H[a]]; (* Proper Hubble Radius by a *) RH[t_] := c/Abs[Ж[t]*a[t]]; (* Comoving Hubble Radius by t *) Rh[a_] := c/Abs[H[a]*a]; (* Comoving Hubble Radius by a *) t0 = Re[t/.FindRoot[a[t]-1, {t, 10 Gyr}]]; tmid = Sqrt[2]/H0; ti = t Gyr; τi = τ Gyr; "t0"->N[t0/Gyr] "Gyr" (* Current Time *) tMax = 2*tmid; tmax = tMax/Gyr; rpN = π/2/H0/Gyr; ηmax = 2 rpN; Ť[η_] := Quiet[FindRoot[RP[τ Gyr]/Glyr-η, (* t by η *) {τ, ((Sin[η π/ηmax-π/2]+1)/2) ptmax}, WorkingPrecision->wp, MaxIterations->1000][[1, 2]]] "PROPER DISTANCES, f(t)" pt = Quiet[Plot[ {rH[τi]/Glyr, rP[τi]/Glyr, π rK a[τi]/Glyr}, {τ, 0, ptmax}, Frame->True, AspectRatio->prmax/ptmax, FrameTicks->None, PlotRange->{{0, ptmax}, {0, prmax}}, 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[ {rL[ti, τi]/Glyr, -rL[ti, τi]/Glyr}, {τ, 0, ptmax}, Frame->True, AspectRatio->prmax/ptmax, FrameTicks->None, PlotRange->{{0, ptmax}, {0, prmax}}, 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, {ptmax/2}}] plot2 = Rasterize[Grid[{{Rotate[Quiet[Plot[ Join[{0}, Table[a[τ Gyr] n/Sqrt[2], {n, 10, 70, 10}], Table[a[τ Gyr] n/Sqrt[2], {n, {100, 160, 340}}]], {τ, 0, ptmax}, Frame->True, AspectRatio->prmax/ptmax, FrameTicks->None, PlotRange->{{0, ptmax}, {0, prmax}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Gray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] "COMOVING DISTANCES, f(t)" ct = Quiet[Plot[ {RH[τi]/Glyr, RP[τi]/Glyr, π rK/Glyr}, {τ, 0, ptmax}, Frame->True, AspectRatio->prmax/ptmax, FrameTicks->None, PlotRange->{{0, ptmax}, {0, prmax}}, 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[ {rL[ti, τi]/(a[τi]Glyr), -rL[ti, τi]/(a[τi]Glyr)}, {τ, 0, ptmax}, Frame->True, AspectRatio->prmax/ptmax, FrameTicks->None, PlotRange->{{0, ptmax}, {0, prmax}}, 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, {ptmax/2}}] plot4 = Rasterize[Grid[{{Rotate[Quiet[Plot[ Join[{0}, Table[n, {n, 10, 100, 10}]], {τ, 0, ptmax}, Frame->True, AspectRatio->prmax/ptmax, FrameTicks->None, PlotRange->{{0, ptmax}, {0, prmax}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Gray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] "CONFORMAL DIAGRAM, f(η)" cη = Quiet[Plot[ {RH[Ť[Ct] Gyr]/Glyr, Ct, π rK/Glyr}, {Ct, 0, ηmax}, Frame->True, AspectRatio->prmax/ηmax, FrameTicks->None, PlotRange->{{0, ηmax}, {0, prmax}}, 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[ {η-Ct, Ct-η}, {Ct, 0, ηmax}, Frame->True, AspectRatio->prmax/ηmax, FrameTicks->None, PlotRange->{{0, ηmax}, {0, prmax}}, 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[plot9[η]], {η, {rpN}}] plot10 = Rasterize[Grid[{{Rotate[Quiet[Plot[ Join[{0}, Table[n, {n, 10, 100, 10}]], {Ct, 0, ηmax}, Frame->True, AspectRatio->prmax/ηmax, FrameTicks->None, PlotRange->{{0, ηmax}, {0, prmax}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Gray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] s[text_] := Style[text, FontFamily->"Lucida Console", FontSize->36]
(* Matter Dominated Big Crunch Universe | Simon Tyran | http://org.yukterez.net *) set = {"GlobalAdaptive", "MaxErrorIncreases"->1000}; (* Integration Rule *) n = 1000; (* Recursion Depth *) int[f_, {x_, xmin_, xmax_}] := (* Integral *) Quiet[NIntegrate[f, {x, xmin, xmax}, Method->set, MaxRecursion->n, WorkingPrecision->wp]]; wp = 20; (* Working Precision *) im = 180; (* Image Size *) prmax = 140; ptmax = tmax; (* Regular 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 = 0; (* Radiation Proportion including Neutrinos *) ΩM = 2; (* Matter Proportion including Dark Matter *) ΩΛ = 0; (* Dark Energy Proportion *) ΩT = ΩR+ΩM+ΩΛ; (* Total Density over Critical Density *) ΩK = 1-ΩT; (* Curvature Density *) rK = c/H0/Sqrt[-ΩK]; (* Curvature Radius *) tmid = π/H0; tMax = 2*tmid; tmax = tMax/Gyr; a0 = 1*^-15; (* Initial Scale Factor *) H0 = 67150 m/Mpc/sek; (* Hubble Constant *) f1[x_]:=(ΩM Log[-Sqrt[1-ΩM] Sqrt[x]+Sqrt[ΩM+x-ΩM x]] Sqrt[ΩM+x-ΩM x]+ Sqrt[1-ΩM] Sqrt[x] (ΩM+x-ΩM x))/((1-ΩM)^(3/2) x^(3/2) Sqrt[(ΩM+x-ΩM x)/x^3]); x1[t_]:=(-Sqrt[a0 (1-ΩM)] (-ΩM+a0 (-1+ΩM) (1+H0 t Sqrt[(a0+ΩM-a0 ΩM)/a0^3]))+ ΩM Sqrt[a0+ΩM-a0 ΩM] Log[-Sqrt[a0 (1-ΩM)]+Sqrt[a0+ΩM- a0 ΩM]])/(a0^(3/2) (1-ΩM)^(3/2) Sqrt[(a0+ΩM-a0 ΩM)/a0^3]); sol = Quiet[NDSolve[{A'[t]/A[t] == H[A[t]], A[0] == a0}, A, {t, 0, tmid}, MaxSteps->∞, WorkingPrecision->wp]]; H[a_] := H0 Sqrt[ΩR/a^4+ΩM/a^3+ΩK/a^2+ΩΛ]; Ж[t_] := H[a[t]]; (* Hubbleparameter *) Â[t_] := InverseFunction[f1[#1] &][x1[t]]; (*Scale Factor a by Time t*) Å[t_] := Evaluate[(A[t]/.sol)[[1]]]; a[t_] := If[t<tmid, Re[Å[t]], Re[Å[tMax-t]]]; (* a = Quiet[Interpolation[Join[{{0, a0}}, ParallelTable[{tmid (Sin[π t/tmid/2])^4, Re[ã[tmid (Sin[π t/tmid/2])^4]]}, {t, tmid/im, tmid-tmid/im, tmid/im}], {{tmid, 2}}, ParallelTable[{t, Re[ã[t]]}, {t, tmid+tmid/im, tMax-tmid/im, tmid/im}], {{tMax, a0}}]]]; *) rP[t_] := a[t] int[c/a[т], {т, 0, t}]; (* Proper Particle Horizon by t *) rp[a_] := a int[c/A^2/H[A], {A, a0, a}]; (* Proper Particle Horizon by a *) RP[t_] := int[c/a[т], {т, 0, t}]; (* Comoving Particle Horizon by t *) Rp[a_] := int[c/A^2/H[A], {A, a0, a}]; (* Comoving Particle Horizon by a *) rE[t_] := Nothing; (* Proper Event Horizon by t *) re[a_] := Nothing; (* Proper Event Horizon by a *) RE[t_] := Nothing; (* Comoving Event Horizon by t *) Rε[a_] := Nothing; (* Comoving Event Horizon by a *) rL[t0_, t_] := a[t] int[c/a[т], {т, t, t0}]; (* Proper Light Cone by t *) rl[a0_, a_] := a int[c/A^2/H[A], {A, a, a0}]; (* Proper Light Cone by a *) RL[t0_, t_] := int[c/a[т], {т, t, t0}]; (* Comoving Light Cone by t *) Rl[a0_, a_] := int[c/A^2/H[A], {A, a, a0}]; (* Comoving Light Cone by a *) rH[t_] := c/Abs[Ж[t]]; (* Proper Hubble Radius by t *) rh[a_] := If[1.99<a, Infinity, c/Abs[H[a]]]; (* Proper Hubble Radius by a *) RH[t_] := c/Abs[Ж[t]a[t]]; (* Comoving Hubble Radius by t *) Rh[a_] := If[a>1.99, Infinity, c/Abs[H[a]a]]; (* Comoving Hubble Radius by a *) t0 = Quiet[Re[t/.FindRoot[a[t]-1, {t, 10 Gyr}, WorkingPrecision->wp]]]; ti = t Gyr; τi = τ Gyr; "t0"->t0/Gyr "Gyr" (* Current Time *) ηmax = 2 π/H0/Gyr; rpN = ηmax/2; fx = 100/101; ηH = Quiet[Interpolation[Join[{{0, 0}}, (* Hubble Radius by Conformal Time *) ParallelTable[ {RP[tMax (Sin[π t/tmax/2])^4]/Glyr, Rh[a[tMax (Sin[π t/tmax/2])^4]]/Glyr}, {t, tmax/im, tmax-tmax/im, tmax/im}], {{ηmax, 0}}]]]; "PROPER DISTANCES, f(t)" pt = Quiet[Plot[ {rh[a[τi]]/Glyr, rP[τi]/Glyr, π rK a[τi]/Glyr}, {τ, 0, ptmax}, Frame->True, AspectRatio->prmax/ptmax, FrameTicks->None, PlotRange->{{0, ptmax}, {0, prmax}}, 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[ {rL[ti, τi]/Glyr, -rL[ti, τi]/Glyr}, {τ, 0, ptmax}, Frame->True, AspectRatio->prmax/ptmax, FrameTicks->None, PlotRange->{{0, ptmax}, {0, prmax}}, 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, {tmax/2}}] plot2 = Rasterize[Grid[{{Rotate[Quiet[Plot[ Join[{0}, Table[a[τ Gyr] n/2, {n, 20, 140, 20}], Table[a[τ Gyr] n/2, {n, {180, 260, 420}}]], {τ, 0, ptmax}, Frame->True, AspectRatio->prmax/ptmax, FrameTicks->None, PlotRange->{{0, ptmax}, {0, prmax}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Gray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] "COMOVING DISTANCES, f(t)" ct = Quiet[Plot[ {Rh[a[τi]]/Glyr, RP[τi]/Glyr, π rK/Glyr}, {τ, 0, ptmax}, Frame->True, AspectRatio->prmax/ptmax, FrameTicks->None, PlotRange->{{0, ptmax}, {0, prmax}}, 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[ {RL[ti, τi]/Glyr, -RL[ti, τi]/Glyr}, {τ, 0, ptmax}, Frame->True, AspectRatio->prmax/ptmax, FrameTicks->None, PlotRange->{{0, ptmax}, {0, prmax}}, 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, {tmax/2}}] plot4 = Rasterize[Grid[{{Rotate[Quiet[Plot[ Join[{0}, Table[n, {n, 10, prmax, 10}]], {τ, 0, ptmax}, Frame->True, AspectRatio->prmax/ptmax, FrameTicks->None, PlotRange->{{0, ptmax}, {0, prmax}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Gray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] "CONFORMAL DIAGRAM, f(η)" cη = Quiet[Plot[ {Piecewise[{ {Piecewise[{{ηH[Ct], Ct<fx rpN}, {1*^6, Ct>fx rpN}}], Ct<rpN}, {Piecewise[{{1*^6, Ct<rpN/fx}, {ηH[ηmax-Ct], Ct>rpN/fx}}], Ct>rpN}}], Ct, π rK/Glyr}, {Ct, 0, ηmax}, Frame->True, AspectRatio->prmax/ηmax, FrameTicks->None, PlotRange->{{0, ηmax}, {0, prmax}}, 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[ {η-Ct, Ct-η}, {Ct, 0, ηmax}, Frame->True, AspectRatio->prmax/ηmax, FrameTicks->None, PlotRange->{{0, ηmax}, {0, prmax}}, 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[plot9[η]], {η, {rpN}}] plot10 = Rasterize[Grid[{{Rotate[Quiet[Plot[ Join[{0}, Table[n, {n, 20, prmax, 20}]], {Ct, 0, ηmax}, Frame->True, AspectRatio->prmax/ηmax, FrameTicks->None, PlotRange->{{0, ηmax}, {0, prmax}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Gray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] s[text_] := Style[text, FontFamily->"Lucida Console", FontSize->36]