# -*- coding: utf-8 -*-

"""
Ugo Hincelin - Mai 2017
Portrait de phase d'un pendule simple avec/sans amortissement.
"""

from __future__ import division
from scipy import *
from pylab import *
from scipy.integrate import odeint      # Module de résolution des équations différentielles
import matplotlib.pyplot as plt

k = 1       # pulsation^2 (= omega_0^2)
k1 = 0.1    # amortissement (= omega_0/Q ; Q>0.5 régime pseudo-périodique ; Q=0.5 critique ; Q<0.5 apériodique)

def deriv(syst, t):
    X = syst[0]                       # Variable1 NBi(t)
    Y = syst[1]                       # Variable2 NPo(t)
    dXdt = Y                       # Equation différentielle 1
    dYdt = -k*sin(X)-k1*Y                  # Equation différentielle 2
    return [dXdt,dYdt]       # Dérivées des variables

# Paramètres d'intégration
start = 0
end = 1000
numsteps = 10000
t = linspace(start,end,numsteps)

# Conditions initiales et résolution
X0 = pi # pi
Y0 = 3 # 0
syst_CI=array([X0,Y0])         # Tableau des CI
Sols=odeint(deriv,syst_CI,t)            # Résolution numérique des équations différentielles

# Récupération des solutions
X = Sols[:, 0]
Y = Sols[:, 1]

# Graphiques des solutions
# Solutions numériques
#plot(t, X, 'o', ms=6, mfc='w', mec='b',label=u"A")
#plot(t, Y, 'o', ms=6, mfc='w', mec='r',label=u"B")
plt.plot(X, Y,ms=6, mfc='w', mec='b')
plt.plot(X[0], Y[0],'x', ms=6, mfc='w', mec='r',label=u"Depart")
plt.plot(X[numsteps-1], Y[numsteps-1],'x', ms=6, mfc='w', mec='g',label=u"Arrivee")
#plt.axis([-1,1,-1,1])

plt.xlabel("$X \, (unité \ arbitraire)$", fontsize=16)     # Label de l'axe des abscisses
plt.ylabel("$Y \, (unité \ arbitraire)$", fontsize=16) # Label de l'axe des ordonnées

legend()                                                # Appel de la légende
show()