#!/home/cbp/Documents/Python-3.4.0/python
# -*- coding: utf-8 -*-
"""
Créé le 19 juillet 2014

@author: Matthieu
"""
#from __future__ import print_function
#from __future__ import division

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


# ceci est l'astuce pour transformer une ED2 en 2 ED1
def pendule(var,t):
    return [acc(var[0],var[1],t),var[0]]



#################################
### pendule non amorti, régime libre
#################################

if True:
    # là est l'équation différentielle sous la accélération = fonction(v,x,t)
    def acc(v,x,t):
        return(-np.sin(x))
    
    # plage temporelle
    t=np.linspace(0,20,200)
    
    # les tracés des états liés
    for v0 in np.linspace(.4,1.99,5):
        simulation = odeint(pendule,[v0,0],t)
        plt.plot(simulation[:,1],simulation[:,0],color='b',lw='2')
        simulation = odeint(pendule,[v0,-2*np.pi],t)
        plt.plot(simulation[:,1],simulation[:,0],color='r',lw='2')
        simulation = odeint(pendule,[v0,2*np.pi],t)
        plt.plot(simulation[:,1],simulation[:,0],color='g',lw='2')
    
    # les tracés des états de diffusion
    for v0 in np.linspace(2.4,4,5):
        simulation = odeint(pendule,[v0,-2*np.pi],t)
        plt.plot(simulation[0:900,1],simulation[0:900,0],color='m',lw='2')
        simulation = odeint(pendule,[-v0,2*np.pi],t)
        plt.plot(simulation[0:900,1],simulation[0:900,0],color='c',lw='2')
    
        # esthétisme de l'affichage
    plt.title(u'portrait de phase')
    plt.axis([-2*np.pi,2*np.pi,-4.2,4.2])
    #plt.adjust(left=0.05,top=0.95,bottom=0.1,right=0.95)
    plt.show()



#################################
### pendule faiblement amorti, régime libre
#################################

if True:
    # là est l'équation différentielle sous la accélération = fonction(v,x,t)
    def acc(v,x,t):
        return(-np.sin(x)-h*v)
    
    #coeff amortissement
    h =0.1
    
    # plage temporelle
    t=np.linspace(0,20,200)
    
    # les tracés des états liés
    for v0 in np.linspace(.4,1.99,3):
        simulation = odeint(pendule,[v0,0],t)
        plt.plot(simulation[:,1],simulation[:,0],color='b',lw='2')
        simulation = odeint(pendule,[v0,-2*np.pi],t)
        plt.plot(simulation[:,1],simulation[:,0],color='r',lw='2')
        simulation = odeint(pendule,[v0,2*np.pi],t)
        plt.plot(simulation[:,1],simulation[:,0],color='g',lw='2')
    
    # les tracés des états de diffusion
    for v0 in np.linspace(2.4,4,3):
        simulation = odeint(pendule,[v0,-2*np.pi],t)
        plt.plot(simulation[0:900,1],simulation[0:900,0],color='m',lw='2')
        simulation = odeint(pendule,[-v0,2*np.pi],t)
        plt.plot(simulation[0:900,1],simulation[0:900,0],color='c',lw='2')
    
        # esthétisme de l'affichage
    plt.title(u'portrait de phase')
    plt.axis([-2*np.pi,2*np.pi,-4.2,4.2])
    #plt.adjust(left=0.05,top=0.95,bottom=0.1,right=0.95)
    plt.show()



#################################
### pendule amorti, régime libre
#################################

if True:
    # là est l'équation différentielle sous la accélération = fonction(v,x,t)
    def acc(v,x,t):
        return(-np.sin(x)-h*v)
    
    #coeff amortissement
    h =0.5
    
    # plage temporelle
    t=np.linspace(0,20,200)
    
    # les tracés des états liés
    for v0 in np.linspace(.4,1.99,3):
        simulation = odeint(pendule,[v0,0],t)
        plt.plot(simulation[:,1],simulation[:,0],color='b',lw='2')
        simulation = odeint(pendule,[v0,-2*np.pi],t)
        plt.plot(simulation[:,1],simulation[:,0],color='r',lw='2')
        simulation = odeint(pendule,[v0,2*np.pi],t)
        plt.plot(simulation[:,1],simulation[:,0],color='g',lw='2')
    
    # les tracés des états de diffusion
    for v0 in np.linspace(2.4,4,3):
        simulation = odeint(pendule,[v0,-2*np.pi],t)
        plt.plot(simulation[0:900,1],simulation[0:900,0],color='m',lw='2')
        simulation = odeint(pendule,[-v0,2*np.pi],t)
        plt.plot(simulation[0:900,1],simulation[0:900,0],color='c',lw='2')
    
        # esthétisme de l'affichage
    plt.title(u'portrait de phase')
    plt.axis([-2*np.pi,2*np.pi,-4.2,4.2])
    #plt.adjust(left=0.05,top=0.95,bottom=0.1,right=0.95)
    plt.show()



