function ydot=RHSinfct3a(t,y)
%  RHSabalone3a
%
%    define the RHS of the model equations. 
%        variables are 
%  y(1): S, susceptible, uninfected individuals [number]
%  y(2): I, infected individuals [number]
%  y(3): D, dead animals from the infected population [number]
%  y(4): P, infectious particles close 
%            to the susceptible population [number]
%   array indexes of variables
iS=1;iI=2;iD=3;iP=4;Nvar=4;
global PAR 
%          most rates are 1/day. Area format
Aydot=zeros(PAR.Ny,PAR.Nx,Nvar);

%  convert column to space and variable structure
Ay=reshape(y,PAR.Ny,PAR.Nx,Nvar);

for ix=1:PAR.Nx
  for iy=1:PAR.Ny
%   time changes of susceptable animals
    Aydot(iy,ix,iS)= - PAR.IPinfect(iy,ix)*PAR.depth(iy,ix)*Ay(iy,ix,iP)*Ay(iy,ix,iS)...
                     - PAR.Iinfect(iy,ix) * Ay(iy,ix,iI) * Ay(iy,ix,iS) ...
                     - PAR.Dinfect(iy,ix) * Ay(iy,ix,iD) * Ay(iy,ix,iS) ...
                     - PAR.Bmort(iy,ix) * Ay(iy,ix,iS);

%   time changes of infected animals
    Aydot(iy,ix,iI)= PAR.IPinfect(iy,ix)*PAR.depth(iy,ix)*Ay(iy,ix,iP)*Ay(iy,ix,iS)... 
                  + PAR.Iinfect(iy,ix) * Ay(iy,ix,iI) * Ay(iy,ix,iS)... 
                  + PAR.Dinfect(iy,ix) * Ay(iy,ix,iD) * Ay(iy,ix,iS)... 
                  - PAR.Imort(iy,ix) * Ay(iy,ix,iI);

%   time changes of dead infected animals
   Aydot(iy,ix,iD)= PAR.Imort(iy,ix) * Ay(iy,ix,iI)...
                   - PAR.DeadDecay(iy,ix) * Ay(iy,ix,iD);

%   time changes of infectious particles
   Aydot(iy,ix,iP)= PAR.Irelease(iy,ix) * Ay(iy,ix,iI)...
                   + PAR.Drelease(iy,ix) * Ay(iy,ix,iD) ...
                   - PAR.IPremove(iy,ix) * Ay(iy,ix,iP);
  end
end
%             calculate the exchange of IP between grid cells
%    Uex and Vex represent exchange rates (1/day)
%    Uex(i,j) > 0 means particles move from cell(i,j) to cell(i+1,j)
%    Uex(i,j) < 0 means particles move from cell(i+1,j) to cell(i,j)
%    and similarly for Vex. 
%        this will be Uex/Vex if it is positive and zero if negative.
   PUex=.5*(abs(PAR.Uex) + PAR.Uex); 
   PVex=.5*(abs(PAR.Vex) + PAR.Vex); 
%        this will be Uex/Vex if it is positive and zero if negative.
   NUex=.5*(abs(PAR.Uex) - PAR.Uex); 
   NVex=.5*(abs(PAR.Vex) - PAR.Vex); 

%    E/W transfer
for ix=1:PAR.Nx-1
  for iy=1:PAR.Ny
   Aydot(iy,ix,iP)=  Aydot(iy,ix,iP) ...
    - PUex(iy,ix)*Ay(iy,ix,iP)  + NUex(iy,ix)*Ay(iy,ix+1,iP);
   Aydot(iy,ix+1,iP)=  Aydot(iy,ix+1,iP) ...
    + PUex(iy,ix)*Ay(iy,ix,iP)  - NUex(iy,ix)*Ay(iy,ix+1,iP);
  end
end 
%   N/S transfer
for ix=1:PAR.Nx
  for iy=1:PAR.Ny-1
   Aydot(iy,ix,iP)=  Aydot(iy,ix,iP)...
    - PVex(iy,ix)*Ay(iy,ix,iP)  + NVex(iy,ix)*Ay(iy+1,ix,iP);
   Aydot(iy+1,ix,iP)=  Aydot(iy+1,ix,iP)...
    + PVex(iy,ix)*Ay(iy,ix,iP)  - NVex(iy,ix)*Ay(iy+1,ix,iP);
  end
end 

clear PUex PVex NUex NVex
%          put in column format
ydot=Aydot(:);