%friedmann
% Function which plots the cosmological scale factor a(t) vs time (in
% billions of years) for various solutions to the Friedmann equations.
%
% LAST UPDATED by Andy French April 2017

function friedmann

%Hubble constant at current Epoch /s^-1
%1Mpc = 3.086e22 metres
H0 = 71.9e3/3.086e22;

%A billion years in seconds
Gy = 1e9*365*24*3600;

%Age of Universe (2017 estimate) /billions of years
t_U = 13.799;

%Cosmological constant /m^-2
A = 1.19e-52;

%Speed of light
c = 2.998e8;

%Define maximum time for model /billions of years
tmax = 2*t_U;

%Define number of time values
N = 500;

%Define fontsize for graphs
fsize = 18;

%

%Define time array in billions of years
t = linspace(0, tmax, N );

%Define universe scale factor a = R/R0 where R0 is a spatial distance at
%the current epoch. Note this means a(t0) where t0 is the age of the
%Universe. Define a for different models - note they will predict different
%ages of the Universe given the modern value of H0.

%Linear expansion (Hubble's constant is constant)
t0_linear = (1/H0)/Gy;
a_linear = t/t0_linear;

%Matter dominated
t0_matter = (1/H0)/Gy;
a_matter = ( t/t0_matter ).^(2/3);

%Radiation dominated
t0_rad = (1/H0)/Gy;
a_radiation =  ( t/t0_rad ).^(1/2);
 
% Exponential expansion
a_exp = exp( c*sqrt( A/3 )*(t-t_U)*Gy );

% Plot scale factors against time since Big Bang
figure('name','friedmann','color',[1 1 1])
plot(t,a_linear,t,a_radiation,t,a_matter,t,a_exp);
legend({'Linear','Radiation','Matter','Exponential'},'fontsize',fsize,...
    'location','best');
xlabel('t /billion years','fontsize',fsize);
ylabel('scale factor a','fontsize',fsize);
title('Universe expansion scale factor a','fontsize',fsize)
set(gca,'fontsize',fsize)
grid on
print(gcf,'friedmann.png','-dpng','-r300')
close(gcf);

% Zoom in to current age of Universe to see differences
t = linspace(4, 15, N );

%Define universe scale factor a = R/R0 where R0 is a spatial distance at
%the current epoch. Note this means a(t0) where t0 is the age of the
%Universe. Define a for different models - note they will predict different
%ages of the Universe given the modern value of H0.

%Linear expansion (Hubble's constant is constant)
t0_linear = (1/H0)/Gy;
a_linear = t/t0_linear;

%Matter dominated
t0_matter = (2/3)*(1/H0)/Gy;
a_matter = ( t/t0_matter ).^(2/3);

%Radiation dominated
t0_rad = (1/2)*(1/H0)/Gy;
a_radiation =  ( t/t0_rad ).^(1/2);
 
% Exponential expansion
a_exp = exp( c*sqrt( A/3 )*(t-t_U)*Gy );

% Plot scale factors against time since Big Bang
figure('name','friedmann','color',[1 1 1])
plot(t,a_linear,t,a_radiation,t,a_matter,t,a_exp);
legend({'Linear','Radiation','Matter','Exponential'},'fontsize',fsize,...
    'location','best');
xlabel('t /billion years','fontsize',fsize);
ylabel('scale factor a','fontsize',fsize);
title('Universe expansion scale factor a','fontsize',fsize)
set(gca,'fontsize',fsize)
xlim([4,15])
grid on
print(gcf,'friedmann_zoom.png','-dpng','-r300')
close(gcf);

%End of code