forked from satellogic/orbit-predictor
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsatellite_nearby.py
executable file
·133 lines (106 loc) · 4.52 KB
/
satellite_nearby.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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
#!/usr/bin/env python3
from glob import glob
import datetime
import numpy as np
import matplotlib.pyplot as plt
from orbit_predictor.sources import EtcTLESource
from orbit_predictor.locations import Location
from orbit_predictor.predictors import Position
BNL_LOC = Location("BNL", latitude_deg=40.878180, longitude_deg=-72.856640, elevation_m=79)
def ll_distance(loc1, loc2):
if hasattr(loc1, "position_llh"):
loc1 = loc1.position_llh
if hasattr(loc2, "position_llh"):
loc2 = loc2.position_llh
return ((loc1[0] - loc2[0])**2 + (loc1[1] - loc2[1])**2)**0.5
def llh_to_altaz(loc1, loc2, radian=False):
'''Return the altaz looking from loc2'''
loc1_llh = loc1.position_llh
loc2_llh = loc2.position_llh
loc1_xyz = loc1.position_ecef
loc2_xyz = loc2.position_ecef
if radian:
coeff = 1
else:
coeff = 180 / np.pi
dx = loc1_llh[1] - loc2_llh[1]
dy = loc1_llh[0] - loc2_llh[0]
az = np.arctan(np.float64(dx) / np.float64(dy)) * coeff
if dy < 0:
az += np.pi * coeff
if az > np.pi * coeff:
az -= 2 * np.pi * coeff
# Earth ellipsoid parameters
a = 6378.1370
b = 6356.752314
# earth radius
n1 = np.sqrt(loc1_xyz[0]**2 + loc1_xyz[1]**2 + loc1_xyz[2]**2)
n2 = np.sqrt(loc2_xyz[0]**2 + loc2_xyz[1]**2 + loc2_xyz[2]**2)
dist = np.sqrt((loc1_xyz[0] - loc2_xyz[0])**2 + (loc1_xyz[1] - loc2_xyz[1])**2 + (loc1_xyz[2] - loc2_xyz[2])**2)\
# cosA = (b^2 + c^2 - a^2) / 2bc
cosalt_center = (n2**2 + dist**2 - n1**2) / (2 * n2 * dist)
#print(cosalt, n1+loc1_llh[2], n2+loc2_llh[2], dist)
alt_center = np.pi - np.arccos(cosalt_center)
lat2_center = np.arctan(loc2_xyz[2] / np.sqrt(loc2_xyz[0]**2 + loc2_xyz[1]**2))
lat2_corr = loc2_llh[0] / 180 * np.pi - lat2_center
alt = (0.5*np.pi - (alt_center + lat2_corr)) * coeff
print(alt, lat2_corr)
return alt, az
def nearest_distance(predictor, loc, start_time=datetime.datetime.now(), delta=datetime.timedelta(seconds=300)):
min_distance = float("inf")
for i in range(86500 // delta.seconds):
curr_time = start_time + delta * i
position = predictor.get_position(curr_time)
distance = ll_distance(position, loc)
if min_distance > distance:
min_distance = distance
min_time = curr_time
return min_distance, min_time
def plot_orbit(ax, predictor, start_time, delta=datetime.timedelta(seconds=60), count=60, polar=False):
x = []
y = []
for i in range(count):
curr_time = start_time + delta * i
position = predictor.get_position(curr_time)
loc = position.position_llh[0:2]
if polar:
#altaz = llh_to_altaz(position, BNL_LOC, radian=True)
az, alt = BNL_LOC.get_azimuth_elev(position)
x.append(az)
y.append((0.5*np.pi - alt) * 180 / np.pi)
else:
x.append(loc[0])
y.append(loc[1])
ax.plot(x, y, label=predictor.sate_id)
sources = EtcTLESource()
for i in glob("gps/*.txt"):
sources.load_file(i)
ax = plt.subplot(111)
print("# Satellite ID, min distance, min time, (satellite lat, lon, height)")
for sate_id in sources.satellite_list():
predictor = sources.get_predictor(sate_id, True)
min_distance, min_time = nearest_distance(predictor, BNL_LOC, start_time=datetime.datetime(2018, 4, 2, 0, 0, 0))
if min_distance < 10:
print(sate_id, ",", min_distance, ",", min_time)
plot_orbit(ax, predictor, min_time - datetime.timedelta(seconds=1800))
plt.plot(BNL_LOC.position_llh[0], BNL_LOC.position_llh[1], marker='o', markersize=5, color='red', label="BMX")
plt.grid()
plt.title("GPS satellites passed over BMX")
plt.xlabel("Latitude, degree")
plt.ylabel("Longitude, degree")
lgd = plt.legend(bbox_to_anchor=(1, 1.2))
plt.savefig("satellite_plot.png", bbox_extra_artists=(lgd,), bbox_inches='tight')
plt.clf()
ax = plt.subplot(111, projection='polar')
ax.set_theta_zero_location("N")
ax.set_theta_direction(-1)
for sate_id in sources.satellite_list():
predictor = sources.get_predictor(sate_id)
min_distance, min_time = nearest_distance(predictor, BNL_LOC)
if min_distance < 10:
#print(sate_id, min_distance, predictor.get_position(min_time).position_llh)
plot_orbit(ax, predictor, min_time - datetime.timedelta(seconds=1800), polar=True)
ax.plot(0, 0, marker='o', markersize=5, color='red', label="BMX")
plt.title("GPS satellites passed over BMX")
lgd = plt.legend(bbox_to_anchor=(1, 1.2))
plt.savefig("satellite_polar.png", bbox_extra_artists=(lgd,), bbox_inches='tight')