# -*- coding: utf-8 -*-
"""
solution of a population growth model
plot numerical solution and analytic solution (if known)
"""

from scipy.integrate import odeint
from numpy import array,linspace
import matplotlib.pylab as plt                # for plotting commands
import numbers

# - - - global variables
aa = 0.1     # parameter of population model


# - - - population model - ODE
def deriv(y,t): # return derivative of y
    return aa*y


# - - - Main program - - -

# - - - prepare graph
plt.figure(figsize=(6,6))
plt.xlabel('Time',fontsize=20,labelpad=0)
plt.ylabel('Population',fontsize=20,labelpad=0)
#plt.title("Population model")

xmin = 0.0; xmax = 10.0
ymin = 0.0; ymax = 5.0
plt.axis([xmin,xmax,ymin,ymax])

dx = 1.0; dy = 1.0
plt.xticks(plt.arange(xmin,xmax+dx,dx),fontsize=14)
plt.yticks(plt.arange(ymin,ymax+dy,dy),fontsize=14)


# - - - solve the initial value problem numerically
# input initial condition
yinit = input('Give initial condition N0: ')
if not isinstance(yinit,numbers.Number):
    print yinit, " is not a number"
    sys.exit()

# construct time interval to solve equation
#time_interval = input('Give time interval t0,t1 (separate by comma): ')
#if time_interval[0] == time_interval[1]: sys.exit()
time_interval = (xmin,xmax)
ntime = 100
time = linspace(time_interval[0],time_interval[1],ntime)

# integrate ODE
popul = odeint(deriv,yinit,time)
print "Final point: ", popul[ntime-1]
                
# plot line
plt.plot(time,popul,'-b')
plt.show()

# save figure
#plt.savefig('plot1.png', dpi=200, transparent=True)