function [T,S,I] = Program_8_4(beta,gamma,mu,p,S,I, tV, MaxTime)
%
% 
%
  % Program_8_4( beta, gamma, mu, p, S I, tV, MaxTime)
%      This is the MATLAB version of program 8.4 from page 306 of 
% "Modeling Infectious Disease in humans and animals" 
% by Keeling & Rohani.
% 
% It is the SIR model with two different risk-groups, each group ha
% an associated birth and vaccination rate.
% 
% Note, gamma, mu, p, S and I are all vectors. beta is a matrix.
% tV is the time at which vaccination starts.


% Sets up default parameters if necessary.
if nargin == 0
   beta=[1 0.01; 0.01 0.1];
   gamma=[0.1 0.1];
   mu=[0.2 0.8]*5e-5;
   p=[0.4 0.1];
   tV=50*365;
   S=[0.1 0.7];
   I=[1e-5 1e-5];
   MaxTime=100*365;
end

% Check all parameters are the right size
beta=CheckSize(beta,2,2,'beta');
gamma=CheckSize(gamma,1,2,'gamma');
mu=CheckSize(mu,1,2,'mu');
p=CheckSize(p,1,2,'p');
S=CheckSize(S,1,2,'S');
I=CheckSize(I,1,2,'I');

% Checks all the parameters are valid
CheckGreaterOrEqual(beta,0,'beta');
CheckGreaterOrEqual(gamma,0,'gamma');
CheckGreaterOrEqual(mu,0,'mu');
m=find(p<0 | p>1);
if length(m)>0
    error('Vaccination probability p(%g) (%g) is not between zero and one',m(1),p(m(1)));
end
CheckGreaterOrEqual(S,0,'Initial S');
CheckGreaterOrEqual(I,0,'Initial I');
CheckGreaterOrEqual(MaxTime,0,'Max Time');

if tV<0 | tV>MaxTime
    error('Vaccination start time tV (%g) is not between zero and MaxTime (%g)',tV,MaxTime);
end


% The main iteration 
options = odeset('RelTol', 1e-5);
[t, pop]=ode45(@Diff_8_4,[0 tV],[S I],options,[reshape(beta,1,4) gamma mu [0 0]]);
T=t; S=pop(:,1:2); I=pop(:,3:4);
[t, pop]=ode45(@Diff_8_4,[tV MaxTime],pop(end,:),options,[reshape(beta,1,4) gamma mu p]);
T=[T; t]; S=[S; pop(:,1:2)]; I=[I; pop(:,3:4)];

% plots the graphs with scaled colours
subplot(2,1,1)
h=plot(T/365,S(:,1),'-g',T/365,S(:,2),'--g'); 
legend(h,'High Risk','Low Risk');
xlabel 'Time (years)';
ylabel 'Susceptibles'

subplot(2,1,2) 
h=plot(T/365,I(:,1),'-r',T/365,I(:,2),'--r'); 
legend(h,'High Risk','Low Risk');
xlabel 'Time (years)';
ylabel 'Infectious'



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

beta=reshape(parameter(1:4),2,2);  gamma=parameter(5:6); mu=parameter(7:8); MU=sum(mu);
p=parameter(9:10);

S=pop(1:2); I=pop(3:4);

dPop=zeros(4,1);

for i=1:2
    dPop(i)= mu(i)*(1-p(i)) - (beta(i,1)*I(1)+beta(i,2)*I(2))*S(i) - MU*S(i);
    dPop(i+2)=(beta(i,1)*I(1)+beta(i,2)*I(2))*S(i) - gamma(i)*I(i) - MU*I(i);
end


% 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
    if length(Parameter)>1
        error('Parameter %s(%g) (=%g) is less than %g',str,m(1),Parameter(m(1)),Value);
    else
        error('Parameter %s (=%g) is less than %g',str,Parameter(1),Value);
    end
end
