Artificial satellite code in C/C++

Last revised 24 Sep 2009

After some hunting for C/C++ code to implement the assorted satellite theories SGP, SGP4, SGP8, and the "deep-space" variants SDP4 and SDP4, I found the source code posted by Neoklis Kyriazis. As he mentions in the readme.txt file, the C source code is derived from the NORAD FORTRAN routines in the Spacetrack #3 report, and from T S Kelso's Pascal routines (which also came from Spacetrack #3.) In my case, all I've done is to take Neoklis Kyriazis' source and revise it somewhat, in several directions.

Copyright: The code is derived from Neoklis Kyriazis' code, and you may want to read the COPYRIGH.TXT file from the source code on his site. But essentially, his code is in turn derived from, and therefore is itself, public domain materials, so that his code (and the code posted here) is in the public domain, usable for all purposes commercial and otherwise (I'm using it in my own commercial software; that's part of what is driving my interest in it.)

OSes/compilers supported thus far: Four 'make' files are provided: linmake for gcc on Linux, watcom.mak for OpenWATCOM , sat_code.mak for a statically-linked version in MS Visual C/C++, and sat_dll.mak for a version that creates a Windoze DLL. This last may prove useful for cross-language support; there has been at least one inquiry about using the library with Visual Basic, which can make use of DLLs.

Changes: Here are the differences between my code "du jour" and the original code provided by Neoklis Kyriazis:

         deep_arg->ssh = sh / deep_arg->sinio;

to:

         deep_arg->ssh = (deep_arg->sinio ? sh / deep_arg->sinio : 0.);

fixes the first problem. (The problem was already fixed in the Dundee and Center for Space versions described below, using slightly different logic.) The same problem occured about thirty lines further down: numbers are divided by deep_arg->sinio without checking to see if that quantity might be zero. I made a similar fix.

To check all this out, I added some test cases to test.tle with inclinations of 0 and 90 degrees. All appears to work.

Also, cubic polynomials are evaluated frequently in deep.cpp. I added a trivial function that takes the coefficients of the polynomial, mostly just to make the code a little neater.

Also, I made a small improvement so that when the code is in Dundee-compliant mode, it can save integration states from one call to the next. Previously, when in Dundee-compliant mode, this code always had to do the full integration from t=0 to t=tsince.

Also, I figured out what some of the 'mystery constants' in deep.cpp were, and did some documenting.

In the process of doing this, I made many changes simply to make it easier to follow what's going on. (For example, the code for Deep_dpper() is much more "logical" now.)

Note that to get agreement with the Dundee code, one must set a flag in the code for "Dundee compliance". An example of this is given in test2.cpp, with the line

         sxpx_set_implementation_param( SXPX_DUNDEE_COMPLIANCE, 1);

In the default, non-Dundee-compliant mode, the results from the SDP4() function may depend on what calls have been made to the function previously. For example, if you asked to know where a satellite was at times 1000 and 900 minutes after epoch, the code would integrate forward 1000 minutes, then integrate backward 100 minutes. If you asked for a position at tsince=900 minutes, the code would simply integrate forward, and the result would be slightly different.

Similarly, in the default mode, the effects of lunar/solar perturbations are only recomputed if tsince changes by more than 30 minutes. This is a generally safe thing to do, since those perturbations are slowly varying. But it does mean that if you asked for positions at tsince=1000 and tsince=1020, the perturbations as of 1000 minutes would be applied on the second call. So again, results are not perfectly replicable; they depend on previous calls to SDP4().

In Dundee-compliant mode, however, the integration always runs from epoch to tsince, and the perturbations are always recomputed. Those are the only differences from the default mode. The default mode is good for all purposes except benchmarking against the Dundee implementation; for that, I wanted to have exact agreement in results.

By default, this code still uses only the second-order terms. But I've found that if one goes to higher orders, you can increase the integration step size without a penalty in precision, and it runs faster. (Not usually an issue, unless you're using elderly TLEs. Unfortunately, due to restrictions in the distribution of TLEs, some of us are having to do just that.)

As a side effect of this, the tle_t structure had to be expanded a bit. It now contains the NORAD number and international designation, bulletin number, revolution number, classification (almost always "U", or "Unclassified"), and ephemeris type... in other words, all the bits and pieces of data given in a TLE.

(20 October 2002) There's a new dynamic.cpp file for Windows use only, for C/C++ programmers looking to bind to the sat_code.dll routines at run-time. (In other words, probably nobody will have much interest in it.) I intend to use this in my existing planetarium/star charting software. This program currently supports only SGP4; using run-time linking, it'll be able to use SDP4 and SxP8 functions if people download the sat_code.dll, without bulking up the (already bloated) main executable. People who don't care about artificial satellites (i.e., most users of the program) can simply not download the DLL.

Also included is a test3.cpp program that makes use of the new 'dynamic' functions, showing how they are used in comparison to the 'usual', statically-linked functions.

(25 September 2002) I also realized that the example output files were out of date. I've fixed this. You should be able to run test_sat and test2 and get the results given in result and test2.out, respectively.

However, the Spacetrack examples and the six TLEs still don't exercise all of the various paths through which the code can go. The behavior can depend on previous calls to the functions and so forth. test2.cpp, therefore, reads in test.tle and tries out a long list of test cases. The output of test2 in 32-bit Windows, 16-bit DOS, 32-bit DOS extender, and Linux forms agrees down to the last decimal place (about 10 microns).

 http://www.satobs.org/seesat/Oct-1997/0440.html 
 http://www.satobs.org/seesat/Oct-1997/0497.html 

