% version using data up to n March (where n is in the file name), i.e. n+1 observations

close all
clear
clc

% data  
Nsick = [66, 127, 152, 196, 262, 400, 639, 795, 902, 1139, 1296 ]; % from Robert Koch Institut
NObs = numel(Nsick);
day = linspace(0,NObs-1,NObs);
% date = ["29.2", "1.3","2.3.", ...]

% growth rate number of sick individuals
Nsick_advanced = Nsick(2:NObs); % first observation left out
newly_sick = Nsick_advanced-Nsick(1:NObs-1); 
gSick = newly_sick./Nsick(1:NObs-1)*100; % in percent
gSick_av = (gSick(end)+gSick(end-1)+gSick(end-2)+gSick(end-3)+gSick(end-4)+gSick(end-5)+gSick(end-6))/7;
Nsick_7_days=Nsick(end)*(1+gSick_av/100)^7;
Nsick_30_days=Nsick(end)*(1+gSick_av/100)^30;
Nsick_40_days=Nsick(end)*(1+gSick_av/100)^40;


% individual infection rate
p = (Nsick_advanced-Nsick(1:NObs-1))./(83100000-Nsick(1:NObs-1))*1000; % probability per 1000 "Promille"
groupSize = 1000./p; % chance to be infected is 1 in groupsize

% prob individual healthy after n days with most recent probability
prob = p(end)/1000; 
healthy_7_days = (1-prob)^7 * 100;
healthy_30_days = (1-prob)^30 * 100;

% average over last 7 days
probAverage7days = (p(end)+p(end-1)+p(end-2)+p(end-3)+p(end-4)+p(end-5)+p(end-6))/1000/7; 
healthy_av_7_days = (1-probAverage7days)^7 * 100;
healthy_av_30_days = (1-probAverage7days)^30 * 100;

% overall average growth rate
gSick_av_all = mean(gSick);
% this computes the number n of days until everybody is sick in Germany: 83*10^6= Nsick(end) * (1+gSick_av/100)^n; 
n_all_sick = log(83*10^6/Nsick(end))/log(1+gSick_av_all/100);

% display results
disp(' ')
disp('Background - newly sick since 1 March 2020')
disp(newly_sick)

disp(' ')
disp('most current number - number of the day')
fprintf('indi rate to get sick is %g\n',prob)
fprintf('This is a probability of %g %%\n',prob*100)
fprintf('In other words 1 out of %g people got sick\n',1/prob)

disp(' ')
fprintf('With the most recent sickness rate of %g %%, the probability\n',prob*100)
fprintf('to be healthy is %g%% for 1 week and %g%% for 1 month (30 days)\n',healthy_7_days, healthy_30_days)

disp(' ')
disp('average last 7 days')
fprintf('indi prob over last seven days to get sick is %g\n',probAverage7days)
fprintf('In other words 1 out of %g people got sick\n',1/probAverage7days)
fprintf('To stay healthy with this average rate is %g%% for 1 week and %g%% for 1 month (30 days)\n',healthy_av_7_days, healthy_av_30_days)

disp(' ')
disp('growth rates and projection')
fprintf('most recent growth rate number of sick is %g%%\n', gSick(end));
fprintf('average growth rate number of sick over last 7 days is %g%%\n', gSick_av);
fprintf('number of sick after 7 days %g\n', Nsick_7_days);
fprintf('number of sick after 30 days %g\n', Nsick_30_days);
fprintf('number of sick after 40 days %g\n', Nsick_40_days);


disp(' ')
fprintf('average growth rate since beginning is %g%%\n', gSick_av_all);
fprintf('with this average growth rate, all of Germany would be sick after %g days\n', n_all_sick);

% prob groups healthy after n days
healthy_7_100_indis = (healthy_7_days/100)^100;
healthy_7_1000_indis = (healthy_7_days/100)^1000;

healthy_30_100_indis = (healthy_30_days/100)^100;
healthy_30_1000_indis = (healthy_30_days/100)^1000;

healthy_7_1000000_indis = (healthy_7_days/100)^1000000;
healthy_30_1000000_indis = (healthy_30_days/100)^1000000;

disp(' ')
fprintf('The probability that a group of 100 (1000) individuals is healthy after 7 days\n')
fprintf('is given by %g %% (%g%%)\n', healthy_7_100_indis*100, healthy_7_1000_indis*100)

fprintf('The probability that a group of 100 (1000) individuals is healthy after 30 days\n')
fprintf('is given by %g %% (%g%%)\n', healthy_30_100_indis*100, healthy_30_1000_indis*100)

fprintf('With a million indis, prob after 7 days is %g %% \n', healthy_7_1000000_indis*100)
fprintf('With a million indis, prob after 30 days is %g %% \n', healthy_30_1000000_indis*100)


% figures
figure()
subplot(1,2,1)
plot(day,Nsick,'o')
ylabel('Fallzahlen BRD (RKI)','FontSize',20)
xlim([0,NObs])
xlabel('Tag (März 2020)','FontSize',20)
set(gca,'FontSize',14) % Achsenbeschriftung und Legende

subplot(1,2,2)
plot(day(2:end),gSick,'x')
ylabel('Zuwachs pro Tag (in %)','FontSize',20)
xlabel('Tag (März 2020)','FontSize',20)
xlim([0,NObs])
set(gca,'FontSize',14) % Achsenbeschriftung und Legende

figure()
plot(day(2:end), newly_sick,'-o')
ylabel('newly sick')

figure()
subplot(1,2,1)
plot(day(1:end-1),p,'-o')
ylabel('p (in Promille)','FontSize',20)
xlim([0,NObs-1])
xlabel('Tag (März 2020)','FontSize',20)
set(gca,'FontSize',14) % Achsenbeschriftung und Legende

subplot(1,2,2)
plot(day(1:end-1),groupSize,'-x')
ylabel('Gruppengröße','FontSize',20)
xlim([0,NObs-1])
xlabel('Tag (März 2020)','FontSize',20)
set(gca,'FontSize',14) % Achsenbeschriftung und Legende