Artificial satellite code in C/C++

Artificial satellite code in C/C++

Last revised 26 Jun 2007

Click here to download the source code and SAT_CODE.DLL (128728 bytes)

Other source code implementations of SGP4/SDP4

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 expect to be 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 WATCOM C/C++ 10, 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:

(26 & 13 Jun 2007) Aleks ???, author of the Heavensat satellite tracking software, located some divide-by-zero problem with satellites with zero inclination, in the Deep_dpinit function. Changing this line:

         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.

(20 Oct 2006) I've been swapping an e-mail or two with Paul Crawford, one of the authors of a paper entitled "Revisiting Spacetrack Report #3." I can highly recommend this paper. These researchers caught several problems I'd overlooked, and made some good improvements. They fixed issues that can occur for eccentricities near zero and one; smoothed out the handling of the Lyddane problem; and corrected a few constants and coefficients. Furthermore, they described the fixes in a manner that allowed me to duplicate them in my code. On the test suite of TLEs mentioned in the above paper, we now match to within a few millimeters, worst case.

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.

(20 Oct 2006) The integration scheme in Deep_dpsec() has been changed extensively (again). Instead of splitting up the integration range into equal steps of 720 minutes or less, steps of 720 minutes are taken until the final step. I did that in order to match the results of the Dundee code. Separately, I realized that the second-order Taylor series integration method used in this part of the code could almost trivially be extended to n-order integration; the higher derivatives fall out with very little math required, once the first two are computed.

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

(28 Aug 2006) From the beginning, there was code in get_el.cpp to take the two input lines of a TLE, and parse them out into a tle_t structure (the function parse_elements() does this.) Now, there's a function in the new file tle_out.cpp that does the reverse of this: write_elements_in_tle_format() will put your tle_t structure back out into the two-line format.

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.

(14 May 2005) Caught a small bug in the Deep_dpinit() function in deep.cpp. The code checked to see if the TLE epoch matched that of the previous TLE; if it did, some computations could be skipped, because the variables wouldn't change. That was true back when the variables in question were stored statically (and would still be true in certain conditions now), but the variables are now stored in the deep_arg_t structure; I was running into uninitialized variable problems when feeding the code a series of TLEs with identical epochs.

(15 April 2005) I'm now looking at generating TLEs from my Find_Orb orbit determination software. (Actually, it already generates them, but they're not very accurate. I just convert a state vector to classical orbital elements, without taking into account the fact that SxPx elements don't quite match up to classical ones.) Anyway, doing this has caused me to add some comments to norad.h concerning units, and to clean up get_el.h a little. No actual changes in functionality are involved; I just needed to make sure that I (and others) could see what was going on.

(16 February 2005) After almost two years of no changes, I regarded this code as rather stable. But today, I received a set of TLEs wherein the decade character was a space, rather than a zero. A small fix to get_el.cpp allowed the code to function despite this oddity.

(22 April 2003) The dpsec function in deep.cpp has been a long-time source of confusion. It contained a very bizarre numerical integration routine which turns out to have numerous Bad Things in it, resulting in unnecessary complexity and performance degredation. (But not, as best I can tell, actual bugs.) I've fixed these; the resulting code is a lot easier to follow. Explaining how the code used to work and what I changed is a long story; click here for details.

(22 April 2003) Earlier, I recognized that the SDP4 and SGP4 code had large chunks of common code, and moved them to common.cpp. It's not so obvious, but their initialization routines, SGP4_init and SDP4_init, also had a lot of redundant code. I've extracted that into a new sxpx_init function, eliminating a lot of waste.

(29 March 2003) Peter Birtwhistle ran into a problem with sat_id that was traceable to a very high-drag satellite. Projected a few days into the future, this satellite ended up with an "orbit" a few kilometers above the center of the earth, and generated a math error (square root of a negative number). I evaded this by adding some error checks in common.cpp. If certain parameters exceed their proper limits, the state vector is now set to all zeroes and the function returns. This is at least more "graceful" than having it crash and/or return utter garbage.

