Astronomy C/C++ source code
Last updated 6 June 2002
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; click here for info on source code in other languages, or here for info on other C/C++ source code for astronomy. 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. However, Mark Huss has produced a C++ implementation of much of this code, and you can click here to download the ZIPped C++ version (about 128 KBytes). If that link fails, try this one.
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.
Changes: Assorted changes have been made in the past with the regrettable absence of a change log. In the future, the .ZIP file will contain a CHANGES.TXT file, explaining what's different about the newly-posted files. (In the case of the 6 June 2002 update, the changes are mostly cosmetic, with a couple of small bug fixes and no really new functions introduced. A 'make' file for producing a 32-bit Windows DLL version of the library was added.)
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. A make file for Linux is provided, and the code works there (but I've only used it to create a Linux version of Find_Orb, and not much else.) Mark Huss' C++ version also compiles in Linux. 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. Some source has been posted for Find_Orb and Charon. 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 are two exceptions to this: the RealSky/DSS image extraction software and the FITS image compression software are entirely in the public domain. (Certain parts of the RealSky/DSS code came from the Space Telescope Science Institute (STScI), and the FITS compression software is a very slightly modified version of an STScI package.)
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.
Astronomy source code in other languages (so far, only Java): Personally, I haven't written much of anything in a language other than C or C++ since 1988. (Not a good thing, I realize.) I get occasional questions about source code in Java and Pascal (surprisingly, not in BASIC... though you can click here for a list of BASIC programs from Sky & Telescope.)
In addition to his C++ astronomy source code, Mark Huss has written a lot of astronomy source code in Java, and has updated and improved it a bit recently (March 2002) as described in the following e-mail:
hi Bill, I finally got things organized - the updated library, source, and info is now available on my website: http://mhuss.com/AstroLib.html. javadoc is an amazing thing - take a glance at http://mhuss.com/AstroLib/docs/api/index.html. regards, --mark huss
(And from a previous e-mail...)
This is what's currently running the lunar phase info and 'darkest hours' page on the DVAA website (http://dvaa.org). You should know, however, that I used Meeus 'mystic formulas' instead of applying your 'witness the truth thereof' philosophy in a couple of places for expediency ;-). regards, --mark
Some other implementation of assorted astronomical algorithms in C/C++: Mark Huss has written quite a bit of astronomy-related code in C++. P. J. Naughter has implemented much of Meeus' Astronomical Algorithms (second edition) in C++. You can click here for his code and info about it.
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 accessing JPL's DE ephemerides: This has grown to the point where it needed a page of its own. Click here for C/C++ source code for accessing JPL DE ephemerides.
Source code for RealSky/DSS image extraction: This is the same library used in Guide for image extraction; it's also used in version of SkyMap Pro, and, I suspect, in some other commercial astronomy software. You can read about it and download it from this page.
Source code for FITS image compression: A very slightly modified version of an STScI program. Click here for information and source code.
Source code and a DOS executable for tables of Galilean satellite events: This is the code I used to generate the lists of occultations, eclipses, transits, and shadow events shown in Guide. Click here for an example of the output of this program. The code is quite ugly, and I don't recommend its use very highly. Click here to download the source code and .EXE files (about 69 KBytes).
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...
At present, it has many of these, but not all. Here's a summary of the purpose of most of the .CPP files:
Click here to download the ZIPped code (about 93 KBytes). Here's some discussion of the meaning of the various source files supplied:ALT_AZ.CPP: Code to convert RA/dec to alt/az, and galactic coordinates to RA/dec.
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.
BIG_VSOP.CPP: Code to compute planetary positions from the full VSOP series, which is stored in binary form on the Guide CD-ROM. The logic is very similar to that used for the "truncated" VSOP (see VSOPSON.CPP), except that the function works by reading in pieces of a file; VSOPSON.CPP assumes the data file has been read into memory.
CLASSEL.CPP: Takes a state vector (position and velocity) and computes the "classical elements" (semimajor axis, eccentricity, mean anomaly, argument of perihelion, etc.) Currently, only the FIND_ORB orbit determination software really makes use of this; after it's computed an orbit in vector format, it uses this routine to produce something in a form suitable for human consumption.
COLORS.CPP and COLORS2.CPP: Given, say, a B-V color value, it's possible to compute a "pretty good" V-R value. Similar conversions are possible between most other color indices, with varying degrees of accuracy. The general form is a polynomial, say,
index2 = a0 + a1 * index + a2 * index2 + a3 * index3 +....
These resemble Taylor series, though I suspect they are actually curve-fit polynomials. I got them from some Fortran snippets supplied by Brian Skiff, and they have been converted to C here.
Certain conversions are handled by inverting the polynomial. The need to do this indicates a very unreliable color conversion, to be used only in cases of dire emergency.
Also, code to convert Tycho B and V to Johnson B and V is given, along with a little piece of test code. (That Tycho-to-Johnson algorithm is one based the book that was distributed with the Hipparcos data, "Introduction and Guide to the Data." COLORS2.CPP is a more recent, and probably better, Tycho-to-Johnson piece of source code.)
COM_FILE.CPP: This is pretty much Guide-specific code. Guide goes through some odd hoops to draw comets as speedily as possible. The tricks involved were absolutely necessary back when computers were slower than they are today, but are probably not all that helpful on modern PCs.
COSPAR.CPP: Code to compute the "planet orientation" matrix for the Sun, all planets, and quite a few satellites. A few years ago, the IAU set up a committee to define pole positions and longitude systems for these objects. The pole position is usually defined to be a linear function in RA/dec, with the meridian intersecting the J2000 plane as another linear function. Often, some small trig terms are added to those linear functions.
At present, COSPAR.CPP includes many, but not all, of the satellites shown by Guide. Some (such as the "captured asteroids" of gas giants) have no defined orientations. Others (small inner objects such as Amalthea) just haven't been tackled by me yet (no big job to do, but not much interest, either).
DATE.CPP: Source code for the following calendrical systems: Julian, Gregorian, Hebrew, Persian, Islamic, and Chinese, and French Republican. Basically, you can do conversions between JD and the calendar systems.
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.
DE_PLAN.CPP: Code to compute very precise planetary positions using the PS-1996 series of Chapront. Includes Mercury through Pluto. To use this, you'll need to download this file of coefficients.
DIST_PA.CPP: Code to compute the distance and position angle between two spherical coordinates. Theoretically, this ought to be a straightforward task. In practice, there are an amazing number of ways to get roundoff/truncation errors and divides by zero, mostly involving how you deal with cases where the distance between the two points is small. I am reasonably sure now that this code evades all such cases.
EART2000.CPP: Code to compute the Earth's location relative to the Sun, in J2000 coordinates, using VSOP.
EASTER.CPP: Code to figure out the date of Easter for a given (Gregorian) year.
ELP82DAT.CPP: Code to compute the position of the moon using the ELP-82 (Ephemerides Lunaire Parisienne 1982) theory of Chapront and Chapront-Touzé. To make use of this code, you'll need to download this file of coefficients.
GETPLANE.CPP: Code to compute the location of any planet, or the Moon, relative to the Sun. The resulting vector is given in five different systems.
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.)
JSATS.CPP: Code to compute the positions of the Galilean satellites, using Lieske's theory as given in Meeus' Astronomical Algorithms. Quite a few smaller terms are omitted.
LUNAR2.CPP: Code to compute the position of the moon, using the truncated series from Meeus' Astronomical Algorithms. This uses the same data file (VSOP.BIN) as the VSOPSON.CPP code discussed below.
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.)
NUTATION.CPP: Computes the nutation of the earth's equator.
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.
PLUTO.CPP: Code to compute the position of Pluto, using the series given in Meeus' Astronomical Algorithms. This uses the same data file (VSOP.BIN) as the VSOPSON.CPP code discussed below.
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.
REFRACT.CPP: Contains functions to convert an observed altitude (one affected by refraction) to a true altitude (one that would be seen on an airless planet), or vice versa. Two different methods are provided. (See next paragraph for a third method.)
REFRACT4.CPP: Contains a more complex/sophisticated method of computing refraction, as given in the Explanatory Supplement to the Astronomical Almanac . This one does a numerical integration through the atmosphere, including effects due to observer altitude, temperature, barometric pressure, humidity, and wavelength of light being refracted. (But the end result is usually not all that different from that given by the two methods in the previous paragraph, except for very low altitudes or extreme conditions.)
RELATIVI.CPP: Code to test out the precession of Mercury's orbit due to general relativity, and also to test a quick and easy way to include first-order general relativity in orbit computations. The code may also be of interest because it integrates using the method of Encke, with a pretty good Runge-Kutta numerical integrator. Unfortunately, it could stand some documenting...
ROCKS.CPP: Code to compute positions for 25 faint, inner satellites of the gas giants. These are all modelled as precessing ellipses.
SSATS.CPP: Code to compute the positions of eight of Saturn's satellites (Mimas, Enceladus, Tethys, Dione, Rhea, Titan, Hyperion, Japetus), using the theory due to G. Dourneau.
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".
VSOPSON.CPP: Code to compute planetary positions from a truncated VSOP series. That series is stored in binary form on the Guide CD-ROM, and you have to load it into a buffer to use this function. The logic is very similar to that used for the full VSOP (see BIGVSOP.CPP), except that that function works by reading in the file of coefficients on an "as needed" basis. This function (and several others listed on this page) require that you download this file of coefficient data.
Tycho-2 CD access: Given the delays in getting Guide 8 out the door, I've had requests to provide some means of accessing the Tycho-2 CD-ROMs. (These disks have the additional advantage of being free, provided by ESA and USNO. Click here for details. Though I wouldn't be surprised if they're running low on disks.)
While I haven't done much with this data (except process it for use in Guide 8, of course), I have written a small utility to access Tycho-2 data for a given star. Click here to download the DOS software, source code, and documentation. Given this, you can extract positions, proper motion, and magnitudes (both Tycho and Johnson), with error estimates, for a given Tycho-2 star. (Both the transformation to the Johnson system and the error estimates are important. "Raw" Tycho mags match Johnson magnitudes only in a loose manner, and the Tycho-2 magnitude data runs the full range from extremely accurate to moderately accurate... you _must_ check error estimates to determine which is which.)
Artificial satellites: This has grown to the point where it deserved a page of its own. Click here for details on C/C++ code to compute artificial satellite ephemerides. (The currently-posted version is a nice improvement over some older artificial satellite code I posted previously. I'm leaving the old stuff because it's still getting a little use here and there.)
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 David'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 (or decreased) over billions of years, we may assume that the whole thing would mysteriously explode with stars flying everywhere (or collapse down to a black hole). 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) uses a method from
PAGE2.CPP: This code
can take an RA/dec position and figure out the corresponding pages to
be used in the Millennium Star Atlas, Uranometria,
and Sky Atlas 2000. There's also a little snippet that takes
a (lunar) lat/lon and figures out the corresponding page to be used
from Antonín Rükl's
Atlas of the
Moon. Source code for Charon: I've started the process of turning
the Charon astrometry software into an open-source
program. This begins by posting the source code
for Charon (about 84 KBytes). You will also need the
basic astronomical source code discussed
above. Why the source for Charon is being posted:
I do not now, nor have I ever, regarded astrometry software as a
commercially profitable idea. By a slight margin, I'm able to keep
myself fed and indoors through sales of Guide; trying to sell Charon
is not apt to help very much. By posting it, I hope to accomplish the
following: -- People will be able to improve/modify the source code to accomplish
tasks of their own. (The basic idea behind almost all open source,
of course.) With any luck, some of that code will return to me and
I'll post it here for others to use. I know that some people have
written Charon-like code and functions for their own project; I'm
hoping that we can stop re-inventing wheels. -- The number of users of Charon has risen to the point where I'm
not doing too well at keeping up with requests for improvement. I've
become a bottleneck to progress. -- Posting the source code will probably encourage somebody to start
work on making a version of Charon that does not require a Guide CD.
If it can access catalog data from the Internet, then it would require
no CD-ROM at all. In this case, people could do astrometry with
Totally Free Software. This might encourage some people to enter
the hobby. -- Verifying the fact that the code does what I tell people it does
becomes easier. Some people build their own astrometry software
because they would rather not have to place trust in me and my
software. Such trust would be well-placed, especially since the
results of the software have been examined by dozens if not hundreds of
people. (The worst risk in 'rolling your own' astrometry software is
that it will usually be tested only by you, and not by dozens
if not hundreds of people.) But being able to check the source code
ought to relieve some anxieties. Current state of the Charon source code:
All of the rest of the source posted on this site has been cleaned up
nicely and documented well enough to make it useful to somebody. Right
now, this is not in the least true of Charon. You can use the source
to satisfy curiosity about parts of Charon, and maybe to grab bits
and pieces from the code. But there are some things I need to do to
make it usable. Charon is compiled with WATCOM C/C++, the only compiler I've got
that can produce 32-bit DOS code that can access graphics. (Micro$oft's
compiler can produce 32-bit DOS console applications, but none that
access graphics.) I'm in the process of getting the make file into a
usable form; right now, it runs across assorted directories and is not
'ready to go out of the box'. There are parts of the code that run for hundreds of lines without
the faintest trace of the residue of a comment. The code isn't
particularly ugly (except for the enormous main() and parts
of matcher.cpp), but figuring out what does what is a
frightening prospect right now. Overall structure of Charon:
The program is overdue to be broken up into certain components. To
a minor extent, this breakup is present now, but it's incomplete and
Charon is really one big monolith at heart. This is really a shame,
because astrometry software lends itself naturally to being broken
into a set of component tasks. Breaking Charon into components would make following the structure of
the code much easier. It would let us rip out certain parts,
replace them, or use them by themselves. The components (or, if you
prefer, "near-components") are: Load an image from a file:
A piece that loads an image, be it FITS, SBIG, or other, into a
standardized in-memory format. This is currently done using the
load_image() function in load_img.cpp. A side note:
it's possible that, in preference to writing one function that can load
a zillion image formats, there should be a separate routine that can
convert a zillion image formats into FITS. If I had it to do over
again, I'd do it that way. Find and build a list of stars in an image:
Takes an image in the aforementioned standardized in-memory format and
return a list of stars in it, with their pixel coordinates and
brightnesses. The code in findstar.cpp does this right now.
(I'm not totally thrilled with the job it does on noisy images; it
often finds hordes of spurious stars in noise. Ideally, the code I
wrote might get pitched in favor of SExtractor. I started
playing with SExtractor and gave up, though; I had no real success in
puzzling out what was supposed to do what.) Get the RA/dec of an object at a given time:
Given an object ID such as "V4191 Sgr" or "P/Linear S4" or "1997 XF11",
and a date/time and observer position, this figures out the RA/dec of
that object. The code in find_obj.cpp does this right now.
If somebody tells Charon, for example, "Here's an image of asteroid
4179", Charon can figure out a rough RA/dec for 4179 at the time the
image was taken, so it knows where to look for a pattern match. In some ways, this is just a convenience. People can and do feed
Charon images and just say, "Here's an image taken at about 14h51m12s,
+10 43' 56"... go process it." But the 'convenience' is a major one. Get stars from a catalog within a given radius of a given RA/dec:
Given a catalog name and an RA/dec (normally provided by the preceding
piece, but not always) and a radius, this piece grabs catalog data
for that area and returns it in a standard way. Right now, assorted
functions in grab_gsc.cpp do that. They can extract ACT,
GSC, or GSC-ACT data from the Guide disks, or USNO Ax.0 or SAx.0
data from the CDs distributed by USNO, or GSC 1.2 data from the files
that come from the STScI server. All the functions feed through
the grab_catalog_stars() function in charon.cpp.
For Guide users, this is pretty nice; it means they can get all the
data they need (comet/asteroid orbits, plus the GSC, plus other handy
catalogs) all on one disk. For an open-source project, though, it'll
be nice if the 'catalog grabbing' code can 'grab' from the two GSC CD-ROMs
distributed by STScI, or via Internet. The first, at least, ought to
be very straightforward. The second will be slightly trickier, but
opens the program to anyone connected to the Internet... this would
be a huge plus. I note that VizieR, for example, provides
a way to ask for a given list of data from a catalog, specified by
RA/dec and radius, using a single URL... almost exactly what we'd
need for a stunt like this. Right now, people have to shell out some money to buy the Guide
disk, or the two GSC disks, to get into astrometry. Neither is
really pricey, as astronomy products go ($89 and $69, respectively).
But if trying out astrometry software was totally free, it might draw
in some fresh blood. Pattern-match the list of image stars to the list of catalog stars:
Given a list of image stars (probably from the
'get stars from image' component), and
a list of catalog stars (probably from the
'get stars from catalog' component), this part does the automated
pattern-matching required to say, "Use the following transformation to
go from image to RA/dec space and vice versa." The code for this is
currently in matcher.cpp, and is a horror to behold... it's the
sort of neighborhood where programmers only go in pairs.
Display the results and interact with user:
This piece takes the image, lists of image stars and catalog
stars, and transformation, and display the whole mess on the device
of your choice. At present, this is done partly with code in the
main( ) of charon.cpp, with a big chunk of it
in dispimg.cpp and smaller pieces throughout the code.
When applied to multiple images, this could allow for blinking.
(The current blinking scheme in Charon is considerably less logical.
I'd ignore it if possible.) This piece is intimately hooked into
the entire user interface problem of how people can click on the
image and manipulate it on-screen.
The previous components could be written in totally
'clean', generic ANSI source. Porting them to other OSes would be
relatively straightforward. Unfortunately, the 'display/interact'
component will be messy, heavily OS-dependent code. Assorted minor functions:
A whole slew of small pieces to do things such as write out MPC headers
and reports, figure out which image star or catalog star is nearest to
a given cursor position, provide ways to adjust settings, handle the
list of MPC stations... (FUTURE, not yet implemented) Automated object detection:
Once the pattern matching code has done
its work for a particular image, we can get a list of the
RA/dec/magnitude values in that image. (Such lists would have value in
and of themselves, of course.) Then we can take two or three such lists
at different instants, and look for moving objects.