# brandon-rhodes/pyephem

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

# Pyephem seems to compute wrong range_velocity #34

Closed
opened this Issue Feb 3, 2014 · 32 comments

Projects
None yet
6 participants

### Box88 commented Feb 3, 2014

 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?
Contributor

### drbitboy commented Feb 3, 2014

 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 commented Feb 3, 2014

 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) ``````
Contributor

### drbitboy commented Feb 3, 2014

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

### drbitboy commented Feb 3, 2014

 i proabably meant geostationary

### Box88 commented Feb 3, 2014

 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
Contributor

### drbitboy commented Feb 3, 2014

 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
Owner

### brandon-rhodes commented Feb 3, 2014

 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!

Contributor

### drbitboy commented Feb 3, 2014

 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.
Owner

### brandon-rhodes commented Feb 3, 2014

 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.
Contributor

### drbitboy commented Feb 3, 2014

 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 .
Contributor

### drbitboy commented Feb 3, 2014

 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 .
Contributor

### drbitboy commented Feb 3, 2014

 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)
Contributor

### drbitboy commented Feb 3, 2014

 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)
Contributor

### drbitboy commented Feb 3, 2014

 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 commented Feb 3, 2014

 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 .
Contributor

### drbitboy commented Feb 3, 2014

 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, [...]
Contributor

### drbitboy commented Feb 3, 2014

 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 .
Contributor

### drbitboy commented Feb 3, 2014

 D'Oh sent that twice, sorry Anyway, great discussion About minutae, but great nonetheless
Owner

### brandon-rhodes commented Feb 3, 2014

 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 commented Feb 3, 2014

 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 commented Feb 3, 2014

 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 .
Owner

### brandon-rhodes commented Feb 3, 2014

 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 :)
Contributor

### drbitboy commented Feb 7, 2014

 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.
Contributor

### drbitboy commented Feb 7, 2014

 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
Contributor

### drbitboy commented Feb 8, 2014

 [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.
Contributor

### drbitboy commented Feb 8, 2014

 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.
Contributor

### drbitboy commented Feb 8, 2014

 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.
Contributor

### drbitboy commented Feb 11, 2014

 http://xkcd.com/1318/ Actually, don't miss the tooltip.
Owner

### brandon-rhodes commented May 21, 2014

 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 commented Nov 11, 2014

 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. :)
Owner

### brandon-rhodes commented Nov 18, 2014

 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
Contributor

### alexschultze commented Mar 30, 2015

 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
to join this conversation on GitHub. Already have an account? Sign in to comment