# How Epidemics Spread 4, Exercise 1: testing the coverage formula by simulation
# Exercise 1: test the formula Vc = (1 - 1/R0) / E by simulation (all values assumed).
r0, E, gamma, population, dt = 4.0, 0.9, 0.2, 10_000, 0.05
beta = r0 * gamma
critical = (1 - 1 / r0) / E
print(f"formula: coverage needed = {critical:.1%}")

def final_share(coverage):
    protected = coverage * E * population        # all-or-nothing vaccine
    S, I = population - protected - 1.0, 1.0
    peak = I
    while I >= 0.01:
        new_inf = beta * S * I / population * dt
        S -= new_inf
        I += new_inf - gamma * I * dt
        peak = max(peak, I)
    return peak, (population - protected - S) / population

print("coverage   peak infectious   share of population infected")
for coverage in [critical - 0.10, critical - 0.05, critical, critical + 0.05]:
    peak, share = final_share(coverage)
    print(f"{coverage:8.1%} {peak:17.0f} {share:20.1%}")

# Output when run:
# formula: coverage needed = 83.3%
# coverage   peak infectious   share of population infected
#    73.3%               132                16.3%
#    78.3%                37                 8.6%
#    83.3%                 1                 0.7%
#    88.3%                 1                 0.1%