1# -*- coding: utf-8 -*-
2"""Test Program for spktype01.
3
4This program computes differences of the result of spktype01 from the
5ephemeris provided by JPL HORIZONS system with three celestial bodies:
6dwarf planet Ceres, and asteroid Vesta.
7
8This program reads two groups of data files.
9First group is SPK data files; Ceres_Vesta.bsp.  this file was
10created by HORIZONS system upon user requests through telnet.
11Second group is ephemeris data files; Ceres_jpl.csv, and Vesta_jpl.csv.
12Each file contains sun centered positions and velocities of respective
13celestial body.  These data are copied from outputs of HORIZONS system.
14
15This program computes 'ranges' for position and velocity, and compare them to
16limits, dposlim and dvellim.  The limits are arbitrary, and physical meanings
17are unclear.
18
19@author: Shushi Uetsuki (whiskie14142)
20"""
21
22from spktype01 import SPKType01
23import csv
24import numpy as np
25
26bspfiles = 'Ceres_Vesta.bsp'
27csvfiles = ['Ceres_jpl.csv', 'Vesta_jpl.csv']
28testnames = ['Ceres', 'Vesta']
29bodyids = [2000001, 2000004]
30
31dposlim = 1.0       # position difference limit = 1.0 kilometer
32dvellim = 1e-7      # velocity difference limit = 0.0001 m/s
33
34kernel = SPKType01.open(bspfiles)
35print(kernel)
36
37for testcase in range(2):
38    csvfile = open(csvfiles[testcase])
39    print()
40    print(testnames[testcase])
41    print()
42
43    err = 0
44    count = 0
45    for row in csv.reader(csvfile):
46        count += 1
47        refjd = float(row[0])
48        refpos = np.array([float(row[2]), float(row[3]), float(row[4])])
49        refvel = np.array([float(row[5]), float(row[6]), float(row[7])])
50
51        spkpos, spkvel = kernel.compute_type01(10, bodyids[testcase], refjd)
52
53        dpos = refpos - spkpos
54        dvel = refvel - spkvel
55        prange = np.sqrt(np.dot(dpos, dpos))
56        vrange = np.sqrt(np.dot(dvel, dvel))
57        if prange > dposlim or vrange > dvellim:
58            print(count, refjd, prange, vrange, ' *')
59            err += 1
60        else:
61            print(count, refjd, prange, vrange)
62    print('Checked count = ', count, ',  Error count = ', err)
63    csvfile.close()
64kernel.close()
65