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

a=.2
b=1.
c=.2
omega=1.5
T=5000.
NPT=50000
dt=T/NPT
x=np.linspace(1,2,NPT)
y=np.linspace(1,2,NPT)
r=np.linspace(1,2,NPT)
phi=np.linspace(1,2,NPT)

x[0]=.0
y[0]=.001
r[0]=math.sqrt(x[0]**2+y[0]**2)
phi[0]=math.pi/2



for i in range(1,NPT):
	r[i]=r[i-1]+dt*(a*r[i-1]-b*r[i-1]**3)
	phi[i]=phi[i-1]+dt*(omega+c*r[i]**2)
	x[i]=r[i]*math.cos(phi[i])
	y[i]=-r[i]*math.sin(phi[i])
	


fig1=plt.figure()            
plt.plot(x,y)
plt.title('Solution de l equation resonante d ordre 3 avec cycle limite')
plt.xlabel ('x')
plt.ylabel('dx/dt')

plt.show(fig1)