Pyephem seems to compute wrong range_velocity

[Box88]

I am comparing the results from satellite trackers (e.g. Gpredict, Orbitron) and pyephem results. The range is the same but the range_velocity is different. Any idea?

[drbitboy]

That is not a lot of info for debugging, so with a lot of assumptions, just the obvious:

The satellite-observer range is correct, so the basic position calculation is probably correct, which means the satellite position is probably correct in whatever base reffrm (reference frame (e.g. J2000 inertial?)) it is stored, and the observer position is also probably correct in whatever base reffrm it is stored. Assuming the base satellite position is correct, I would be surprised if the satellite velocity in that base reffrm were wrong, so that leaves the observer velocity. You don't say what the observer is, but if it's a point on the surface of the earth, are all calculations (Gpredct, Orbitron, pyephem) using the same velocity calculation to account for the earth's rotation, adjusted for latitude? A diagnostic for this would be a range_velocity error of up to 900nm/h at the equator for a satellite on the horizon to zero for a satellite overhead, assume range_velocity = dRange/dt.

How about a test case? Define the observer (earth lat, lonWest, R) and the satellite. If you give me a case with the moon then I can do a quick check using NAIF/SPICE.

[Box88]

This is the simple python script I am using to compute the range and the range-rate. I tested it with different cubesats and ground station locations.

import ephem
import time
import datetime
import urllib
from itertools import izip, count
#TLE
myfile = urllib.urlopen('http://www.celestrak.com/NORAD/elements/cubesat.txt')
lines = myfile.readlines()
data = []
line1 = 'HUMSAT-D'
for num, line in izip(count(),lines) :
    if line1 in line:
        line2 = lines[num+1]
        line3 = lines[num+2]
satellite = ephem.readtle(line1, line2, line3) # create ephem object from tle information
while True:
    city = ephem.Observer() # recreate Observer with current time
    city.lon, city.lat, city.elevation = '-8.6881' , '42.1702' , 440
    city.date = datetime.datetime.utcnow() # set the current UTC time
    satellite.compute(city)
    RangeRate = satellite.range_velocity/1000 # get RangeRate in km/sec
    print ("RangeRate: " + str(RangeRate))
    Range = satellite.range/1000 # get Range in km
    print ("Range: " + str(Range))
    time.sleep(1)

[drbitboy]

I just thought of a simple diagnostic: which one gives near-zero values for geosynchronous satellites?

[drbitboy]

i proabably meant geostationary

[Box88]

Good idea! I tried it using intelsat 22 (an operative geosynchronous satellite). The results from Gpredict are: Range = 41903 km Range Rate = 0.000 km/s The results from pyephem are: Range = 41892 km Range Rate = -0.020122 km/s

[drbitboy]

yeah i did the same thing with GOES-2 and -6, and put the ground station under GOES-2; both range_vel's are order .01

eta: GOES-2 are about 70deg apart, so errors in ground station velocity should show up.

eta2: also, CELESTIS 03 appears to be polar and just flew nearly over that sub-GOES-2 ground station; range velocity went from -6 to -0.06 (two orders of magnitude) and is now climbing into positive range. So that looks about right

[brandon-rhodes]

To diagnose a problem with range_velocity I will need (a) a Python script that builds an observer, builds a satellite, and prints one single range at a specific date and time; and (b) a cut-and-paste from another prediction service for that satellite, observer, and time showing roughly the same position (because if the position itself is off too far then there are other problems we track down!) but showing a different range velocity.

As for the numbers for IntelSat 22 that you pasted in above: the numbers from Gpredict look clearly wrong, because any position on Earth will be moving and thus will experience a velocity relative to the geostationary satellite. Unless, of course, your Observer is at the North or South pole — but that is not what it looks like from your example code. :)

The reason that I need you to paste in a self-contained script that contains inline orbital elements in a string is that satellite orbital elements go out of date so quickly — after just one or two days, they are often several kilometers off, and so celestrak replaces them with new elements. So unless you include the TLE in your code, my attempt to re-run your code will probably be with a new set of elements that I grab from Celestrak, and will give different numbers that I can't directly compare to understand what is going on. Thanks!

