# -*- 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
import os  

plt.close()

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

x=np.linspace(1,3,NPT)
xdot=np.linspace(1,2,NPT)
xdotm=np.linspace(1,2,NPT)

for i in range(NPT):
    s=k*(.5+2*x[i]-2-x[i]**2/2)/m
    xdot[i]=math.sqrt(s)
    xdotm[i]=-xdot[i]
    

fig1=plt.figure()
plt.plot(x,xdot,x,xdotm)# trace de H
plt.plot([-1.0, -1.0], [-2.0, 2], '-', lw=2) 
plt.plot([1.0, 1.0], [-2.0, 2], '-', lw=2) 
plt.plot([-2, 4.0], [0.0, 0.0], '-', lw=2) 
plt.plot([0.0, 0.0], [-2.0, 2], '-', lw=2) 
plt.xlabel('x')    
plt.ylabel('dx/dt')
plt.title('Courbes R- et H ')
plt.show(False)

raw_input("Press enter to continue")

x[0]=.0
x[1]=0.
xdot[0]=0.
xdot[1]=0.

for j in range(2,NPT):
	x[j]=2*x[j-1]-x[j-2]-k*(x[j-1]-2)*dt**2+dt*(1-x[j-1]**2)*(x[j-1]-x[j-2])
	xdot[j]=(x[j]-x[j-1])/dt

fig2=plt.figure()            
plt.plot(x,xdot)
plt.title('Solution de VDP decalee avec point limite')
plt.xlabel('x')
plt.ylabel('dx/dt')

plt.show(fig2)