import matplotlib.pyplot as plt
import math as m
import scipy.optimize as op 
#déclaration des listes
t,M,u,theta,R,X,Y = [],[],[],[],[],[],[]

 # Données d’astronomie de Mercure
T_rev = 0.240 #période de révolution (an) 
e = 0.206 # excentricité
a = 0.387# demi grand axe en (UA)
N = 40# Nombre de positions

# résolution des équations de Kepler 
# détermination de la position de l'astre autour de son orbite
for i in range(N):
	t.append(i*T_rev /N)
	M.append(2*m.pi/T_rev *t[i]) 
	u.append(float(op.fsolve (lambda x:x-e*m.sin(x)-M[i],0) ))
# Calcul des coordonnées polaires
	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]))) 
# calcul des corrdonnées cartésiennes
	X.append(R[i]*m.cos(theta[i])) 
	Y.append(R[i]*m.sin(theta[i]))

# Aires des triangles
# t1,t2 dates première aire et t2, t3 date de la seconde aire 
t1,t2 = 0,2 # intervalle entre la position 0 et 2
t3,t4 = 10,12 # intervalle entre la position 10 et 12 de l'orbite
# les positions peuvent changer mais l'intervalle entre deux dates doit rester le même


AIRE1,AIRE2 = 0,0 # initialisation des aires
i1,i2 = 0,0 

#calcul de l’aire balayée entre t1 et t2	
Delta_t1 =t2-t1 # calcul de l'intevalle de temps
for i1 in range(Delta_t1):
# Calcul des longueur des cotés des triangles
   long1 = m.sqrt((X[t1+i1])**2+(Y[t1+i1])**2)
   long2 = m.sqrt((X[t2+i1])**2+(Y[t2+i1])**2)
   long3 = m.sqrt((X[t2+i1]-X[t1+i1])**2+(Y[t2+i1]-Y[t1+i1])**2) 
# calcul du demi périmètre
   S_1 = 1/2*(long1+long2+long3)
# Calcul de l'aire par la formule de Héron
   AIRE1 = m.sqrt(S_1*(S_1-long1)*(S_1-long2)*(S_1-long3))+AIRE1
	
#calcul de l’aire balayée entre t2 et t3 
Delta_t2 =t4-t3 # calcul de l'intevalle de temps

for i2 in range(Delta_t2):
# Calcul des longueur des cotés des triangles
   long1b = m.sqrt((X[t3+i2])**2+(Y[t3+i2])**2)
   long2b = m.sqrt((X[t4+i2])**2+(Y[t4+i2])**2)
   long3b = m.sqrt((X[t4+i2]-X[t3+i2])**2+(Y[t4+i2]-Y[t3+i2])**2)

# calcul du demi périmètre
   S_1b = 1/2*(long1b+long2b+long3b)

# Calcul de l'aire par la formule de Héron
   AIRE2 = m.sqrt(S_1b*(S_1b-long1b)*(S_1b-long2b)*(S_1b-long3b)) +AIRE2

#Affichage de l'aire calculée (unité : UA au carré)
print('aire balayée entre t1 et t2 --> ' + str(AIRE1)) 
print('aire balayée entre t3 et t4 --> ' + str(AIRE2))

#affichage de l’orbite 
plt.grid(True) 
plt.xlabel("distance (U.A)") 
plt.ylabel("distance (U.A)") 
plt.axis('equal')
plt.plot(X,Y,"bo") 
plt.plot(0,0,"go")
plt.show()