# DEVIATION VERS L'EST, CALCUL NUMÉRIQUE.

import numpy as np
import matplotlib.pyplot as plt

g=9.806                # g locale à Freiberg
h=158                  # altitude initiale, z=0 au fond du puits
tmax=((2*h/g)**0.5)    # temps de chute libre approximatif
n=1000
dt=tmax/n
w=2*np.pi/86164        # vitesse angulaire de la Terre
l=(50+54/60)/180*np.pi # latitude 50° 54' de Freiberg en radians

t=0
vx=0
vy=0
vz=0
x=0
y=0
z=h

lx=[x]
ly=[y]

# calcul par Euler simple...
while z>0:                     # tant que la masse n'est pas au fond
    ax=2*w*np.sin(l)*vy
    ay=-2*w*(np.sin(l)*vx+np.cos(l)*vz)
    az=-g+2*w*np.cos(l)*vy

    vx=vx+ax*dt
    vy=vy+ay*dt
    vz=vz+az*dt

    x=x+vx*dt
    y=y+vy*dt
    z=z+vz*dt

    t=t+dt

    lx.append(x*1000) # "enregistré" en millimètres
    ly.append(y*1000)

# vue de dessus...
plt.figure(figsize=(15,7))
plt.xlim([-50,50])
plt.ylim([-5e-2,5e-2])
plt.plot(ly,lx)
plt.scatter(y*1000,x*1000,marker='+',c='r',s=200)

titre="Déviation en millimètres vue de dessus - temps de chute : "\
+str(round(t,2))+" s - vers l'Est : "+str(round(y*1000,1))\
+" mm, vers le Sud : "+str(round(x*1e6,1))+" micromètres !"

plt.title(titre)
plt.ylabel("vers le Sud zoomé 1000 fois !")
plt.xlabel("vers l'Est")
plt.grid()
plt.show()