Inaccurate motion of planets in orbit animations

[DraconianOrder]

Functions that animate the motion of a planet in orbit do not appear to follow Kepler's Second Law, where the planet should travel faster when closer to the Sun.

Animated orbit of Pluto, saved from planet.animate_orbit() An orbit with a higher eccentricity clearly shows a higher velocity when further away from the Sun, and a lower velocity when near it, the opposite of Kepler's Second Law.

Maths on slideshow When comparing to Kepler's own work, this appears to be an oversimplification, which could be the cause of the inaccuracy.

Code sample from planet.animate_orbit(self)

frames = int(1000 * 50 / self.period)
a = self.sm_axis
e = self.eccentricity
time = np.linspace(0, 1000, frames + 1)
...
...
theta = 2 * np.pi * time[frame] / self.period
r = a * (1 - e ** 2) / (1 - e * np.cos(theta))
x = r * np.cos(theta)
y = r * np.sin(theta)

Possible causes:

• Slideshow's inaccurate maths (highly likely)
• Incorrect interpretation/implementation of slideshow's maths
• Matplotlib/Python lag/unexpected lapses from low intervals (20 ms)

Tasks affected: 2, 3, 4, 6, 7 unconfirmed and hypothetical, dependent on cause

[DraconianOrder]

Slideshow oversimplified Kepler's calculations; this has been fixed with the following function.

# returns true anomaly and heliocentric distance as a function of time
# solves kepler's equation using kepler's 1621 fixed-point iteration
def kepler_eq(time, sm_axis, period, eccentricity):
    n = 2 * np.pi / period  # mean motion
    M = n * time  # mean anomaly

    # kepler's equation: E = M + eccentricity * sin(E)
    e = eccentricity
    E = M
    for _ in range(10):  # 10 iterations balances accuracy with speed
        E = M + e * np.sin(E)  # eccentric anomaly

    # true anomaly theta
    theta = 2 * np.arctan(np.sqrt((1 + e) / (1 - e)) * np.tan(E / 2))
    a = sm_axis
    r = a * (1 - e * np.cos(E))  # heliocentric distance

    return theta, r

In non-animations, r has now been written like this r = a * (1 - e ** 2) / (1 + e * np.cos(theta)) instead of using the slideshow's equations: r = a * (1 - e ** 2) / (1 - e * np.cos(theta))

Planetary orbits now appear to follow Kepler's Second Law, like in this example