Displaying orbit with vpython using kepler's equation but the planet won't orbit

Question

I've been trying to make some sort of simulation of the solar system, and I came across Kepler's equation in order to give me the position of the planet at any given time. I used the Newton-Raphson method in order to calculate the eccentric anomaly. However, I can't seem to make the planet actually orbit properly, and I am not sure what is wrong. The planet simply goes back and forth in a very tiny spot.

I appreciate the help in advance.

from datetime import datetime, timedelta
from math import *
from vpython import *

 scene = canvas(title='Solar system simulation',
 width=1024, height=720,
 center=vector(0,0,0), background = color.black)

 sun = sphere(pos = vector(0,0,0), radius = 5, color = color.yellow)
 earth = sphere(pos=vector(1,1,0), radius = 1, color = color.blue, make_trail = True)
 earth.trail_color = color.white

 d1 = datetime.now()

while True:

    d1 = d1 + timedelta(minutes = 1)
    d2 = datetime(2000, 1, 1, 00, 00)

    SecondsFrom2000 = (d1 - d2).total_seconds()  #seconds between today and 01.01.2000 at midnight
    CenturiesFrom2000 = SecondsFrom2000/(60*60*24*365.25*100) #centuries between today and 2000
    e0Earth = 0.01671123 #eccentricity
    edotEarth = -0.00004392
    eEarth = e0Earth + edotEarth*CenturiesFrom2000

    a0Earth = 1.00000018 #semi-major axis[au], from 3000 BC - 3000 AD
    adotEarth = -0.00000003 #[au/century]
    aEarth = a0Earth + adotEarth*CenturiesFrom2000

    L0Earth = 100.46691572 #mean longitude [deg]
    LdotEarth = 35999.37306329 #[deg/century]
    LEarth = (L0Earth + LdotEarth*CenturiesFrom2000)*(pi/180) #[rad/century]

    Pi0Earth = 102.93005885 #longitude of the perihelion [deg]
    PidotEarth = 0.31795260 #[deg/century]
    PiEarth = (Pi0Earth + PidotEarth*CenturiesFrom2000)*(pi/180)

    W0Earth = -5.11260389 #longitude of the ascending node [deg]
    WdotEarth = -0.24123856
    WEarth = (W0Earth + WdotEarth*CenturiesFrom2000)*(pi/180)

    I0Earth = -0.00054346 #inclination [deg]
    IdotEarth = -0.01337178
    IEarth = (I0Earth + IdotEarth*CenturiesFrom2000)*(pi/180)

    MEarth = LEarth - PiEarth #mean anomaly
    wEarth = PiEarth - WEarth #argument of perihelion

    E0 = 0.1
    E1 = MEarth #E = eccentric anomaly
    error = 1

    while error > 0.000001:
        E2 = E1 - (E1 - eEarth*sin(E1)/(1 - eEarth*cos(E1)))
        E1 = E2
        erreur = abs(E2 - E1)

    E_Earth = E2

    PEarth = aEarth * (cos(E_Earth) - eEarth) #[au]
    QEarth = aEarth * sin(E_Earth)*sqrt(1 - eEarth*eEarth) #[au]

    xEarth = cos(wEarth)*PEarth - sin(wEarth)*QEarth
    yEarth = sin(wEarth)*PEarth + cos(wEarth)*QEarth
    zEarth = sin(IEarth)*xEarth
    xEarth = cos(IEarth)*xEarth
    xtempEarth = xEarth
    xEarth = cos(WEarth)*xtempEarth - sin(WEarth)*yEarth
    yEarth = sin(WEarth)*xtempEarth + cos(WEarth)*yEarth

    earth.pos = vector(xEarth*10,yEarth*10,zEarth*10)

Answer 1

Your program is too complicated for me to follow, but I would recommend inserting some print statements in the hope that they would uncover where your code is doing something different than want you wanted.