function [t,IH,IL] = Program_3_1(beta,gamma,nH,IH,IL, MaxTime)
%
% 
%
% Program_3_1( beta, gamma, nH, IH, IL, MaxTime)
%      This is the MATLAB version of program 3.1 from page 58 of 
% "Modeling Infectious Disease in humans and animals" 
% by Keeling & Rohani.
% 
% It is the SIS model with two different risk-groups.
% Note we no-longer explicitly model the susceptible class.


% Sets up default parameters if necessary.
if nargin == 0
   beta=[10 0.1; 0.1 1];
   gamma=1;
   nH=0.2;
   IH=1e-5;
   IL=1e-3;
   MaxTime=15;
end

% Checks all the parameters are valid
if IH<0 
    error('Initial level of infection in high risk group (%g) is less than zero',IH);
end

if IH>nH 
    error('Initial level of infection in high risk group (%g) is great than the high risk population size (%g)',IH,nH);
end

if IL<0 
    error('Initial level of infection in low risk group (%g) is less than zero',IH);
end

if IL>(1-nH)
    error('Initial level of infection in low risk group (%g) is great than the low risk population size (%g)',IH,nH);
end

if IH+IL<=0 
    error('Initial level of total infection (%g) is less than or equal to zero',IH+IL);
end

if nH<=0 
    error('High risk population (%g) is less than or equal to zero',nH);
end

if nH>=1 
    error('Low risk population (%g) is less than or equal to zero',1-nH);
end

if length(find(beta<0))>0 
    error('Transmission rate beta ([%g %g; %g %g]) is less than zero',beta(1,1),beta(1,2),beta(2,1),beta(2,2));
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
    
if beta(1,2)~=beta(2,1)
    warning('Transmission rate beta ([%g %g; %g %g]) is not symmetric',beta(1,1),beta(1,2),beta(2,1),beta(2,2));

end

% The main iteration 
options = odeset('RelTol', 1e-5);
[t, pop]=ode45(@Diff_3_1,[0 MaxTime],[IH IL],options,[reshape(beta,1,4) gamma nH]);

IH=pop(:,1); IL=pop(:,2);

% plots the graphs with scaled colours
subplot(2,1,1)
h=plot(t,IH,'-r',t,IL,'--r'); 
legend(h,'High Risk','Low Risk');
xlabel 'Time';
ylabel 'Infectious'

subplot(2,1,2) 
h=semilogy(t,IH,'-r',t,IL,'--r'); 
legend(h,'High Risk','Low Risk');
xlabel 'Time';
ylabel 'Infectious'


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

beta=reshape(parameter(1:4),2,2);  gamma=parameter(5); nH=parameter(6);
IH=pop(1); IL=pop(2); SH=nH-IH; SL=(1-nH)-IL;

dPop=zeros(2,1);

dPop(1)=(beta(1,1)*IH+beta(1,2)*IL)*SH - gamma*IH;
dPop(2)=(beta(2,1)*IH+beta(2,2)*IL)*SL - gamma*IL;



