import numpy as np
import matplotlib.pyplot as plt

# Parametres du problème et conditions initiales en USI
w0 = 1e3 # en rad/s
Q = 2 # Q > 0.5 => transitoire pseudo-périodique
taupseudo = 2*Q/w0  # tps caractéristique du RT en s
xeq = 0.05  # position d'équilibre en m
x0 = xeq + 0.01 # position initiale : écart de 1 cm / à la position d'équilibre
v0 = 0 # vitesse initiale

# Nombre d’iterations, instant final et pas de temps en USI
n = 1000
tf = 5*2*np.pi/w0  # 5 périodes propres
dt = tf/(n-1)

# Listes des instants 
t = [i*dt for i in range(n)]
# Création d’une liste pour les positions et les vitesses
x = [0 for i in range(n)]
v = [0 for i in range(n)]


# Initialisation
x[0] = x0
v[0] = v0
# Resolution par la methode d’Euler
for i in range(n-1):
    x[i+1] = x[i]+dt*v[i]
    v[i+1] = v[i]*(1-dt*w0/Q)+dt*w0**2*(xeq-x[i])

# Solution analytique
delta=w0**2/Q**2-4*w0**2
wp=np.sqrt(-delta)/2
sol=[xeq+(x0-xeq)*np.exp(-ti/taupseudo)*np.cos(wp*ti) for ti in t]

# Representation graphique
plt.plot(t,x,"--",t,sol)
plt.xlabel("t (s)")
plt.ylabel("x (m)")
plt.title("Euler ED2 oscillateur amorti")
plt.grid()
plt.show()

# Affichage des differentes durees caracteristiques du probleme
print("tau en s:",taupseudo)
print("période propre en s:",2*np.pi/w0)
print("dt en s:",dt)