Self-contained test of range_velocity values in pyephem

testpyephem.py

"""
Usage:  python testephem.py [ExportTLEs|exporttles]
"""
import re
import sys
import time
import math
import ephem
import urllib
import datetime
from itertools import izip, count

#Sample TLEs ca. 2014-02-03
lines = [ line.rstrip() for line in """
CELESTIS-02
1 25160U 98007D   14033.23774194  .00000074  00000-0  91275-4 0  3911
2 25160 107.9536 241.0825 0064627 311.1050 144.6915 14.21167072828655
CELESTIS-03
1 26034U 99070C   14033.51828286  .00008891  00000-0  94575-3 0    14
2 26034  98.2209 133.6899 0014527 296.6954  63.2775 14.89134992760942
GOES 2 [-]
1 10061U 77048A   14033.09888698 -.00000285  00000-0  00000+0 0  9402
2 10061  14.4473 340.9297 0012689 287.2582 231.6870  0.99400760 78600
GOES 6 [-]
1 14050U 83041A   14033.40063398 -.00000113  00000-0  00000+0 0  9995
2 14050  14.7558   3.0962 0004166 103.2528  61.8160  1.00327243170765
""".strip().split('\n')]

### Output TLEs for use in another program
if sys.argv[1:] and sys.argv[1].lower()=='exporttles':
  tleOutputFilename = 'goes_other.txt'
  f = open(tleOutputFilename,'wb')
  for line in lines: f.write( '%s\n' % (line,))
  f.close()
  sys.stderr.write("Saved TLEs to file %s\n" % (tleOutputFilename,))

### degrees/radian factor
dpr=180.0/math.pi

### Select and read GOES and CELESTIS satellites
satellites = []
rgx = re.compile('^(GOES|CELESTIS)')
for num, line in izip(count(),lines):
  if rgx.match(line):
    ### create ephem object from TLE information
    satellites.append( ephem.readtle(line, lines[num+1], lines[num+2]) )

### Create Observer under GOES-2
city = ephem.Observer()
city.lon, city.lat, city.elevation = '-32.0' , '-2.0' , 0

### Timestep
deltaSeconds = 10
deltaTime=datetime.timedelta(0,deltaSeconds)

### Run from noon to 4PM on 2014-02-03
t0=datetime.datetime(2014,2,3,12,0,0,500000)
tend=datetime.datetime(2014,2,3,16,0,1)

### List to save ranges and range_velocities for each satellite;
### Used to calculate delta on next pass through loop
rs = [None]*len(satellites)
rvs = rs[:]

### Lists to save delta ranges, and rate percentage errors
drs = [[],[],[],[]]
fwderrs = [[],[],[],[]]
miderrs = [[],[],[],[]]
bckerrs = [[],[],[],[]]

while t0 < tend:

  city.date = t0

  for iSat,satellite in enumerate(satellites):

    ### Save last state for satellite
    lastr,lastrv = rs[iSat],rvs[iSat]

    ### Compute state of satellite for city
    satellite.compute(city)

    ### convert m/s and m to km/s and km, respectively
    rs[iSat],rvs[iSat] = satellite.range/1000., satellite.range_velocity/1000.

    ### Output data
    print ("%s:  RangeRate=%10.3fkm/s; Range=%12.3fkm; Sublat,Sublong=%7.2f,%7.2f Elev.=%10.3fkm; Name=%s; datetime=%s"
           % ( city.date
             , rvs[iSat], rs[iSat]
             , satellite.sublat*dpr, satellite.sublong*dpr, satellite.elevation/1000.
             , satellite.name
             , t0.__repr__()
             ,)
          )

    ####################################################################
    ### Calculate percentage error of rates vs. difference/deltaTime

    ### Do nothing if last Range is None
    if lastr is None: continue

    ### Calculate deltaRange three ways:  by difference; by two rates (before and after)
    dr, fwddr, bckdr = rs[iSat]-lastr, lastrv*deltaSeconds, rvs[iSat]*deltaSeconds

    ### If delta Range is zero, skip this line and the next
    if dr == 0.0:
      rs[iSat] = None
      continue

    ### Percentage error function
    pct = lambda num,den : 100.0 * (num-den) / den

    ### Append delta range to appropriate list
    drs[iSat].append(dr)

    ### Append percentage errors to appropriate list
    fwderrs[iSat].append(pct(fwddr,dr))
    bckerrs[iSat].append(pct(bckdr,dr))

    ### Use mean of before and after for midpoint rate-based delta Range
    miderrs[iSat].append(pct((fwddr+bckdr) / 2.,dr))

  ### Increment time and print a blank line
  t0 += deltaTime
  print( '' )


### Plot error results
for iSat,satellite in enumerate(satellites):
  import matplotlib.pyplot as pl
  pl.plot( fwderrs[iSat], label='Fwd ' + len(fwderrs[iSat]).__str__() )
  pl.plot( bckerrs[iSat], label='Bck ' + len(bckerrs[iSat]).__str__() )
  pl.plot( miderrs[iSat], label='Mid ' + len(miderrs[iSat]).__str__() )
  pl.legend()
  pl.title(satellite.name)
  pl.xlabel('obs @ 0.1Hz')
  pl.ylabel('Error, percent')
  pl.show()