# -*- coding: utf-8 -*-
"""
Created on Tue Jul 19 11:42:04 2022

@author: starons
"""

"""
Ce programme résoud numériquement l'équation différentielle de l'oscillateur 
harmonique non amorti en utilisant la fonction odeint et compare à
la solution analytique. 
"""

# On importe les bibliothèques utiles
# -----------------------------------
import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import odeint

# Valeurs numériques caractéristiques du problème
# -----------------------------------------------
omega0 = 2                                  # pulsation propre
x0 = 4                                      # conditions initiales 
xp0 = 0

# Paramètres de la résolution numérique
# ------------------------------------- 
t_fin = 15/omega0                           # durée
N = 100                                     # nombre de points
t = np.linspace(0,t_fin,N)                  # découpage du temps
W0 = [x0,xp0]                               # conditions initiales

# Définition du système différentiel d'ordre 1 à résoudre 
# ------------------------------------------------------- 
def OH (W,t):
    Wp = np.zeros(2)                        # tableau vide pour stocker les dérivées
    Wp[0] = W[1]                            # première équation
    Wp[1] = -omega0**2*W[0]                 # deuxième équation 
    return Wp
    
# Intégration numérique à l'aide d'odeint
# ---------------------------------------
Wsol=odeint(OH,W0,t)                        # la ligne qui fait le job ! 
x=Wsol[:,0]                                 # extraction de la position
xp=Wsol[:,1]                                # extraction de la vitesse

# On connait une solution analytique
# ---------------------------------- 
t_analytique = np.linspace(0,t_fin,1000)# pour disposer d'un tableau plus dense de valeurs de temps   
x_analytique = x0*np.cos(omega0*t_analytique)+xp0/omega0*np.sin(omega0*t_analytique)
xp_analytique = -x0*omega0*np.sin(omega0*t_analytique)+xp0*np.cos(omega0*t_analytique)

# Représentations graphiques
# --------------------------
plt.subplot(1,2,1)
plt.plot(t,x,color='black',marker='.',label="solution numérique")
plt.plot(t_analytique,x_analytique,color='red',linestyle='dashed',label="solution analytique")
plt.xlabel("temps"),plt.ylabel("position")
plt.title("Oscillateur harmonique non amorti")
plt.legend(loc="lower left"),plt.grid()

plt.subplot(1,2,2)
plt.plot(x,xp/omega0,color='black',marker='.',label="solution numérique")
plt.plot(x_analytique,xp_analytique/omega0,color='red',linestyle='dashed',label="solution analytique")
plt.axis("equal")
plt.xlabel("position"),plt.ylabel("vitesse")
plt.title("Trajectoire de phase")
plt.legend(),plt.grid()