% Unit 5 (MATLAB Course @ UAM)
% (c) Alexander Knebe

%=========================================
% Unit 5:     the cosmic expansion
%=========================================

%-----------------------------------------
% some general variables
%-----------------------------------------
clear all
close all
clc

kmsecMpc2Gyr = 3.1536/3.086 * 1E-3;

Nsteps = 10000;
t0     = 0;
a0     = 1;


%-----------------------------------------
% backwards integration
%-----------------------------------------

tN = -1.5;


figure(1)
hold on
xlim([tN 0])
ylim([0 1])
title('past expansion')
xlabel('(time-today)*H0')
ylabel('a')



% LCDM model
%-----------------------------------------
Omega0  = 0.28;
OmegaL0 = 0.72;

g  = @(a,t)(sqrt(Omega0./a+OmegaL0*a.^2+(1-Omega0-OmegaL0)));
[a,t]  = rk2(g,t0,a0,tN,Nsteps);

plot(t,a,'r')

%
% SCDM model
%-----------------------------------------
Omega0  = 1.0;
OmegaL0 = 0.0;

g  = @(a,t)(sqrt(Omega0./a+OmegaL0*a.^2+(1-Omega0-OmegaL0)));
[a,t]  = rk2(g,t0,a0,tN,Nsteps);

plot(t,a,'b')


% OCDM model
%-----------------------------------------
Omega0  = 0.24;
OmegaL0 = 0.0;

g  = @(a,t)(sqrt(Omega0./a+OmegaL0*a.^2+(1-Omega0-OmegaL0)));
[a,t]  = rk2(g,t0,a0,tN,Nsteps);

plot(t,a,'g')


% closed CDM model
%-----------------------------------------
Omega0  = 5.0;
OmegaL0 = 0.0;

g  = @(a,t)(sqrt(Omega0./a+OmegaL0*a.^2+(1-Omega0-OmegaL0)));
[a,t]  = rk2(g,t0,a0,tN,Nsteps);

plot(t,a,'c')

legend('LCDM','SCDM','OCDM','closed CDM','Location','NorthWest')


%-----------------------------------------
% forward integration
%-----------------------------------------
tN = 1.5;

figure(2)
hold on
xlim([0 tN])
ylim([1 4])
title('future expansion')
xlabel('(time-today)*H0')
ylabel('a')


% LCDM model
%-----------------------------------------
Omega0  = 0.24;
OmegaL0 = 0.76;

g  = @(a,t)(sqrt(Omega0./a+OmegaL0*a.^2+(1-Omega0-OmegaL0)));
[a,t]  = rk2(g,t0,a0,tN,Nsteps);

plot(t,a,'r')

t(1:10)
a(1:10)
t(Nsteps-10:Nsteps)
a(Nsteps-10:Nsteps)

% SCDM model
%-----------------------------------------
Omega0  = 1.0;
OmegaL0 = 0.0;

g  = @(a,t)(sqrt(Omega0./a+OmegaL0*a.^2+(1-Omega0-OmegaL0)));
[a,t]  = rk2(g,t0,a0,tN,Nsteps);

plot(t,a,'b')


% OCDM model
%-----------------------------------------
Omega0  = 0.24;
OmegaL0 = 0.0;

g  = @(a,t)(sqrt(Omega0./a+OmegaL0*a.^2+(1-Omega0-OmegaL0)));
[a,t]  = rk2(g,t0,a0,tN,Nsteps);

plot(t,a,'g')


% closed CDM model
% ATTENTION: model does not close as \dot{a}=0 appears...
%-----------------------------------------
Omega0  = 5.0;
OmegaL0 = 0.0;

g  = @(a,t)(sqrt(Omega0./a+OmegaL0*a.^2+(1-Omega0-OmegaL0)));
[a,t]  = rk2(g,t0,a0,tN,Nsteps);

plot(t,a,'c')

legend('LCDM','SCDM','OCDM','closed CDM','Location','NorthWest')
