function ydot=RHSinfct4(t,y)
%  RHSabalone4
%
%    define the RHS of the model equations. 
%      for each species, x, y, 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
global PAR 
iS=PAR.iS;iI=PAR.iI;iD=PAR.iD;iP=PAR.iP;
%          most rates are 1/day
Aydot=zeros(PAR.Nvar,PAR.Nspecies,PAR.Ny,PAR.Nx);
Ay=reshape(y,PAR.Nvar,PAR.Nspecies,PAR.Ny,PAR.Nx);

KK=ones(PAR.Nspecies,PAR.Ny,PAR.Nx);
%   calculate the carrying capacity term for each area
for is=1:PAR.Nspecies
  for iy=1:PAR.Ny
   for ix=1:PAR.Nx
    KK(is,iy,ix)=1-(Ay(iS,is,iy,ix)+Ay(iI,is,iy,ix))/PAR.Carry(is);
   end
  end
end 
for s=1:PAR.Nspecies
 for iy=1:PAR.Ny
  for ix=1:PAR.Nx
%   time changes of susceptable animals
    Aydot(iS,s,iy,ix)= - PAR.Bmort(s) * Ay(iS,s,iy,ix) ...
        + KK(s,iy,ix)*(PAR.Srepro(s)*Ay(iS,s,iy,ix) ...
                      +PAR.Irepro(s)*Ay(iI,s,iy,ix));
%   time changes of infected animals
    Aydot(iI,s,iy,ix)= - PAR.Imort(s) * Ay(iI,s,iy,ix);
%   add infection by different species
%   convention is Infect(source, infected)
    for i=1:PAR.Nspecies
     Aydot(iS,s,iy,ix)=Aydot(iS,s,iy,ix) ...
            - PAR.IPinfect(i,s)* PAR.depth(iy,ix)*Ay(iP,i,iy,ix) * Ay(iS,s,iy,ix)...
            - PAR.Iinfect(i,s) * Ay(iI,i,iy,ix) * Ay(iS,s,iy,ix) ...
            - PAR.Dinfect(i,s) * Ay(iD,i,iy,ix) * Ay(iS,s,iy,ix);
     Aydot(iI,s,iy,ix)=Aydot(iI,s,iy,ix) ...
            + PAR.IPinfect(i,s)* PAR.depth(iy,ix)*Ay(iP,i,iy,ix) * Ay(iS,s,iy,ix)...
            + PAR.Iinfect(i,s) * Ay(iI,i,iy,ix) * Ay(iS,s,iy,ix) ...
            + PAR.Dinfect(i,s) * Ay(iD,i,iy,ix) * Ay(iS,s,iy,ix);
    end

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

%   time changes of infectious particles
    Aydot(iP,s,iy,ix)= PAR.Irelease(s) * Ay(iI,s,iy,ix) ...
            + PAR.Drelease(s) * Ay(iD,s,iy,ix) ...
            - PAR.IPremove(s) * Ay(iP,s,iy,ix);
  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(iP,s,iy,ix)=  Aydot(iP,s,iy,ix) ...
    - PUex(iy,ix)*Ay(iP,s,iy,ix)  + NUex(iy,ix)*Ay(iP,s,iy,ix+1);
   Aydot(iP,s,iy,ix+1)=  Aydot(iP,s,iy,ix+1) ...
    + PUex(iy,ix)*Ay(iP,s,iy,ix)  - NUex(iy,ix)*Ay(iP,s,iy,ix+1);
  end
 end 
%   N/S transfer
 for ix=1:PAR.Nx
  for iy=1:PAR.Ny-1
   Aydot(iP,s,iy,ix)=  Aydot(iP,s,iy,ix)...
    - PVex(iy,ix)*Ay(iP,s,iy,ix)  + NVex(iy,ix)*Ay(iP,s,iy+1,ix);
   Aydot(iP,s,iy+1,ix)=  Aydot(iP,s,iy+1,ix)...
    + PVex(iy,ix)*Ay(iP,s,iy,ix)  - NVex(iy,ix)*Ay(iP,s,iy+1,ix);
  end
 end 
end

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