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

## Definition des constantes du probleme
g=9.81                  # en m.s^-2
l=67                    # en m
omega=7.3e-5            # en rad.s^-1
lambd=48.52*np.pi/180   # en rad
a=4                     # en m

## Definition de l'ensemble des valeurs t (en s) pour lesquelles
## on cherche la solution
t=np.linspace(0,6*3600,10000)  # 10000 points entre 0 et 6 h

## Definition du systeme differentiel
def systDiff(pos,t):
    # pos est le vecteur position inconnu de dimension 4
    # (pos[0] : x, pos[1] : y, pos[2] : u, pos[3] : v)
    # t est le vecteur contenant les instants t de résolution
    # La fonction renvoit la derivée 1ere temporelle du vecteur pos
    dx=pos[2]
    dy=pos[3]
    du=-g/l*pos[0]+2*omega*np.sin(lambd)*pos[3]
    dv=-g/l*pos[1]-2*omega*np.sin(lambd)*pos[2]
    return([dx,dy,du,dv])

## Definition des conditions initiales
CI=[a,0,0,0]

## Resolution
sol=odeint(systDiff, CI, t)
x=sol[:,0]      # Extraction des valeurs x pour chaque instant t
y=sol[:,1]      # Extraction des valeurs y pour chaque instant t

## Trace graphique
plt.figure(figsize=(8,8))
plt.plot(x,y)
plt.title('Trajectoire du pendule de Foucault',fontsize=15)
plt.xlabel('x (en m)',fontsize=15)
plt.ylabel('y (en m)',fontsize=15)
plt.xticks(size=15)
plt.yticks(size=15)
plt.axis('equal')
plt.grid()
plt.tight_layout()
plt.show()

# print(np.max(np.abs(sol[:,2])))