%  abalone3c
%
%   SI model of abalone in an array of subpopulations (nx,ny)
%
%    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
%
%        variables are 
%  y(1): S, susceptible, uninfected individuals [number/m^2]
%  y(2): I, infected individuals [number/m^2]
%  y(3): D, dead infected animals [number/m^2]
%  y(4): P, infectious particles close 
%            to the susceptible population [number/m^3]
%   PAR is a structure containing the parameters for the model
%
%   PARabalone3c.m sets the values of the various model parameters
%   RHSabalone3c.m defines the equations and interaction
%
global PAR
%   array indexes of variables
Nvar=5;  %  number of variables in the model
iS=1;iE=2;iI=3;iD=4;iP=5;

PAR=PARabalone3c;  %   define model parameters

y0=zeros(PAR.Ny,PAR.Nx,Nvar);  %     initial conditions
y0(:,:,iS)=[0 0 10; 0 0 50; 0 0 100; 0 0 100];
%y0(:,:,iS)=[100 100 100; 100 100 100; 100 100 100; 100 100 100];
y0(1,3,iI)=1;
y0(1,3,iP)=10;

%  put initial conditions in a column
Cy0=y0(:);

tspan=[0 300];  %   time span

[t,Cy]=ode45(@RHSabalone3c,tspan,Cy0);
%   rename model results and recover the 2d grid
Nt=length(t);
y=reshape(Cy,Nt,PAR.Ny,PAR.Nx,Nvar);
S=y(:,:,:,iS);
E=y(:,:,:,iE);
I=y(:,:,:,iI);
D=y(:,:,:,iD);
P=y(:,:,:,iP);