在下面所示的时间测试中，我发现Skyfield http://rhodesmill.org/skyfield/需要几百微秒到一毫秒才能返回obj.at(jd).position.km对于单个时间值jd，但较长时间的增量成本JulianDate对象（时间点列表）每个点大约只有一微秒。我看到类似的速度使用Jplephem https://pypi.python.org/pypi/jplephem并有两个不同的星历表。
我的问题是：如果我想随机访问时间点，例如作为使用其自己的变量步长的外部 Runge-Kutta 例程的从属，有没有一种方法可以在 python 中更快地完成此操作（无需学习编译代码） ？
我明白这根本不是 Skyfield 的典型使用方式。通常我们会加载一个JulianDate对象具有一长串时间点，然后立即计算它们，并且可能会这样做几次，而不是像轨道积分器那样数千次（或更多）。
解决方法：我可以想象一种解决方法，我可以自己构建NumPy通过运行 Skyfield 一次使用JulianDate具有精细时间粒度的对象，然后编写我自己的 Runge-Kutta 例程，该例程以离散量上下更改步长，以便时间步始终直接对应于 NumPy 数组的跨步。
或者我什至可以尝试重新插值。我没有进行高精度计算，因此简单的 NumPy 或 SciPy 二阶可能就可以了。
最终我想尝试整合太阳系重力场影响下的物体（例如深空卫星、彗星、小行星）的路径。在寻找轨道解时，人们可能会在 6D 相空间中尝试数百万个起始状态向量。我知道我应该使用类似的东西ob.at(jd).observe(large_body).position.km方法，因为引力像其他东西一样以光速传播。这似乎花费了大量时间，因为（我猜）这是一个迭代计算（“让我们看看……木星会在哪里，让我现在感觉到它的引力”）。但让我们一次剥一层宇宙洋葱。
图1。Skyfield 和 JPLephem 在我的笔记本电脑上不同长度的性能JulianDate对象，适用于 de405 和 de421。它们都大致相同 - 第一个点（非常）大约为半毫秒，每个附加点为一微秒。还有第一点脚本运行时计算（地球（蓝色）与len(jd) = 1）有一个额外的毫秒伪影。
地球和月球速度较慢，因为它内部是两步计算（地月重心加上绕重心的各个轨道）。水星可能会更慢，因为与星历时间步长相比，它移动得如此之快，以至于在（昂贵的）切比雪夫插值中需要更多系数？
天场数据脚本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()
