% populations.m
% Brad Poon & Les Fletcher
% Math 164 Spring 2003
% Professor de Pillis
% Last Modified: April 30, 2002
%
% DESCRIPTION:
% A model of population based on 4 age groups:  young(0-9), adolescent(10-19),
% adult(20-64), and elderly (65-death).
%

function yprime = populations(t,y)

% Birth and Death Rates

a = 0.1144;     % Birth rate of 'Young' age group
b = 0.005134;   % Rate at which population leaves 'Young' age group (through death)
c = .1;         % Rate at which population leaves 'Young' age group (through acension)

d = 0.001026;   % Rate at which population leaves 'Adolescent' age group (through death)
e = .1;         % Rate at which population leaves 'Adolescent' age group (through acension)

f = 0.03938;    % Rate at which population leaves 'Adults' age group (through death)
g = .022;       % Rate at which population leaves 'Adults' age group (through acension)


h = 0.32649;    % Rate at which population leaves 'Elderly' age group (through death)

baby_boom_rate = 1;   % Multiplier for the birth rate to simulate a baby boom

if (t > 100 & t < 105) a = a*baby_boom_rate; end


% Production and Consumption Rates (per person)

ry1 = 77;         % Resource production rate of 'Young' and 'Adolescent" age group
ry2 = 103.6;      % Resource consumption rate of 'Young' and 'Adolescent" age group

rz1 = 820.1;      % Resource production rate of 'Adults" age group
rz2 = 774.3;      % Resource consumption rate of 'Adults" age group

re1 = 102.1;      % Resource production rate of 'Elderly" age group
re2 = 122.1;      % Resource consumption rate of 'Elderly" age group


% Critical resource values, amount of resources needed, amount resources available

rc = ry2*(y(1)+y(2)) + rz2*y(3) + re2*y(4);   % Resources consumed
ra = y(5);                                    % Resources available


% Modified birth / death rates due to lack of (or abundance of) resources available

if (rc/ra < 1)
    sensitivity_death = 1;
    max_change_death = .3183;
else
    sensitivity_death = .1;
    max_change_death = 10;
end

if (ra/rc < 1)
    max_change_birth = 1.2732;   % -1 / atan(-1)
else
    max_change_birth = 1/(2*pi);
end

birth_modifier = max_change_birth*atan((ra/rc)-1)+1;
death_modifier = max_change_death*atan(sensitivity_death*((rc/ra)-1))+1;

% Modify the birth and death rates

a = a*(birth_modifier);
b = b*(death_modifier);
d = d*(death_modifier);
f = f*(death_modifier);
h = h*(death_modifier);


% Population model with resource consumption      
 
yprime = [a*y(3) - (b + (1-b)*c)*y(1); 
          (1-b)*c*y(1)-(d+(1-d)*e)*y(2); 
          (1-d)*e*y(2)-(f+(1-f)*g)*y(3); 
          (1-f)*g*y(3)-h*y(4);
          ((ry1-ry2)*(y(2)+y(1)))+(rz1-rz2)*y(3)+(re1-re2)*y(4)];