% BAM16
% Modified version of BAM15 -- removed 1-1/mu factors
% Also fixed the phi problem (no longer mult by tstep)

%s(t) = total population susceptible at time t
%e(t) = total population exposed at time t (incubation state)
%i(t) = total population infectious at time t
%q1(t) = total population who are non-symptomatic & in quarantine at time t
%q2(t) = total population e who are symptomatic (infectious) & in quarantine at time t
%d(t) = total population dead at time t
%m(t) = total population disabled/maimed at time t
%r(t) = total population recovered at time t

%   where
%       s = w(1)
%       e = w(2)
%       i = w(3)
%       q1 = w(4)
%       q2 = w(5)
%       d  = w(6)
%       m  = w(7)
%       r  = w(8)



function wBAM = BAM17(t,w,params)
% Right hand sides of system of O.D.E.'s as per BAM_ModelEquations_v10
%global alpha beta gamma mu1 mu2 mu3 x phi 

persistent tprev  tstep steps

if (nargin < 3)
    params = [5 3 0.3 0.3 13 20 0.05 0.005];
end     

if steps == 0
    tstep = 0;
    tprev = 0;
else   
    tstep  = t - tprev;  % compute the time period in this step
end

steps = steps + 1; % iteration number
tprev = t;         % store the current time

if or(tstep == 0 , w(1) < params(8))
   phi     = 0;
else phi   = params(8);
end 

alpha   = params(1);   
beta    = params(2); 
gamma   = params(3); 
dead    = params(4); 
mu1     = params(5); 
mu2     = params(6); 
maimed  = params(7); 

SE = beta*w(3)*w(1);
SQ1 = alpha*gamma*w(3)*w(1);
SR = phi;
EI = w(2)/mu1;
EQ1 = w(3)*w(2);
ID = dead*w(3)/mu2;
IM = maimed*w(3)/mu2;
IR = (1-dead-maimed)*w(3)/mu2;
IQ2 = gamma*w(3);
Q1R = w(4)*w(1)/(mu1*(w(1)+w(2)));
Q1Q2 = w(4)*w(2)/(mu1*(w(1)+w(2)));
Q2D = dead*w(5)/mu2;
Q2M = maimed*w(5)/mu2;
Q2R = (1-dead-maimed)*w(5)/mu2;

 s  =  -SQ1-SE-SR;
 e  =  SE-EI-EQ1;
 i  =  EI-ID-IM-IR-IQ2;
 q1 =  SQ1+EQ1-Q1R-Q1Q2;   
 q2 =  IQ2-Q2D-Q2M-Q2R+Q1Q2;
 d  =  ID+Q2D;  
 m  =  IM+Q2M;
 r  =  IR+Q1R+Q2R+SR;
 
 wBAM = [ s ; e ; i ; q1 ;q2; d ;m ;r; SE; SQ1; SR; EI; EQ1; ID; IM; IR; IQ2; Q1R; Q1Q2; Q2D; Q2M; Q2R];

