Wrong Range Rate with Pyephem

Question

I am trying to calculate a satellites Range Rate using Python and pyephem. Unfortunately pyephems result seems to be wrong.

After comparing the value with calculations made by other satellite tracking programs like GPredict or Ham Radio Deluxe the the difference goes up to 2km/sec.The calculated values for the Azemuth and Elevation ankle are almost the same thought. TLE's are new and the system clock is the same.

Do you see any mistake I made in my code or do you have an idea what else could cause the error?

Thank you very much!

Here is my Code:

import ephem
import time
#TLE Kepler elements
line1 = "ESTCUBE 1"       
line2 = "1 39161U 13021C   13255.21187718  .00000558  00000-0  10331-3 0  3586"
line3 = "2 39161  98.1264 332.9982 0009258 190.0328 170.0700 14.69100578 18774"
satellite = ephem.readtle(line1, line2, line3) # create ephem object from tle information
while True:
    city = ephem.Observer() # recreate Oberserver with current time
    city.lon, city.lat, city.elevation = '52.5186' , '13.4080' , 100

    satellite.compute(city)
    RangeRate = satellite.range_velocity/1000 # get RangeRate in km/sec
    print ("RangeRate: " + str(RangeRate))
    time.sleep(1)

I recorded some Range Rate values from the script and from GPRedict to make the error reproducibly:

ESTCUBE 1               
1 39161U 13021C   13255.96108453  .00000546  00000-0  10138-3 0  3602
2 39161  98.1264 333.7428 0009246 187.4393 172.6674 14.69101320 18883

date: 2013-09-13  
time       pyephem-Script  Gpredict          
14:07:02   -1.636          -3.204  
14:12:59   -2.154          -4.355  
14:15:15   -2.277          -4.747  
14:18:48   -2.368          -5.291  

And I added some lines to calculate the satellites elevation and coordinates:

elevation = satellite.elevation
sat_latitude = satellite.sublat
sat_longitude = satellite.sublong

The results with time stamp are:

2013-09-13 14:58:13  
RangeRate: 2.15717797852 km/s  
Range: 9199834.0  
Sat Elevation: 660743.6875  
Sat_Latitude: -2:22:27.3  
Sat_Longitude: -33:15:15.4    

2013-09-13 14:58:14  
RangeRate: 2.15695092773 km/s  
Range: 9202106.0  
Sat Elevation: 660750.9375  
Sat_Latitude: -2:26:05.8  
Sat_Longitude: -33:16:01.7  

Another important information might be that I am trying to calculate the Doppler Frequency for a satellite pass. And therefore I need the Range Rate:

f_Doppler_corrected = (c0/(c0 + RangeRate))*f0

Range Rate describes the velocity of a moving object on the visual axis to the observer. Maybe the range_velocity is something different?

Answer 1

It seems pyephem (libastro as a backend) and gpredict (predict) as a backend use different ways to calculate the satellite velocity. I am attaching detailed output of comparison for an actual reference observation. It can be seen that both output the correct position, while only gpredict outputs reasonable range_rate values. The error seems to occur in the satellite velocity vector. I would say that the reasons from gpredict are more reasonable (and the similar code is with question marks in libastro ..) therefore I will propose a fix in libastro to handle it as in gpredict, however maybe someone who understands the math behind it can add to this.

I added another tool, PyPredict (also predict based), to get some calculations here. However these values are off, so must be something else.

Pyephem: 3.7.5.3
Gpredict: 1.3
PyPredict 1.1 (Git: 10/02/2015)


OS: Ubuntu x64
Python 2.7.6

Time:
Epoch timestamp: 1420086600
Timestamp in milliseconds: 1420086600000
Human time (GMT): Thu, 01 Jan 2015 04:30:00 GMT

ISS (ZARYA)             
1 25544U 98067A   15096.52834639  .00016216  00000-0  24016-3 0  9993
2 25544  51.6469  82.0200 0006014 185.1879 274.8446 15.55408008936880

observation point: N0 E0 alt=0

Test 1:

Gpredict: (Time, Az, El, Slant Range, Range Velocity)

2015 01 01 04:30:00 202.31 -21.46 5638 -5.646
2015 01 01 04:40:00 157.31 -2.35 2618 -3.107
2015 01 01 04:50:00 72.68 -10.26 3731 5.262

Pyephem 3.7.5.3 (default atmospheric refraction)

