-
Notifications
You must be signed in to change notification settings - Fork 0
/
main.py
70 lines (52 loc) · 2 KB
/
main.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#!/usr/bin/python3
import sys
import numpy as np
from earthposition import EarthPosition
from icrsposition import IcrsPosition
from observations import Observations
from orbitalelements import OrbitalElements
if len(sys.argv) != 5:
print("Simple orbit determination program.")
print("Usage:")
print(" python3 main.py [obs_file] i j k")
print("")
print("where [obs_file] is an Observer Table ephemeris generated by ")
print("NASA JPL HORIZONS system containing minimum three observer records ")
print("with: UT date, RA(ICRF), DEC(ICRF) each,")
print("and i j and k are three observation numbers from [obs_file]")
print("which will be used for orbit determination.")
print()
print("The program has built in Earth ephemeris for dates 01.01.2018 - 01.01.2021.")
print("Only dates after 31.12.2016 are properly calculated (see: timestamp.py line 42).")
quit(0)
observations_fname = sys.argv[1]
i = int(sys.argv[2])
j = int(sys.argv[3])
k = int(sys.argv[4])
ijk = (i, j, k)
earthPosition = EarthPosition("earth.csv")
observations = Observations(observations_fname)
icrsposition = IcrsPosition(earthPosition, observations)
phi = icrsposition.phi_angle(ijk)
positions = icrsposition.positions(ijk)
radtodeg = 360/(2*np.pi)
i = 0
for p in positions:
i += 1
print(f'Solution: {i}:')
orbitalElements = OrbitalElements(p)
elem = orbitalElements.elements()
inclination = elem.inclination * radtodeg
if inclination < 0: inclination += 180
Omega = elem.Omega * radtodeg
if Omega < 0: Omega += 360
omega = elem.omega * radtodeg
if omega < 0: omega += 360
print(f'i | Inclination: {inclination} deg')
print(f'Omega | Longitude of ascending node: {Omega} deg')
print(f'omega | Argument of periapsis: {omega} deg')
print(f'a | Semi-major axis: {elem.a} au')
print(f'e | Eccentricity: {elem.e}')
print(f'tp | Time of periapsis: {elem.tp} JDTT')
print()
pass