FLRW evolution, Milne model

(* 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 *)
   
   t0 = Re[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, 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->{{}, {}}]];
   
   plot1[t_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[Evaluate[
   {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]], 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[Evaluate[
   {rH[τi]/(a[τi]Glyr), rP[τi]/(a[τi]Glyr), rE[τi]/(a[τ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->{{}, {}}]];
   
   plot3[t_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[Evaluate[
   {rL[ti, τi]/(a[τi]Glyr), -rL[ti, τi]/(a[τ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]], 90 Degree]}}]];
   
   Do[Print[plot3[t]], {t, {t0/Gyr}}]
   
   plot4 = Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[
   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]}}]]
   
   "COMOVING DISTANCES, f(η)"
   
   cη = Quiet[Plot[Evaluate[
   {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->{{}, {}}]];
   
   plot5[η_] := Rasterize[Grid[{{Rotate[Quiet[Show[Plot[Evaluate[
   {η-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η]], 90 Degree]}}]];
   
   Do[Print[plot5[η]], {η, {ηmax/2}}]
   
   plot6 = Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[
   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[Evaluate[Join[Table[Sqrt[τ^2-n^2], {n, 7, 70, 7}], {τ}]], 
   {τ, 0, pmax}, Frame->True, AspectRatio->pmax/pmax, FrameTicks->None, 
   PlotRange->{{0, pmax}, {0, pmax}}, 
   PlotStyle->Table[{Black, Thickness[0.0025]}, {n, 1, 10, 1}], 
   ImageSize->im, ImagePadding->1];
   
   plot8 = Plot[V x, {x, 0, pmax}, 
   Frame -> True, AspectRatio -> pmax/pmax, FrameTicks -> None, 
   PlotRange -> {{0, pmax}, {0, pmax}}, 
   PlotStyle -> Table[{Blue, Thickness[0.005]}, {n, 1, 10, 1}], 
   ImageSize -> im, ImagePadding -> 1];
   
   plot9[T_] := Rasterize[Grid[{{Rotate[Quiet[Show[
   Plot[Evaluate[{T-Ct, Ct-T}], {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->{{}, {}}],
   plot7, plot8]], 90 Degree]}}]];
   
   Do[Print[plot9[T]], {T, {t0/Gyr}}]
   
   plot10 = Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[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]}}]]