kg=1; m=1; sek=1;  K=1;(*Units *) 
set={"GlobalAdaptive","MaxErrorIncreases"->100,
Method->"GaussKronrodRule"}; (* Integration Rule *)
int[f_, {x_, xmin_, xmax_}]:=Quiet[NIntegrate[f,{x, xmin, xmax}, Method->set, MaxRecursion->100]];

c=299792458m/sek; (* Lightspeed *) 
G=667384*^-16(m^3)(kg^(-1))(sek^(-2));  (* Newton Constant *) 
Gyr=(10^7)*36525*24*3600*sek; (* Billion Year Scale *) 
Glyr=Gyr*c; (* Billion Lightyear Scale *) 
Mpc=30856775777948584200000m; (* Megaparsec *) 
kB=13806488*^-30kg (m^2)/(sek^2)/K; (* Boltzmann *) 
h=662606957*^-42kg (m^2)/sek; (* Planck *) 
ρ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 *)
ΩR=ρR/ρc[H0]; ΩM=315/1000; ΩΛ=1-ΩM-ΩR; ΩK=0; (* Density Parameters *)

F[a_]:=(Sqrt[ΩR (a^(-4))+ΩM (a^(-3))+ΩK (a^(-2))+ΩΛ]);  (* Density Function *)
H[a_]:=H0 F[a]; (* Hubble Parameter *)
ta[A_]:=int[1/a/ H[a], {a,0,A}]; (* Time by Scalefactor *) 
α[τ_]:=Quiet[A/.FindRoot[ta[A]-τ, {A,1}]]; (* Scalefactor by Time *)

V[v_,τ_, t_]:=((α[t] c v )/(Sqrt[((α[τ])^2) (c^2)+((α[t])^2) (v^2)-((α[τ])^2) (v^2)]));  (* Peculiar Target Velocity *)
Pc[τ_, t_]:=int[((α[τ])/((a^2)H[a])) V[v,τ, t],{a,α[t],α[τ]}]; (* Particle Trajectory *) 
Lc[τ_, t_]:=int[-c α[τ]/(a^2)/H[a],{a,α[t],α[τ]}]; (* Light Cone *) 
lC[τ_]:=int[-c α[τ]/(a^2)/H[a],{a,1,α[τ]}]; (* Light Cone of t0 *) 
rH[τ_]:=c/H[α[τ]]; (* Hubble Radius *)
eH[τ_]:= α[τ]int[c/(α[κ]), {κ, τ, tE}]; (* Event Horizon *) 
pH[τ_]:=int[-α[τ] c/(a^2)/H[a],{a,α[τ],0}]; (* Particle Horizon *) 
g[τ_]:=tc/.Quiet[FindRoot[pH[tc]/c-τ, {tc,τ}]]; (* Conformal Time *)

ωR[τ_]:=ΩR ((α[τ])^(-4))/ρc[H[α[τ]]]; (* Radiation Evolution *) 
ωM[τ_]:=ΩM ((α[τ])^(-3))/ρc[H[α[τ]]]; (* Matter Evolution *) 
ωK[τ_]:=ΩK ((α[τ])^(-2))/ρc[H[α[τ]]]; (* Curvature Evolution *)

t0=NIntegrate[1/a/ F[a]/H0, {a,0,1}]; 
tE=50 t0;

d0=10 Glyr; (* Comoving Distance *)
Quiet[т/.FindRoot[-Lc[t0+т, t0]/α[t0+т]-d0, {т,t0}]]/Gyr  "Gyr" (* Light Travel Time from today until d0 *)

т=10 Gyr; (* Particle Travel Time *)
v=c/2; (* Initial Velocity *)

V[v,t0+т,t0]/c "c"  (* Peculiar Velocity after т *)

(H[α[t0+т]]Pc[t0+т,t0]+V[v,t0+т,t0])/c  "c"(* Total of Recessional and Peculiar Velocities *)

Pc[t0+т,t0]/Glyr "Glyr" (* Distance travelled *)

c=1; a=2;
Plot[((v(Sqrt[(a^2) (c^2)+ (v^2)-(a^2) (v^2)]))/(c v)), {v,0,c}, Frame->True]

Quit[]