# -*- coding: utf-8 -*-
"""
Created on Wed Jan  1 14:30:33 2015
ok
@author: joseorellana
"""
# import basic modules
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
import pylab as pyl
import math
import cmath 

plt.close()

NPT=10000
T=100.
dt=T/NPT
m=1.
k=1.

x=np.linspace(1.32,2.88,NPT)-2.1
y=np.linspace(1.32,2.88,NPT)
xdot=np.linspace(1,2,NPT)
xdotm=np.linspace(1,2,NPT)
c0=.6
c1=.66

for i in range(NPT):
    s3=(x[i]-x[i]**3/3-c1)**2/m+k*x[i]**2-c0
    s2=x[i]-x[i]**3/3-c1
    s1=m
    delta=s2**2-s1*s3
    
    if(delta>0):
		xdot[i]=(-s2+math.sqrt(delta))/m
		xdotm[i]=(-s2-math.sqrt(delta))/m
    else:
        xdot[i]=-s2
        xdotm[i]=-s2
    
    y[i]=x[i]**2-x[i]**3/3-c1


fig1=plt.figure()
plt.plot(x,y)
plt.plot([-1, 2],[0,0])
plt.xlabel('x')
plt.ylabel('Fonction x Q(x) ici p=0')
plt.title(' Representation de x Q(x) ici p=0')

#plt.show(False)

#raw_input("Press enter to continue")

fig2=plt.figure()   
plt.plot(x,xdot,x,xdotm)
plt.xlabel('x')
plt.ylabel('dx/dt')
plt.title('Courbe J  definissant un cycle limite')
plt.show(fig2)