import numpy as np
import matplotlib.pyplot as plt

# Parametres
g = 9.8 ; omega = 7.3e-5 ; R = 6370. ; h = 1. # A MODIFIER
l = 10 #latitude en degres, A MODIFIER
cosl = np.cos(l*np.pi/180)
sinl = np.sin(l*np.pi/180)

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

# Calcul du pas de temps Delta_t en fonction d'un entier N
# Pour cela on se base sur une chute libre classique, en ref gal
N = 2000
Temps_chute = np.sqrt(2*h/g)
Delta_t=Temps_chute / N


def euler(t0,x0,y0,z0,vx0,vy0,vz0,Delta_t) :

    t=[t0]
    x=[x0] ; y=[y0] ; z=[z0]
    vx=[vx0] ; vy=[vy0] ; vz=[vz0]
    i=0
    while z[i] > 0 :
        t.append( t[i] + Delta_t )
        x.append( x[i]  )     # A MODIFIER
        y.append( y[i] + Delta_t*vy[i] )
        z.append( z[i] + Delta_t*vz[i] )
        vx.append( vx[i] )    # A MODIFIER
        vy.append( vy[i] + Delta_t * ( - 2 * omega * sinl * vx[i]) )
        vz.append( vz[i] + Delta_t * ( - g + 2 * omega * cosl * vx[i]) )
        i += 1
    
    return (t, x, y, z, vx, vy, vz)


sol=euler(t0,x0,y0,z0,vx0,vy0,vz0,Delta_t)

print('La derniere valeur de x est', sol[1][-1], 'm')
print('La derniere valeur de y est', sol[2][-1], 'm')

# Pour la question plus difficile (*) :
print('La derniere valeur de vx est', sol[4][-1], 'm/s')
print('La derniere valeur de vy est', sol[5][-1], 'm/s')
print('La derniere valeur de vz est', sol[6][-1], 'm/s')

plt.plot(sol[0],sol[1],'--')
plt.plot(sol[0],sol[2],':')

plt.grid()
plt.legend(('x(t)','y(t)'),loc=4)
plt.show()