clear figure

N = 4*365;

S = 410000000;
I = 55595;
Q = 0;
R = 0;

coeffs = options;
coeffs.InfectionRate = 2e-2;
coeffs.ImmunisationRate = 0.029;
coeffs.ImmunisationLossRate = 2.7e-4;
coeffs.DeathRate = 1.3e-2;
coeffs.NaturalDeathRate = 2e-4;
coeffs.NatalityRate = 5.3e-4;
coeffs.QuarantineRate = 4.4e-4;
VaccinationVector = build_vacc_rate(0, 4.4e-4, 0, N);
% VaccinationVector = build_vacc_rate(0, 0, 50, N);

m = model(S, I, Q, R, 12, N);

t = 1:N;

for i=t
    m.step(i, coeffs, VaccinationVector);
end

area([m.Slist.', m.Wlist.', m.Ilist.', m.Qlist.', m.Rlist.']);

legend('Susceptible', 'Waiting', 'Infected', 'In quarantine', 'Immunized');
title('Evolution of virus in population');
xlabel('Time');
ylabel('Population');
axis([1, N, 0, inf]);

% utility functions

function vacc_rate=build_vacc_rate(start_rate, end_rate, change_step, N)
    vacc_rate = zeros(1, N);
    for i=1:change_step
        vacc_rate(i) = start_rate;
    end
    for i=change_step+1:N
        vacc_rate(i) = end_rate;
    end
end