(2015/1/1 04:30:00, 202:18:45.3, -21:27:43.0, 5638.0685, -5.3014228515625)
(2015/1/1 04:40:00, 157:19:08.3, -1:21:28.6, 2617.9915, -2.934402099609375)
(2015/1/1 04:50:00, 72:40:59.9, -10:15:15.1, 3730.78375, 4.92381201171875)
No atmospheric refraction
(2015/1/1 04:30:00, 202:18:45.3, -21:27:43.0, 5638.0685, -5.3014228515625)
(2015/1/1 04:40:00, 157:19:08.3, -1:21:28.6, 2617.9915, -2.934402099609375)
(2015/1/1 04:50:00, 72:40:59.9, -10:15:15.1, 3730.78375, 4.92381201171875)

Pypredict

1420086600.0
{'decayed': 0, 'elevation': -19.608647085869123, 'name': 'ISS (ZARYA)', 'norad_id': 25544, 'altitude': 426.45804846615556, 'orbit': 92208, 'longitude': 335.2203454719759, 'sunlit': 1, 'geostationary': 0, 'footprint': 4540.173580837984, 'epoch': 1420086600.0, 'doppler': 1635.3621339278857, 'visibility': 'D', 'azimuth': 194.02436209048014, 'latitude': -45.784314563471646, 'orbital_model': 'SGP4', 'orbital_phase': 73.46488929141783, 'eclipse_depth': -8.890253049060693, 'slant_range': 5311.3721164183535, 'has_aos': 1, 'orbital_velocity': 27556.552465256085}
1420087200.0
{'decayed': 0, 'elevation': -6.757496200551716, 'name': 'ISS (ZARYA)', 'norad_id': 25544, 'altitude': 419.11153234752874, 'orbit': 92208, 'longitude': 9.137628905963876, 'sunlit': 1, 'geostationary': 0, 'footprint': 4502.939901708917, 'epoch': 1420087200.0, 'doppler': 270.6901377419433, 'visibility': 'D', 'azimuth': 139.21315598291235, 'latitude': -20.925997669236732, 'orbital_model': 'SGP4', 'orbital_phase': 101.06301876416072, 'eclipse_depth': -18.410968838249545, 'slant_range': 3209.8444916123644, 'has_aos': 1, 'orbital_velocity': 27568.150821416708}
1420087800.0
{'decayed': 0, 'elevation': -16.546383900323555, 'name': 'ISS (ZARYA)', 'norad_id': 25544, 'altitude': 414.1342802649042, 'orbit': 92208, 'longitude': 31.52356804788407, 'sunlit': 1, 'geostationary': 0, 'footprint': 4477.499436144489, 'epoch': 1420087800.0000002, 'doppler': -1597.032808834609, 'visibility': 'D', 'azimuth': 76.1840387294104, 'latitude': 9.316828913183791, 'orbital_model': 'SGP4', 'orbital_phase': 128.66115193399546, 'eclipse_depth': -28.67721196244149, 'slant_range': 4773.838774518728, 'has_aos': 1, 'orbital_velocity': 27583.591664378775}


Test 2 (short time):
Gpredict: (Slant Range, Range Velocity)
2015 01 01 04:30:00 5638 -5.646
2015 01 01 04:30:10 5581 -5.648
->5.7 km/s avg

(2015/1/1 04:30:00, 5638.0685, -5.3014228515625)
(2015/1/1 04:30:10, 5581.596, -5.30395361328125)
->5.7 km/s avg

Pyephem

import ephem
import time
#TLE Kepler elements

line1 = "ISS (ZARYA)"             
line2 = "1 25544U 98067A   15096.52834639  .00016216  00000-0  24016-3 0  9993"
line3 = "2 25544  51.6469  82.0200 0006014 185.1879 274.8446 15.55408008936880"

satellite = ephem.readtle(line1, line2, line3) # create ephem object from tle information

obs = ephem.Observer() # recreate Oberserver with current time
obs.lon, obs.lat, obs.elevation = '0' , '0' , 0

print('Pyephem Default (atmospheric refraction)')

obs.date = '2015/1/1 04:30:00'
satellite.compute(obs)
print(obs.date, satellite.az, satellite.alt,satellite.range/1000, satellite.range_velocity/1000)

obs.date = '2015/1/1 04:40:00'
satellite.compute(obs)
print(obs.date, satellite.az, satellite.alt,satellite.range/1000, satellite.range_velocity/1000)

obs.date = '2015/1/1 04:50:00'
satellite.compute(obs)
print(obs.date, satellite.az, satellite.alt,satellite.range/1000, satellite.range_velocity/1000)

obs.pressure = 0 # disable atmospheric refraction
print('Pyephem No atmospheric refraction')

