%Small Pox modeling

%s = susceptable 

%I = infectious

%P = prodromal

%MC = most contagious

%OS = overtly symptomatic

%R = recovered

%D = dead

%t = time step [days] 

clear all clf

fprintf('\nSmallpox Simulator\n\n')

%------------Error Variables-----------------------

USpop = 3e8; 

%---------------DAY #1-----------------------------

t = 1; %day 1

DaysOfOutbreak = 500; %length of time the outbreak is to be simulated

%Dependent variables 

initProdromal = 0;

initMostContagious = 0;

initOvertlySymptomatic = 0;

initRecovered = 0;

initDead = 0;

%---------initial conditions hard coded-------------------

initPop = 3e6; %assumed

initInfected = 10; %assumed 

TransRate = .75; %assumed 

initialVaccinated = 30000; %assumed

interrupt = 60+1; %the time where new vaccination program begins, specific day 

vacc = 111100; %number vaccinated in the new program 

% %---------initial conditions User chosen-------------------------------------------------

% %DaysOfOutbreak = input('How many days would you like to simulate the outbreak?');

% initPop = input('What is your initial population? '); 

% while initPop > USpop

%     fprintf(1,'\nYour initial population cannot exceed the population\n of the United States, which is

%     about %9.0f!\n',USpop)

%     initPop = input('What is your new initial population? '); 

% end

% 

% initInfected = input('How many people are initially infected? '); 

% while initInfected >= initPop

%     fprintf(1,'\nYour initial infected would result in collapse of social order!\n')

%     initInfected = input('What is your new amount of initial infected? '); 

% end

% 

% TransRate = input('What is the rate of transmission? '); 

% while TransRate > 10

%     fprintf(1,'\nYour initial transmission rate may not exceed 10 people\nwithin per day!')

%     TransRate = input('Please enter a transmission rate less than or equal to 10:  \n'); 

% end

% 

% initialVaccinated = input('What is the initial vaccinated population? '); 

% while initialVaccinated >= initPop/2

%     fprintf(1,'\nYour initial vaccinated may not exceed half of your initial population.\n',1)

%     initialVaccinated = input('What is your new amount of initially vaccinated? '); 

% end

% 

% interrupt = input('At what day of the outbreak would you like\n to begin supplemental vaccinations: Day =

%')+1;

% 

% vacc = input('How many people in your population\n do you wish to vaccinate supplementally? ');

% while vacc >= (initPop - initInfected -initialVaccinated)

%     fprintf(1,'\nYou may not vaccinate more people than are still able to be vaccinated \nwhich is about

%9.0f',...

%         (initPop - initInfected -initialVaccinated))

%     vacc = input('Please enter a new amount of people \nto receive supplemental vaccinations '); 

% end

s(t) = initPop - initInfected - initialVaccinated -

(TransRate*initMostContagious)/initPop; 

I(t) = initInfected*(12/13) +

TransRate*initMostContagious;  

P(t) = initProdromal + (1/13)*initInfected;

MC(t) = initMostContagious;

OS(t) = initOvertlySymptomatic;

R(t) = initRecovered;

D(t) = initDead;

%--------------Day #2 and on...

-------------------------

for t = 2:DaysOfOutbreak

    if t == interrupt

        s(t) = s(t-1) -

(TransRate*MC(t-1))*(s(t-1)/initPop) - vacc; 

        I(t) = (12/13)*I(t-1) +

(TransRate*MC(t-1))*(s(t-1)/initPop);

        P(t) = (2/3)*P(t-1) + (1/13)*I(t-1);

        MC(t) = (3/4)*MC(t-1) + (1/3)*P(t-1);

        OS(t) = (15/16)*OS(t-1) + (1/4)*MC(t-1);

        R(t) = R(t-1) + .7*(1/16)*OS(t-1) + vacc;

        D(t) = D(t-1) + .3*(1/16)*OS(t-1);

    else

        s(t) = s(t-1) -

(TransRate*MC(t-1))*(s(t-1)/initPop); 

        I(t) = (12/13)*I(t-1) +

(TransRate*MC(t-1))*(s(t-1)/initPop);

        P(t) = (2/3)*P(t-1) + (1/13)*I(t-1);

        MC(t) = (3/4)*MC(t-1) + (1/3)*P(t-1);

        OS(t) = (15/16)*OS(t-1) + (1/4)*MC(t-1);

        R(t) = R(t-1) + .7*(1/16)*OS(t-1);

        D(t) = D(t-1) + .3*(1/16)*OS(t-1);

    end

end

for timmy = 1:length(s);

    if s(timmy) < 0

        s(timmy) = 0;

    end

end

time = [1:length(s)];

totSick = P+MC+OS;

plot(time,s,'r-',time,D,'b-',time,R,'g-',time,I,'m',time,totSick,'k',[interrupt-1

interrupt-1],[s(interrupt) 0],'k')

legend('Susceptable','Dead','Recovered','Infected','Total

Sick','Time Vaccionations Occured',0)

grid