#My model to describe the number of Corona cases in Germany, from data in March 2020
#Gerald Schuller, March 2020

import numpy as np
import matplotlib.pyplot as plt

def logistic(K, valueondate, factorincrease, t):
   #logistic function for exponential growth with saturation at population size K
   #starting from a value on a given date "valueondate"
   #factor for the increase from day to day is "factorincrease"
   #prediction for the number of days "t" after the given day,
   #prediction before the given date is obtained by using negative values for t
   return K/(1+ ((K-valueondate)/valueondate)*np.exp(-np.log(factorincrease)*t ))

#Model for Corona cases:
#639 cases on March 6th, 25% increase per day
#639*(1.25**3) #Days since March 6

factorincrease=1.25  #increase factor day to day
valueondate=639  #value on a given day
dayofvalue=6 #the given day of the month of the value
K=83e6 #Population of Germany
timetorecovery=16 #days from dicovering to recovery

#Plot range in days:
t=np.arange(1,100) #prediction for this range of days of the month. For the next month, use day numbers larger that 30.
plt.plot(t, logistic(K, valueondate, factorincrease, t-dayofvalue) )
plt.plot(t, logistic(K, valueondate, factorincrease, t-timetorecovery-dayofvalue) ) #estimated recovered
#plt.plot(t, valueondate*(factorincrease**(t-dayofvalue)))
plt.xlabel('Day after beginning of March')
plt.ylabel('Corona Cases')
plt.legend(('Cases', 'Recovered Cases'))
plt.grid()
plt.title('Corona Cases in Germany, from logistic model')
plt.show()

#Plot of currently sick patients
#Difference cases minus recovered
plt.plot(t, logistic(K, valueondate, factorincrease, t-dayofvalue) -logistic(K, valueondate, factorincrease, t-timetorecovery-dayofvalue) )
plt.xlabel('Day after beginning of March')
plt.ylabel('Current Corona Cases')
plt.grid()
plt.title('Currently Sick Cases in Germany, from logistic model')
plt.show()

#Death cases: 2 on march 9, 5 on March 13. 
#20 days back for cases: 18% death rate, compensation for the undiscovered cases
#2*(1.25**4) #todesfaelle, tage seit 9.3.
#valueondate=2  #value on a given day
#dayofvalue=9 #the given day of the month of the value

#plt.plot(t, valueondate*(factorincrease**(t-dayofvalue)))
plt.plot(t,logistic(K, valueondate, factorincrease, t-timetorecovery-dayofvalue)*0.02) #0.07 early, 0.02 long term
plt.xlabel('Day of March')
plt.ylabel('Accumulated Corona Death Cases')
plt.grid()
plt.title('Corona Death Cases in Germany, from logistic model')
plt.show()