Astronomy C/C++ source code

Last updated 12 Jan 1999

From time to time, I'm asked to provide source code for doing some sort of astronomical calculation, or for providing direct access to the numerous compressed datasets on the Guide CD-ROM. In the future, I intend to gather information about and links to all such source code on this page.

The source code in question is written in C. I use an extension of .CPP just to force the compiler to do certain error checking (type casts, etc.) that C++ compilers do better than C compilers. The only use I make of C++ is in Windows user interface (MFC) code, and I haven't posted any of that yet.

I will post remaining code here from time to time. If you have any priorities (bits you would really like to see soon), please let me know; I can probably rearrange the order to get your code uploaded faster.

Porting: As much as possible, I've stuck with ANSI C. Most of the code has already been used in 16-bit DOS and Windows, and in 32-bit DOS and Windows. I suspect that much of it would port to Linux and similar systems without difficulty. The only warning I'd offer in this regard is that I haven't dealt with a "wrong-end" byte order platform yet, and I am quite certain that much of the code will crash in grisly fashion on such systems.

Documentation: Almost all the code has quite thorough descriptions of how the innards work. (This requirement is the main one slowing posting of code. Eventually, I expect that most of the underlying code in Guide, Find_Orb, and Charon will be posted on this site. But most of it is complex enough that posting it without explanation would not be very helpful.)

Copyright restrictions: All this code is free for use in non-commercial applications. If you wish to use any of it in a commercial application, please let me know; I'm not averse to a little bartering in such cases.

There is one exception to this: the RealSky/DSS image extraction software is entirely in the public domain. (Certain parts of the code came from the Space Telescope Science Institute, so I don't really have a lot of choice in the matter.)

For all the source code, I would ask that you inform me of any bugs you find, and if you make interesting improvements, I'd very much like to hear about them.

Source code for accessing Guide datasets and orbital elements: There is rather a lot of code in this category, and it has its own separate page. Click here for information on this subject.

Source code for RealSky/DSS image extraction: This is the same library used in Guide for image extraction; it will also be used in the next version of SkyMap Pro. You can read about it and download it from this page.

Source code for basic astronomical calculations: This file will eventually encompass all the low-level routines I use in Guide: nutation, aberration, computing planet positions from VSOP and PS_1996, refraction, lunar positions, the orientations of planets, coordinate systems, sidereal time...

However, at present, it just contains functions for computing comet/asteroid ephemerides, for extracting orbital data from the Guide datasets, and a few odds and ends. As of 21 Jan 1999, it also contains code for handling several calendrical systems: Julian, Gregorian, Hebrew, Persian, Islamic, and Chinese, and French Republican calendars. Click here to download the ZIPped code (about 35 KBytes). Here's some discussion of the meaning of the various source files supplied:

ASTEPHEM.CPP: The main function for the code to compute asteroid ephemerides. It demonstrates accessing orbital elements from the Guide disk, and then shows how to compute an RA/dec, magnitude, elongation, and other data based on using those elements to compute the asteroid position, and VSOP to compute the Earth's position. The most significant thing omitted right now is the topocentric offset.

ASTFUNCS.CPP: Some basic functions to set up asteroid/comet elements and to compute their locations (solving Kepler's equation and so forth.)

The Kepler's equation solution draws quite a bit of interest, so much so that there is a separate page discussing Kepler's Equation.

DATE.CPP: Code to handle calendar conversions in assorted systems (Hebrew, Islamic, Persian, Gregorian, etc.)

DELTA_T.CPP: Code to compute Delta_T = TD - UT for any date, using a lookup table for years 1620 to 2002 and extrapolations beyond those ranges.

EART2000.CPP: Code to compute the Earth's location relative to the Sun, in J2000 coordinates, using VSOP.

JD.CPP: Produces an "example" program showing the use of the functions in DATE.CPP. You can type in, say, "jd 10 7 1999", and be told what day 10 July 1999 corresponds to in assorted calendars. (You can also get the actual JD corresponding to that date, or type in, say, "JD 2451324.879" and get the calendar date... which was the original reason I wrote the program.)

MISCELL.CPP: Code for basic matrix functions (inverting an orthogonal matrix, rotation, setting identity matrices, etc.) Also code to handle the odd system used for variable star designations, and to provide a version of the C-language ctime() function that takes a JD and can work over a full time range (ctime only works between 1970 and 2028.)

OBLIQUIT.CPP: Computes the obliquity of the earth's equator for dates from -8000 to +12000.

PERSIAN.CPP: Very special-purpose code, written to demonstrate how the Persian calendar is computed. Most people will either skip this, or just use the functions in DATE.CPP without really needing to know how the functions were put together.

PRECESS.CPP: Code to build the matrix used to convert coordinates between epochs, plus code to use those matrices to convert vectors, RA/dec coordinates, and ecliptic coordinates between epochs.

VISLIMIT.CPP: This code computes the limiting visual magnitude, sky brightness, and extinction coefficients. It's basically lifted/ported from Brad Schaefer's article and code on pages 57-60, May 1998 Sky & Telescope, "To the Visual Limits".

Symplectic integration: This single source file is derived from code posted by David Whysong on the newsgroup sci.astro.amateur. I simplified the code a bit and added a "test code" main() function. You can now get the results of our collaboration at Dave's Web site. It's a very interesting way of doing numerical integration, on several counts.

The major use for it is in cases where one wants to maintain certain constants of motion (total energy, angular momentum, etc.) Usual methods of numerical integration (Runge-Kutta, Adams-Bashforth, etc.) don't necessarily do a good job of this over really long integration spans. This was the reason for David's interest in the method. He is a graduate student working on globular cluster simulations. If, say, the total energy of his simulated globular increased over billions of years, we may assume that the whole thing would mysteriously explode with stars flying everywhere. More details are given in the source code.

DIST.CPP: This code resulted from a discussion on the newsgroup sci.astro.amateur about how to compute the geodetic distance between two lat/lon points on the earth's surface. It incorporates two methods: one assumes a spherical earth, the other (much more complex) accounts, to first order, for the fact that the earth is an oblate spheroid.