# -*- coding: utf-8 -*-
"""
Created on Thu Feb  9 11:36:40 2023

@author: vleprince
"""
import numpy as np
import matplotlib.pyplot as plt 


"""1.1 première méthode"""

def Euler_pendule(l,T,theta_0,omega_0,n):
    g = 9.81
    h = T/n
    list_theta = [0 for i in range(n+1)]
    list_theta[0] = theta_0
    list_omega = [0 for i in range(n+1)]
    list_omega[0] = omega_0
    for i in range(0,n):
        list_theta[i+1] = list_theta[i] + h * list_omega[i]
        list_omega[i+1] = list_omega[i] + \
                              h * (-g/l)*np.sin(list_theta[i])
        
    return [list_theta , list_omega]


# tracé avec les valeurs données :
l = 1.
T = 5.
theta_0 = 0.5
omega_0 = 0
n = 5000  # on peut faire varier n pour observer les conséquences



Sol_theta , Sol_omega = Euler_pendule(l,T,theta_0,omega_0,n)

h = T/n
temps = [h * i for i in range(n+1) ]

plt.plot(temps,Sol_theta)



"""1.2 deuxième méthode : avec numpy"""

# g et l sont ici définis de manière globale ;
# la fonction F en dépend
l = 1.
g = 9.81

def F(Y,t):
    theta , omega = Y
    return np.array([omega , -(g/l)*np.sin(theta)])
    
def Euler_pendule2(T,theta_0,omega_0,n):
    h = T/n
    Sol = np.zeros((2,n+1))
    Sol[:,0] = [theta_0,omega_0]
    for i in range(0,n):
        Sol[:,i+1] = Sol[:,i] + h*F(Sol[:,i], h*i)
        
    return Sol

T = 5.
theta_0 = 0.5
omega_0 = 0
n = 500

SolEuler = Euler_pendule2(T,theta_0,omega_0,n)
temps = np.linspace(0,T,n+1)
plt.plot(temps, SolEuler[0,:])

#portrait de phase :
plt.plot(SolEuler[0,:],SolEuler[1,:])