129 lines
3.8 KiB
Python
129 lines
3.8 KiB
Python
"""Test harness for checking ``jplephem`` against actual JPL results.
|
|
|
|
This test can be invoked with a simple::
|
|
|
|
python -m jplephem.jpltest
|
|
|
|
"""
|
|
import numpy as np
|
|
from .spk import SPK
|
|
|
|
AU = eval('0.149597870700000000D+09'.replace('D', 'E')) # km
|
|
meter = 0.001 / AU
|
|
epsilon = 0.01 * meter
|
|
field_names = {1: 'x', 2: 'y', 3: 'z', 4: 'vx', 5: 'vy', 6: 'vz'}
|
|
|
|
def run_testpo(spk, testpo_file):
|
|
"""Compare the positions we compute against those computed by the JPL."""
|
|
|
|
lines = iter(testpo_file)
|
|
|
|
while next(lines).strip() != 'EOT':
|
|
continue
|
|
|
|
skips = 0
|
|
successes = 0
|
|
targets = set([segment.target for segment in spk.segments])
|
|
|
|
# Old codes, special-cased in _position():
|
|
targets.add(11)
|
|
targets.add(12)
|
|
targets.add(13)
|
|
|
|
for line in lines:
|
|
de, date, jed, target, center, number, value = [f(v) for f, v
|
|
in zip((str, str, float, int, int, int, float), line.split())]
|
|
|
|
if (target not in targets) or (center not in targets):
|
|
assert target == 14 or target == 15
|
|
assert center == 0
|
|
skips += 1
|
|
continue
|
|
|
|
if 14 <= target <= 15:
|
|
r = _position(spk, jed, target)
|
|
else:
|
|
tpos = _position(spk, jed, target)
|
|
cpos = _position(spk, jed, center)
|
|
r = (tpos - cpos) / AU
|
|
|
|
delta = r[number - 1] - value
|
|
# if (target == 15 and number == 3):
|
|
# delta = delta / (0.23 * (jed - 2451545.0))
|
|
# elif (target == 15 and number == 6):
|
|
# delta = delta * 0.01 / (1.0 + (jed - 2451545.0) / 365.25)
|
|
|
|
if abs(delta) >= epsilon:
|
|
print('%s %s %s->%s field %d (%s)'
|
|
% (date, jed, center, target, number, field_names[number]))
|
|
print(' JPL result: %.15f' % value)
|
|
print(' Our result: %.15f' % r[number - 1])
|
|
print(' ERROR: difference = %s' % (delta,))
|
|
exit(1)
|
|
|
|
successes += 1
|
|
|
|
print(' {0} tests successful, {1} skipped'.format(successes, skips))
|
|
|
|
|
|
def _position(spk, jed, target):
|
|
"""Compute position given a JPL test file target integer identifier."""
|
|
|
|
if target == 3:
|
|
p1, v1 = spk[0,3].compute_and_differentiate(jed)
|
|
p2, v2 = spk[3,399].compute_and_differentiate(jed)
|
|
p = p1 + p2
|
|
v = v1 + v2
|
|
elif target == 10:
|
|
p1, v1 = spk[0,3].compute_and_differentiate(jed)
|
|
p2, v2 = spk[3,301].compute_and_differentiate(jed)
|
|
p = p1 + p2
|
|
v = v1 + v2
|
|
elif target == 11:
|
|
p, v = spk[0,10].compute_and_differentiate(jed)
|
|
elif target == 12:
|
|
return np.zeros((6,)) # solar system barycenter is the origin
|
|
elif target == 13:
|
|
p, v = spk[0,3].compute_and_differentiate(jed)
|
|
else:
|
|
p, v = spk[0,target].compute_and_differentiate(jed)
|
|
|
|
return np.concatenate((p, v))
|
|
|
|
|
|
class MissingFile(Exception):
|
|
pass
|
|
|
|
|
|
def test_all():
|
|
number = 430
|
|
spk_path = 'de%d_test_excerpt.bsp' % number
|
|
testpo_path = 'testpo.%d' % number
|
|
try:
|
|
spk = SPK.open(spk_path)
|
|
except IOError:
|
|
raise MissingFile('cannot open: %s' % spk_path)
|
|
try:
|
|
testpo_file = open(testpo_path)
|
|
except IOError:
|
|
raise MissingFile('cannot open: %s' % testpo_path)
|
|
run_testpo(spk, testpo_file)
|
|
|
|
|
|
if __name__ == '__main__':
|
|
try:
|
|
test_all()
|
|
except MissingFile as e:
|
|
print("""
|
|
Cannot find the JPL SPK files or test files against which this test
|
|
suite validates that the positions it generates are correct. To fetch
|
|
them, run these commands in your current working directory:
|
|
|
|
wget http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de430.bsp
|
|
wget ftp://ssd.jpl.nasa.gov/pub/eph/planets/test-data/430/testpo.430
|
|
|
|
When you are done running the tests, you can remove the files.
|
|
""")
|
|
print(str(e))
|
|
exit(1)
|