function radius_universe

c=3*10^(5);
H0=71;

Omega1_m=0.9;
Omega1_red=0;
Omega1_k=1-Omega1_m-Omega1_red;
Omega2_m=0.3;
Omega2_red=0.7;
Omega2_k=1-Omega2_m-Omega2_red;
Omega3_m=0.3;
Omega3_red=0.8;
Omega3_k=1-Omega3_m-Omega3_red;

z_begin=0;
z_final=1100;
inter=1;

j=1;
for z=z_begin:inter:z_final
dc1(j)=c/(H0*sqrt(abs(Omega1_k)))*sinh(sqrt(abs(Omega1_k))*quad(@(x)myfunc(x,Omega1_m,Omega1_red,Omega1_k),z_begin,z));
j=j+1;
end

j=1;
for z=z_begin:inter:z_final
dc2(j)=c/H0*(quad(@(x)myfunc(x,Omega2_m,Omega2_red,Omega2_k),z_begin,z));
j=j+1;
end

j=1;
for z=z_begin:inter:z_final
dc3(j)=c/(H0*sqrt(abs(Omega3_k)))*sin(sqrt(abs(Omega3_k))*quad(@(x)myfunc(x,Omega3_m,Omega3_red,Omega3_k),z_begin,z));
j=j+1;
end

figure(1);
plot(z_begin:inter:z_final,dc1,'b',z_begin:inter:z_final,dc2,'r',z_begin:inter:z_final,dc3,'k');
box on;
grid on;

xlabel('Redshift z');
ylabel('Cosmic Horizon (Mpc)');


%-------------------%
function y = myfunc(x,Omega_m,Omega_red,Omega_k)
y=(Omega_m*(1+x).^(3)+Omega_red+Omega_k*(1+x).^(2)).^(-1/2);
end
%-------------------%

end