import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as scin

# Parametres du problème et conditions initiales en USI
g = 9.81
l = 0.2
T0 = 2*np.pi*np.sqrt(l/g)
theta0 = 45*np.pi/180
w0 = 0

# Nombre d’iterations, instant final et pas de temps en USI
n = 3*1000+1
tf = 3*T0  # 3 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 angles et les vitesses angulaires
theta = [0 for i in range(n)]
w = [0 for i in range(n)]

# Initialisation
theta[0] = theta0
w[0] = w0
# Resolution par la methode d’Euler
for i in range(n-1):
    theta[i+1] = theta[i]+dt*w[i]
    w[i+1] = w[i]-dt*g/l*np.sin(theta[i])


# Définition de la fonction F(y,t)
def F(y,t):
    theta,omega=y
    yp=[omega,-g/l*np.sin(theta)]
    return(yp)

# Resolution par la fonction odeint
theta2=scin.odeint(F,[theta0,w0],t)



# Solution analytique dans le cas des petits angles
theta3 = [theta0*np.cos(2*np.pi/T0*ti) for ti in t]


# Representation graphique
plt.plot(t,theta,t,theta2[:,0],t,theta3)
plt.xlabel("t (s)")
plt.ylabel("theta (rad)")
plt.legend(["Euler","odeint","OH"],loc="best")
plt.grid()
plt.show()

# Affichage des differentes durees caracteristiques du probleme
print("période propre en s:",T0)
print("dt en s:",dt)