import matplotlib.pyplot as plt
import math as m
import scipy.optimize as op

#declaration des listes de données 
t,M,u,theta,R,X,Y,VX,VY,AX,AY = [],[],[],[],[],[],[],[],[],[],[]
# Données d’astronomie
T_revolution = 0.240 #période de révolution en année 
e = 0.206 # excentricité
a = 3.66e12 # demi grand axe en kilomètre (km)
N = 40 # Nombre d’intervalles de temps 

# résolution des équations de Kepler 
for i in range(N+5):
   t.append(i*T_revolution/N)	
   M.append(2*m.pi/T_revolution *t[i]) 
   u.append(float(op.fsolve(lambda x:x-e*m.sin(x)-M[i],0) ))

# Calculs des coordonnées 
   theta.append(2*m.atan((m.sqrt((1+e)/(1-e)) *m.tan(u[i]/2))))
   R.append(a*(1-e**2)/(1+e*m.cos(theta[i]))) 
   X.append(R[i]*m.cos(theta[i])) 
   Y.append(R[i]*m.sin(theta[i]))

# calcul du vecteurs vitesses 
for i in range (1,N+4) :
   VX.append((X[i+1]-X[i-1])/(2*T_revolution/N))
   VY.append((Y[i+1]-Y[i-1])/(2*T_revolution/N)) 
# calcul des vecteurs accélérations
for i in range (1,N+2) :
   AX.append((VX[i+1]-VX[i-1])/(2*T_revolution/N)) 
   AY.append((VY[i+1]-VY[i-1])/(2*T_revolution/N))
#Affichage
plt.figure('Mouvement de Mercure')
plt.axis('equal')
plt.text(0,0 ," Soleil", fontsize = 15,color ='g') 
plt.plot(X,Y,"bo")
plt.plot(0,0,"go")
#Affichage vecteurs accélérations
for i in range (N+1) :
   plt.arrow(X[i+2],Y[i+2],5e-4*AX[i],5e-4*AY[i], head_width=1e11, head_length =1e11, length_includes_head = True)
plt.show()