How to plot an orbit with pyephem/matplotlib

Question

With pyephem I can take orbital elements from the Minor Planet Center and add them into a ephem.EllipticalBody().

To plot the orbit I've been getting the computing the sun distance, heliocentric latitiude and longitude over a the period of the orbit and plotting that:

dt = body._epoch
period = (sqrt(a**3))

timespace  = np.linspace(dt-((period*365)/2),
                         dt+((period*365)/2), 720)

theta_e = []
r_e = []
phi_e = []
for t in timespace:
    body.compute(t)
    theta_e.append(body.hlon)
    phi_e.append(body.hlat)
    r_e.append(body.sun_distance)

subplot = figure(figsize=(20, 20)).add_subplot(111, polar=True)
subplot.scatter(theta_e, r_e*cos(phi_e), s=0.5)

This approch works ok for minor planets that have short periods (a few years). However when I go to the extreme and plot say Sedna I get:

https://i.stack.imgur.com/wxRCW.png

So it looks as though pyephem is being clever and taking apsidal precession into account. Which isn't something I need for this.

Any suggestions so on how I can plot the orbit from the orbital elements or what I'm doing wrong here? Ideally I'll like to be able to do plots from different 'views' to show the inclination as well as eccentricity or perhaps a live mayavi scene.

Use Kepler's equation look at simulation part of this answer of mine [Is it possible to make realistic n-body solar system simulation in matter of size and mass?](http://stackoverflow.com/a/28020934/2521214) and all its sub-links. From that loop `M` through full circle and compute `E` then compute coordinate on the ellipsis and rotate by node,inclination,... then just plot curve along these points ... — Spektre, Sep 12 '15 at 07:45

score 0 · Answer 1 · answered Sep 14 '15 at 05:12

I suspect that the phenomenon here is not apsidal precession, because the libastro library beneath PyEphem uses simple Keplerian orbits when given elliptical coordinates.

Instead, my guess is that you are seeing the Earth’s polar precession, which is moving the heliocentric latitude and longitude system out from under you as you are trying to plot. For most objects the effect is negligible, but an 11,000-year orbit is enough time for the Earth’s pole to make it a good bit of the way around the sky.

The only coordinates that libastro produces that are likely to be immune to the effect (assuming that a compute(…, epoch=J2000) does not make the heliocentric longitude stable?) are the a_ra and a_dec coordinates, so you might have to pull those for both the Sun and Sedna, convert them into vectors, and subtract.

How to plot an orbit with pyephem/matplotlib

1 Answers1