function ydot=RHSinfct3b(t,y)
%  RHSabalone3b
%
%    define the RHS of the model equations. 
%      for each species, variables are 
%  y(is,1): S, susceptible, uninfected individuals [number]
%  y(is,2): I, infected individuals [number]
%  y(is,3): D, dead animals from the infected population [number]
%  y(is,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.Nspecies,PAR.Nvar);
Ay=reshape(y,PAR.Nspecies,PAR.Nvar);

%   calculate the carrying capacity term
KK=1-(Ay(:,iS)+Ay(:,iI))./PAR.Carry(:);
i=find(KK < 0); KK(i)=0;

for s=1:PAR.Nspecies
%   time changes of susceptable animals
 Aydot(s,iS)= - PAR.Bmort(s) * Ay(s,iS) ...
       + KK(s)*(PAR.Srepro(s)*Ay(s,iS) ...
               +PAR.Irepro(s)*Ay(s,iI));
%   time changes of infected animals
 Aydot(s,iI)= - PAR.Imort(s) * Ay(s,iI);
%   add infection by different species
%   convention is Infect(source, infected)
 for i=1:PAR.Nspecies
   Aydot(s,iS)=Aydot(s,iS) ...
            - PAR.IPinfect(i,s)* Ay(i,iP) * Ay(s,iS)...
            - PAR.Iinfect(i,s) * Ay(i,iI) * Ay(s,iS) ...
            - PAR.Dinfect(i,s) * Ay(i,iD) * Ay(s,iS);
   Aydot(s,iI)=Aydot(s,iI) ...
            + PAR.IPinfect(i,s)* Ay(i,iP) * Ay(s,iS)...
            + PAR.Iinfect(i,s) * Ay(i,iI) * Ay(s,iS) ...
            + PAR.Dinfect(i,s) * Ay(i,iD) * Ay(s,iS);
 end

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

%   time changes of infectious particles
 Aydot(s,iP)= PAR.Irelease(s) * Ay(s,iI) ...
            + PAR.Drelease(s) * Ay(s,iD) ...
            - PAR.IPremove(s) * Ay(s,iP);
end
%  put in column format
ydot=Aydot(:);