# -*- coding: utf-8 -*-
"""
Resolution du pendule pesant par methode d'Euler et odeint'
"""

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

## Parametres et CI
w0 = 2           # pulsation propre en s**-1
theta0 = 1       # angle initial en rad
thetapoint0 = 0  # vitesse angulaire initiale en s**-1
t0 = 0                 # instant initial
tf = 30                # instant final
N = 10000              # nombre d'iterations
t = np.linspace(t0, tf, N+1)       # instants
dt = (tf-t0)/N         # pas de temps

## Resolution par methode d'Euler
theta = [theta0]       
thetapoint = [thetapoint0]  
for i in range(N):
    thetapoint.append(thetapoint[i] - dt*np.sin(theta[i])*w0**2)
    theta.append(theta[i] + dt*thetapoint[i])

## Resolution par odeint, X = [thetapoint, theta]
CI = [thetapoint0,theta0]
def F(X,t):
    return [-np.sin(X[1])*w0**2,X[0]]

solution = odeint(F, CI, t)

plt.plot(t,theta, label = 'Euler') 
plt.plot(t,solution[:,1], label = 'odeint')
plt.title('Résolution par Euler et odeint')
plt.xlabel('t(s)')
plt.ylabel('theta')
plt.legend(loc = "upper left")
plt.show()