(8 March 2003) There is a new example program, sat_id.cpp. This is similar to obs_tes2.cpp (in fact, large chunks of the code are identical). The only real change is that it reads through a file containing astrometric observations in the MPC format. As it comes across them, it looks for satellites within a given radius of those coordinates as seen from that location at that time. This is both an "example" program for programmers and a tool for non-programming astronomers, who routinely get images of satellites, but don't know if they might actually be Earth-passing asteroids. There is a writeup describing the program from their viewpoint here.

(12 December 2002) There is a new example program, obs_tes2.cpp, showing how to use the routines to get a list of satellites within a given search radius of a given RA/dec, as seen from a given lat/lon/altitude.

(20 October 2002) Mostly, a bunch of cosmetic changes. I realized that large chunks of the SGP4 and SDP4 code are in common. With some more work, I can probably eliminate a lot of redundant code (and, thereby, opportunities for bugs.) So far, all I've done is to combine a lot of code in the new common.cpp code. This still leaves a lot of redundancy in the init routines. All of this did require me to rearrange some of the items in the params array and in the deep_arg_t structures. (Nothing added or changed; just swapped around so that the order in SGP4 would match that in SDP4.)

(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 realized that in the DLL version, the functions exported were not declared as extern "C". They would therefore not be callable from a Visual Basic (or, probably, a Pascal or FORTRAN) app. I've fixed this. (Only the DLL version is affected by this; the changes compile to "no change" in the other versions.)

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

(23 September 2002) I've added a few functions to support topocentric ephemerides. The method for computing the satellite position hasn't changed, but there are now some functions for computing an "observer" position, and for converting the difference between the two into an RA/dec and distance. Click here for details.

(17 August 2002) Jim Lewis pointed out a page on the STK site discussing some of the changes AGI had to make to the SxP4 code. One of those changes was to avoid problems with zero-eccentricity satellites. We did a little hunting and found two cases where this would lead to a division by zero, and made two small changes (both in sgp4.cpp and marked with '17 Aug 2002') to repair them.

(16 August 2002) Jim Lewis pointed out a slight blunder in the dpper code in norad.cpp: two variables, 'pgh' and 'ph', were being modified between calls in ways that ought not to have been happening. In most test cases, this didn't matter; in others, it could lead to weird results. A couple of lines have been changed to kill this bug.

(2 August 2002) Jim Lewis, of Ten Sigma Research, pointed out that the sources in the sat_code.zip file were a bit confused; one of the older sources had been left in, while several new ones were omitted! This was due to a mistake I made in the batch file used to create the ZIP file. If you downloaded sat_code.zip before 2 August, you should definitely click here to download the corrected version.

I renamed *.c to *.cpp, just to get extra compiler checking, and tried it under MS Visual C/C++ 5.0, WATCOM C 10, MS Visual C/C++ 1.0 (to get a 16-bit DOS app), and Linux gcc. No problems under any of these, no alterations required. Renaming back to .c should also be okay.

There are now two different test programs. test_sat.cpp the "test data" values from Spacetrack 3 stored in it, and if you run test_sat with a command line argument, the differences between the "test data" values and the ones actually being computed is shown. Any errors stand out as non-zero values. At the suggestion of Thierry Marais, I added six other TLEs to exercise assorted possible execution paths for SDPx.

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

I added get_el.c, which contains the function needed to parse TLE data into the tle_t structures.

After hearing about assorted errors in the original Spacetrack Fortran code, I contacted Rob Matson after finding these posts:

 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.

The AcTan() function is essentially the C atan2( ) function rearranged so that the output runs from 0 to 2π (instead of -π to +π). So it can be reduced to a couple of lines of code.

Similarly, the Modulus() and FMod2p() functions are slight variants of the fmod() library functions. I wiped Modulus() out completely.

In a few places, items were not explicitly cast from double to integer types. That caused my compiler to belch some warnings. Just for cosmetics, I added "(int)" and "(long)" in those places.

Calling pow() to square, cube, or raise a number to the fourth power is a Bad Idea in terms of speed, and occurred several times in the code. pow() has to bring in the heavy machinery of taking a logarithm and then an exp()-type call. In such cases, I revised the code slightly to just multiply the number by itself. (Not that this optimization is immensely helpful. Neoklis Kyriazis, who did the original port of this code from Fortran, already did the major thing I planned to do, namely, computing certain "constants" that are independent of 'tsince' during an initialization stage, so they don't get recomputed at each call. That helps a lot.)

You may sometimes want to get a position and not care about the velocity. I revised the code so that, if the 'vel' argument to any of the five satellite functions is NULL, the velocity goes uncomputed. (Saves peanuts in computation time, I admit, but spares you having to declare a velocity vector, compute it, and then ignore it.)

I moved certain #include statements from norad.h to the C files. norad.c really only needed math.h, test_sat.c only needed stdio.h and stdlib.h.

I broke up norad.c into components, for several reasons. It looks to me as if the SxP8 methods get almost no use, and I may not actually include them in my commercial software. Modularizing code is a Good Idea in any case. With the breakup, a person could include only SGP, or only SGP4, etc., thereby drastically reducing the amount of code brought into play.

I separated the 'initialization' step from the rest of the process for all five models. Thus, you now call (for example) SDP4_init() once after loading a TLE. After that, you can call SDP4 as desired. (This is essentially a cosmetic change. Previously, the SDP4 function would use a flag to determine if it should initialize.)

In the interests of further "encapsulating" the code and letting people avoid looking at the internals of the code, I arranged it so that most structures were kept internally, with the "public interface" declared in norad.h and "private routines" declared in norad_in.h. I also realized that in many cases, variables were declared with much greater scope than was absolutely necessary. Most of that sort of thing has been cleaned up now.

I eliminated the 'vector_t' type. The position and velocity vectors are handled simply as arrays of three doubles. This just reflects a personal preference; I like being able to write code such as

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.

Deep() had a slew of static variables that, it appeared to me, more logically belonged as part of the deep_arg_t structure. So I moved them there.

All five SxPx functions make no changes to tsince or the input TLE data, so I made those arguments 'const'.

Deep() logically ought to be three separate functions: one to initialize (Deep_dpinit), one for secular effects (Deep_dpsec), and one for periodic effects (Deep_dpper). I've made that division.

At one point, I got slightly differing results from test_sat, down at the millimeter level, but I couldn't see why there should be any difference, and suspected an uninitialized variable that might, at a future date, cause kilometer or teraparsec-level errors. I eventually traced it to π being given to a mere thirteen places in norad.c. That was overkill in the Bad Old Days of single-precision, but double-precision values are good to sixteen or seventeen decimal places. I added a few more digits and made sure that constants such as π/180 would be accurately rendered.

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

When compiled with full optimizations with Microsoft Visual C, the select_ephemeris() function misbehaved to the extent of setting the DEEP_SPACE_EPHEM_FLAG incorrectly on occasion. I traced it to this line (in basics.cpp):

   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.")

The code set the satellite "epoch" to be one read directly from a TLE. Day 41 of year 2002 would be 02041, day 333 of year 1997 would be 97333. This is very TLE-specific; almost all my code (and most astronomical code) works in terms of JD (and internally, norad.c has to convert this sort of "epoch" to a JD before it can do anything with it.) Fixing this involved some modest changes to ThetaG() and allowed the complete removal of Julian_Date_of_Year().

(25 Jun 2002) After a brief bug scare (I thought the code might be producing wrong values), I went through deep.cpp and did a lot of housekeeping. The actual functionality of the code is unchanged, but the logic is much easier to follow, constant variables are set to be 'const', and most variables have been moved to be within their proper scope. The 'flags' turn out to be almost totally unnecessary, except for those indicating if the object is synchronous or in a resonant (12-hour, high-eccentricity) orbit.

(27 Jun 2002) Careful checking against Thierry Marais' code revealed some small differences. Some of these were due to the constants 'ck4', 's', 'qoms2t', and 'xke' being given to insufficient precision (see norad_in.h for details). Also, the ThetaG() function in deep.cpp had miscellaneous very strange problems, including a very subtle loss-of-precision one causing differences at the 2-centimeter level. (After assorted fixes, my code and Thierry's quite different code agree to the 10-micron level.)

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

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. He is one of the authors of this useful paper on SGP4/SDP4.

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.