%  abalone4
%
%   multi-species, multi-population 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 and area, the  variables are 
%  y(1): S, susceptible, uninfected individuals [number/m^2]
%  y(2): I, infected individuals [number/m^2]
%  y(3): DI, dead infected animals [number/m^2]
%  y(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;Nspecies=2;Nx=3;Ny=4;
PAR=PARabalone4(Nvar,Nspecies,Nx,Ny);  %   define model parameters
iS=PAR.iS;iI=PAR.iI;iD=PAR.iD;iP=PAR.iP;

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

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

tspan=[0 150];  %   time span

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