% File runBAM17.m
% Modified version of BAM15 --- removed the 1-1/mu factors 
%   (as recommended by Dr. Liu)
%    Also fixed the phi problem (no longer mult by tstep)

% Runs the BAM SIRQ disease & quarantine model
%

% Definitions used by this implementation
% 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
% Right hand sides of system of  O.D.E.'s
%  dw(1)/dt=...
%  dw(2)/dt=...
%  dw(3)/dt=...   
%  .... 
%
%   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 RunBAM17(filenam,plist)

totalpopulation = 1000000

initExposed = 300

Sinit   = (totalpopulation-initExposed)/totalpopulation  % initial number of susceptible a very large city 
Einit   = initExposed/totalpopulation   % initial number of exposed due to the attack
Iinit   = 0.0000   % initial number of infected at the time of the attack
Q1init  = 0       % initial number of people in Q1
Q2init  = 0       % initial number of people in Q2
xinit   = 0       % initial number of Dead because of attack
minit   = 0       % initial number of Maimed because of attack
Rinit   = 0       % initial number of Recovered at the time of attack

SEinit = 0.0;
SQ1init = 0.0;
SRinit = 0.0;
EIinit = 0.0;
EQ1init = 0.0;
IDinit = 0.0;
IMinit = 0.0;
IRinit = 0.0;
IQ2init = 0.0;
Q1Rinit = 0.0;
Q1Q2init = 0.0;
Q2Dinit = 0.0;
Q2Minit = 0.0;
Q2Rinit = 0.0;

tmax    = 200 % amount of time to run the model

% Setup Intial values for all states

%w0 = [Sinit; Einit; Iinit; Q1init; Q2init; xinit; minit ;Rinit];

w0 = [Sinit;Einit;Iinit;Q1init;Q2init;xinit;minit;Rinit;SEinit;SQ1init;SRinit;EIinit;EQ1init;IDinit;IMinit;IRinit;IQ2init;Q1Rinit;Q1Q2init;Q2Dinit;Q2Minit;Q2Rinit];

%Setup the time range for the simulation
tlim = [0 tmax];

options = odeset('RelTol',.0000000000001,'AbsTol',.00000000000001,'Stats','on');

% run the model
[t, w, stats ] = ode45(@BAM17, tlim, w0,[], plist);
 


wcopy = w;
tempw = w;

tempw (:,[1 2 3 4 5 6 7 8]) =[];    % Discard columnes 1-8

temptransit = tempw;         

temptransit = [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0; temptransit];  % Add a row to the beginining

temptransit(end,:) = [];                                    % Delete the last row   

transit = tempw - temptransit;

w(:,[9 10 11 12 13 14 15 16 17 18 19 20 21 22]) = [];         % Delete columns from w

% write output to excel

par = {'alpha','beta','gamma','dead','mu1','mu2','maimed','phi' };

filenam = strcat(filenam, '.xls');

delete(filenam);

hdr  = {'t' , 'S' , 'E' , 'I' , 'Q1' ,'Q2', 'd','m', 'r','SE', 'SQ1', 'SR', 'EI', 'EQ1', 'ID', 'IM', 'IR', 'IQ2', 'Q1R', 'Q1Q2', 'Q2D', 'Q2M', 'Q2R'};

%wtemp = [t wcopy];
%xlswrite(filenam, hdr,   'Debug','A3'); 
%xlswrite(filenam, wtemp, 'Debug', 'A4'); 
%wtemp = []; % Deallocate

wtemp = [t w  transit];
xlswrite(filenam,par,'Percentage','A1');
xlswrite(filenam,plist,'Percentage','A2');
xlswrite(filenam, hdr,   'Percentage','A3');           % Creates a sheet called Output and writes a header row for all data items
xlswrite(filenam, wtemp, 'Percentage','A4');           % Writes the output of the wtemp vector to sheet outputs starting at cell A2
wtemp   = []; % Deallocate


w       = w*totalpopulation;
transit = transit* totalpopulation;
wtemp = [t w transit];

xlswrite(filenam,par,'Counts','A1');
xlswrite(filenam,plist,'Counts','A2');
xlswrite(filenam, hdr,   'Counts','A3');           % Creates a sheet called Counts and writes a header row for all data items
xlswrite(filenam, wtemp, 'Counts','A4');
wtemp = []; % Deallocate


 %make a plot
plot(t, w, 'LineWidth', 2)
title(filenam)
xlabel('time')
ylabel('Individuals in Each State')
legend('S', 'E', 'I','Q1', 'Q2', 'D','M','R')

w =[];