Simulation de l'orbite terrestre

Bonjour,

Je souhaite modéliser l'orbite terrestre sur python. J'ai fais une simulation sur python; le problème est que la terre dévie un peu de son orbite et au bout de quelque période la terre, l'orbite se décale d'une 100 de millions de kilomètres (ça fais bcp :) )
Pour la simulation je n'ai utilisé que la mécanique newtonienne. Je n'ai pas pris en compte les autres planètes. Le soleil est fixe. J'ai utilisé la méthode d'Euler pour intégrer.

Je précise que je suis d'un niveau d'un bon élève de prépa math sup en python.

Voici mon code, si vous pouviez jeter un œil pour me dire comment le corriger. Je suis pas sur mais je crois que l'indentation n'arrive pas à passer sur le forum, sorry. J'ai fais comme j'ai pu j'espère que ça n'a pas rajouté des erreurs d'indentations.

# initialisation des constantes
#ra=rayon aphélie
#rb=rayon périphélie

ra=152093407000
rb=147093407000
a=(ra+rb)/2

#masse de la terre
mt=5.92*10**24

#masse du soleil
ms=1.9884*10**30

#constante gravitationnelle
G=6.6742*10**(-11)


va=28851 #vitesse aphélie

ax=[]
ay=[]


vx=[0]
vy=[va]


x=[ra]
y=[0]


d=1

F=G*ms/(d**2) 

theta=0



import numpy as np
import matplotlib.pyplot as plt


dt=3600
for k in range(300000):




         ---d=(x[k]**2+y[k]**2)**0.5

         ---F=G*ms/(d**2)

 
         ---theta=np.abs(np.arctan(y[k]/x[k]))
 
         ---if x[k]<=0:
             ------if y[k]>=0:
                 ---------ax.append(F*np.cos(theta))
                 ---------ay.append(-F*np.sin(theta))

             ------if y[k]<0:
                 ---------ax.append(F*np.cos(theta))
                 ---------ay.append(F*np.sin(theta))
         ---if x[k]>0:
             ------if y[k]>=0:
                 ---------ax.append(-F*np.cos(theta))
                 ---------ay.append(-F*np.sin(theta))
        ------if y[k]<0:
            ---------ax.append(-F*np.cos(theta))
            ---------ay.append(F*np.sin(theta))


 
         ---vx.append(vx[k]+dt*ax[k])
         ---vy.append(vy[k]+dt*ay[k])



         ---x.append(x[k]+dt*vx[k])
         ---y.append(y[k]+dt*vy[k])



         ---plt.scatter(0,0,s=50,c='black',marker='+')

         ---plt.axis([-2*10**11,2*10**11,-2*10**11,2*10**11])
         ---plt.plot(x,y)

         ---plt.pause(0.01)

plt.show()