# CAPACITE 3 : FORCE DE CORIOLIS - METHODE D'EULER
###############################################################################
# Déviation vers l'est d'une chute libre 
# Ox vers le sud, Oy vers l'est, Oz vertical du lieu
# Système d'équations différentielles x''= F(x,x',y,y',z,z',t)
#                                     y''= G(x,x',y,y',z,z',t)
#                                     z''= H(x,x',y,y',z,z',t)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Text3D
###############################################################################
# Définition de la fonction "Méthode d'Euler"

# Les paramètres d'entrée sont :
#   - Les trois fonctions F,G,H du système d'équations différentielles
#   - L'instant initial t0
#   - L'instant final tf
#   - Le pas de temps dt 
#   - Les conditions initiales x0, vx0, y0, vy0, z0, vz0

def Euler(F,G,H,t0,tf,dt,x0,vx0,y0,vy0,z0,vz0):   
    
    t = np.arange(t0,tf,dt)      # Tableau des instants ti : t0 , t0+dt , t0+2dt ,....la dernière valeur appartient à l'intervalle [tf-dt,tf[   
    N = len(t)                   # Nombre de valeurs de ti
    
    # Initialisation des variables de sortie
    x = np.zeros(N)              # Mise à zéro du tableau des N valeurs xi de x
    vx = np.zeros(N)             # Mise à zéro du tableau des N valeurs vxi de vx
    y = np.zeros(N)              # Mise à zéro du tableau des N valeurs yi de y
    vy = np.zeros(N)             # Mise à zéro du tableau des N valeurs vyi de vy
    z = np.zeros(N)              # Mise à zéro du tableau des N valeurs zi de z
    vz = np.zeros(N)             # Mise à zéro du tableau des N valeurs vzi de vz
    
    x[0] = x0                    # Condition initiale : le premier élément du tableau x[0] est égal à x0 (c'est à dire x(t=0)) 
    vx[0] = vx0                  # Condition initiale : le premier élément du tableau vx[0] est égal à vx0 (c'est à dire x'(t=0))
    y[0] = y0                    # Condition initiale : le premier élément du tableau y[0] est égal à y0 (c'est à dire y(t=0)) 
    vy[0] = vy0                  # Condition initiale : le premier élément du tableau vy[0] est égal à vy0 (c'est à dire y'(t=0))
    z[0] = z0                    # Condition initiale : le premier élément du tableau z[0] est égal à z0 (c'est à dire z(t=0)) 
    vz[0] = vz0                  # Condition initiale : le premier élément du tableau vz[0] est égal à vz0 (c'est à dire z'(t=0))
    
    # Boucle permettant le calcul des (xi,vxi,yi,vyi,zi,vzi) par récurrence
    
    for i in range(0,N-1):       # i prend N valeurs de 0 à N-1
    
        x[i+1] = x[i] + vx[i]*dt
        vx[i+1] = vx[i] + F(x[i],vx[i],y[i],vy[i],z[i],vz[i],t[i])*dt
        
        y[i+1] = y[i] + vy[i]*dt
        vy[i+1] = vy[i] + G(x[i],vx[i],y[i],vy[i],z[i],vz[i],t[i])*dt
        
        z[i+1] = z[i] + vz[i]*dt
        vz[i+1] = vz[i] + H(x[i],vx[i],y[i],vy[i],z[i],vz[i],t[i])*dt
    
    return t,x,vx,y,vy,z,vz
################################################################################
def F(x,vx,y,vy,z,vz,t):                     # Définition de la fonction F (première ligne du système)
    return 2*omega*np.sin(latitude)*vy
    
def G(x,vx,y,vy,z,vz,t):                     # Définition de la fonction G (Deuxième ligne du système)
    return -2*omega*(np.sin(latitude)*vx+np.cos(latitude)*vz)

def H(x,vx,y,vy,z,vz,t):                     # Définition de la fonction H (Troisième ligne du système)
    return -g+2*omega*np.cos(latitude)*vy

# Constantes

g = 9.81                      # Champ de pesanteur                          (m.s-2)              
omega = 7.27e-5               # Vitesse angulaire de rotation de la Terre   (rad.s-1)
latitude = 0.79               # Latitude                                    (rad)
h = 150                       # Hauteur initiale                            (m)

dt = 0.01                     # Pas de temps                                                            (s)
t0 = 0                        # Instant initial                                                         (s)
tf = np.sqrt(2*h/g)           # Instant final calculé en prenant la cas d'une chute libre simple        (s)

x0 = 0                        # Conditions initiales
vx0 = 0
y0 = 0
vy0 = 0
z0 = h
vz0 = 0

Solution = Euler(F,G,H,t0,tf,dt,x0,vx0,y0,vy0,z0,vz0)               # Résolution

t = Solution[0]
x = Solution[1]                   # Tableau des abscisses en fonction du temps
y = Solution[3]                   # Tableau des ordonnées en fonction du temps
z = Solution[5]                   # Tableau des cotes en fonction du temps

print("La déviation vers l'est vaut :",y[len(t)-1]*100,'cm')
################################################################################
# Tracés de x(t),y(t) et z(t)

plt.subplot(1,3,1)
plt.xlabel('t')
plt.ylabel('x(t)')
plt.grid()
plt.plot(t,x)

plt.subplot(1,3,2)
plt.xlabel('t')
plt.ylabel('y(t)')
plt.grid()
plt.plot(t,y)

plt.subplot(1,3,3)
plt.xlabel('t')
plt.ylabel('z(t)')
plt.grid()
plt.plot(t,z)

plt.show()
################################################################################
# Tracés de z(x) et z(y)

plt.subplot(1,2,1)
plt.xlabel('x')
plt.ylabel('z')
plt.grid()
plt.plot(x,z)

plt.subplot(1,2,2)
plt.xlabel('y')
plt.ylabel('z')
plt.grid()
plt.plot(y,z)

plt.show()
################################################################################
# Tracé de la trajectoire en trois dimensions

fig = plt.figure()
axe = Axes3D(fig)

plt.ticklabel_format(axis='x', style='sci', scilimits=(0, 0))
axe.set_xlabel('x')
axe.set_ylabel('y')
axe.set_zlabel('z')

axe.plot(x,y,z)

plt.show()
################################################################################