%  abalone3b
%
%   multi-species SIP model
%
%    The model considers a number of benthic animals in a volume of water
%    Infectious particles can be released which spread the infection
%  
%    Units:
%      time is in days
%      populations are in number of animals/m^2
%      concentrations are in number/m^3
%
%     for each species (is) the  variables are 
%  y(is,1): S, susceptible, uninfected individuals [number/m^2]
%  y(is,2): I, infected individuals [number/m^2]
%  y(is,3): DI, dead infected animals [number/m^2]
%  y(is,4): IP, infectious particles close 
%            to the susceptible population [number/m^3]
%   PAR is a structure containing the parameters for the model
%
%   PARabalone3b.m sets the values of the various model parameters
%   RHSabalone3b.m defines the equations and interaction
%
global PAR
%   array indexes of variables
Nvar=4;  %  number of variables in the model
Nspecies=2;  %  number of species in the model

PAR=PARabalone3b(Nspecies,Nvar);  %   define model parameters
iS=PAR.iS;iI=PAR.iI;iD=PAR.iD;iP=PAR.iP;

y0=zeros(Nspecies,Nvar);  %     initial conditions
y0(:,iS)=[100; 100];
y0(1,iI)=1;
y0(2,iI)=0;
y0(1,iP)=0;
y0(2,iP)=0;

%  put initial conditions in a column
Cy0=y0(:);
%  recover the area format (test)
%Ry0=reshape(Cy0,Nspecies,Nvar);

%tspan=[0 150];  %   time span
tspan=[0 300];  %   time span

[t,Cy]=ode45(@RHSabalone3b,tspan,Cy0);
%   rename model results and recover array structure
Nt=length(t);
y=reshape(Cy,Nt,Nspecies,Nvar);
S=y(:,:,iS);
I=y(:,:,iI);
D=y(:,:,iD);
P=y(:,:,iP);