obs.date = '2015/1/1 04:30:00'
satellite.compute(obs)
print(obs.date, satellite.az, satellite.alt,satellite.range/1000, satellite.range_velocity/1000)

obs.date = '2015/1/1 04:40:00'
satellite.compute(obs)
print(obs.date, satellite.az, satellite.alt,satellite.range/1000, satellite.range_velocity/1000)

obs.date = '2015/1/1 04:50:00'
satellite.compute(obs)
print(obs.date, satellite.az, satellite.alt,satellite.range/1000, satellite.range_velocity/1000)

print('10 s timing')
obs.date = '2015/1/1 04:30:00'
satellite.compute(obs)
print(obs.date, satellite.range/1000, satellite.range_velocity/1000)
obs.date = '2015/1/1 04:30:10'
satellite.compute(obs)
print(obs.date, satellite.range/1000, satellite.range_velocity/1000)

Pypredict

import predict
import datetime
import time

format = '%Y/%m/%d %H:%M:%S'

tle = """ISS (ZARYA)
1 25544U 98067A   15096.52834639  .00016216  00000-0  24016-3 0  9993
2 25544  51.6469  82.0200 0006014 185.1879 274.8446 15.55408008936880"""
qth = (0, 10, 0)  # lat (N), long (W), alt (meters)

#expect time as epoch time float

time= (datetime.datetime.strptime('2015/1/1 04:30:00', format)  -datetime.datetime(1970,1,1)).total_seconds()
result = predict.observe(tle, qth, time)
print time
print result

time= (datetime.datetime.strptime('2015/1/1 04:40:00', format)  -datetime.datetime(1970,1,1)).total_seconds()
result = predict.observe(tle, qth, time)
print time
print result

time= (datetime.datetime.strptime('2015/1/1 04:50:00', format)  -datetime.datetime(1970,1,1)).total_seconds()
result = predict.observe(tle, qth, time)
print time
print result

Debug output of Gpredict and PyEphem

PyPredict

Name = ISS (ZARYA)
current jd = 2457023.68750
current mjd = 42003.7
satellite jd = 2457119.02835
satellite mjd = 42099
SiteLat = 0
SiteLong = 6.28319
SiteAltitude = 0
se_EPOCH  :    115096.52834638999775052071
se_XNO    :         0.06786747737871574870
se_XINCL  :         0.90140843391418457031
se_XNODEO :         1.43151903152465820312
se_EO     :         0.00060139998095110059
se_OMEGAO :         3.23213863372802734375
se_XMO    :         4.79694318771362304688
se_BSTAR  :         0.00024016000679694116
se_XNDT20 :         0.00000000049135865048
se_orbit  :                          93688
dt        :   -137290.81880159676074981689
CrntTime = 42004.2
SatX = -3807.5
SatY = 2844.85
SatZ = -4854.26
Radius = 6793.68
SatVX = -5.72752
SatVY = -3.69533
SatVZ = 2.32194
SiteX = -6239.11
SiteY = 1324.55
SiteZ = 0
SiteVX = -0.0965879
SiteVY = -0.454963
Height = 426.426
SSPLat = -0.795946
SSPLong = 0.432494
Azimuth = 3.53102
Elevation = -0.374582
Range = 5638.07
RangeRate = -5.30142
(2015/1/1 04:30:00, 5638.0685, -5.3014228515625)

Gpredict

time: 2457023,687500
pos obs: -6239,093574, 1324,506494, 0,000000
pos sat: -3807,793748, 2844,641722, -4854,112635
vel obs: -0,096585, -0,454962, 0,000000
vel sat: -6,088242, -3,928388, 2,468585

Gpredict (sgp_math.h)

/------------------------------------------------------------------/

/* Converts the satellite's position and velocity  */
/* vectors from normalised values to km and km/sec */ 
void
Convert_Sat_State( vector_t *pos, vector_t *vel )
{
      Scale_Vector( xkmper, pos );
      Scale_Vector( xkmper*xmnpda/secday, vel );

} /* Procedure Convert_Sat_State */

Ephem (Libastro)

*SatX = ERAD*posvec.x/1000; /* earth radii to km */
*SatY = ERAD*posvec.y/1000;
*SatZ = ERAD*posvec.z/1000;
*SatVX = 100*velvec.x;      /* ?? */
*SatVY = 100*velvec.y;
*SatVZ = 100*velvec.z;

Answer 2

Updating to the most recent release of pyephem (I tried V3.7.6.0) seems to solve the problem. The range rate now agrees closely with the values given by other commonly used tracking software.