#!/usr/bin/python
# # -*- coding: utf_8 -*-

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import os
from pylab import *
from  scipy  import *
from  scipy.integrate  import  odeint
from math import*

#--------------------------------------------------------------------
#                       Oscillateur Lotka-Volterra
#                         xp = x*(1-y)
#                         yp = -y*(1-x)
#--------------------------------------------------------------------

#Ce programme permet de tracer les trajectoires de phase de l'oscillateur
#harmonique pour plusieurs valeurs de vitesse initiale (l'angle initial est nul)
#et de visualiser l'enrichissment spectral causé par les non-linéarités en
#traçant l'évolution temporelle d'une part et la TF des oscillations à 
#grande amplitude d'aure part.

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

# Conditions initiales
Ninit=5 # nbr de conditions initiales
x0=np.ones(Ninit) # angle initial fixé à 0
v0=np.linspace(0.1,0.9,Ninit) #vitesse initiale en rad/s
CI=array([x0,v0])

# Paramètres d'intégration
start = 0                        # debut
end = 15                         # fin
numsteps = 750                  # nombre de pas d'integration
t = np.linspace(start,end,numsteps)



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

# Definition des deux équa diff 
def lotka(y, t):
    x, v = y                    # Vecteur variable
    dydt = [x*(1-v),-v*(1-x)]
    return dydt                 # Solutions
                   
# Résolution numérique des équations différentielles
# Fait un tableau où les lignes correspondent au temps
# La premiere colone contient la solution pour la variable y[0]=x
# La deuxième colone contient la solution pour la variable y[1]=v
for i in range(Ninit):
    sol=odeint(lotka,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((0,4,0,4))
title('Systeme proie-predateur : Portait de phase')#, fontsize=20)
xlabel('x', fontsize=20)
ylabel('y', fontsize=20,rotation=90)
grid(True)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.savefig('lotkavolterra_phase.png')
plt.show()


# Evolution temporelle
fig = plt.figure(figsize=(16,6))
title('Evolution temporelle d un systeme proie-predateur',fontsize=22)
plt.plot(t, X[:,2], color='blue', label='proie') #on choisit la solution d'amplitude moyenne
plt.plot(t, V[:,2], color='green', label='predateur')
plt.axis((0,15,0,2))
plt.legend()
xlabel('temps', fontsize=20)
ylabel('densites de populations', fontsize=20, rotation=90)
grid(True)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.savefig('lotka_volterra_temp.png')
plt.show()
#Spectre
#plt.subplot(1,2,2)
#freq=np.fft.fftfreq(numsteps,(end-start)/numsteps) #construction du vecteur des fréquences, contient autant de points que de pas de temps
#plt.plot(freq,abs(np.fft.fft(X[:,14],numsteps))) #tracé de la partie réelle de la FFT
#plt.axis((0,4,0,5000))
#plt.xticks(fontsize=20)
#plt.yticks(fontsize=20)
#plt.grid()
#xlabel('frequence (Hz)', fontsize=20)
#plt.ylabel('amplitude',fontsize=20)
#plt.title('Spectre de Fourier',fontsize=20)
#plt.legend(fontsize=10)
#plt.savefig('harm_temp.png')
#plt.show()