In [1]:
%matplotlib inline
import sgp4
import sgp4.io
import sgp4.earth_gravity
import urllib.request
import http.cookiejar
import numpy as np
import matplotlib.pyplot as plt
# Larger figure size
fig_size = [14, 10]
plt.rcParams['figure.figsize'] = fig_size
In [2]:
def opener_spacetrack(username, password):
# todo: clean login info from ipynb file
cj = http.cookiejar.CookieJar()
opener = urllib.request.build_opener(urllib.request.HTTPCookieProcessor(cj))
auth_url = 'https://www.space-track.org/ajaxauth/login/'
auth_data = urllib.parse.urlencode({'identity' : username, 'password' : password}).encode('utf-8')
auth_req = urllib.request.Request(auth_url, auth_data)
r = opener.open(auth_req)
return opener
In [3]:
def get_tles(norad_id, opener):
url = 'https://www.space-track.org/basicspacedata/query/class/tle/EPOCH/%3E2017-01-01/NORAD_CAT_ID/{}/orderby/EPOCH%20ASC/format/tle'.format(norad_id)
r = opener.open(url)
tle_lines = r.read().decode('ascii').split('\r\n')
return [sgp4.io.twoline2rv(*x, whichconst=sgp4.earth_gravity.wgs72) for x in zip(tle_lines[::2], tle_lines[1::2])]
In [4]:
seconds_in_day = 24*3600
def prediction_error(tle1, tle2):
x1, x2 = (t.propagate(year = tle2.epochyr, day = tle2.epochdays)[0] for t in (tle1, tle2))
return 1e3*np.sqrt(sum((x1[j] - x2[j])**2 for j in range(3)))
def time_difference(tle1, tle2): # in days
return (tle2.jdsatepoch - tle1.jdsatepoch)
def tle_variation(tle1, tle2):
time = time_difference(tle1, tle2)*seconds_in_day
if time == 0:
return np.nan
return prediction_error(tle1, tle2)/time
def compute_variation(tles):
variations = np.array([tle_variation(*x) for x in zip(tles[:-1], tles[1:])])
times = np.array([t.epochdays for t in tles[1:]])
filt = ~np.isnan(variations)
return (times[filt], variations[filt])
def compute_errors(tles):
l = len(tles)
times = [time_difference(tles[j], tles[k]) for j in range(l) for k in range(j+1,l)]
errors = [prediction_error(tles[j], tles[k]) for j in range(l) for k in range(j+1,l)]
return (times, errors)
def plot_variation(var, **kwargs):
return plt.semilogy(var[0], var[1], **kwargs)
def plot_errors(err, **kwargs):
return plt.loglog(err[0], err[1], ',', **kwargs)
In [5]:
with open('spacetrack_auth', 'r') as f:
username, password = f.read().split('\n')[:2]
opener = opener_spacetrack(username, password)
In [6]:
satellites = { 7530 : 'AO-7', 24278 : 'FO-29', 27607 : 'SO-50', 30776 : 'FALCONSAT-3', 39444 : 'AO-73', 40074 : 'UKUBE-1', 40903 : 'XW-2A', 40911 : 'XW-2B', 40906 : 'XW-2C', 40907 : 'XW-2D', 40910 : 'XW-2F', 40908 : 'LilacSat-2', 42017 : 'NAYIF-1', 42725 : 'LilacSat-1', 42761 : 'CAS-4A', 42759 : 'CAS-4B'}
In [7]:
tles = {k : get_tles(k, opener) for k in satellites}
In [8]:
variations = {k : compute_variation(tles[k]) for k in satellites}
In [9]:
errors = {k : compute_errors(tles[k]) for k in satellites}
In [10]:
def plot_graph(ids, **kwargs):
[plot_variation(variations[k], **kwargs) for k in ids]
plt.legend([satellites[k] for k in ids])
plt.ylabel('TLE variation (m/s)')
plt.xlabel('Time since 2017-01-01 (days)')
def plot_graph_errs(ids):
[plot_errors(errors[k], alpha = 0.1) for k in ids]
leg = plt.legend([satellites[k] for k in ids])
for l in leg.get_lines():
l.set_alpha(1)
l.set_marker('.')
plt.ylabel('Error in prediction (m)')
plt.xlabel('TLE age (days)')
In [11]:
xws = [k for (k,v) in satellites.items() if v.startswith('XW-2')]
plot_graph(xws, alpha = 0.4)
plt.title('TLE variation: XW-2 satellites');
In [12]:
plot_graph([39444, 40074, 42017], alpha = 0.7)
plt.title('TLE variation: FUNcube satellites');
In [13]:
plot_graph([7530, 24278, 27607, 30776], alpha = 0.7)
plt.title('TLE variation: older satellites');
In [14]:
lilac = [k for (k,v) in satellites.items() if v.startswith('LilacSat-')]
plot_graph(lilac, alpha = 0.7)
plt.title('TLE variation: LilacSat\'s');
In [15]:
cas = [k for (k,v) in satellites.items() if v.startswith('CAS-4')]
plot_graph(cas, alpha = 0.7)
plt.title('TLE variation: CAS-4 satellites');
In [16]:
leo_orbital_vel = 7.8e3