function [t,NSS,NIS,NSI,NRS,NSR,NRI,NIR,NRR] = Program_4_1(beta,gamma,mu,alpha,a,MaxTime)
%
% 
%
% PROGRAM_4_1( beta, gamma, mu, alpha a, MaxTime)
%      This is the MATLAB version of program 4.1 from page 118 of 
% "Modeling Infectious Disease in humans and animals" 
% by Keeling & Rohani.
% 
% It is a 2 strain SIR disease model, where parameter alpha governs partial
% susceptibility and parameter a governs partial transmissibility. Note
% that in theory alpha or a could be greater than one to imply enhanced
% susceptibility or transmissibility.
%

% Sets up default parameters if necessary.
if nargin == 0
   beta=[260 520]/365.0;
   gamma=[1 1]/7.0;
   mu=1/(70*365.0);
   alpha=[0.5 0.4];
   a=[0.4 0.5];
   MaxTime=100*365;
end

% Check sizes for terms that are 'vectors'
beta=CheckSize(beta,2,1,'beta');
gamma=CheckSize(gamma,2,1,'gamma');
alpha=CheckSize(alpha,2,1,'alpha');
a=CheckSize(a,2,1,'a');


% Checks all the parameters are valid
CheckGreaterOrEqual(beta,0,'beta');
CheckGreaterOrEqual(gamma,0,'gamma');
CheckGreaterOrEqual(alpha,0,'alpha');
CheckGreaterOrEqual(a,0,'a');
CheckGreaterOrEqual(mu,0,'mu');
CheckGreaterOrEqual(MaxTime,0,'MaxTime');

% initial conditions are set to suit default parameters
NSS=0.1; NSI=1e-4; NIS=1e-4; NSR=0.5; NRS=0.02; NRI=0; NIR=0;

% The main iteration 
options = odeset('RelTol', 1e-5);
[t, pop]=ode45(@Diff_4_1,[0 MaxTime],[NSS NIS NRS NRI NSI NSR NIR],options,beta,gamma,mu,alpha,a);

NSS=pop(:,1); NIS=pop(:,2); NRS=pop(:,3); NRI=pop(:,4);
NSI=pop(:,5); NSR=pop(:,6); NIR=pop(:,7); NRR=1-sum(pop')';

% plots the graphs with scaled colours
subplot(3,1,1)
h=plot(t,[NSS NRS NSR NRR],'-'); 
xlabel 'Time';
ylabel 'Uninfected'
legend(h,'N_{SS}','N_{RS}','N_{SR}','N_{RR}');

subplot(3,1,2) 
h=plot(t,[NIS NIR (NIS+a(1)*NIR)],'-');
xlabel 'Time';
ylabel 'Infectious 1'
legend(h,'N_{IS}','N_{IR}','I_1');

subplot(3,1,3) 
h=plot(t,[NSI NRI (NSI+a(2)*NRI)],'-');
xlabel 'Time';
ylabel 'Infectious 2'
legend(h,'N_{SI}','N_{RI}','I_2');


% Calculates the differential rates used in the integration.
function dPop=Diff_4_1(t, pop, beta,gamma,mu,alpha,a)

NSS=pop(1); NIS=pop(2); NRS=pop(3); NRI=pop(4);
NSI=pop(5); NSR=pop(6); NIR=pop(7); NRR=1-sum(pop);

lambda1=beta(1)*(NIS+a(1)*NIR);  lambda2=beta(2)*(NSI+a(2)*NRI);

dPop=zeros(7,1);

dPop(1)= mu - NSS*(lambda1 + lambda2)  - mu*NSS;

dPop(2)= NSS*lambda1 - gamma(1)*NIS - mu*NIS;
dPop(3)= gamma(1)*NIS - alpha(2)*NRS*lambda2 - mu*NRS;
dPop(4)= alpha(2)*NRS*lambda2 - gamma(2)*NRI - mu*NRI;

dPop(5)= NSS*lambda2 - gamma(2)*NSI - mu*NSI;
dPop(6)= gamma(1)*NSI - alpha(1)*NSR*lambda1 - mu*NSR;
dPop(7)= alpha(1)*NSR*lambda1 - gamma(1)*NIR - mu*NIR;


% Does a simple check on the size and possible transpose
function Parameter=CheckSize(Parameter, L, W, str)

[l w]=size(Parameter);

if(l==W & w==L)
    Parameter=Parameter'; [l w]=size(Parameter);
end

if(l==1 & w==1)
    warning('Parameter %s is a scaler value, expanding to size %d x %d',str,L,W);
    Parameter=Parameter*ones(L,W);
    [l w]=size(Parameter);
end

if(l~=L | w~=W)
    error('Parameter %s is of size %d x %d and not %d x %d',str,l,w,L,W);
end


% Does a simple check on the value
function []=CheckGreaterOrEqual(Parameter, Value, str)

m=find(Parameter<Value);
if length(m)>0
    error('Parameter %s(%g) (=%g) is less than %g',str,m(1),Parameter(m(1)),Value);
end