import numpy as np                 #appel de la bibliothèque numpy
from scipy.integrate import odeint #appel de la bibliothèque scipy
import matplotlib.pyplot as plt    #module matplotlib.pyplot renommé plt

tmax=5.676             #pour une hauteur de chute de 158 m
                       #expérience dans un puits à Freiberg (Allemagne)
N=5000                 #N points de la discrétisation 
lambda1=51*np.pi/180   #51° latitude nord de Freiberg
omega=2*np.pi/86164    #vitesse angulaire de rotation de la terre
g=9.81                 #champ de pesanteur terrestre

def F(Yn, tn):       
    #définition de la fonction F pour odeint
    #Yn[0]=x à l'instant tn, Yn[1]=y et Yn[2]=z
    #Yn[3]=dx/dt à l'instant tn, Yn[4]=dy/dt et Yn[5]=dz/dt
    L=[Yn[3], Yn[4], Yn[5]]
    L=L+[-2*Yn[5]*omega*np.cos(lambda1)+2*Yn[4]*omega*np.sin(lambda1)]
    L=L+[-2*Yn[3]*omega*np.sin(lambda1)]
    L=L+[2*Yn[3]*omega*np.cos(omega)-g]
    return L

Y0=[0, 0, 0, 0, 0, 0]      #conditions initiales
t=np.linspace(0, tmax, N)  #N valeurs comprises entre 0 inclus et tmax inclus
sol=odeint(F, Y0, t)       #tableau à N lignes et 6 colonnes
                           #première colonne : valeurs xn
                           #deuxième colonne : valeurs yn
x=sol[:, 0]        #tableau des x
y=sol[:, 1]        #tableau des y
z=sol[:, 2]        #tableau des z
der_x=sol[:, 3]    #tableau des vx
der_x=sol[:, 4]    #tableau des vy
der_x=sol[:, 5]    #tableau des vz
 
plt.figure()  #nouvelle fenêtre graphique
plt.grid()    #affiche la grille
plt.title("x(t)")
plt.xlabel("t")
plt.ylabel("x")
plt.plot(t, x)
plt.show()  #affiche la figure à l'écran

plt.figure()  #nouvelle fenêtre graphique
plt.grid()    #affiche la grille
plt.title("y(t)")
plt.xlabel("t")
plt.ylabel("y")
plt.plot(t, y)
plt.show()  #affiche la figure à l'écran


plt.figure()  #nouvelle fenêtre graphique
plt.grid()    #affiche la grille
plt.title("z(t)")
plt.xlabel("t")
plt.ylabel("z")
plt.plot(t, z)
plt.show()  #affiche la figure à l'écran

print(x[-1])  #on obtient 2,7 mm vers l'est
print(y[-1])  #on obtient 4,3 µm vers le sud
print(z[-1])  #hauteur de chute 158 m
print()

#Approximation du premier ordre Wikipédia déviation vers l'est
print(2/3*omega*tmax*158*np.cos(lambda1))
#Approximation du deuxième ordre Wikipédia déviation vers l'est
print(2*158**2*omega**2/3/g*np.sin(lambda1)*np.cos(lambda1))