(* FLRW/ΛCDM cosmology ||| lcdm.yukterez.net ||| Simon Tyran, Vienna *)

set={"GlobalAdaptive","MaxErrorIncreases"->100,
Method->"GaussKronrodRule"}; (* Integration Rule *)
n=100; (* Recursion Depth *)
int[f_, {x_, xmin_, xmax_}]:=NIntegrate[f,{x, xmin, xmax}, Method->set, MaxRecursion->n]; 
im=640;

c=299792458m/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=30856775777948584200000m; (* Megaparsec *) 
kB=13806488*^-30kg (m^2)/(sek^2)/K; (* Boltzmann Constant *) 
h=662606957*^-42kg (m^2)/sek; (* Planck Constant *) 
ρc[H_]:=3 (H^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/1000K; (* CMB Temperature *) 
H0=67150m/Mpc/sek; (* Hubble Constant *)
kg=m=sek=1; (* SI Units *)

ΩR=92136/1000000000;
ΩM=315/1000; 
ΩΛ=1-ΩM-ΩR;
ΩT=ΩR+ΩM+ΩΛ; 
ΩK=1-ΩT; (* Density Parameters *)

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 *)
rP[t_]:=a[t]int[c/a[т] ,{т,0,t}]; (* Particle Horizon *)
rp[a_]:=a int[c/(A^2)/H[A],{A,0,a}]; (* Particle Horizon *)
rE[t_]:=a[t]int[c/a[т] ,{т,t,∞}]; (* Event Horizon *)
re[a_]:=a int[c/(A^2)/H[A],{A,a,∞}]; (* Event Horizon *)
rL[t0_,t_]:=a[t]int[c/a[т] ,{т,t,t0}]; (* Event Horizon *)
rl[a0_,a_]:=a int[c/(A^2)/H[A],{A,a,a0}]; (* Event Horizon *)
rH[t_]:=c/Evaluate[(A'[t]/A[t]/.sol)[[1]]]; (* Hubble Radius *)
rh[a_]:=c/H[a]; (* Hubble Radius *)

t0=t/.FindRoot[a[t]-1,{t,10 Gyr}];
"t0"->t0/Gyr "Gyr" (* Current time *)

"PROPER DISTANCES, f(t)"

plot1[t_]:=Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[
{((rH[τ Gyr])/Glyr),((rP[τ Gyr])/Glyr),((rE[τ Gyr])/Glyr),((rL[t Gyr,τ Gyr])/Glyr),((-rL[t Gyr,τ Gyr])/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]},{Orange,Thickness[0.005]},{{Orange,Thickness[0.005]},Dashed}},
ImageSize->im,Filling->Top,FillingStyle->Opacity[0.1],ImagePadding->1,GridLines->{{t0/Gyr},{}}]],90 Degree]}}]];
Do[Print[plot1[t]],{t,1t0/Gyr,1t0/Gyr,1t0/Gyr}]

plot2[t_]:=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],
Gray},{n,1,100}],ImageSize->im,ImagePadding->1]],90 Degree]}}]];
Do[Print[plot2[t]],{t,1t0/Gyr,1t0/Gyr,1t0/Gyr}]

"COMOVING DISTANCES, f(t)"

plot3[t_]:=Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[
{((rH[τ Gyr])/(a[τ Gyr]Glyr)),((rP[τ Gyr])/(a[τ Gyr]Glyr)),((rE[τ Gyr])/(a[τ Gyr]Glyr)),((rL[t Gyr,τ Gyr])/(a[τ Gyr]Glyr)),((-rL[t Gyr,τ Gyr])/(a[τ Gyr]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]},{Orange,Thickness[0.005]},{{Orange,Thickness[0.005]},Dashed}},
ImageSize->im,Filling->Top,FillingStyle->Opacity[0.1],ImagePadding->1,GridLines->{{t0/Gyr},{}}]],90 Degree]}}]];
Do[Print[plot3[t]],{t,1t0/Gyr,1t0/Gyr,1t0/Gyr}]

plot4[t_]:=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],
Gray},{n,1,100}],ImageSize->im,ImagePadding->1]],90 Degree]}}]];
Do[Print[plot4[t]],{t,1t0/Gyr,1t0/Gyr,1t0/Gyr}]

"PROPER DISTANCES, f(a)"

plot5[å_]:=Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[
{((rh[α])/Glyr),((rp[α])/Glyr),((re[α])/Glyr),((rl[å,α])/Glyr),((-rl[å,α])/Glyr)}],
{α,0,a[3t0]},Frame->True,AspectRatio->1,PlotRange->{{0,a[3t0]},{0,70}},PlotStyle->{{Thickness[0.005]},
{Darker[Green],Thickness[0.005]},{Purple,Thickness[0.005]},{Orange,Thickness[0.005]},{{Orange,Thickness[0.005]},Dashed}},
ImageSize->im,Filling->Top,FillingStyle->Opacity[0.1],ImagePadding->1,GridLines->{{1},{}}]],90 Degree]}}]];
Do[Print[plot5[å]],{å,1,1,1}]

plot6[å_]:=Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[
Join[{0},Table[α n^(7/2)/250,{n,4,55}]]],
{α,0,a[3t0]},Frame->True,AspectRatio->1,PlotRange->{{0,a[3t0]},{0,70}},PlotStyle->Table[{Dashing->Large,Thickness[0.005],
Gray},{n,1,100}],ImageSize->im,ImagePadding->1]],90 Degree]}}]];
Do[Print[plot6[å]],{å,1,1,1}]

"COMOVING DISTANCES, f(a)"

plot7[å_]:=Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[
{((rh[α])/(α Glyr)),((rp[α])/(α Glyr)),((re[α])/(α Glyr)),((rl[å,α])/(α Glyr)),((-rl[å,α])/(α Glyr))}],
{α,0,a[3t0]},Frame->True,AspectRatio->1,PlotRange->{{0,a[3t0]},{0,70}},PlotStyle->{{Thickness[0.005]},
{Darker[Green],Thickness[0.005]},{Purple,Thickness[0.005]},{Orange,Thickness[0.005]},{{Orange,Thickness[0.005]},Dashed}},
ImageSize->im,Filling->Top,FillingStyle->Opacity[0.1],ImagePadding->1,GridLines->{{1},{}}]],90 Degree]}}]];
Do[Print[plot7[å]],{å,1,1,1}]

plot8[å_]:=Rasterize[Grid[{{Rotate[Quiet[Plot[Evaluate[
Join[{0},Table[n^(7/2)/250,{n,4,55}]]],
{α,0,a[3 t0]},Frame->True,AspectRatio->1,PlotRange->{{0,a[3t0]},{0,70}},PlotStyle->Table[{Dashing->Large,Thickness[0.005],
Gray},{n,1,100}],ImageSize->im,ImagePadding->1]],90 Degree]}}]];
Do[Print[plot8[å]],{å,1,1,1}]

Export["p1.png", plot1[t0/Gyr]]
Export["p2.png", plot2[t0/Gyr]]
Export["p3.png", plot3[t0/Gyr]]
Export["p4.png", plot4[t0/Gyr]]

Export["p5.png", plot5[1]]
Export["p6.png", plot6[1]]
Export["p7.png", plot7[1]]
Export["p8.png", plot8[1]]