Calculate eccentricity vector of planetary ring particles

Question

I wish to set-up an initially circular (e=0) system of planetary rings which I can later perturb over time and see how the eccentricity changes. However, my calculation of the eccentricity vector returns -1 as the value of my initial ring, rather than zero.

The eccentricity vector equation takes this form

import numpy as np
import matplotlib.pyplot as plt

G = 6.674e-20 # km^3 kg^-1 s^-2
day = 60.0 * 60.0 * 24.0
dt = day / 10.0 
Mass = 5.683e26

N = 30000
delta = np.random.random(1) * 2.0 * np.pi / N
angles = np.linspace(0.0, 2.0 * np.pi, N) + delta

radius = np.random.uniform(low = 1e6, high = 2e6, size = (N)) # ring radius
xrings, yrings = radius * np.cos(angles), radius * np.sin(angles) # positions
vxrings, vyrings = np.sqrt((G * Mass) / radius) * -np.sin(angles), np.sqrt((G * Mass) / radius) * np.cos(angles) # velocities

dist = np.hypot(xrings, yrings) # distance between particles

# update positions
xrings += (vxrings * dt)
yrings += (vyrings * dt)

#update velocities
vxrings -= (G * Mass * xrings / (dist ** 3.0 + 1.0e6) * dt)
vyrings -= (G * Mass * yrings / (dist ** 3.0 + 1.0e6) * dt)

v = np.hypot(vxrings, vyrings) 
mu = G*Mass
e = (((abs(v)**2) / mu) - (1/abs(dist)))*radius - (((radius*v) / mu)*v)


plt.plot(e, radius)
plt.show()

I have tried interchanging dist and radius in various ways within the equation as I know the radius needs to be with respect to the central mass, but to no avail.

Answer 1

I think the problem is arising due to the fact that it is a vector equation and when you have implemented it, you've used the magnitudes of radius and velocity when you should have used their vectors. Implementing either equation from the wiki (with the vectors for r and v) gives the expected result of e being 0 when dt is 0:

import numpy as np
import matplotlib.pyplot as plt

G = 6.674e-20 # km^3 kg^-1 s^-2
day = 60.0 * 60.0 * 24.0
dt = day / 10.0
Mass = 5.683e26
mu = G*Mass
dt = 0

N = 30000
delta = np.random.random(1) * 2.0 * np.pi / N
angles = np.linspace(0.0, 2.0 * np.pi, N) + delta

radius = np.random.uniform(low = 1e6, high = 2e6, size = (N)) # ring radius
xrings, yrings = radius * np.cos(angles), radius * np.sin(angles) # positions

vxrings, vyrings = np.sqrt((G * Mass) / radius) * -np.sin(angles), np.sqrt((G * Mass) / radius) * np.cos(angles) # velocities

dist = np.hypot(xrings, yrings) # distance between particles

# update positions
xrings += (vxrings * dt)
yrings += (vyrings * dt)

# #update velocities
vxrings -= (G * Mass * xrings / (dist ** 3.0 + 1.0e6) * dt)
vyrings -= (G * Mass * yrings / (dist ** 3.0 + 1.0e6) * dt)

# Convert to array of vectors assuming there is no motion out of the plane
r_vector = np.array([[i, j, 0 ] for i, j in zip(xrings, yrings)])
v_vector = np.array([[i, j, 0] for i, j in zip(vxrings, vyrings)])

# Implement the equation as given in the wiki page
# Cross product method
h = [np.cross(i, j) for i, j in zip(r_vector, v_vector)] # r cross v
v_h = [np.cross(i, j)/mu for i, j in zip(v_vector, h)] # v cross h over mu
r_normalised = [i/np.linalg.norm(i) for i in r_vector]
e_vector_cross = [i-j for i,j in zip(v_h, r_normalised)]
absolute_e_cross = [np.linalg.norm(i) for i in e_vector_cross]

plt.figure()
plt.title('Cross product method')
plt.xlabel('Eccentricity')
plt.ylabel('Radius')
plt.plot(absolute_e_cross, radius)

# Dot product method
first_factor = [np.linalg.norm(i)**2/mu -1/np.linalg.norm(j) for i, j in zip(v_vector, r_vector)]
first = [i*j for i, j in zip(first_factor, r_vector)]

second_factor = [np.dot(i, j)/mu for i, j in zip(r_vector, v_vector)]
second = [i*j for i, j in zip(second_factor, v_vector)]
e_vector_dot = [i-j for i, j in zip(first, second)]
absolute_e_dot = [np.linalg.norm(i) for i in e_vector_dot]
plt.figure()
plt.title('Dot product method')
plt.xlabel('Eccentricity')
plt.ylabel('Radius')
plt.plot(absolute_e_dot, radius)

Output: