function [t,S,P0,lambda] = Program_4_4(r,T,gamma,nu,mu,X0,Y0,MaxTime)
%
% 
%
% PROGRAM_4_4( r, T, gamma, nu, mu, X0, Y0, MaxTime)
%      This is the MATLAB version of program 4.2 from page 123 of 
% "Modeling Infectious Disease in humans and animals" 
% by Keeling & Rohani.
% 
% It is a simple model of vector borne transmission of infection. We
% implicitly assume a vector-notation, which humans being the first entry
% and mosquitoes being the second.
% In the default parameterisation of the model mosquitoes never recover 
% from infection, although this is allowed by the general formulation.
%

% Sets up default parameters if necessary.
if nargin == 0
   r=0.5/1e3;
   T=[0 0.5; 0.8 0];
   gamma=[0.033 0];
   mu=[5.5e-5 0.143];
   nu=[5.5e-2 1.443e3];
   X0=[1e3 1e4];
   Y0=[1 1];
   MaxTime=1000;
end

% Check sizes for terms that are 'vectors'
X0=CheckSize(X0,2,1,'Susceptibles');
Y0=CheckSize(Y0,2,1,'Infecteds');
mu=CheckSize(mu,2,1,'Per capita death rate');
nu=CheckSize(nu,2,1,'Total birth rate');

A=size(T);
if A(1)==1 & A(2)==1
    T=[0 T; T 0];
end
if ((A(1)==2 & A(2)==1) | (A(2)==1 & A(1)==2))
    T=[0 T(1); T(2) 0];
end
T=CheckSize(T,2,2,'Transmission probability');


% Checks all the parameters are valid
CheckGreaterOrEqual(r,0,'r');
CheckGreaterOrEqual(gamma,0,'gamma');
CheckGreaterOrEqual(nu,0,'nu');
CheckGreaterOrEqual(mu,0,'m');
CheckGreaterOrEqual(X0,0,'X0');
CheckGreaterOrEqual(Y0,0,'Y0');
CheckGreaterOrEqual(MaxTime,0,'MaxTime');
if T(1,1)~=0 | T(2,2)~=0
    warning('Transmission probability within species (the diagonal of T) should be zero');
end
if T(1,2)<0 | T(1,2)>1
    error('Transmission probability to humans from mosquitoes (%g) is not between zero and one',T(1,2));
end
if T(2,1)<0 | T(2,1)>1
    error('Transmission probability to mosquitoes from humans (%g) is not between zero and one',T(2,1));
end


% initial conditions are set to suit default parameters

% The main iteration 
options = odeset('RelTol', 1e-5);
[t, pop]=ode45(@Diff_4_4,[0 MaxTime],[X0 Y0],options,r,T,gamma,nu,mu);

X=pop(:,1:2); Y=pop(:,3:4); 

% plots the graphs with scaled colours
subplot(2,2,1)
h=semilogy(t,X(:,1),'-g'); 
xlabel 'Time';
ylabel 'Susceptible Humans'

subplot(2,2,2) 
h=semilogy(t,X(:,2),'-g');
xlabel 'Time';
ylabel 'Suceptible Mosquitoes'

subplot(2,2,3)
h=semilogy(t,Y(:,1),'-r'); 
xlabel 'Time';
ylabel 'Infected Humans'

subplot(2,2,4) 
h=semilogy(t,Y(:,2),'-r');
xlabel 'Time';
ylabel 'Infected Mosquitoes'


% Calculates the differential rates used in the integration.
function dPop=Diff_4_4(t, pop, r, T,gamma,nu,mu)

X=pop(1:2); Y=pop(3:4); 

dPop=zeros(4,1);

% note formulated to allow human-to-human and mosquito-to-mosquito transmission
% we also potentially let infected mosquitoes recover.
dPop(1)=nu(1) - r*(T(1,2)*Y(2) + T(1,1)*Y(1))*X(1) - mu(1)*X(1);  %X_H
dPop(2)=nu(2) - r*(T(2,2)*Y(2) + T(2,1)*Y(1))*X(2) - mu(2)*X(2);  %X_M

dPop(3)=r*(T(1,2)*Y(2) + T(1,1)*Y(1))*X(1) - gamma(1)*Y(1) - mu(1)*Y(1);  %Y_H
dPop(4)=r*(T(2,2)*Y(2) + T(2,1)*Y(1))*X(2) - gamma(2)*Y(2) - mu(2)*Y(2);  %Y_M



% 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