[drbitboy]

https://gist.github.com/drbitboy/8786954

That's my code; it doesn't do the fixed elements yet (will have to make it read ~/.config/Gpredict/satdata/.

From first-order eyeballing, the change in Range over 10s of running is != ten times the Range_velocity; I wouldn't expect it to be exact, but it should be close.

eta: Landing at DEN shortly; will look at this later tonight.

[brandon-rhodes]

Thanks so much for the script! Here are two of the ten-second-separated readouts, from my machine:

2014/2/3 16:38:43:  RangeRate=     2.052m/s; Range=     10428.1m; Name=CELESTIS-02
2014/2/3 16:38:43:  RangeRate=     0.496m/s; Range=      4218.4m; Name=CELESTIS-03
2014/2/3 16:38:43:  RangeRate=     0.014m/s; Range=     36048.2m; Name=GOES 2 [-]
2014/2/3 16:38:43:  RangeRate=     0.012m/s; Range=     41112.0m; Name=GOES 6 [-]

2014/2/3 16:38:53:  RangeRate=     2.041m/s; Range=     10449.6m; Name=CELESTIS-02
2014/2/3 16:38:53:  RangeRate=     0.592m/s; Range=      4223.9m; Name=CELESTIS-03
2014/2/3 16:38:53:  RangeRate=     0.014m/s; Range=     36048.3m; Name=GOES 2 [-]
2014/2/3 16:38:53:  RangeRate=     0.012m/s; Range=     41111.8m; Name=GOES 6 [-]

Everything here is exactly as I would expect. The difference between the first and second range for each satellite is almost exactly 10 times the rate per second.

[drbitboy]

Not equal On Feb 3, 2014 9:36 AM, "Brandon Rhodes" notifications@github.com wrote:

Do you mean "about equal" or "does not equal" when you type "!="? :)

Reply to this email directly or view it on GitHubhttps://github.com/brandon-rhodes/pyephem/issues/34#issuecomment-33972366 .

[drbitboy]

Close but off about 5%; I would expect better than that.

This is all done in double precision, right? On Feb 3, 2014 9:41 AM, "Brandon Rhodes" notifications@github.com wrote:

Thanks so much for the script! Here are two of the ten-second-separated readouts, from my machine:

2014/2/3 16:38:43: RangeRate= 2.052m/s; Range= 10428.1m; Name=CELESTIS-02 2014/2/3 16:38:43: RangeRate= 0.496m/s; Range= 4218.4m; Name=CELESTIS-03 2014/2/3 16:38:43: RangeRate= 0.014m/s; Range= 36048.2m; Name=GOES 2 [-] 2014/2/3 16:38:43: RangeRate= 0.012m/s; Range= 41112.0m; Name=GOES 6 [-]

2014/2/3 16:38:53: RangeRate= 2.041m/s; Range= 10449.6m; Name=CELESTIS-02 2014/2/3 16:38:53: RangeRate= 0.592m/s; Range= 4223.9m; Name=CELESTIS-03 2014/2/3 16:38:53: RangeRate= 0.014m/s; Range= 36048.3m; Name=GOES 2 [-] 2014/2/3 16:38:53: RangeRate= 0.012m/s; Range= 41111.8m; Name=GOES 6 [-]

Everything here is exactly as I would expect. The difference between the first and second range for each satellite is almost exactly 10 times the rate per second.

Reply to this email directly or view it on GitHubhttps://github.com/brandon-rhodes/pyephem/issues/34#issuecomment-33973017 .

[drbitboy]

Another trick:

Look at 3 groups over 20s (0, 10,20)

Middle group is avg velocity, double and shift decimal (x10) to get expected approx 20s delta.

End groups are start and end range to calculate actual delta over 20s

Hmm, switch to 5s delta an wliminate doubling

