(* Simon Tyran, Vienna. Yukerez 2015 *)

kg=1; m=1; sek=1;  K=1;(*Units *) 
set={"GlobalAdaptive","MaxErrorIncreases"->100,
Method->"GaussKronrodRule"}; (* Integration Rule *)
n=100; (* Recursion Depth *)
tE=300 Gyr;  (* Eventhorizon Limit *) 

c=299792458m/sek; (* Lightspeed *) 
ca=1; (* Perturbation Velocity *)
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=1.69ρR/ρc[H0]; ΩM=317/1000; ΩΛ=683/1000-ΩR; ΩT=ΩR+ΩM+ΩΛ; ΩK=1-ΩT; (* Density Parameters *) 

aE[t_]:=(((((Sqrt[((ΩM)/(ΩΛ))])Sinh[((3 H0 (Sqrt[ΩΛ]))/2)t]))^2)^(1/3)); (* Solving Region *) 
w[a_, w0_]:=(1+w0)((((Sqrt[1+(((ΩΛ)^(-1))-1)(a^(-3))])-(((ΩΛ)^(-1))-1)(a^(-3))((Tanh[(1/(Sqrt[1+(((ΩΛ)^(-1))-1)(a^(-3))]))])^(-1))))^2)((((1/(Sqrt[ΩΛ]))-(((ΩΛ)^(-1))-1)((Tanh[(Sqrt[ΩΛ])])^(-1))))^(-2))-1;

F[a_, w0_]:=(Sqrt[ΩR (a^(-4))+ΩM (a^(-3))+ΩK (a^(-2))+ΩΛ (a^(-3(w[a, w0]+1)))]); (* Density Function by Scalefactor*) 
F[a_, w0_]:=(Sqrt[ΩR (a^(-4))+ΩM (a^(-3))+ΩK (a^(-2))+ΩΛ ]); (* Density Function by Scalefactor, Λ=constant*) 
H[a_, w0_]:=H0 F[a, w0]; (* Hubble Parameter by Scalefactor *) 


int[f_, {x_, xmin_, xmax_}]:=Quiet[NIntegrate[f,{x, xmin, xmax}, Method->set, MaxRecursion->n]];

ta[A_, w0_]:=int[1/a/ H[a, w0], {a,0,A}]; (* Time by Scalefactor *) 
α[τ_, w0_]:=Quiet[A/.FindRoot[ta[A, w0]-τ, {A,1}]] (* Scalefactor by Time *)

rH[τ_, w0_]:=c/H[α[τ, w0], w0]; (* Hubble Radius *)
lC[τ_, w0_]:=int[-c α[τ, w0]/(a^2)/H[a, w0],{a,1,α[τ, w0]}]; (* Light Cone of t0 *) 
Lc[τ_, t_, w0_]:=int[-c α[τ, w0]/(a^2)/H[a, w0],{a,α[t, w0],α[τ, w0]}]; (* Light Cone of t *) 
eH[τ_, w0_]:= α[τ, w0]int[c/(α[time, w0]), {time, τ, tE}]; (* Event Horizon *) 
pH[τ_, w0_]:=int[-α[τ, w0] c/(a^2)/H[a, w0],{a,α[τ, w0],0}]; (* Particle Horizon *) 
g[τ_, w0_]:=tc/.Quiet[FindRoot[pH[tc, w0]/c-τ, {tc,τ}]]; (* Conformal Time *)

ωR[τ_, w0_]:=ΩR ((α[τ, w0])^(-4))/ρc[H[α[τ, w0]]]; (* Radiation Evolution *) 
ωM[τ_, w0_]:=ΩM ((α[τ, w0])^(-3))/ρc[H[α[τ, w0]]]; (* Matter Evolution *) 
ωK[τ_, w0_]:=ΩK ((α[τ, w0])^(-2))/ρc[H[α[τ, w0]]]; (* Curvature Evolution *) 
ωΛ[τ_, w0_]:=ΩΛ ((α[τ, w0])^(-3(w[α[τ, w0], w0]+1)))/ρc[H[α[τ, w]]]; (* Dark Energy Evolution *) 

t0[w0_]:=ta[1, w0]/Gyr; (* Age of the Universe, now *)
"t0 in Gyr"->t0[-1]

plot[w0_]:={
ImageRotate[Quiet[Plot[{((rH[τ Gyr, w0])/Glyr),((lC[τ Gyr, w0])/Glyr),((pH[τ Gyr, w0])/Glyr),((eH[τ Gyr, w0])/Glyr),((lC[τ Gyr, w0])/(-Glyr)),
α[τ Gyr, w0] NestList[#*1.5&,1.5,11]},{τ,0,4 t0[-1]},GridLines->{{t0[w0]},{}},Frame->True,AspectRatio->1,PlotRange->{{0,4 t0[-1]},{0,4 t0[-1]}},PlotStyle->{{Thick},{Orange,Thick},{Thick},{Thick},{{Orange,Thick},Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed}}, ImageSize->500]],Pi/2],

ImageRotate[Quiet[Plot[{((rH[τ Gyr, w0])/(Glyr α[τ Gyr, w0])),((lC[τ Gyr, w0])/(Glyr α[τ Gyr, w0])),((pH[τ Gyr, w0])/(Glyr α[τ Gyr, w0])),((eH[τ Gyr, w0])/(Glyr α[τ Gyr, w0])),((lC[τ Gyr, w0])/(-Glyr α[τ Gyr, w0])), 5, 10, 15, 20, 25, 30, 35, 40, 45, 50,55,60,65,70,75,80,85,90,95,100,105},{τ,0,4t0[-1]},GridLines->{{t0[w0]},{}},Frame->True, AspectRatio->1, PlotRange->{{0,4 t0[-1]},{0,4 t0[-1]}},PlotStyle->{{Thick},{Orange,Thick},{Thick},{Thick},{{Orange,Thick},Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed}, {LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed},{LightGray,Dashed}},ImageSize->500]],Pi/2]
};

Do[Print[plot[w0]],{w0,-1,-1,1}]