Orbit integration with gravitational fields of the solar system from Skyfield - speed issues

Question

In the time tests shown below, I found that Skyfield takes several hundred microseconds to a millisecond to return obj.at(jd).position.kmfor one time value in jd, but the extra cost for longer objects JulianDate(list of points in time) is only one microsecond per point. I see similar speeds using Jplephem and with two different ephemeris.

My question is: if I want to randomly get points in time, for example, as a subordinate to an external Runge-Kutta subroutine that uses its own stepize variable, is there a way I can do this faster in python (without having to learn to compile code)?

I understand that this is not a typical way to use Skyfield. Usually we loaded an object JulianDatewith a long list of time points, and then calculated them immediately and probably did it several times, and not thousands of times (or more), as the orbit integrator can do.

Workaround: . I can imagine where I create my own database NumPyby launching Skyfield once using an object JulianDatewith precise time detail, and then creating my own Runge-Kutta that resizes steps up and down by discrete amounts, so that timestamps always correspond directly to a step NumPy array.

Or I could even repeat the interpolation. I don't do very accurate calculations, so a simple NumPy or SciPy 2nd order can be great.

(, , , ). 6D. , , ob.at(jd).observe(large_body).position.km, , . , , , ( "... , " ). .

1. Skyfield JPLephem JulianDate de405 de421. - () . , , script ( () len(jd) = 1), .

, ( - ). , , () ?

SCRIPT SKYFIELD JPLephem script

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 JPLEPHEM Skyfield script

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