import matplotlib.pyplot as plt
import math as m
import numpy as np
import scipy.optimize as op 
# déclaration des listes
t, M, u, theta, R, X, Y = [], [], [], [], [], [], [] 
# données d'astronomie
T_rev = 0.240   #période de révolution (an) 
e = 0.206       #excentricité
a = 0.387       #demi-grand axe (UA) 
N = 40          #nombre de positions
#résolution des équations de Kepler
for i in range(N):
    t.append(i*T_rev/N)
    M.append(2*m.pi/T_rev*t[i])
    u.append(op.fsolve(lambda x: x - e*np.sin(x) -M[i] , 0)[0])
    # calcul 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]))
#affichage de l'orbite
plt.grid(True)
plt.xlabel("distance (UA)")
plt.ylabel("distance (UA)")
plt.axis('equal')
plt.plot(X, Y, "bo")
plt.plot(0, 0, "go")
plt.show()