% Lecture Optimal Control and Estimation 2015
% Exercise 1
% Script by Greg Horn and Jonas Koenemann

clear all

% setup constants
n = 100; % max age
ages = (1:n)';
N = 150; % years to simulate

[beta, mu] = birth_mortality_rates();
 
figure(1)
clf
subplot(121)
bar(beta*100)
ylabel('birth rate [%]')
xlabel('age [years]')
subplot(122)
bar(mu*100)
ylabel('mortality rate [%]')
xlabel('age [years]')


A = zeros(n,n);
A(1,:) = beta';
for k=1:n
    A(2:n,1:n-1) = diag(1 - mu);
end

% eigenvalues of A
[v,lam] =eig(A);
lam = diag(lam);

% plot eigenvalues of A
figure(2)
clf
subplot(121)
plot(real(lam), imag(lam),'.')
hold on
plot(sin(0:0.01:2*pi), cos(0:0.01:2*pi),'r')
title('all eigenvalues')

% plot eigenvectors of A
[~, lam_order] = sort(abs(lam),'descend');
subplot(122)
hold on
nbig = 2;
plot(repmat(1:n,nbig,1)',real(v(:,lam_order(1:nbig))))
xlabel('age')
legend('biggest eigenvector','second biggest eigenvector')

% simulation with impulse response
% simulate!
x0 = zeros(n,1);
x0(5) = 100;
X = zeros(n,N);
X(:,1) = x0;
for k=1:N-1
    X(:,k+1) = A*X(:,k);
end
% plot!
figure(3)
clf
hold on
for k=1:N
  plot3(k*ones(n,1), (1:n)', X(:,k))
end
grid on
xlabel('year')
ylabel('age')
zlabel('population')

% simulation of eigenvectors w/ biggest eigenvalues
% simulate!
for eigk=[1,2,4]
    x0 = real(v(:,lam_order(eigk)));
    X = zeros(n,N);
    X(:,1) = x0;
    for k=1:N-1
        X(:,k+1) = A*X(:,k);
    end
    % plot!
    figure(3+eigk)
    clf
    hold on
    for k=1:N
      plot3(k*ones(n,1), (1:n)', X(:,k))
    end
    xlabel('year')
    ylabel('age')
    zlabel('population')
    title(sprintf('eigenvector %d',eigk))
    grid on
end