The "Lyddane modifications" had already been made by Neoklis Kyriazis, but the lunisolar perturbation modifications (a matter of retaining some offsets as to what the perturbations were at the epoch) were not made. I made those modifications, and promptly got values a few km different from the original Spacetrack values. At present, if you compile norad.c with RETAIN_PERTURBATION_VALUES_AT_EPOCH #defined, you'll get that modification. But by default, it will (as apparently was done for the Spacetrack examples) leave those values at zero (by omitting a few lines of code.)

If you do include this modification, be aware that you will get mismatches between the code and the Spacetrack 3 examples of order 20 km. Omitting this modification (as is done by default) gets a match within about 6 m for SDP4, and about 100 meters for SDP8.

for( i = 0; i < 3; i++)
   do something with vector element i; 

as opposed to

do something with vector element x;
do it again,  but this time to vector element y;
do it again,  but this time to vector element z;

Also, I changed the units from Earth radii to kilometers. I've never been able to come up with a case where you would actually want the former.

   "...The primary purpose of the Data statement is to give names
   to constants; instead of referring to π as 3.141592653589793
   at every appearance, the variable Pi can be given that value
   with a Data statement and used instead of the longer form of
   the constant. This also simplifies modifying the program,
   should the value of π change."

    --  Fortran manual for Xerox Computers

   temp = ck2 * 1.5 * (r1*r1*3.0-1.0) / pow( 1.0-tle->eo*tle->eo, 1.5);

which sets 'temp' to INF, despite the fact that ck2, r1, and tle->eo had reasonable values. Changing it to read as follows:

   temp = ck2 * 1.5 * (r1*r1*3.0-1.0) * pow( 1.0-tle->eo*tle->eo, -1.5);

should, in theory, have changed nothing. In practice, it fixes the bug, so I assume an optimization error occured. (I mention it because in the past, I've sometimes seen "optimization errors" turn out to be "errors that were there all along, but optimizing made it more likely that they'd pop up at you.")

"Observer" functions: An example program showing how to compute topocentric ephemerides is given in obs_test.cpp. The topocentric observer's viewpoint is defined in terms of a longitude, latitude, and height in meters above the ellipsoid. The lat_alt_to_parallax() function converts these last two to "parallax constant" form, rho_sin_phi and rho_cos_phi. Next, given a particular time, the observer_cartesian_coords() function converts all this to a Cartesian, earth-centered vector, just like that returned by the satellite functions.

Now that you have the observer's coordinates and the satellite's coordinates, you can hand both off to the get_satellite_ra_dec_delta() to get an RA/dec and distance (in kilometers) to the satellite. Please note that, so far, everything has been done in the mean ecliptic of date; we'll usually want to use the epoch_of_date_to_j2000() function to convert the RA/dec to J2000 coordinates.

Acknowledgments: Thanks go first to Neoklis Kyriazis for supplying the source on which these efforts are based. I got some helpful data on test cases for SD* routines from Thierry Marais (the last six TLEs computed by test_sat are examples given by him that should exercise all the paths of those functions.) I was quite baffled by the assorted "versions" of SD* out there, and Rob Matson gave me a wonderful rundown of the tweaks I would have to make to the code to get the solar-lunar perturbations to work in a reasonable manner. Paul Gabriel and Mike McCants gave me some pushes in the right directions when I initially began looking into the problem of getting this sort of library of routines together. Jim Lewis tested the code on a Sun box and found several small errors in the process.

Other source code implementations of SGP4/SDP4: The code on this page was derived from Neoklis Kyriazis' source code, which in turn was derived from Dr. TS Kelso's Pascal implementation and the FORTRAN routines in NORAD's Spacetrack report #3. (The preceding link is in PDF, but you can cut-and-paste FORTRAN source from it.) Another C/C++ adaptation of Kelso's code is provided with GPS 2.4, a program intended primarily for tracking GPS satellites, but which is not limited to them.

Dr Kelso's code implements only SGP4 and SDP4. Dominik Brodowski has posted Pascal source for SGP and SxP8.

The implementation of SGP4/SDP4 at "Revisiting Spacetrack Report #3", by David A. Vallado, Paul Crawford, Richard Hujsak, and T. S. Kelso, is well-thought of by many, including me. (I swapped a series of e-mails with the authors, resulting in some bug fixes, almost all of them in my code and not theirs.) Source in C++, FORTRAN, Java, MATLAB, and Pascal is provided.

An advantage of this implementation is that the authors were at pains to make the reasoning behind their changes clear. That removes much of the mystery behind what can be a very mysterious chunk of code. (Of course, I'd like to think that I've achieved similar clarity!)

Paul Crawford has posted a C implementation of SDP4 and SGP4 here. It provided several bug fixes; in fact, as of 15 October 2006, my code and Crawford's code produce identical results on a variety of test cases.

You can also download C source code for the PREDICT software. This program runs on Linux, DOS, and Sharp Zaurus PDA. Looks as if SGP4 and SDP4 are implemented.

There is also an implementation of SGP4 in Java.

Claudio Pizzillo has announced an implementation of SGP, SGP4, SDP4, SGP8, and SDP8 in VB.NET.

Note that most of the bug fixes described for my code will also apply to these sources (except the Dundee code). In particular, almost all of them lack the change in dpper to get certain variables initialized at the right time; search through deep.cpp for SPACETRACK_3 for details. If you do decide to make use of a different implementation, I suggest you keep the URL for this page around; the bug comments will still be useful to you.