# import des librairies
import numpy as np
from scipy.integrate import odeint
from scipy import fft
import matplotlib.pyplot as plt

# Données propres au système stable
m = 1.      # masse du système
g = 9.81 # accélération de pesanteur supposée uniforme
T = 86140 # période de rotation propre de la Terre
omega_T = 2*np.pi/T # vitesse angulaire propre de la Terre
lambda_T = 49*np.pi/180 # latitude de Paris
R_T =6.3672e6 # rayon de la Terre

Tmax = 10 # Durée maximale de la chute étudiée
Num = 1000000 # Nombre d'intervalles retenus pour l'étude numérique
t = np.linspace(0,Tmax,Num) # vecteur instants t

xv0=[0,0,160,0,0,0] # conditions initiales sous la forme [x0,y0,z0,vx0,vy0,vz0] selon schéma

def dev_est(xv,t):
    return [xv[3], xv[4], xv[5], 0+2*omega_T*np.sin(lambda_T)*xv[4], 0-2*omega_T*(xv[3]*np.sin(lambda_T)+xv[5]*np.cos(lambda_T)),-g + 2*omega_T*np.cos(lambda_T)*xv[4]]


xv = odeint(dev_est,xv0,t) # résolution de l'équation différentielle

x,y,z = xv[:,0], xv[:,1], xv[:,2] # extraction des résultats de position
dotx,doty,dotz = xv[:,3], xv[:,4], xv[:,5] # extraction des résultats de vitesse

ind=np.argmin(abs(z)) # recherche de l'indice ou z=0, en recherchant le min de |z|
print("Le point d'impact retenu a pour hauteur ", format(z[ind],"#.2e"), "m")

##

#Les lois horaires calculées, il ne reste qu'à exploiter les données précédentes. On commence par tracer la loi horaire de la chute libre jusqu'à l'instant d'impact approché, calculé à l'étape précédente.

plt.figure('Loi horaire de la chute libre')
plt.plot(t[:ind],z[:ind], label=r"$z(t)$")
plt.legend()
plt.xlabel(r"t en $s$")
plt.ylabel(r"z en $m$")
plt.grid(True)
plt.show()

print("La durée calculée de la chute est de ", format(t[ind], "#.3e"), "s")
print("La durée théorique d'une chute libre sans prise en compte des forces de Coriolis est ", format(np.sqrt(2*xv0[2]/g), "#.3e"),"s")

##

#On trace ensuite la projection de la trajectoire sur le plan  pour évaluer la valeur de la déviation vers l'Est que l'on compare avec une valeur approchée (cours)

plt.figure("Déviation vers l'Est lors de la chute libre")
plt.plot(y[:ind],z[:ind], label=r"$z(y)$")
plt.legend()
plt.xlabel(r"y en $m$")
plt.ylabel(r"z en $m$")
plt.grid(True)
plt.show()

print("La déviation vers l'Est calculée est de ", format(y[ind], "#.3e"), "m")
print("La déviation vers l'Est théorique approchée est de ", format(omega_T*np.cos(lambda_T)/3*np.sqrt(8*xv0[2]**3/g), "#.3e"),"m")