# -*- coding: utf-8 -*-
"""
Created on Wed Jan 21 13:53:19 2015

listing 11.5SP : mouvement d'un pendule à l'aide d'un
sous programme de scipy
u[0]: position; u[1]: vitesse
@author: grivet
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import *

def pdl(u,t):
    return [u[1],-np.sin(u[0])- r*u[1] + g*np.cos(omg*t)]
#
#attention à l'ordre des arguments !!!
#
tmax = float(input("duree de la simulation: "))
u00 = float(input("position initiale: "))
u10 = float(input("vitesse initiale: "))
u0 = np.array([u00,u10])
r = float(input("coefficient d amortissement: "))
omg = float(input("pulsation du couple exterieur: "))
g = float(input("amplitude du couple exterieur: "))


t0 = 0; t = np.arange(t0,tmax,0.05)

u = odeint(pdl, u0,t);
G = g*np.cos(omg*t);

plt.plot(t,u[:,0],t,u[:,1], t, G)
plt.title("mouvement d'un pendule: position, vitesse et couple excitateur")
plt.xlabel("temps")
Ep = 1 - np.cos(u[:,0])
Ec = 0.5*u[:,1]*u[:,1]
Et = Ec + Ep
plt.figure(2)
plt.plot(t,Ec,label = "E_cin")
plt.plot(t, Ep, label="E_pot")
plt.plot(t,Et, label = "E_tot")
plt.title('énergies pour le pendule')
plt.xlabel("temps")
plt.legend()
plt.figure(3)
plt.title("plan de phase")
plt.plot(u[:,0],u[:,1])
plt.plot(u0[0],u0[1],'ok')
plt.plot(u[-1,0],u[-1,1],'*k')
plt.xlabel("position"); plt.ylabel("vitesse")