I have been writing a Python script (GitHub LINK) for visualizing asteroid/comet/meteoroid orbits. The script also plots the position of planets and their orbits.
It works just right for orbits with small semi-major axis (i.e. "smaller" orbits). But when I have an orbit that goes way beyond Neptune (e.g. of a Halley-type comet), and from certain perspectives, there is a weird "wraparound" (for lack of a better word) effect.
Let me show you what I mean:
Image compilation: https://i.stack.imgur.com/1WgX0.png
This image shows the plot from a perspective where it does not break.
When you rotate the same plot a bit to the right, it is as if the orbit folded in half and reversed its direction!
And if you look at the plot from a great distance, you can see that the elipse is plotted as it should be.
And here is a minimal version of the code with which the issue can be reproduced. The "wraparound" occurs only when the perspective of the camera is closely parallel with the large orbit.
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
def orbitalElements2Cartesian(a, e, I, peri, node, E):
""" Convert orbital elements to Cartesian coordinates in the Solar System.
Args:
a (float): semi-major axis (AU)
e (float): eccentricity
I (float): inclination (degrees)
peri (float): longitude of perihelion (degrees)
node (float): longitude of ascending node (degrees)
E (float): eccentric anomaly (radians)
"""
# The source of equations used:
# http://farside.ph.utexas.edu/teaching/celestial/Celestialhtml/node34.html
# Check if the orbit is parabolic or hyperbolic
if e >=1:
e = 0.99999999
# Convert degrees to radians
I, peri, node = map(np.radians, [I, peri, node])
# True anomaly
theta = 2*np.arctan(np.sqrt((1.0 + e)/(1.0 - e))*np.tan(E/2.0))
# Distance from the Sun to the poin on orbit
r = a*(1.0 - e*np.cos(E))
# Cartesian coordinates
x = r*(np.cos(node)*np.cos(peri + theta) - np.sin(node)*np.sin(peri + theta)*np.cos(I))
y = r*(np.sin(node)*np.cos(peri + theta) + np.cos(node)*np.sin(peri + theta)*np.cos(I))
z = r*np.sin(peri + theta)*np.sin(I)
return x, y, z
if __name__ == '__main__':
# Example orbital elements
# a, e, incl, peri, node
orb_elements = np.array([
[2.363, 0.515, 4.0, 205.0, 346.1],
[0.989, 0.089, 3.1, 55.6, 21.2],
[0.898, 0.460, 1.3, 77.1, 331.2],
[104.585332285, 0.994914, 89.3950, 130.8767, 282.4633]
])
# Setup the plot
fig = plt.figure()
ax = fig.gca(projection='3d')
# Eccentric anomaly (full range)
E = np.linspace(-np.pi, np.pi, 100)
# Plot the given orbits
for i, orbit in enumerate(orb_elements):
a, e, I, peri, node = orbit
# Take extra steps in E if the orbit is very large
if a > 50:
E = np.linspace(-np.pi, np.pi, (a/20.0)*100)
# Get the orbit in the cartesian space
x, y, z = orbitalElements2Cartesian(a, e, I, peri, node, E)
# Plot orbits
ax.plot(x, y, z, c='#32CD32')
# Add limits (in AU)
ax.set_xlim3d(-5,5)
ax.set_ylim3d(-5,5)
ax.set_zlim3d(-5,5)
plt.tight_layout()
plt.show()
I am a bit dumbfounded by this and cannot seem to find a proper solution. I would greatly appreciate some help!
