performance - スカイフィールドからの太陽系重力場と軌道の統合 - 速度の問題

Question

以下に示す時間テストでは、Skyfieldobj.at(jd).position.kmがの単一の時間値を返すのに数百マイクロ秒から 1 ミリ秒かかることがわかりましたjdが、より長いオブジェクト (時間内のポイントのリスト) の増分コストはJulianDate、ポイントごとに約 1 マイクロ秒に過ぎません。 . Jplephemと 2 つの異なるエフェメリスを使用した場合、同様の速度が見られます。

ここでの私の質問は次のとおりです。たとえば、独自の可変ステップサイズを使用する外部ルンゲクッタルーチンのスレーブとして、特定の時点にランダムアクセスしたい場合、Python内でこれをより高速に実行できる方法はありますか(コードのコンパイル方法を学ぶ)?

これは、Skyfield の一般的な使用方法ではないことを理解しています。通常、JulianDateオブジェクトを長いリストの時点でロードし、それらを一度に計算します。軌道積分器が行う方法のように、おそらく数千回 (またはそれ以上) ではなく、数回実行します。

回避策:時間粒度の細かいオブジェクトをNumPy使用して Skyfield を 1 回実行し、タイムステップが常にJulianDateNumPy 配列のストライディングに直接対応します。

または、再補間を試すこともできます。私は非常に正確な計算を行っていないので、単純な NumPy または SciPy の 2 次で問題ないかもしれません。

最終的には、太陽系の重力場の影響下にあるオブジェクト (深宇宙衛星、彗星、小惑星など) の経路を統合してみたいと思います。軌道解を探すとき、6D 位相空間で何百万もの開始状態ベクトルを試すかもしれません。ob.at(jd).observe(large_body).position.km重力は他のすべてのものと同じように光速で移動するため、メソッドなどを使用する必要があることはわかっています。これは反復計算であるため（推測では）かなりの時間がかかるようです（「見てみましょう...木星はどこにあったので、今すぐ重力だと感じます」）。しかし、コズミックオニオンを一度に1層ずつ剥がしましょう.

図 1. de405 と de421 のさまざまな長さJulianDateのオブジェクトに対する私のラップトップでの Skyfield と JPLephem のパフォーマンス。それらはすべてほぼ同じです-（非常に）最初のポイントで約0.5ミリ秒、追加のポイントごとに1マイクロ秒です。また、スクリプトの実行時に計算される最初のポイント(Earth (blue) with len(jd) = 1) には、追加のミリ秒のアーティファクトがあります。

地球と月は、内部で 2 段階の計算 (地球と月の重心と重心に関する個々の軌道) であるため、処理が遅くなります。水星は、エフェメリスの時間ステップと比較して非常に速く移動するため、（コストのかかる）チェビシェフ補間でより多くの係数が必要になるため、遅くなる可能性がありますか?

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

SCRIPT FOR JPLEPHEM DATA 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()