#Script to support Exercise 2.2.8, modeling yeast population growth.

#Times at which yeast population was measured
times = [i for i in range(0,18)];

#Yeast population (millions) at each time
data = [9.6, 18.3, 29, 47.2, 71.1, 119.1, 174.6, 257.3, 350.7, 441, 513.3, 559.7, 594.8, 629.4, 640.8, 651.1, 655.9,
659.6];

#Plot the raw data versus time. Call the plot "plot1".
pdata = list(zip(times,data))
plt1 = scatter_plot(pdata)
show(plt1)

#Given that u(0) = 9.6, the solution to the logistic equation with intrinsic growth rate "r" and carrying capacity "K"
is
def u(t,r,K):
    return K/(1+exp(-r*t)*(K/9.6-1))

#Take a guess r = 1 and K = 600, plot, compare to the data
var('t');
plt2 = plot(u(t,1.0,600.0),(t,0,17),color='red')
pp = plt1+plt2;
pp.axes_labels(['Time (hours)','Population'])
show(pp)

#Perhaps we can do better...