在下面显示的时间测试中，我发现Skyfield需要几百微秒到一毫秒才能返回obj.at(jd).position.km单个时间值jd，但较长对象（时间点列表）的增量成本JulianDate仅为每点大约 1 微秒. 我使用Jplephem和两种不同的星历表看到了相似的速度。
我的问题是：如果我想随机访问时间点，例如作为使用自己的变量步长的外部 Runge-Kutta 例程的奴隶，有没有办法可以在 python 中更快地做到这一点（无需学习编译代码）？
我知道这根本不是 Skyfield 的典型使用方式。通常我们会加载一个JulianDate包含一长串时间点的对象，然后一次计算它们，并且可能会像轨道积分器那样做几次，而不是数千次（或更多）。
解决方法：我可以想象一种解决方法，我NumPy通过使用具有精细时间粒度的对象运行 Skyfield 一次来构建自己的数据库JulianDate，然后编写我自己的 Runge-Kutta 例程，该例程通过离散量上下改变步长，这样时间步长总是直接对应 NumPy 数组的跨步。
或者我什至可以尝试重新插值。我没有进行高度精确的计算，所以简单的 NumPy 或 SciPy 2 阶可能没问题。
最终我想尝试在太阳系重力场的影响下整合物体的路径（例如深空卫星、彗星、小行星）。在寻找轨道解决方案时，可能会在 6D 相空间中尝试数百万个起始状态向量。我知道我应该使用ob.at(jd).observe(large_body).position.km方法之类的东西，因为重力像其他一切事物一样以光速传播。这似乎花费了大量时间，因为（我猜）它是一个迭代计算（“让我们看看......木星会在哪里让我现在感觉到它的引力”）。但是，让我们一次剥一层宇宙洋葱。
图 1.我的笔记本电脑上针对 de405 和 de421 的不同长度JulianDate对象的 Skyfield 和 JPLephem 性能。它们都差不多——（非常）第一个点大约是半毫秒，每个附加点大约是一微秒。此外，脚本运行时要计算的第一个点len(jd) = 1（地球（蓝色），带有）还有一个额外的毫秒伪影。
地球和月球较慢，因为它是内部的两步计算（地球-月球质心加上围绕质心的各个轨道）。水星可能更慢，因为它与星历时间步长相比移动得如此之快，以至于它在（昂贵的）切比雪夫插值中需要更多的系数？
SKYFIELD DATA 脚本 JPLephem 脚本位于更下方
import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import load, JulianDate
import time
ephem = 'de421.bsp'
ephem = 'de405.bsp'
de = load(ephem)
earth = de['earth']
moon = de['moon']
earth_barycenter = de['earth barycenter']
mercury = de['mercury']
jupiter = de['jupiter barycenter']
pluto = de['pluto barycenter']
things = [ earth, moon, earth_barycenter, mercury, jupiter, pluto ]
names = ['earth', 'moon', 'earth barycenter', 'mercury', 'jupiter', 'pluto']
ntimes = [i*10**n for n in range(5) for i in [1, 2, 5]]
years = [np.zeros(1)] + [np.linspace(0, 100, n) for n in ntimes[1:]] # 100 years
microsecs = []
for y in years:
jd = JulianDate(utc=(1900 + y, 1, 1))
mics = []
for thing in things:
tstart = time.clock()
answer = thing.at(jd).position.km
mics.append(1E+06 * (time.clock() - tstart))
microsecs.append(mics)
microsecs = np.array(microsecs).T
many = [len(y) for y in years]
fig = plt.figure()
ax = plt.subplot(111, xlabel='length of JD object',
ylabel='microseconds',
title='time for thing.at(jd).position.km with ' + ephem )
for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] +
ax.get_xticklabels() + ax.get_yticklabels()):
item.set_fontsize(item.get_fontsize() + 4) # http://stackoverflow.com/a/14971193/3904031
for name, mics in zip(names, microsecs):
ax.plot(many, mics, lw=2, label=name)
plt.legend(loc='upper left', shadow=False, fontsize='x-large')
plt.xscale('log')
plt.yscale('log')
plt.savefig("skyfield speed test " + ephem.split('.')[0])
plt.show()
JPLEPHEM 数据 的脚本 Skyfield 脚本在上面
import numpy as np
import matplotlib.pyplot as plt
from jplephem.spk import SPK
import time
ephem = 'de421.bsp'
ephem = 'de405.bsp'
kernel = SPK.open(ephem)
jd_1900_01_01 = 2415020.5004882407
ntimes = [i*10**n for n in range(5) for i in [1, 2, 5]]
years = [np.zeros(1)] + [np.linspace(0, 100, n) for n in ntimes[1:]] # 100 years
barytup = (0, 3)
earthtup = (3, 399)
# moontup = (3, 301)
microsecs = []
for y in years:
mics = []
#for thing in things:
jd = jd_1900_01_01 + y * 365.25 # roughly, it doesn't matter here
tstart = time.clock()
answer = kernel[earthtup].compute(jd) + kernel[barytup].compute(jd)
mics.append(1E+06 * (time.clock() - tstart))
microsecs.append(mics)
microsecs = np.array(microsecs)
many = [len(y) for y in years]
fig = plt.figure()
ax = plt.subplot(111, xlabel='length of JD object',
ylabel='microseconds',
title='time for jplephem [0,3] and [3,399] with ' + ephem )
# from here: http://stackoverflow.com/a/14971193/3904031
for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] +
ax.get_xticklabels() + ax.get_yticklabels()):
item.set_fontsize(item.get_fontsize() + 4)
#for name, mics in zip(names, microsecs):
ax.plot(many, microsecs, lw=2, label='earth')
plt.legend(loc='upper left', shadow=False, fontsize='x-large')
plt.xscale('log')
plt.yscale('log')
plt.ylim(1E+02, 1E+06)
plt.savefig("jplephem speed test " + ephem.split('.')[0])
plt.show()