问题描述
请原谅标题的长度，但这是一个非常具体的问题。我目前正在模拟在 2022 年发射窗口向火星发射火箭，我注意到我的火箭离火星很远，即使它朝着正确的方向行进。在简化我的代码以缩小问题的范围后，我简单地绘制了地球和火星的轨道（使用来自 NASA 的 SPICE 库的数据）并传播了我实现的 LAmbert 求解器（通用变量）给我的位置和速度来绘制火箭的最终轨道。
我只是让太阳的引力影响火箭，而不是地球或火星，以尽量减少我的问题空间。然而，尽管到目前为止我已经简化了我的问题，但火星和我的火箭轨道之间的交叉点发生在飞行时间之前，一路模拟，两个物体之间的最小距离超过一百万公里任何时候。
话虽如此，一定是出了什么问题，但我找不到问题所在。我通过将它与 Dario Izzo 的方法进行比较来确保我复制的 LAmbert 求解器代码是正确的，并且两者都给出了相同的结果。此外，我还通过传播火星和地球的轨道并将这些椭圆与来自 SPICE 的数据进行比较来检查我的轨道传播器是否工作。
总而言之，我认为这一定是我在某个地方犯的一个愚蠢的小错误，但由于我缺乏该领域的经验而找不到。感谢您的任何帮助！ :)
这是我使用的 JupyterLab 笔记本：
import numpy as np
import matplotlib.pyplot as plt
import json
import math
import spiceypy as spice
# Physics
G = 6.6741e-11
class Entity:
def __init__(self,x,v,m,do_gravity):
self.x = x
self.v = v
self.a = np.array([0,0])
self.m = m
self.do_gravity = do_gravity
def apply_force(self,f):
self.a = np.add(self.a,f / self.m);
def step(self,dt):
self.v = np.add(self.v,self.a * dt)
self.x = np.add(self.x,self.v * dt)
self.a = np.array([0,0])
class StaticEntity(Entity):
def __init__(self,body,t,do_gravity):
super().__init__(self.get_state(body,t)[:3],self.get_state(body,t)[3:],self.get_mass(body),do_gravity)
self.body = body
self.t = t
def step(self,dt):
self.t += dt
state = self.get_state(self.body,self.t)
self.x = state[:3]
self.v = state[3:]
@classmethod
def get_state(self,t):
[state,lt] = spice.spkezr(body,"J2000","NONE","SSB")
return state * 1000
@classmethod
def get_mass(self,body):
[dim,gm] = spice.bodvrd(body,"GM",1)
return gm * 1e9 / G
def get_position(self,t):
return self.get_state(self.body,t)[:3]
def get_veLocity(self,t)[3:]
class Propagator:
def __init__(self,entities):
self.entities = entities
def step(self,dt):
for e1 in self.entities:
for e2 in self.entities:
if (e1 is e2) or (not e1.do_gravity) or isinstance(e2,StaticEntity):
continue
diff = np.subtract(e1.x,e2.x)
fg = G * e1.m * e2.m / np.dot(diff,diff)
force = fg * diff / np.linalg.norm(diff)
e2.apply_force(force)
for entity in self.entities:
entity.step(dt)
# LAmbert solver
def C2(psi):
if psi >= 0.0:
sp = math.sqrt(psi)
return (1 - math.cos(sp)) / psi
else:
sp = math.sqrt(-psi)
return (1 - math.cosh(sp)) / psi
def C3(psi):
if psi >= 0.0:
sp = math.sqrt(psi)
return (sp - math.sin(sp)) / (psi * sp)
else:
sp = math.sqrt(-psi)
return (sp - math.sinh(sp)) / (psi * sp)
def lAmbert_solve(r1,r2,tof,mu,iterations,tolerance):
R1 = np.linalg.norm(r1)
R2 = np.linalg.norm(r2)
cos_a = np.dot(r1,r2) / (R1 * R2)
A = math.sqrt(R1 * R2 * (1.0 + cos_a))
sqrt_mu = math.sqrt(mu)
if A == 0.0:
return None
psi = 0.0
psi_lower = -4.0 * math.pi * math.pi
psi_upper = 4.0 * math.pi * math.pi
c2 = 1.0 / 2.0
c3 = 1.0 / 6.0
for i in range(iterations):
B = R1 + R2 + A * (psi * c3 - 1.0) / math.sqrt(c2)
if A > 0.0 and B < 0.0:
psi_lower += math.pi
B = -B
chi = math.sqrt(B / c2)
chi3 = chi * chi * chi
tof_new = (chi3 * c3 + A * math.sqrt(B)) / sqrt_mu
if math.fabs(tof_new - tof) < tolerance:
f = 1.0 - B / R1
g = A * math.sqrt(B / mu)
g_dot = 1.0 - B / R2
v1 = (r2 - f * r1) / g
v2 = (g_dot * r2 - r1) / g
return (v1,v2)
if tof_new <= tof:
psi_lower = psi
else:
psi_upper = psi
psi = (psi_lower + psi_upper) * 0.5
c2 = C2(psi)
c3 = C3(psi)
return None
# Set up solar system
spice.furnsh('solar_system.tm')
inject_time = spice.str2et('2022 Sep 28 00:00:00')
exit_time = spice.str2et('2023 Jun 1 00:00:00')
sun = StaticEntity("Sun",inject_time,True)
earth = StaticEntity("Earth",False)
mars = StaticEntity("Mars Barycenter",False)
(v1,v2) = lAmbert_solve(earth.get_position(inject_time),mars.get_position(exit_time),exit_time - inject_time,G * sun.m,1000,1e-4)
rocket = Entity(earth.x,v1,100000,False)
propagator = Propagator([sun,earth,mars,rocket])
# Generate data
earth_pos = [[],[],[]]
mars_pos = [[],[]]
rocket_pos = [[],[]]
t = inject_time
dt = 3600 # seconds
while t < exit_time:
propagator.step(dt)
earth_pos[0].append(earth.x[0])
earth_pos[1].append(earth.x[1])
earth_pos[2].append(earth.x[2])
mars_pos[0].append(mars.x[0])
mars_pos[1].append(mars.x[1])
mars_pos[2].append(mars.x[2])
rocket_pos[0].append(rocket.x[0])
rocket_pos[1].append(rocket.x[1])
rocket_pos[2].append(rocket.x[2])
t += dt
# Plot data
plt.figure()
plt.title('Transfer orbit')
plt.xlabel('x-axis')
plt.ylabel('y-axis')
plt.plot(earth_pos[0],earth_pos[1],color='blue')
plt.plot(mars_pos[0],mars_pos[1],color='orange')
plt.plot(rocket_pos[0],rocket_pos[1],color='green')
编辑：
我最近重构了我的代码，以便它使用轨道类来表示实体。这实际上给了我可以接受的结果，即使代码在理论上没有做任何不同的事情（据我所知；显然有些事情必须有所不同）
def norm(a):
return np.dot(a,a)**0.5
def fabs(a):
return -a if a < 0 else a
def newton_raphson(f,f_dot,x0,n):
res = x0
for i in range(n):
res = res - f(res) / f_dot(res)
return res
def get_ephemeris(body,time):
state,_ = sp.spkezr(body,time,"SSB")
return np.array(state[:3]) * ap.units.km,np.array(state[3:]) * ap.units.km / ap.units.s
def get_mu(body):
_,mu = sp.bodvrd(body,1)
return mu * ap.units.km**3 / ap.units.s**2
class orbit:
def __init__(self,position,veLocity,mu):
self.position = position
self.veLocity = veLocity
self.mu = mu
@staticmethod
def from_body(name,center,time):
return static_orbit(name,time)
def get_ephemerides(self,dt):
time = 0
positions = []
veLocities = []
#M = self.M
position = self.position
veLocity = self.veLocity
delta_t = dt * ap.units.s
t1 = t * ap.units.s
while time < t1:
g = self.mu / np.dot(position,position)
g_vec = g * -position / norm(position)
veLocity = np.add(veLocity,g_vec * delta_t)
position = np.add(position,veLocity * delta_t)
positions.append(position)
veLocities.append(veLocity)
time = time + delta_t
return positions,veLocities
class static_orbit(orbit):
def __init__(self,name,time):
p,v = get_ephemeris(name,time)
pc,vc = get_ephemeris(center,time)
super().__init__(p - pc,v - vc,get_mu(center))
self.name = name
self.center = center
self.time = time
def get_ephemerides(self,dt):
time = 0
positions = []
veLocities = []
while time < t:
p,v = get_ephemeris(self.name,time + self.time)
pc,vc = get_ephemeris(self.center,time + self.time)
positions.append(p - pc)
veLocities.append(v - vc)
time += dt
return positions,veLocities
sp.furnsh('solar_system.tm')
t1 = sp.str2et('2022 Sep 28 00:00:00')
t2 = sp.str2et('2023 Jun 10 00:00:00')
eo = orbit.from_body("Earth","Sun",t1)
mo = orbit.from_body("Mars Barycenter",t1)
earth_x,earth_v = eo.get_ephemerides(t2 - t1,3600)
mars_x,mars_v = mo.get_ephemerides(t2 - t1,3600)
l = lAmbert(earth_x[0],mars_x[-1],t2 - t1,get_mu("Sun"),1e-6)
ro = orbit(earth_x[0],l.v1,get_mu("Sun"))
rocket_x,rocket_v = ro.get_ephemerides(t2 - t1,60)
earth_x = np.array(earth_x)
mars_x = np.array(mars_x)
rocket_x = np.array(rocket_x)
fig = go.figure()
fig.add_trace(go.Scatter3d(x=earth_x[:,0],y=earth_x[:,1],z=earth_x[:,2],marker_size=1,marker_color='blue'))
fig.add_trace(go.Scatter3d(x=mars_x[:,y=mars_x[:,z=mars_x[:,marker_color='orange'))
fig.add_trace(go.Scatter3d(x=rocket_x[:,y=rocket_x[:,z=rocket_x[:,marker_color='green'))
fig.show()
另外，在再次提到这一点之前，我已经改变了我的积分步长和朗伯求解器容差，但结果却完全不同。
解决方法
所以，经过一番折腾，我设法弄清楚问题出在哪里。我只是没有考虑到太阳不在我的坐标系中的 (0,0) 处。我认为这是微不足道的，但这就是区别所在。最后，我简单地将地球与火星的位置矢量和太阳的位置矢量之间的差值传递给兰伯特求解器。这终于给了我想要的结果。
错误最终变得如此“小”的原因（起初这似乎不是一个明显的错误）是因为我的坐标以距离太阳几百万公里的太阳系重心为中心，正如人们所期望的那样。
感谢您的评论！