#pour la division euclidienne
from __future__ import division
import numpy as np
# Pour les figures
import matplotlib.pyplot as plt
# Pour lire des images
import os
from pylab import *
# Pour les equa diff
from  scipy  import *
from  scipy.integrate  import  odeint
# Pour toutes les fonctions mathematiques
from math import*

#--------------------------------------------------------------------
#                       Oscillateur amorti
#                       xpp+2*eps*w0*xp+(w0^2)*x=0
#--------------------------------------------------------------------
# Parametres du probleme
w0 = 1
# pulsation propre en rad/s
eps = 1.1
# eps=1/2Q
# eps = 0 oscillateur harmonique
# valeur critique : eps =1

# Conditions initiales
Ninit=8 # nbr de conditions initiales
x0=0*np.ones(Ninit) # angle initial fixe a 0
v0=np.linspace(-3,3,Ninit) #vitesse initiale en rad/s
CI=array([x0,v0])


#--------------------------------------------------------------------
#       Obtention des solutions
#--------------------------------------------------------------------

# -------------------------------------------------------------------
# Definition de l'equa. diff.

def amort(y, t):
    x, v = y                            # Vecteur variable
    dvdt = -w0*2*eps*v-w0**2*x              # Equa. diff. 2
    dydt = [v,dvdt]                     # Vecteur solution
    return dydt                         

# -------------------------------------------------------------------


# Parametres d'integration
start = 0                        # debut
end = 60                         # fin
numsteps = 500                 # nombre de pas d'integration
t = np.linspace(start,end,numsteps)


# Tableau des solutions
X=np.zeros((numsteps, Ninit))
V=np.zeros((numsteps, Ninit))


# Resolution numerique des equations differentielles
# Fait un tableau ou les lignes correspondent au temps
# La premiere colone contient la solution pour la variable y[0]=x
# La deuxieme colone contient la solution pour la variable y[1]=v
for i in range(Ninit):
    sol=odeint(amort,CI[:,i],t)            
    X[:,i]=sol[:,0]
    V[:,i]=sol[:,1]



# Portrait de phase
fig = plt.figure(figsize=(10,10))
for i in np.arange(Ninit):
    plot(X[:, i], V[:, i], '-', color='blue') #label='v_0={}'.format(v0[i])
    plt.axis((-3.2,3.2,-3.1,3.1))
title('Oscillateur amorti : Portait de phase')#, fontsize=20)
xlabel(r'$x$', fontsize=20)
ylabel(r'$\dot{x} / {\omega_0}$', fontsize=20,rotation=0)
grid(True)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)

#plt.savefig('criti_phase.png')
plt.show()



##--------------------------------------------------------------------
## Description
## routine amorti.py
## function amort(y,t,eps,w0)
##     arguments: y: solution, t:temps d'integration, Q,w0: variables du probleme.
## commentaires:
## On resoud numeriquement l'equa. diff. de l'oscillateur amorti pour plusieurs valeurs du facteur de qualite (boucle sur i avec i le nombre de Q testes differents.
## La resolution numerique des equations differentielles donne un tableau ou les lignes correspondent au temps, la premiere colone  a la solution pour la variable y[0]=x et la deuxieme colone a la solution pour la variable y[1]=v
##--------------------------------------------------------------------