Also add decimal to range (%12.2f)

[drbitboy]

Will save goes.txt and other.txt as local files, and inhibit tle update in Gpredict.

That will ensure same inputs are used. On Feb 3, 2014 10:06 AM, "Brian Carcich" drbitboy@gmail.com wrote:

Another trick:

Look at 3 groups over 20s (0, 10,20)

Middle group is avg velocity, double and shift decimal (x10) to get expected approx 20s delta.

End groups are start and end range to calculate actual delta over 20s

Hmm, switch to 5s delta an wliminate doubling

Also add decimal to range (%12.2f)

[drbitboy]

I am not sure I can agree with what you say here.

Any lat-lon observer ("city") is geostationary and so should have zero velocity wrt any other geostationary objects including satellites,

by definition!

(We're not modeling plate tectonics ;-)

On Feb 3, 2014 9:07 AM, "Brandon Rhodes" notifications@github.com wrote:

[...] As for the numbers for IntelSat 22 that you pasted in above: the numbers from Gpredict look clearly wrong, because any position on Earth will be moving and thus will experience a velocity relative to the geostationary satellite.

[dcajacob]

Geostationary is an academic term. No satellite is in such a perfectly synchronous and equatorial orbit so as to be perfectly geostationary wrt to any observer. There will be some (small) relative velocity component and doppler. Occasionally, you will hit a zero point, like in a LEO pass when the vectors align just right, but this is not going to be the case in general. Plate tectonics have nothing to do with this, it's orbital mechanics.

Very Respectfully,

Dan CaJacob

On Mon, Feb 3, 2014 at 12:47 PM, Brian Carcich notifications@github.comwrote:

I am not sure I can agree with what you say here.

Any lat-lon observer ("city") is geostationary and so should have zero velocity wrt any other geostationary objects including satellites,

by definition!

(We're not modeling plate tectonics ;-)

On Feb 3, 2014 9:07 AM, "Brandon Rhodes" notifications@github.com wrote:

[...]

As for the numbers for IntelSat 22 that you pasted in above: the numbers from Gpredict look clearly wrong, because any position on Earth will be moving and thus will experience a velocity relative to the geostationary satellite.

—

Reply to this email directly or view it on GitHubhttps://github.com/brandon-rhodes/pyephem/issues/34#issuecomment-33980001 .

[drbitboy]

Of course, eppur si muovi (summat like that - it's a West Wing episode title;-), and in my earlier comment I qualified the velocity as "near-zero" IIRC.

But I am afk and getting tired of typing on phone so I left it off here.

But it will be and remain very near zero

Equally respectfully,

Brian Carcich

On Feb 3, 2014 10:51 AM, "dcajacob" notifications@github.com wrote:

Geostationary is an academic term. No satellite is in such a perfectly synchronous and equatorial orbit so as to be perfectly geostationary wrt to any observer. There will be some (small) relative velocity component and doppler. Occasionally, you will hit a zero point, like in a LEO pass when the vectors align just right, but this is not going to be the case in general. Plate tectonics have nothing to do with this, it's orbital mechanics.

Very Respectfully, [...]

[drbitboy]

Of course, eppur suo movi (summat like that - it's a West Wing episode title;-), and in my earlier comment I qualified the velocity as "near-zero" IIRC On Feb 3, 2014 10:51 AM, "dcajacob" notifications@github.com wrote:

Geostationary is an academic term. No satellite is in such a perfectly synchronous and equatorial orbit so as to be perfectly geostationary wrt to any observer. There will be some (small) relative velocity component and doppler. Occasionally, you will hit a zero point, like in a LEO pass when the vectors align just right, but this is not going to be the case in general. Plate tectonics have nothing to do with this, it's orbital mechanics.

Very Respectfully,

Dan CaJacob

On Mon, Feb 3, 2014 at 12:47 PM, Brian Carcich notifications@github.comwrote:

I am not sure I can agree with what you say here.

Any lat-lon observer ("city") is geostationary and so should have zero velocity wrt any other geostationary objects including satellites,

by definition!

(We're not modeling plate tectonics ;-)

On Feb 3, 2014 9:07 AM, "Brandon Rhodes" notifications@github.com wrote:

[...]

As for the numbers for IntelSat 22 that you pasted in above: the numbers from Gpredict look clearly wrong, because any position on Earth will be moving and thus will experience a velocity relative to the geostationary satellite.

Reply to this email directly or view it on GitHub< https://github.com/brandon-rhodes/pyephem/issues/34#issuecomment-33980001>

.

Reply to this email directly or view it on GitHubhttps://github.com/brandon-rhodes/pyephem/issues/34#issuecomment-33980488 .

[drbitboy]

D'Oh sent that twice, sorry

Anyway, great discussion

About minutae, but great nonetheless

[brandon-rhodes]

Note that your script waits more than 10 seconds between runs, because it takes some fraction of a second to do the calculation, then waits an additional ten full seconds before starting the computation again. For more accuracy, do:

while True:
    t0 = time.time()
    ...
    seconds_remaining = t0 + 10.0 - time.time()
    time.sleep(seconds_remaining)

Even with this improvement to your code, though, I am still seeing the few-percent (maybe 3% or 4%) difference that you are reporting between the range rate, which always seems a bit too low, and the actual rate at which the range itself is changing.

When I have time later this evening, I am going to try this same computation with Skyfield and see whether I can get a rate out or not.

[Box88]

I would say that the doppler shift for geostationary satellite is zero and since the doppler shift = -(range_velocity/speed_of_light)*f_0, I would expect that the range_velocity is zero. Am I wrong?

[dcajacob]

The doppler shift (and therefore range rate) for "geostationary" satellites will be very near zero, but that doesn't mean that it is zero. The doppler is probably negligible for most geo sats, but that is not the issue. The goal is to have an accurate orbit propagator.

I suggest downloading the free version of STK to use as a benchmark against any other propagators. Predict and others like Orbitron probably return a rounded, or reduced-precision doppler/range rate and may even zero it erroneously if they deem that it fits the definition of geostationary.

Very Respectfully,

Dan CaJacob

On Mon, Feb 3, 2014 at 3:28 PM, Box88 notifications@github.com wrote:

I would say that the doppler shift for geostationary satellite is zero and since the doppler shift = -(range_velocity/speed_of_light)*f_0, I would expect that the range_velocity is zero. Am I wrong?

— Reply to this email directly or view it on GitHubhttps://github.com/brandon-rhodes/pyephem/issues/34#issuecomment-33996160 .

[brandon-rhodes]

In this case, I should note, my immediate goal is quite modest: it would be nice of the range velocity given by PyEphem really predicted, instead of kind-of-within-a-few-percent-predicted, the rate at which the range velocity was actually changing :)

[drbitboy]

Also note that a difference in observer altitude, or reference radius (geoid?) with respect to the earth's center of mass, would affect the range_velocity.

[drbitboy]

Also don't forget this, in case one of the methods being compared is FORTRAN-based:

$ cat x.f ; make x ; ./x

  programx
  implicitnone
  doubleprecisionsp,dp

  datasp,dp/1.19459,1.19459D0/

  print'(e20.12,x,a)',dp,'Double',sp,'Single',dp-sp,'Difference'
  end

gfortran x.f -lsqlite3 -L/usr/lib/x86_64-linux-gnu -lcurl -Wl,-Bsymbolic-functions -Wl,-z,relro -o x

0.119459000000E+01 Double 0.119458997250E+01 Single 0.275039673259E-07 Difference

[drbitboy]

[ETA: enhanced and rearranged plot description] [ETA2: correct output line count] [ETA3: fixed one image]

Here is a decent script, with self-contained TLE data:

https://gist.github.com/drbitboy/8886697

The results are below, showing about 5% error between the difference calculated from two Observer.range values 10s apart, and the same difference predicted via Observer.range_velocity (rate) values; there are two plots, one full and one zoomed, for each of the four satellites.

The calculated errors use the rates at the start and end of a time period, as well as a mean rate.

N.B. This is a first-order (linear) approximation only (Trapezoidal Rule) that does not consider second- and higher-order effects, but the time step is small enough (dT ~ 0) that the accuracy is adequate to justify the conclusion.

I think the point here is that the Observer.range and Observer.range_velocity values for pyephem are not self-consistent, whether or not such values are correct in Gpredict, Orbitron, etc.

See further notes at the end of this comment.

Notes:

The Python script produces ~7klines of output, and requires matplotlib to run to completion.

The asymptotes probably occur when dR/dt approaches zero, but even away from those, the error seems to be about 5%.

I suspect the noise in the GOES satellite results is due to roundoff of small numbers, but even so the large errors away from the asymptotes should not occur.

[drbitboy]

Hmm, I just realized that the errors running at a constant level (5-6%) for each TLE means there is something systemic here - my guess is that the issue is somewhere in the evaluation of the inertial velocity of the Observer.

[ETA:] On second thought, maybe not the observer: I put the Observer at different latitudes and it's 5-6% for all cases. This is definitely something systemic.

[drbitboy]

There is also a minor inconsistency between EarthRadius (6378.16) in earthsat.c and XKMPER (6378.135) in sgp4.c and deepconst.h under libastro*/; probably not related to this issue.

[drbitboy]

http://xkcd.com/1318/

Actually, don't miss the tooltip.

[brandon-rhodes]

I agree that the range and range rate seem to disagree systematically. And the graphs above will hopefully suggest to an analyst how to bring the two into closer agreement. But as I am not myself the author of the routines in earthsat.c that are making this systematic error in the derivative, I am going to make the (doubtless disappointing) move of falling back on the mission statement for PyEphem:

PyEphem wraps the astronomy logic of XEphem and makes it available in Python

In both satellite range and range-rates here, I am seeing 100% agreement between PyEphem and the parent project XEphem:

Here is what Python tells me for the same objects and circumstances:

2014/2/3 12:00:00
GOES 2 [-]
39943.08
18.9045314789

2014/2/3 12:00:00
GOES 6 [-]
47234.088
12.3930282593

Given this perfect agreement, and in the absence of any suggestion as to how the derivative code can be fixed, I am going to close this issue temporarily (to keep my actionable to-do list at a manageable length) until we have a suggestion about how the discrepancy might be fixed. At the moment, PyEphem is at least doing what it ought in making XEphem results available to programs.

PyEphem always incorporates whatever improvements happen in XEphem, so if Elwood winds up fixing this in XEphem (should we mention it on the XEphem mailing list?), then the fix will appear in the next version of PyEphem.

I am always hesitant to include fixes in PyEphem that have not been made upstream in XEphem since then I have to field questions from users who know both tools about why their results do not agree. In a case like this — where the output does just look plain wrong — I might bend the rules if someone comes up with a fix, and just hope that it eventually gets made upstream as well if we let them know on the mailing list. So, if someone finds the flaw, please re-open this issue and let us know — thanks!

[johandc]

I just came across the same error in our satellite frequency tracking application. I come to the same conclusion of about 5-6% wrong values when comparing to Gpredict. Have this issue been reported upstream with XEphem? I would hate to go back to the sgp4 module again. :)

[brandon-rhodes]

I do not remember seeing any questions about it on the XEphem mailing list, but I could easily have missed it. Here is their discussion list, in case you want to do your own search or asking the question yourself — good luck!

https://groups.yahoo.com/neo/groups/xephem/info

[alexschultze]

Has there any progress on the discrepancy between xephem and gpredict? Possibly the difference is caused by the fact that Xephem uses atmospheric fraction for its prediction. (So maybe its a feature?)

http://stackoverflow.com/questions/23551019/range-rate-for-doppler-calculation-and-atmospheric-refraction-compensation?rq=1