function cosmological_distances

% Light speed (km/s)
c=3*10^(5);
% Hubble constant (km/s/Mpc)
H0=71;
% Factor for LookBack time with Giga year (Gyr)
H0_Gyr=71000/(3*10^(22))*(3600*24*365*10^(9));
% Cosmological Parameters
% Hyperbolic space
Omega1_m=0.9;
Omega1_r=10^-4;
Omega1_l=0;
Omega1_k=1-Omega1_m-Omega1_r-Omega1_l;
% Flat space
Omega2_m=0.3;
Omega2_r=10^-4;
Omega2_k=0;
Omega2_l=1-Omega2_m-Omega2_r;
% Spherical space
Omega3_m=0.3;
Omega3_r=10^-4;
Omega3_l=0.8;
Omega3_k=1-Omega3_m-Omega3_r-Omega3_l;
% Redshift array
z_begin=0;
z_final=1100;
z_step=0.01;
z=z_begin:z_step:z_final;
% Functions for 3 cases into Cosmological Horizon computation
integralCosmoHorizon1=@(x) c/(H0*sqrt(Omega1_k))*sinh(sqrt(Omega1_k)*integral(@(x)CosmoHorizon(x,Omega1_m,Omega1_r,Omega1_l,Omega1_k),z_begin,x)); 
integralCosmoHorizon2=@(x) c/H0*(integral(@(x)CosmoHorizon(x,Omega2_m,Omega2_r,Omega2_l,Omega2_k),z_begin,x)); 
integralCosmoHorizon3=@(x) c/(H0*sqrt(abs(Omega3_k)))*sin(sqrt(abs(Omega3_k))*integral(@(x)CosmoHorizon(x,Omega3_m,Omega3_r,Omega3_l,Omega3_k),z_begin,x));
% Compute cosmological horizon for 3 cases
dc1=arrayfun(integralCosmoHorizon1,z);
dc2=arrayfun(integralCosmoHorizon2,z);
dc3=arrayfun(integralCosmoHorizon3,z);
% Figure 1 : Cosmological Horizon vs Redshift
figure(1);
plot(log10(z),dc1,'b',log10(z),dc2,'r',log10(z),dc3,'k');
box on;
grid on;
% Set label for each curve 
posX = 1;
posX1=posX;
posX2=posX;
posX3=posX;
posY1=5600;
posY2=4500;
posY3=3400;
strMax=['\Omega^{0}_{m} = 0.9 , \Omega^{0}_{\Lambda} = 0 , \Omega^{0}_{k} > 0']; 
text(posX1, posY1, strMax, 'FontSize', 14, 'Color', 'blue', 'VerticalAlignment', 'top');
strMax=['\Omega^{0}_{m} = 0.3 , \Omega^{0}_{\Lambda} = 0.7 - 10^{-4} , \Omega^{0}_{k} = 0'];
text(posX2, posY2, strMax, 'FontSize', 14, 'Color', 'red', 'VerticalAlignment', 'top');
strMax=['\Omega^{0}_{m} = 0.3 , \Omega^{0}_{\Lambda} = 0.8 , \Omega^{0}_{k} < 0'];
text(posX3, posY3, strMax, 'FontSize', 14, 'Color', 'black', 'VerticalAlignment', 'top');
xlabel('Redshift log10(z)');
ylabel('Cosmological Horizon (Mpc)');
title('Cosmological Horizon vs Redshift');
% Compute angular diameter distance from 0 to z_max
% indexes for z<z_max
z_max=10;
z_index=find(z<=z_max);
add1=dc1(z_index)./(z(z_index)+1);
add2=dc2(z_index)./(z(z_index)+1);
add3=dc3(z_index)./(z(z_index)+1);
% Figure 2 : Angular Diameter Distance vs Redshift
figure(2);
plot(z(z_index),add1,'b',z(z_index),add2,'r',z(z_index),add3,'k');
box on;
grid on;
xlim([z_begin z_max]);
% Searching index of max for 3 cases
indexMax1=find(max(add1) == add1);
indexMax2=find(max(add2) == add2);
indexMax3=find(max(add3) == add3);
% Maximum for 3 cases
maxZ1=z(indexMax1);
maxAdd1=add1(indexMax1);
maxZ2=z(indexMax2);
maxAdd2=add2(indexMax2);
maxZ3=z(indexMax3);
maxAdd3=add3(indexMax3);
% Set maximum on figure
posX=2;
posX1=posX;
posX2=posX;
posX3=posX;
posY1=650;
posY2=490;
posY3=330;
strMax=['\Omega^{0}_{m} = 0.9 , \Omega^{0}_{\Lambda} = 0 , \Omega^{0}_{k} > 0 : z_{max} = ', num2str(maxZ1)]; 
text(posX1, posY1, strMax, 'FontSize', 14, 'Color', 'blue', 'VerticalAlignment', 'top');
strMax=['\Omega^{0}_{m} = 0.3 , \Omega^{0}_{\Lambda} = 0.7 - 10^{-4} , \Omega^{0}_{k} = 0 : z_{max} = ', num2str(maxZ2)];
text(posX2, posY2, strMax, 'FontSize', 14, 'Color', 'red', 'VerticalAlignment','top');
strMax=['\Omega^{0}_{m} = 0.3 , \Omega^{0}_{\Lambda} = 0.8 , \Omega^{0}_{k} < 0 : z_{max} = ', num2str(maxZ3)];
text(posX3, posY3, strMax, 'FontSize', 14, 'Color', 'black', 'VerticalAlignment','top');
% Draw lines for maximums
line([maxZ1 maxZ1], [0 maxAdd1], 'Color', 'blue', 'LineStyle', '-.');
line([maxZ2 maxZ2], [0 maxAdd2], 'Color', 'red', 'LineStyle', '-.');
line([maxZ3 maxZ3], [0 maxAdd3], 'Color', 'black', 'LineStyle', '-.');
% Labels
xlabel('Redshift z');
ylabel('Angular Diameter Distance (Mpc)');
title('Angular Diameter Distance vs Redshift');
% Functions for 3 cases into Lookback Time computation
integralLookBack1=@(x) (1/H0_Gyr*integral(@(x)lookBackTime(x,Omega1_m,Omega1_r,Omega1_l,Omega1_k),z_begin,x)); 
integralLookBack2=@(x) (1/H0_Gyr*integral(@(x)lookBackTime(x,Omega2_m,Omega2_r,Omega2_l,Omega2_k),z_begin,x)); 
integralLookBack3=@(x) (1/H0_Gyr*integral(@(x)lookBackTime(x,Omega3_m,Omega3_r,Omega3_l,Omega3_k),z_begin,x)); 
% Compute LookBack Time only up to z=10
z_final=10;
z=z_begin:z_step:z_final;
% Compute Lookback Time for 3 cases
tlb1=arrayfun(integralLookBack1,z);
tlb2=arrayfun(integralLookBack2,z);
tlb3=arrayfun(integralLookBack3,z);
% Figure 3 : LookBack Time vs Redshift
figure(3);
plot(z,tlb1,'b',z,tlb2,'r',z,tlb3,'k');
box on;
grid on;
% Set label for each curve 
posX=3.5;
posX1=posX;
posX2=posX;
posX3=posX;
posY1=6;
posY2=4.8;
posY3=3.6;
strMax=['\Omega^{0}_{m} = 0.9 , \Omega^{0}_{\Lambda} = 0 , \Omega^{0}_{k} > 0']; 
text(posX1, posY1, strMax, 'FontSize', 14, 'Color', 'blue', 'VerticalAlignment', 'top');
strMax=['\Omega^{0}_{m} = 0.3 , \Omega^{0}_{\Lambda} = 0.7 - 10^{-4} , \Omega^{0}_{k} = 0'];
text(posX2, posY2, strMax, 'FontSize', 14, 'Color', 'red', 'VerticalAlignment', 'top');
strMax=['\Omega^{0}_{m} = 0.3 , \Omega^{0}_{\Lambda} = 0.8 , \Omega^{0}_{k} < 0'];
text(posX3, posY3, strMax, 'FontSize', 14, 'Color', 'black', 'VerticalAlignment', 'top');
xlabel('Redshift z');
ylabel('LookBack Time (Gyr)');
title('LookBack Time vs Redshift');
end
% Function to integrate for Cosmological Horizon 
function y = CosmoHorizon(x,Omega_m,Omega_r,Omega_l,Omega_k)
y=(Omega_m*(1+x).^(3)+Omega_r*(1+x).^(4)+Omega_l+Omega_k*(1+x).^(2)).^(-1/2);
end
% Function to integrate for LookBack Time
function y = lookBackTime(x,Omega_m,Omega_r,Omega_l,Omega_k)
y=1./(1+x).*(Omega_m*(1+x).^(3)+Omega_r*(1+x).^(4)+Omega_l+Omega_k*(1+x).^(2)).^(-1/2);
end