A Correction For Plan13 in C
The short version of this post is as follows: the 'C' versions of Plan13 all share a bug that make the algorithm unable reliably to track HEO and GEO satellites. To fix this bug find the part of the code that reads something like this:
The more in-depth version of the story follows.
Plan13 is a long-standing algorithm for predicting calculating satellite locations. Its original incarnation is found in a typically delightful 1990 article by James Miller, G3RUH, which AMSAT has wisely archive for anyone to read. Even if you don't speak BBC Basic, you'll probably understand what is going on in the code, and you can find BBC Basic interpreters to run the code. (Though in my experience when I compiled the Linux version of brandy it crashed on Miller's Plan13, and I couldn't get the MacOS version to compile. That leaves Windows, and it works well in both XP and Vista.)
In early 2000's Edson Pereira PU1JTE made a C port of the BBC basic code, which has proved to be very useful to many people, myself included. It retains the clarity of Miller's code, while adding some useful features from POSIX-compliant systems.
I have been empirically testing the reliability of this algorithm by comparing its output against that of Predict when using the keplerian elements of AO-07 and other low-earth orbit (LEO) amateur satellites. While doing this, a letter came to the amsat-bb list from Mark K6HX noting that he'd been using the Plan13 algorithm with Python, but had been disappointed to find that it was inaccurate when dealing with high earth orbit (HEO) satellites. This struck me as strange, since Miller's original article uses AO-10, a (now defunct) HEO satellite, as its example!
Nevertheless, I ran AO-10's current keplerian elements through my testing suite, and, sure enough, the results were disastrous. Plan13 in C wasn't even close to Predict, whose results I confirmed with SatPC32 for good measure. Something was up. I returned to the G3RUH code, which has the keplerian elements hard coded into it, and put in some recent AO-10 keplerian elements. Sure enough, the result of this 1993 program, running through a basic emulator, was spot-on with Predict and SatPC32.
Now was the time for a good ol' fashioned bug-hunt. I'm very new to C, but I know how to keep adding print statements until something fishy raises its head, and it was after the loop listed above that things seemed to go bad. The value EA is supposed to be iteratively calculated in that loop, but at the end of the loop the C version had a very different figure from the BBC basic version. One more print statement, and I found out why: the C version never takes more than one trip around the loop; but the BBC basic version was doing three.
The bug is in the use of the abs function. Unless specially called, it takes an integer as its input and returns an integer. EA is a float, so C helpfully casts it as an integer, rounding to 0. The absolute value of 0 is 0, so the loop figures it is good to go. In BBC Basic, of course, abs can take and return floating values, so it properly optimizes the value of EA.
The solution is to call the C function that provides absolute values for floating numbers, fabs. Once I did that, the output of my C code concurred with Predict et al. I'll run it through the test bed soon, but I'm pretty sure all will be well.
It will be interesting to see if this improves the performance of Plan13 with the usual fleet of LEO birds.
(Since doing this, I found that Howard G6LVB has properly implemented a floating point abs function in his port of Plan13 C to the PIC platform for the LVB tracker.)
/* Solve M = EA - EC * sin(EA) for EA given M, by Newton's method */
EA = M; /* Initail solution */
do
{
C = cos(EA);
S = sin(EA);
DNOM = 1.0 - EC * C;
D = (EA - EC * S - M) / DNOM; /* Change EA to better resolution */
EA = EA - D; /* by this amount until converged */
}
while (abs(D) > 1.0E-5);
and replace abs with fabs. (If you are using the original C port by Edson Pereira this is line 501.)The more in-depth version of the story follows.
Plan13 is a long-standing algorithm for predicting calculating satellite locations. Its original incarnation is found in a typically delightful 1990 article by James Miller, G3RUH, which AMSAT has wisely archive for anyone to read. Even if you don't speak BBC Basic, you'll probably understand what is going on in the code, and you can find BBC Basic interpreters to run the code. (Though in my experience when I compiled the Linux version of brandy it crashed on Miller's Plan13, and I couldn't get the MacOS version to compile. That leaves Windows, and it works well in both XP and Vista.)
In early 2000's Edson Pereira PU1JTE made a C port of the BBC basic code, which has proved to be very useful to many people, myself included. It retains the clarity of Miller's code, while adding some useful features from POSIX-compliant systems.
I have been empirically testing the reliability of this algorithm by comparing its output against that of Predict when using the keplerian elements of AO-07 and other low-earth orbit (LEO) amateur satellites. While doing this, a letter came to the amsat-bb list from Mark K6HX noting that he'd been using the Plan13 algorithm with Python, but had been disappointed to find that it was inaccurate when dealing with high earth orbit (HEO) satellites. This struck me as strange, since Miller's original article uses AO-10, a (now defunct) HEO satellite, as its example!
Nevertheless, I ran AO-10's current keplerian elements through my testing suite, and, sure enough, the results were disastrous. Plan13 in C wasn't even close to Predict, whose results I confirmed with SatPC32 for good measure. Something was up. I returned to the G3RUH code, which has the keplerian elements hard coded into it, and put in some recent AO-10 keplerian elements. Sure enough, the result of this 1993 program, running through a basic emulator, was spot-on with Predict and SatPC32.
Now was the time for a good ol' fashioned bug-hunt. I'm very new to C, but I know how to keep adding print statements until something fishy raises its head, and it was after the loop listed above that things seemed to go bad. The value EA is supposed to be iteratively calculated in that loop, but at the end of the loop the C version had a very different figure from the BBC basic version. One more print statement, and I found out why: the C version never takes more than one trip around the loop; but the BBC basic version was doing three.
The bug is in the use of the abs function. Unless specially called, it takes an integer as its input and returns an integer. EA is a float, so C helpfully casts it as an integer, rounding to 0. The absolute value of 0 is 0, so the loop figures it is good to go. In BBC Basic, of course, abs can take and return floating values, so it properly optimizes the value of EA.
The solution is to call the C function that provides absolute values for floating numbers, fabs. Once I did that, the output of my C code concurred with Predict et al. I'll run it through the test bed soon, but I'm pretty sure all will be well.
It will be interesting to see if this improves the performance of Plan13 with the usual fleet of LEO birds.
(Since doing this, I found that Howard G6LVB has properly implemented a floating point abs function in his port of Plan13 C to the PIC platform for the LVB tracker.)
Comments
Post a Comment