%%  Expansion des Universums

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016
% Vers.2.1 15.09.2016

% Parameter
H0 = 71;                           % Hubbelkonstante
omega_m = 0.3;                     % Massendichte
omega_la = 0.7;                    % Vakuumenergiedichte
c = 299792.458;                    % km/s
%% Abstand der Galaxien; heute und zum Zeitpunkt der Lichtemission
Hzint = @(z) c./(H0.*sqrt(omega_m.*(1+z).^3 + omega_la));
z = logspace(-3,3,100);
for n=0:length(z)-1
    n=n+1;
    q(n) = -integral(Hzint,z(n),0)*3.26/1000;
end
semilogx(z,q)
hold on
semilogx(z,q./(1+z));
shg
legend('r_0','r_e')
xlabel('z'), ylabel('r [10^9 Lj]')
%% Abb. 1.43
Haint = @(a) 1./(a.* H0.*sqrt(omega_m./a.^3 + omega_la));
t0 = integral(Haint,0,1)*1000;
for n=0:length(z)-1
    n=n+1;
    ag = 1/(1+z(n));
    t(n) = integral(Haint,0,ag)*1000;
end
figure
plot(q./(1+z),t),shg
hold on
np = 40:10:70;
for nl=np
    plot([q(nl)./(1+z(nl)), q(nl)],[t(nl),t0],'k--','LineWidth',1)
    text(q(nl),t0,num2str(round(z(nl)*100)/100))
end
x0 = q(60)./(1+z(60));
x1 = x0 + t0-t(60);
plot([x0,x1],[t(60),t0]),shg
xlabel('r'),ylabel('t')
% x0 = q(50)./(1+z(50));
% x1 = x0 + t0-t(50);
% plot([x0,x1],[t(50),t0]),shg
hold off


