计算行星环粒子的偏心矢量

Question

我希望建立一个初始为圆形 (e=0) 的行星环系统，稍后我可以随时间对其进行扰动，并查看偏心率如何变化。但是，我计算的偏心向量returns -1 作为我初始环的值，而不是零。

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()

我已经尝试在方程中以各种方式交换 dist 和 radius，因为我知道半径需要相对于中心质量，但无济于事。

Answer 1

我认为问题的出现是因为它是一个向量方程，当你实现它时，你在应该使用它们的向量时使用了半径和速度的大小。从 wiki 实现任一方程（使用 r 和 v 的向量）给出当 dt 为 0 时 e 为 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)

输出：

计算行星环粒子的偏心矢量

Calculate eccentricity vector of planetary ring particles

python

orbital-mechanics