🗎 Code for the
Milne universe spacetime diagrams
(* Milne Model | FLRW coordinates | infinite particle horizon and event horizon *) set = {"GlobalAdaptive", "MaxErrorIncreases"->100, Method->"GaussKronrodRule"}; (* Integration Rule *) if[F_] := If[ΩΛ<=0, Nothing, If[ΩK>0, Nothing, F]] (* Event Horizon Check *) n = 100; (* Recursion Depth *) int[f_, {x_, xmin_, xmax_}] := (* Integral *) NIntegrate[f, {x, xmin, xmax}, Method->set, MaxRecursion->n, WorkingPrecision->wp]; wp = MachinePrecision; (* Working Precision *) im = 200; (* Image Size *) ηmax = 70; pmax = 70; (* 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 = 0; (* Radiation Proportion including Neutrinos *) ΩM = 0; (* Matter Proportion including Dark Matter *) ΩΛ = 0; (* Dark Energy Proportion *) ΩT = 0; (* Total Density over Critical Density *) ΩK = 1-ΩT; (* Curvature Density *) H0 = 67150 m/Mpc/sek; (* Hubble Constant *) H[a_] := H0/a (* Hubble Parameter *) a[t_] := H0 t; (* Scale Factor a by Time t *) т[a_] := a/H0; (* Time t by Scale Factor a *) rP[t_] := Nothing; (* Proper Particle Horizon by t *) rp[a_] := Nothing; (* Proper Particle Horizon by a *) RP[t_] := Nothing; (* Comoving Particle Horizon by t *) Rp[a_] := Nothing; (* 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/H[a[t]]; (* Proper Hubble Radius by t *) rh[a_] := c/H[a]; (* Proper Hubble Radius by a *) RH[t_] := c/H[a[t]]/a[t]; (* Comoving Hubble Radius by t *) Rh[a_] := c/H[a]/a; (* Comoving Hubble Radius by a *) rT[t0_, t_] := c t (ArcCosh[t0/t]); (* Proper T Hypersurface by t *) rt[a0_, a_] := c т[a] (ArcCosh[t0/т[a]]); (* Proper T Hypersurface by a *) RT[t0_, t_] := c t (ArcCosh[t0/t])/a[t]; (* Comoving T Hypersurface by t *) Rt[a0_, a_] := c т[a] (ArcCosh[t0/т[a]])/a; (* Comoving T Hypersurface by a *) t0 = Re[t/.FindRoot[a[t]-1, {t, 10 Gyr}]]; ti = t Gyr; τi = τ Gyr; "t0"->t0/Gyr "Gyr" (* Current Time *) LT = Quiet[Interpolation[ (* conformal Hypersurface of constant T *) Join[{{0, 0}}, ParallelTable[{RL[t0, t Gyr]/Glyr, RT[t0, t Gyr]/Glyr}, {t, t0/Gyr/im, t0/Gyr, t0/Gyr/im}]]]]; "PROPER DISTANCES, f(t)" pt = Quiet[Plot[ {rH[τi]/Glyr, rP[τi]/Glyr, rE[τi]/Glyr}, {τ, 0, pmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->{{Thickness[0.005]}, {Darker[Green], Thickness[0.005]}, {Purple, Thickness[0.005]}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; pT = Quiet[Plot[ {rT[t0, τi]/Glyr}, {τ, 0, pmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->{{Red, Thickness[0.005]}}, ImageSize->im, Filling->None, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; plot1[t_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[ {rL[ti, τi]/Glyr, -rL[ti, τi]/Glyr}, {τ, 0, pmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->{{Orange, Thickness[0.005]}, {{Orange, Thickness[0.005]}, Dashed}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}], pt, pT]], 90 Degree]}}]]; Do[Print[plot1[t]], {t, {t0/Gyr}}] plot2 = Rasterize[Grid[{{Rotate[Quiet[ Plot[Evaluate[ Join[{0}, Table[a[τ Gyr] 10 Tan[n], {n, π/20, π/2, π/20}], {a[τ Gyr] 1000}]], {τ, 0, pmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, 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, RE[τi]/Glyr}, {τ, 0, pmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->{{Thickness[0.005]}, {Darker[Green], Thickness[0.005]}, {Purple, Thickness[0.005]}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; cT = Quiet[Plot[ {RT[t0, τi]/Glyr}, {τ, 0, pmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->{{Red, Thickness[0.005]}}, ImageSize->im, Filling->None, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; plot3[t_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[ {RL[ti, τi]/Glyr, -RL[ti, τi]/Glyr}, {τ, 0, pmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->{{Orange, Thickness[0.005]}, {{Orange, Thickness[0.005]}, Dashed}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}], ct, cT]], 90 Degree]}}]]; Do[Print[plot3[t]], {t, {t0/Gyr}}] plot4 = Rasterize[Grid[{{Rotate[Quiet[Plot[ Join[{0}, Table[n, {n, 10, pmax, 10}]], {τ, 0, pmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Gray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] "CONFORMAL DIAGRAM, f(η)" cη = Quiet[Plot[ {c/H0/Glyr}, {Ct, 0, ηmax}, Frame->True, AspectRatio->pmax/ηmax, FrameTicks->None, PlotRange->{{0, ηmax}, {0, pmax}}, PlotStyle->{{Thickness[0.005]}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; cN = Quiet[Plot[ {LT[ηmax/2-Ct]}, {Ct, 0, ηmax}, Frame->True, AspectRatio->pmax/ηmax, FrameTicks->None, PlotRange->{{0, ηmax}, {0, pmax}}, PlotStyle->{{Red, Thickness[0.005]}}, ImageSize->im, Filling->None, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]; plot5[η_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[ {η-Ct, Ct-η}, {Ct, 0, ηmax}, Frame->True, AspectRatio->pmax/ηmax, FrameTicks->None, PlotRange->{{0, ηmax}, {0, pmax}}, PlotStyle->{{Orange, Thickness[0.005]}, {{Orange, Thickness[0.005]}, Dashed}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}], cη, cN]], 90 Degree]}}]]; Do[Print[plot5[η]], {η, {ηmax/2}}] plot6 = Rasterize[Grid[{{Rotate[Quiet[Plot[ Join[{0}, Table[n, {n, 10, pmax, 10}]], {Ct, 0, ηmax}, Frame->True, AspectRatio->pmax/ηmax, FrameTicks->None, PlotRange->{{0, ηmax}, {0, pmax}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Gray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] "MINKOWSKI COORDINATES, x(T)" int[f_, {x_, xmin_, xmax_}] := Integrate[f, {x, xmin, xmax}]; Z=z/.Solve[H[1/(1+z)]rl[1, 1/(1+z)]/c==1, z][[1]]; V=v/.Solve[1/(Z+1)==Sqrt[(1-v)/(1+v)], v][[1]]; plot7 = Plot[Join[Table[Sqrt[τ^2-n^2], {n, t0/Gyr/2, 7 t0/Gyr, t0/Gyr/2}], {τ}], {τ, 0, pmax}, Frame->True, AspectRatio->pmax/pmax, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->Table[{LightGray, Thickness[0.0025]}, {n, 1, 10, 1}], ImageSize->im, ImagePadding->1]; plot8 = Plot[Sqrt[τ^2-(t0/Gyr)^2], {τ, 0, pmax}, Frame->True, AspectRatio->pmax/pmax, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->{{Red, Thickness[0.005]}}, ImageSize->im, ImagePadding->1]; plot9 = Plot[V x, {x, 0, pmax}, Frame -> True, AspectRatio -> pmax/pmax, FrameTicks -> None, PlotRange -> {{0, pmax}, {0, pmax}}, Filling->Top, FillingStyle->Opacity[0.1], PlotStyle -> Table[{Thickness[0.005]}, {n, 1, 10, 1}], ImageSize -> im, ImagePadding -> 1]; plot10[T_] := Rasterize[Grid[{{Rotate[Quiet[Show[ plot7, plot8, plot9, Plot[{T-Ct, Ct-T, Ct}, {Ct, 0, ηmax}, Frame->True, AspectRatio->pmax/ηmax, FrameTicks->None, PlotRange->{{0, ηmax}, {0, pmax}}, PlotStyle->{{Orange, Thickness[0.005]}, {{Orange, Thickness[0.005]}, Dashed}, {Darker[Green], Thickness[0.005]}}, ImageSize->im, Filling->Top, FillingStyle->Opacity[0.1], ImagePadding->1, GridLines->{{}, {}}]]], 90 Degree]}}]]; Do[Print[plot10[T]], {T, {t0/Gyr}}] plot11 = Rasterize[Grid[{{Rotate[Quiet[Plot[Join[{0}, Table[τ ArcTanh[Tanh[1]n/pmax], {n, 7, 70, 7}]], {τ, 0, pmax}, Frame->True, AspectRatio->pmax/pmax, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, PlotStyle->Table[{Dashing->Large, Thickness[0.005], Gray}, {n, 1, 100}], ImageSize->im, ImagePadding->1]], 90 Degree]}}]] plot12 = Rasterize[Grid[{{Rotate[Quiet[Plot[ Join[{0}, Table[τ Tanh[n], {n, 1/10, 10, 1/10}]], {τ, 0, pmax}, Frame->True, AspectRatio->1, FrameTicks->None, PlotRange->{{0, pmax}, {0, pmax}}, 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]