1-1(* | 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]; 
   im = 640;                                                         (* Image Size *)
   
   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 = NDSolve[{A'[t]/A[t] == H[A[t]], A[0] == 1*^-15}, 
   A, {t, 0, 500 Gyr}, Method->{"ImplicitRungeKutta", "DifferenceOrder"->20}, 
   MaxSteps->Infinity, WorkingPrecision->MachinePrecision]; 
   
   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/Evaluate[(A'[t]/A[t]/.sol)[[1]]];             (* Hubble Radius by t *)
   rh[a_] := c/H[a];                                         (* Hubble Radius by a *) 
   
   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->{{t0/Gyr}, {}}]];
   
   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->{{t0/Gyr}, {}}]];
   
   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->{{1}, {}}]];
   
   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->{{1}, {}}]];
   
   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]}}]]