import numpy as np                       # llama a la libreria numpy (operaciones numericas)
import matplotlib.pyplot as plt          # llama a la libreria grafica
from scipy.integrate import odeint       # llama al integrador numerico 'odeint'

# vamos a definir una funcion matricial que realiza la siguente cuenta
#
# [0       +    dy(t)/dt   ]
# [-w^2*sin(y) - g*dy(t)/dt]
#

def sistema(Y,t,omega):
    y,ydot=Y
    gamma=2.0
    return ydot,-omega**2*np.sin(y)-gamma*ydot   # amortiguamiento (por unidad de masa)

t=np.linspace(0,2*np.pi,101)                     # puntos donde se evalua la funcion

y0=[1,0]                                         # posicion inicial in 1 y velocidad nula
omega=2.0                                        # frecuencia angular del oscilador

f=odeint(sistema,y0,t,args=(omega,))             # integrador (de Euler) 
y,ydot=f[:,0],f[:,1]

plt.figure()                                     # inicio nuevo grafico
plt.plot(y,ydot)                                 # grafico de diagrama de fases
plt.show()                                       # mostrar en pantalla