mirror of
https://github.com/cosinekitty/astronomy.git
synced 2026-01-01 20:20:15 -05:00
98 lines
3.2 KiB
Python
Executable File
98 lines
3.2 KiB
Python
Executable File
#!/usr/bin/env python3
|
|
import sys
|
|
import re
|
|
from math import sqrt, degrees
|
|
|
|
class Record:
|
|
def __init__(self, planet, jd, x, y, z, vx, vy, vz):
|
|
self.planet = planet
|
|
self.jd = jd
|
|
self.x = x
|
|
self.y = y
|
|
self.z = z
|
|
self.vx = vx
|
|
self.vy = vy
|
|
self.vz = vz
|
|
|
|
def JplRead(fn, planet):
|
|
with open(fn) as jplfile:
|
|
seek = True
|
|
jlist = []
|
|
for line in jplfile:
|
|
line = line.strip()
|
|
if seek:
|
|
if line == '$$SOE':
|
|
seek = False
|
|
row = 0
|
|
else:
|
|
if line == '$$EOE':
|
|
break
|
|
if row == 0:
|
|
# 2411545.000000000 = A.D. 1890-Jun-26 12:00:00.0000 TDB
|
|
jd = float(line.split()[0])
|
|
elif row == 1:
|
|
# X = 1.213238742051184E+01 Y = 2.532261522896372E+01 Z = 1.006281191830243E+01
|
|
m = re.match(r'^\s*X\s*=\s*(\S+)\s+Y\s*=\s*(\S+)\s+Z\s*=\s*(\S+)', line)
|
|
if not m:
|
|
raise Exception('Bad data format in [{}]'.format(line))
|
|
x = float(m.group(1))
|
|
y = float(m.group(2))
|
|
z = float(m.group(3))
|
|
elif row == 2:
|
|
# VX=-2.880310222354256E-03 VY= 1.177013903016748E-03 VZ= 5.534782388307750E-04
|
|
m = re.match(r'^\s*VX\s*=\s*(\S+)\s+VY\s*=\s*(\S+)\s+VZ\s*=\s*(\S+)', line)
|
|
if not m:
|
|
raise Exception('Bad data format in [{}]'.format(line))
|
|
vx = float(m.group(1))
|
|
vy = float(m.group(2))
|
|
vz = float(m.group(3))
|
|
else:
|
|
jlist.append(Record(planet, jd, x, y, z, vx, vy, vz))
|
|
row = (row + 1) % 4
|
|
return jlist
|
|
|
|
|
|
def TopRead(fn):
|
|
with open(fn) as topfile:
|
|
row = 0
|
|
tlist = []
|
|
for line in topfile:
|
|
line = line.strip()
|
|
token = line.split()
|
|
if row == 0:
|
|
planet = int(token[0])
|
|
jd = float(token[1])
|
|
elif row == 3:
|
|
x, y, z, vx, vy, vz = [float(t) for t in token]
|
|
tlist.append(Record(planet, jd, x, y, z, vx, vy,vz))
|
|
row = (row + 1) % 4
|
|
return tlist
|
|
|
|
|
|
def Compare(t, j):
|
|
if t.jd != j.jd:
|
|
raise Exception('Mismatching times')
|
|
if t.planet != j.planet:
|
|
raise Exception('Mismatching planets')
|
|
range = sqrt(j.x**2 + j.y**2 + j.z**2)
|
|
diff = sqrt((t.x-j.x)**2 + (t.y-j.y)**2 + (t.z-j.z)**2)
|
|
arcmin = 60 * degrees(diff / range)
|
|
print('planet {:1d} jd {:10.1f} arcmin {:10.6f}'.format(t.planet, t.jd, arcmin))
|
|
if arcmin > 0.09:
|
|
raise Exception('EXCESSIVE ARCMIN ERROR')
|
|
|
|
|
|
def main():
|
|
tlist = TopRead('correct.txt')
|
|
print('TOP2013 record count =', len(tlist))
|
|
jlist = []
|
|
for planet in range(5, 10):
|
|
jlist += JplRead('jplhor_{:d}.txt'.format(planet), planet)
|
|
print('JPLHOR record count =', len(jlist))
|
|
for (t, j) in zip(tlist, jlist):
|
|
Compare(t, j)
|
|
return 0
|
|
|
|
if __name__ == '__main__':
|
|
sys.exit(main())
|