#My model to describe the number of Corona cases in Germany, from data in March and April 2020
#Gerald Schuller, April 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.075  #1.25 increase factor day to day
valueondate=84600  #value on a given day, April 2nd, 2020, from jhu.edu
dayofvalue=2 #the given day of the month of the value
K=83e6 #Population of Germany
timetorecovery=11.5 #16 days from dicovering to recovery

#Plot range in days:
t=np.arange(1,31) #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 April')
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 April')
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.048) #0.07 early, then 0.08 (case based mortality rate), 0.02 long term (infection based mortality)
plt.xlabel('Day after beginning of April')
plt.ylabel('Accumulated Corona Death Cases')
plt.grid()
plt.title('Corona Death Cases in Germany, from logistic model')
plt.show()