function ydot=allometricRHS(t,y)
%  rhsALM
%
%    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);

particles=sum(Ay(:,iP));

for s=1:PAR.Nspecies
%      time changes of susceptable animals
   fecund=PAR.fecundA*PAR.mass(s).^PAR.fecundB;
   mort1=PAR.mort1A*PAR.mass(s).^PAR.mort1B;
   mort2=PAR.mort2A*PAR.mass(s).^PAR.mort2B;
   imort=PAR.ImortA*PAR.mass(s).^PAR.ImortB;
   Ddecay=PAR.DeadDecayA*PAR.mass(s).^PAR.DeadDecayB;
   release=PAR.releaseA*PAR.mass(s).^PAR.releaseB;
   beta=PAR.betaA*PAR.mass(s).^PAR.betaB;
   mort=mort1 + mort2*(Ay(s,iS)+Ay(s,iI));

 Aydot(s,iS)= fecund*(Ay(s,iS)+Ay(s,iI)) ...
             - mort * Ay(s,iS) -beta*particles*Ay(s,iS);
 
%   time changes of infected animals
 Aydot(s,iI)= - (imort+mort) * Ay(s,iI)+beta*particles*Ay(s,iS);

%   time changes of dead infected animals
 Aydot(s,iD)= (imort + mort) * Ay(s,iI) ...
            - Ddecay * Ay(s,iD);
%   time changes of infectious particles
 Aydot(s,iP)= release* Ay(s,iI) - PAR.remove * Ay(s,iP);
end
%  put in column format
ydot=Aydot(:);