#################################
### pendule fortement amorti, régime libre
#################################

if True:
    # là est l'équation différentielle sous la accélération = fonction(v,x,t)
    def acc(v,x,t):
        return(-np.sin(x)-h*v)
    
    # plage temporelle
    t=np.linspace(0,20,200)
    
    #coeff amortissement
    h =1.5
    
    # les tracés des états liés
    for v0 in np.linspace(.4,1.99,3):
        simulation = odeint(pendule,[v0,0],t)
        plt.plot(simulation[:,1],simulation[:,0],color='b',lw='2')
        simulation = odeint(pendule,[v0,-2*np.pi],t)
        plt.plot(simulation[:,1],simulation[:,0],color='r',lw='2')
        simulation = odeint(pendule,[v0,2*np.pi],t)
        plt.plot(simulation[:,1],simulation[:,0],color='g',lw='2')
    
    # les tracés des états de diffusion
    for v0 in np.linspace(2.4,4,3):
        simulation = odeint(pendule,[v0,-2*np.pi],t)
        plt.plot(simulation[0:900,1],simulation[0:900,0],color='m',lw='2')
        simulation = odeint(pendule,[-v0,2*np.pi],t)
        plt.plot(simulation[0:900,1],simulation[0:900,0],color='c',lw='2')
    
    # esthétisme de l'affichage
    plt.title(u'portrait de phase')
    plt.axis([-2*np.pi,2*np.pi,-4.2,4.2])
    #plt.adjust(left=0.05,top=0.95,bottom=0.1,right=0.95)
    plt.show()



#################################
### oscillateur amorti, régime forcé
#################################

if True:
    # là est l'équation différentielle sous la accélération = fonction(v,x,t)
    def acc(v,x,t):
        return(-x-h*v+2*np.cos(2*np.pi*t))
    
    #coeff amortissement
    h =0.5
    
    # conditions initiales sous la forme (v0,x0)
    CI=[0,1]
    
    # plage temporelle
    t=np.linspace(0,40,2000)
    
    # la résolution
    simulation = odeint(pendule,CI,t)
    
    # tracé
    plt.plot(simulation[:,1],simulation[:,0],color='m',lw='1')

    # esthétisme de l'affichage
    plt.title(u'portrait de phase')
    #plt.axis([-2*np.pi,2*np.pi,-4.2,4.2])
    #plt.adjust(left=0.05,top=0.95,bottom=0.1,right=0.95)
    plt.show()


#################################
### oscillateur faiblement amorti, régime forcé
#################################

if True:
    # là est l'équation différentielle sous la accélération = fonction(v,x,t)
    def acc(v,x,t):
        return(-x-h*v+2*np.cos(2*np.pi*t))
    
    #coeff amortissement
    h =0.1
    
    # conditions initiales sous la forme (v0,x0)
    CI=[0,1]
    
    # plage temporelle
    t=np.linspace(0,40,2000)
    
    # la résolution
    simulation = odeint(pendule,CI,t)
    
    # tracé
    plt.plot(simulation[:,1],simulation[:,0],color='m',lw='1')

    # esthétisme de l'affichage
    plt.title(u'portrait de phase')
    #plt.axis([-2*np.pi,2*np.pi,-4.2,4.2])
    #plt.adjust(left=0.05,top=0.95,bottom=0.1,right=0.95)
    plt.show()


#################################
### oscillateur fortement amorti, régime forcé
#################################

if True:
    # là est l'équation différentielle sous la accélération = fonction(v,x,t)
    def acc(v,x,t):
        return(-x-h*v+2*np.cos(2*np.pi*t))
    
    #coeff amortissement
    h =1.5
    
    # conditions initiales sous la forme (v0,x0)
    CI=[0,1]
    
    # plage temporelle
    t=np.linspace(0,40,2000)
    
    # la résolution
    simulation = odeint(pendule,CI,t)
    
    # tracé
    plt.plot(simulation[:,1],simulation[:,0],color='m',lw='1')

    # esthétisme de l'affichage
    plt.title(u'portrait de phase')
    #plt.axis([-2*np.pi,2*np.pi,-4.2,4.2])
    #plt.adjust(left=0.05,top=0.95,bottom=0.1,right=0.95)
    plt.show()