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