function [t,S,I,C,R] = Program_2_7(beta,epsilon,gamma,Gamma,mu,q,S0,I0,C0,MaxTime)
%
% 
%
% RISK_STRUCTURE( beta, epsilon, gamma, Gamma, mu, q, S0, I0, C0, MaxTime)
%      This is the MATLAB version of program 2.7 from page 44 of 
% "Modeling Infectious Disease in humans and animals" 
% by Keeling & Rohani.
% 
% It is the SICR which includes a carrier class.
%

% Sets up default parameters if necessary.
if nargin == 0
   beta=0.2;
   epsilon=0.1;
   gamma=0.01;
   Gamma=0.001;
   mu=1/(50*365.0);
   q=0.4;
   S0=0.1;
   I0=1e-4;
   C0=1e-3;
   MaxTime=60*365;
end

% Checks all the parameters are valid
if S0<=0 
    error('Initial level of susceptibles (%g) is less than or equal to zero',S0);
end

if I0<=0 
    error('Initial level of infecteds (%g) is less than or equal to zero',I0);
end

if C0<=0 
    error('Initial level of carriers (%g) is less than or equal to zero',C0);
end

if beta<=0 
    error('Transmission rate beta (%g) is less than or equal to zero',beta);
end

if epsilon<=0 
    error('Reduced transmission rate epsilon (%g) is less than or equal to zero',epsilon);
end
    
if gamma<=0 
    error('Recovery rate gamma (%g) is less than or equal to zero',gamma);
end

if Gamma<0 
    error('Recovery rate of carriers Gamma (%g) is less than zero',Gamma);
end

if (q<0 | q>1)
       error('Probability of becoming a carrier q (%g) is not a valid probability',q);
end

if mu<0 
    error('Birth / Death rate gamma (%g) is less than zero',mu);
end
    
if MaxTime<=0 
    error('Maximum run time (%g) is less than or equal to zero',MaxTime);
end
    
if S0+I0+C0>1
    warning('Initial level of susceptibles+infecteds+carriers (%g+%g+%g=%g) is greater than one',S0,I0,C0,S0+I0+C0);
end

if beta/(gamma+mu) + beta*q*epsilon/(Gamma+mu) < 1
    warning('Basic reproductive ratio (R_0=%g) is less than one',beta/(gamma+mu) + beta*q*epsilon/(Gamma+mu));
end

S=S0; I=I0; C=C0; R=1-S-I-C;

% The main iteration 
options = odeset('RelTol', 1e-5);
[t, pop]=ode45(@Diff_2_7,[0 MaxTime],[S I C],options,[beta epsilon gamma Gamma mu q]);

S=pop(:,1); I=pop(:,2); C=pop(:,3); R=1-S-I-C;

% plots the graphs with scaled colours
subplot(3,1,1)
h=plot(t,S,'-g'); 
xlabel 'Time';
ylabel 'Susceptibles'

subplot(3,1,2) 
h=plot(t,I,'-r');
xlabel 'Time';
ylabel 'Infectious'

subplot(3,1,3)
h=plot(t,C,'-m'); 
xlabel 'Time';
ylabel 'Carriers'


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

beta=parameter(1); epsilon=parameter(2); 
gamma=parameter(3); Gamma=parameter(4);
mu=parameter(5); q=parameter(6);

S=pop(1); I=pop(2); C=pop(3);

dPop=zeros(3,1);

dPop(1)= mu -beta*S*(I + epsilon*C) - mu*S;
dPop(2)= beta*S*(I+ epsilon*C) - gamma*I - mu*I;
dPop(3)= gamma*q*I - Gamma*C - mu*C;



