##================================================================
## A SOLUTION TO PROJECT 2
##================================================================
## (C) Marek Rychlik (rychlik@alamos.math.arizona.edu), 2003
## All rights reserved
##================================================================
## If x is the population distribution today, and
## y is the pop. distribution a year from now then
## the equations are:
## y(1) = sum f(j) * x(j) over j = 1:100
## y(2) = s(1)*x(1)
## y(3) = s(2)*x(2)
## y(100) = s(100)*x(100)
##================================================================

## Set these values to something reasonable
n = 100;			# Maximum age
m = 25;				# Maximally likely age to have a child
p = 5;				# Power in morbidity

## fertility = exp(-.5 * ( (log(1:n) - log(m)) / .2 ) .^ 2);

## This is the true log-normal distribution
fertility = exp(-.5 * ( (log(1:n) - log(m)) / .2 ) .^ 2 - log(1:n));

c = 1 / sum(fertility);

morbidity = ((1 : n) / n) .^ p;

survival = 1 - morbidity(1:n-1);

trans = [c * fertility;diag(survival), zeros(n - 1, 1)];

[S, D] = eig(trans);

## Find index of eigenvalue 1
tol = 1e-2;
ind = 0;
for i=1:n 
  if( abs(D(i,i)-1) < tol && ind == 0)
    ind = 1;
  endif
endfor;

if(ind == 0)
  '1 is not an eigenvalue';
endif

gset terminal x11;

clearplot;

##title('Age distribution (scaled) depending on fertility and morbidity');
title(strcat('Age distribution - max. age = ', num2str(n), \
		    ', most likely reprod. age = ', num2str(m), \
		    ', morbidity power = ', num2str(p)));

xlabel('Age');

plot(fertility /max(fertility), '-;fertility;',
     morbidity, '-;morbidity;',
     abs(S(:, ind)) / max(abs(S(:, ind))), '-;population;');


## Now this code will send output to a PostScript file
gset terminal postscript color;
gset output 'solution2.ps';

replot;

## ... And this to the PNG file good for the Web
gset terminal png color;
gset output 'solution2.png';

replot;

## Population distribution
population_dist  = S(:,1)/sum(S(:,1));

average_age = (1:n) * population_dist

lifespan_dist = morbidity' .* population_dist; # Distribution of lifespans

average_lifespan = (1:n) * (lifespan_dist / sum(lifespan_dist))
fraction_of_seniors = sum(population_dist(floor(.55 .* n):n))

## ;;; Local Variables: ***
## ;;; mode: octave ***
## ;;; End: ***
