Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
193 lines (162 sloc) 9.09 KB
# -*- coding: utf-8 -*-
import numpy as np
from datetime import datetime, timedelta
import pytz
import graphics
import urllib
pdt = pytz.timezone('US/Pacific')
sqrt = np.sqrt
pi = np.pi
sin = np.sin
cos = np.cos
# Standard Gravitational parameter in km^3 / s^2 of Earth
GM = 398600.4418
def splitElem(tle):
"Splits a two line element set into title and it's two lines with stripped lines"
return map(lambda x: x.strip(), tle.split('\n'))
def checkValid(tle):
"Checks with checksum to make sure element is valid"
title, line1, line2 = splitElem(tle)
return line1[0] == '1' and line2[0] == '2' and \
line1[2:7] == line2[2:7] and \
int(line1[-1]) == doChecksum(line1) and int(line2[-1]) == doChecksum(line2)
def stringScientificNotationToFloat(sn):
"Specific format is 5 digits, a + or -, and 1 digit, ex: 01234-5 which is 0.01234e-5"
return 0.00001*float(sn[5]) * 10**int(sn[6:])
def eccentricAnomalyFromMean(mean_anomaly, eccentricity, initValue,
maxIter=500, maxAccuracy = 0.0001):
"""Approximates Eccentric Anomaly from Mean Anomaly
All input and outputs are in radians"""
mean_anomaly = mean_anomaly
e0 = initValue
for x in xrange(maxIter):
e1 = e0 - (e0 - eccentricity * sin(e0) - mean_anomaly) / (1.0 - eccentricity * cos(e0))
if (abs(e1-e0) < maxAccuracy):
break
return e1
def pretty_print(tle, printInfo = True, labels = True):
"Returns commented information on a two line element"
title, line1, line2 = splitElem(tle)
if not checkValid(tle):
print "Invalid element."
return
satellite_number = int(line1[2:7])
classification = line1[7:8]
international_designator_year = int(line1[9:11])
international_designator_launch_number = int(line1[11:14])
international_designator_piece_of_launch = line1[14:17]
epoch_year = int(line1[18:20])
epoch = float(line1[20:32])
first_time_derivative_of_the_mean_motion_divided_by_two = float(line1[33:43])
second_time_derivative_of_mean_motion_divided_by_six = stringScientificNotationToFloat(line1[44:52])
bstar_drag_term = stringScientificNotationToFloat(line1[53:61])
the_number_0 = float(line1[62:63])
element_number = float(line1[64:68])
checksum1 = float(line1[68:69])
satellite = int(line2[2:7])
inclination = float(line2[8:16])
right_ascension = float(line2[17:25])
eccentricity = float(line2[26:33]) * 0.0000001
argument_perigee = float(line2[34:42])
mean_anomaly = float(line2[43:51])
mean_motion = float(line2[52:63])
revolution = float(line2[63:68])
checksum2 = float(line2[68:69])
# Inferred Epoch date
year = 2000 + epoch_year if epoch_year < 70 else 1900 + epoch_year
epoch_date = datetime(year=year, month=1, day=1, tzinfo=pytz.utc) + timedelta(days=epoch-1) # Have to subtract one day to get correct midnight
# Time difference of now from epoch, offset in radians
diff = datetime.now().replace(tzinfo=pytz.utc) + timedelta(hours=8) - epoch_date # Offset for PDT
diff_seconds = 24*60*60*diff.days + diff.seconds + 1e-6*diff.microseconds # sec
print "Time offset: %s" % diff
motion_per_sec = mean_motion * 2*pi / (24*60*60) # rad/sec
print "Radians per second: %g" % motion_per_sec
offset = diff_seconds * motion_per_sec #rad
print "Offset to apply: %g" % offset
mean_anomaly += offset * 180/pi % 360
# Inferred period
day_seconds = 24*60*60
period = day_seconds * 1./mean_motion
# Inferred semi-major axis (in km)
semi_major_axis = ((period/(2*pi))**2 * GM)**(1./3)
# Inferred true anomaly
eccentric_anomaly = eccentricAnomalyFromMean(mean_anomaly * pi/180, eccentricity, mean_anomaly * pi/180)
true_anomaly = 2*np.arctan2(sqrt(1+eccentricity) * sin(eccentric_anomaly/2.0), sqrt(1-eccentricity) * cos(eccentric_anomaly/2.0))
# Convert to degrees
eccentric_anomaly *= 180/pi
true_anomaly *= 180/pi
if (printInfo):
print "----------------------------------------------------------------------------------------"
print tle
print "---"
print "Satellite Name = %s" % title
print "Satellite number = %g (%s)" % (satellite_number, "Unclassified" if classification == 'U' else "Classified")
print "International Designator = YR: %02d, LAUNCH #%d, PIECE: %s" % (international_designator_year, international_designator_launch_number, international_designator_piece_of_launch)
print "Epoch Date = %s (YR:%02d DAY:%.11g)" % (epoch_date.strftime("%Y-%m-%d %H:%M:%S.%f %Z"), epoch_year, epoch)
print "First Time Derivative of the Mean Motion divided by two = %g" % first_time_derivative_of_the_mean_motion_divided_by_two
print "Second Time Derivative of Mean Motion divided by six = %g" % second_time_derivative_of_mean_motion_divided_by_six
print "BSTAR drag term = %g" % bstar_drag_term
print "The number 0 = %g" % the_number_0
print "Element number = %g" % element_number
print
print "Inclination [Degrees] = %g°" % inclination
print "Right Ascension of the Ascending Node [Degrees] = %g°" % right_ascension
print "Eccentricity = %g" % eccentricity
print "Argument of Perigee [Degrees] = %g°" % argument_perigee
print "Mean Anomaly [Degrees] Anomaly = %g°" % mean_anomaly
print "Eccentric Anomaly = %g°" % eccentric_anomaly
print "True Anomaly = %g°" % true_anomaly
print "Mean Motion [Revs per day] Motion = %g" % mean_motion
print "Period = %s" % timedelta(seconds=period)
print "Revolution number at epoch [Revs] = %g" % revolution
print
print "semi_major_axis = %gkm" % semi_major_axis
print "eccentricity = %g" % eccentricity
print "inclination = %g°" % inclination
print "arg_perigee = %g°" % argument_perigee
print "right_ascension = %g°" % right_ascension
print "true_anomaly = %g°" % true_anomaly
print "----------------------------------------------------------------------------------------"
if labels:
graphics.plotOrbit(semi_major_axis, eccentricity, inclination,
right_ascension, argument_perigee, true_anomaly, title)
else:
graphics.plotOrbit(semi_major_axis, eccentricity, inclination,
right_ascension, argument_perigee, true_anomaly)
def doChecksum(line):
"""The checksums for each line are calculated by adding the all numerical digits on that line, including the
line number. One is added to the checksum for each negative sign (-) on that line. All other non-digit
characters are ignored.
@note this excludes last char for the checksum thats already there."""
return sum(map(int, filter(lambda c: c >= '0' and c <= '9', line[:-1].replace('-','1')))) % 10
# elem1 = """ISS (ZARYA)
# 1 25544U 98067A 08264.51782528 -.00002182 00000-0 -11606-4 0 2927
# 2 25544 51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537"""
# dragon = """DRAGON CRS-2
# 1 39115U 13010A 13062.62492353 .00008823 00000-0 14845-3 0 188
# 2 39115 51.6441 272.5899 0012056 334.2535 68.5574 15.52501943 306"""
graphics.plotEarth()
# Data from NORAD http://www.celestrak.com/NORAD/elements/
# filename = "noaa.txt" # NOAA satellites
# filename = "geo.txt" # Geostationary satellites
# filename = "gps-ops.txt" # GPS sats
# filename = "military.txt" # Some military satellites
# filename = "stations.txt" # Space stations
# filename = "visual.txt" # 100 brightest or so objects
# files = ["noaa.txt", "stations.txt", "military.txt", "gps-ops.txt"]
# files = ["stations.txt"]
names = ["stations"]
for urlname in names:
f = urllib.urlopen("http://www.celestrak.com/NORAD/elements/%s.txt" % urlname)
elem = ""
for line in f:
elem += line
if (line[0] == '2'):
elem = elem.strip()
if elem.startswith("ISS"):
pretty_print(elem, printInfo=True, labels=True)
elem = ""
# pretty_print(elem1)
# pretty_print(dragon)
graphics.doDraw()
# EOF