function [T,S,E,I] = Program_3_4(beta,sigma,gamma,S0,E0,I0, MaxTime)
%
% 
%
% Program_3_4(beta,sigma,gamma,S0,E0,I0, MaxTime)
%      This is the MATLAB version of program 3.4 from page 87 of 
% "Modeling Infectious Disease in humans and animals" 
% by Keeling & Rohani.
% 
% It is the SEIR model with four different age-groups and 
% yearly "movements" between the groups mimicking the school year
%
% S0, E0 and I0 are all vectors.


% Sets up default parameters if necessary.
if nargin == 0
   beta=[2.089 2.089 2.086 2.037; 2.089 9.336 2.086 2.037; 2.086 2.086 2.086 2.037; 2.037 2.037 2.037 2.037];
   sigma=1/8;
   gamma=1/5;
   S0=[0.05 0.01 0.01 0.008];
   E0=[0.0001 0.0001 0.0001 0.0001];
   I0=[0.0001 0.0001 0.0001 0.0001];
   MaxTime=100*365;
end

nu(1)=1/(55*365); nu(2:4)=0;
mu(4)=1/(55*365); mu(1:3)=0;

n=[6 4 10 55]/75;

% Checks all the parameters are valid

for a=1:4
    if I0(a)<0
        error('Initial level of infectious in age-group %d (%g) is less than zero',a,I0(a));
    end
    if E0(a)<0
        error('Initial level of exposed in age-group %d (%g) is less than zero',a,E0(a));
    end
    if S0(a)<0
        error('Initial level of susceptibles in age-group %d (%g) is less than zero',a,S0(a));
    end
    if S0(a)+E0(a)+I0(a)>n(a)
        error('Initial level of infecteds + susceptibles in age-group %d (%g) is greater than calculated group size (%g)',a,I0(a)+E0(a)+S0(a),n(a));
    end
end

if length(find(beta<0))>0 
    error('Transmission rate beta has a term that is less than zero');
end
    
if sigma<=0 
    error('Exposed to Infectious rate sigma (%g) is less than or equal to zero',sigma);
end

if gamma<=0 
    error('Recovery rate gamma (%g) is less than or equal to zero',gamma);
end
    
if MaxTime<=0 
    error('Maximum run time (%g) is less than or equal to zero',MaxTime);
end
    

% The main iteration 
options = odeset('RelTol', 1e-3);
T0=0; S=[]; E=[]; I=[]; T=[];
while T0<MaxTime
  [t, pop]=ode45(@Diff_3_4,[T0 (T0+365)],[S0 E0 I0],options,[reshape(beta,1,16) sigma gamma reshape(nu,1,4) reshape(mu,1,4)]);
  T=[T; t]; S=[S; pop(:,1:4)]; E=[E; pop(:,5:8)]; I=[I; pop(:,9:12)];
  
  % Now do the yearly movement bit
  
  S0(1)=S(end,1)-S(end,1)/6;
  S0(2)=S(end,2)+S(end,1)/6 - S(end,2)/4;
  S0(3)=S(end,3)+S(end,2)/4 - S(end,3)/10;
  S0(4)=S(end,4)+S(end,3)/10;
  
  E0(1)=E(end,1)-E(end,1)/6;
  E0(2)=E(end,2)+E(end,1)/6 - E(end,2)/4;
  E0(3)=E(end,3)+E(end,2)/4 - E(end,3)/10;
  E0(4)=E(end,4)+E(end,3)/10;
  
  I0(1)=I(end,1)-I(end,1)/6;
  I0(2)=I(end,2)+I(end,1)/6 - I(end,2)/4;
  I0(3)=I(end,3)+I(end,2)/4 - I(end,3)/10;
  I0(4)=I(end,4)+I(end,3)/10;
    
  T0=T(end);
end
 
  
% plots the graphs with scaled colours
subplot(3,1,1)
h=plot(T/365,S,'-g'); 
for i=1:4
    set(h(i),'Color',[0 0.2+0.2*i 0]);
end
legend(h,'0-6','6-10','10-20','20+');
xlabel 'Time (years)';
ylabel 'Susceptible'

subplot(3,1,2) 
h=semilogy(T/365,I,'-r');
for i=1:4
    set(h(i),'Color',[0.2+0.2*i 0 0]);
end
legend(h,'0-6','6-10','10-20','20+');
xlabel 'Time (years)';
ylabel 'Infectious'

subplot(3,1,3);
m=find(T>T(end)-365);
R=1-mean(S(m,:))./n;
fill([0 0 6 6  6 6 10 10 10 10 20 20 20 20 75 75],[0 R(1) R(1) 0 0 R(2) R(2) 0 0 R(3) R(3) 0 0 R(4) R(4) 0],'r');
axis([0 25 0 1]);
xlabel 'Age-group'
ylabel 'Proportion Sero-positive'


% Calculates the differential rates used in the integration.
function dPop=Diff_3_4(t, pop, parameter)

beta=reshape(parameter(1:16),4,4);  sigma=parameter(17); gamma=parameter(18); 
nu=reshape(parameter(19:22),1,4); mu=reshape(parameter(23:26),1,4);

S=pop(1:4); E=pop(5:8); I=pop(9:12); 

dPop=zeros(12,1);

for a=1:4
    Inf=(beta(a,:)*I)*S(a);
    dPop(a)=nu(a)*55/75 - Inf - mu(a)*S(a);
    dPop(a+4)=Inf - sigma*E(a) - mu(a)*E(a);
    dPop(a+8)=sigma*E(a) - gamma*I(a) - mu(a)*I(a);
end