
# python code

"""
find solution for a dynamical system (2x2)
and render phase portrait
Plot lines in real-time: ask for initial condition and plot - repeatedly
"""

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


# system of two ODEs
def deriv(y,t): # return derivatives of the array y
    """nonlinear ODEs"""
    return array([ 2+y[0]*(1.0-y[1]), y[0]-y[1] ])


# Main program

# prepare graph
plt.figure(figsize=(6,6))
plt.xlabel(r'$x$',fontsize=20,labelpad=0)
plt.ylabel(r'$y$',fontsize=20,labelpad=0)
xmin = -4.0; xmax = 4.0
ymin = -4.0; ymax = 4.0
dx = 2.0; dy = 2.0
plt.xticks(plt.arange(xmin,xmax+dx,dx),fontsize=14)
plt.yticks(plt.arange(ymin,ymax+dy,dy),fontsize=14)

plt.axis([xmin,xmax,ymin,ymax])
plt.ion()
plt.show()

while True:
    
# input initial condition
    init = input('Give initial condition x,y (0,0 to end): ')
    if init[0] == 0.0 and init[1] == 0.0: break
    yinit = array(init)             # initial values
    #time_interval = input('Give time interval t0,t1 (same to end): ')
    #if time_interval[0] == time_interval[1]: break
    time_interval = (0.0,5.0)
    ntime = 100    # number of points on time interval
    time = linspace(time_interval[0],time_interval[1],ntime)

# integrate ODEs
    y = odeint(deriv,yinit,time)
    print "Final point: ", y[ntime-1,0],",",y[ntime-1,1]

# discard points ourside plotting interval
    yp1 = []
    yp2 = []
    for iline in range(ntime):
        if abs(y[iline,0]) > 1.e4 or abs(y[iline,1]) > 1.e4:
            break
        else:
            yp1.append(y[iline,0])
            yp2.append(y[iline,1])
                
# plot line
    plt.plot(yp1,yp2,'-b')
    plt.pause(0.01)             


#savefig('plot1.png', dpi=96)