Mathcad Worksheets by Astroger

(An Electronic Edition of Topics in Astrodynamics is now available on CD-ROM.)

After many years of developing astronomy & astrodynamics applications in FORTRAN, BASIC, Pascal, and C, I have found a faster and easier way to implement and document the relevant mathematics, using a single vehicle.

The vehicle I am referring to is, of course, Mathcad, originally developed and formerly published by Mathsoft, Inc., of Cambridge, Massachusetts, but now published by Parametric Technology Corporation (PTC) of Needham, Massachusetts (both Cambridge and Needham are in the greater Boston area).

Mathcad does not require students and reviewers of my work to learn a conventional, ASCII-text based programming language, nor to read through a separate mathematical document in order to extract the mathematical content of my software applications. But you do need Mathcad 2001 Professional, or later, to view and work with the Mathcad worksheets that I describe below. (Mathcad 2001 was succeeded by Mathcad 2001i, Mathcad 11, Mathcad 12, Mathcad 13, Mathcad 14, and Mathcad 15.)

How to Purchase Mathcad

Mathcad 15 introduced a new suite of functions for Design of Experiments, and is currently available for purchase directly from PTC. Click on Mathcad Store to order. The purchase price includes a node-locked license and a prepaid year of annual maintenance. I strongly urge past users of Mathcad to upgrade to Mathcad 15, and to suppport the product by renewing the annual maintenance contract each year (I just renewed my own annual maintenance contract).

Annual maintenance entitles the purchaser to download each new maintenance release that occurs during the year. Also, it entitles the user to have a home copy of Mathcad as well as a work copy. If you choose not to renew your annual maintenance contract, your work copy will always work ("in perpetuity"), but your home copy will stop working when your annual maintenance contract expires.

I should note that PTC just introduced on January 10, 2011, a wholly new product called Mathcad Prime 1.0. Mathcad Prime 1.0 comes bundled with the purchase of Mathcad 15, and Mathcad 15 comes bundled with the purchase of Mathcad Prime.

Mathcad Prime is the future of Mathcad. It is PTC's stated goal that existing Mathcad 15 users be able to transition fully to Mathcad Prime by Mathcad Prime 3.0.

In other words, Mathcad Prime 1.0 does not now have all of the capabilities of Mathcad 15. But by releasing it now at Version 1.0, PTC is giving existing Mathcad 15 users the opportunity to work with Prime and to have input to the ongoing development process by joining the Mathcad community at PlanetPTC.

Special Note on Mathcad 14

There may still be at this writing (March 2011) a way to obtain Mathcad 14 for academic use (student or academic faculty/staff). Brent Maxfield's paperback book, Essential Mathcad for Engineering, Science, and Math (2nd. Edition, 2009) bundles a non-expiring, single-user copy of Mathcad 14 on CD-ROM.

The single-user copy of Mathcad bundled with Essential Mathcad is "for educational use only." Commercial use is expressly prohibited, and violates the educational user license agreement. Click on Essential Mathcad to see how to purchase, at a 30% discount, the book-plus-CD-ROM directly from Elsevier. (Note: this offer was valid as of February 2010, but may be withdrawn by Elsevier at any time.)

If you are a student or academic faculty/staff, you may also go to your university bookstore to purchase Mathcad. If you are a student, DO IT NOW while you are still eligible for the academic price. The academic edition consists of the full professional edition, and may include eligibility to receive free electronic books that would cost hundreds of dollars more if purchased separately.

Some Notes on Mathcad Versions 11 and Later

Mathcad 14 had enhanced numerical integrators and some other new features. When it was initially released, and for a long time afterward, it was plagued by the "print preview" crash. PTC's Mathcad 14.0 service release M020 fixes the print preview crash and is currently available for download from PTC. Click Mathcad 14 M020 for further information.

Speaking as a Mathcad 14 beta tester, the two biggest surprises in Mathcad 14 were the new symbolics engine, based upon MuPAD rather than Maple, and the adoption of Unicode. Since I don't use symbolics much, the change from Maple to MuPAD had little effect on my legacy worksheets. But the adoption of Unicode had a serious impact on my prior work, because whenever I edit the text of a legacy worksheet in Mathcad 14, the edited Greek letters don't look the same as the ones that were left unchanged. And the edited Greek letters just look bad.

Mathcad 11, if you are lucky enough to still have it, works with Windows 98SE, Me, NT4, NT5, and 2000, as well as with Windows XP. It is the last "Win32" version of Mathcad, for Mathcad 12 is a Microsoft .NET Framework-based application. Your computer must have Windows 2000 SP4 or Windows XP SP2, at least, to run Mathcad 12, 13, and 14.

Starting with Mathcad 12, worksheets are saved in XML format by default. But you can still save in prior-version formats back to Mathcad 2001 Professional. The new XML format makes the worksheet files externally searchable and traceable. This is a "calculation management" enhancement that is particularly valuable to large-enterprise users of Mathcad.

Mathcad 12 was in some ways a disappointment to this single-user licensee. E.g., it ran slower than Mathcad 11 and didn't appear to be as stable. The move to .NET Framework and XML worksheet format seemed to be at the expense of speed and stability.

But based upon my experience as a beta tester for Mathcad 13 and Mathcad 12, I can say now that Mathcad 13 was a dramatic improvement over Mathcad 12. With my orbital mechanics benchmark, a 498-observation batch filter for orbit determination that I have run with every Mathcad version since Mathcad 7 Professional, Mathcad 13 is about three times faster than Mathcad 11. Also, Mathcad 13 appears to be much more stable than Mathcad 12.

Noteworthy new features of Mathcad 13 were: trace debugging of functions; worksheet autosave; improved genfit and linear algebra functions; 2D plot trace enhancement; documentation enhancements (e.g., new programming and migration guide tutorials; new quicksheets on waterfall plots and ODE solving).

How to Obtain My Worksheets

To date, I have contributed twelve worksheets to Mathsoft's Mathcad Library. They are:

1. Ephemeris of a Comet Via Uniform Path Mechanics (UPM)

2. Herget's Method for an Asteroid

3. Orbit Propagation via State Space Analysis

4. Effect of a Radial Impulse on a Circular Orbit

5. Herget's Method with Cassini's Earth Flyby

6. Sun Altitudes for Sextant Practice

7. Sun-Sight Solutions Without Tables

8. Rectilinear Two-Body Motion ("Earth Falls Into the Sun")

9. Gauss's Angles-Only Method with "Killer Asteroid"

10. Tracking Data Reduction for Galileo's Earth 1 Flyby

11. Calculating the Photoperiod in Plant Physiology

12. Modeling Blackbody Radiation

All twelve worksheets are currently available for downloading from the Mathcad Resource Center. But they have not been kept up to date, e.g., some of the earlier worksheets employ features that were deprecated and then discontinued in subsequent versions of Mathcad.

Here is the link to the Mathcad Resource Center. Note that the worksheet titles in the Mathcad Resource Center are listed in reverse chronological order, whereas the worksheet titles above are listed in strict chronological order.

E.g., you will find that the last worksheet title listed above, "Modeling Blackbody Radiation," is given on the first page in the Mathcad Resource Center with the title "Blackbody Radiation."

If you download the Blackbody Radiation worksheet, be sure to read the comments about it given below: there is a typo in the value of c, the speed of light, that needs to be corrected in the worksheet. It should be assigned as c := 2.99792458 x 10^10 cm/sec in the worksheet. That is, the second place to the right of the decimal should be a 9, not a 7.

All twelve worksheets, with corrections as needed, and verified to work with Mathcad 11-15, are now available in both .mcd and .pdf format on either (a) the CD-ROM that contains the Electronic Edition of Topics in Astrodynamics, or (b) the supplementary CD-ROM that accompanies the print edition of the book. Click Topics in Astrodynamics for further information.

Below I provide additional details and background on my contributed worksheets. I do so for the benefit of my present and former students, and for the benefit of interested colleagues in the fields of astronomy and astrodynamics.

Finally, in Author's Background, at the very end of this webpage, I give a brief outline of my space career and my related interests for the reader who wishes to know more about my professional background and credentials.

1. Ephemeris of a Comet Via Uniform Path Mechanics (UPM)

This worksheet uses the orbit of Comet Hale-Bopp to illustrate UPM, a mathematical procedure for calculating the orbital path of a space object. This procedure has the same form (i.e., is "uniform") for all possible path eccentricities.

  • More about Ephemeris of a Comet

  • The latest version of the "Ephemeris of a Comet" worksheet, with updates and corrections, is now bundled with the Electronic Edition of Topics in Astrodynamics.

    Back | Top

    2. Herget's Method for an Asteroid

    The two worksheets, (the "initiator") and (the "iterator") use astrometric observations of the orbit of the critical-list minor planet 1035 Amata to illustrate Herget's Method of computing a preliminary orbit. The observations were taken by USAF 2Lt Dan Burtz.

  • More about Herget's Method for an Asteroid

  • Determining the Orbits of Comets and Asteroids

  • The latest versions of the two "Herget's Method for an Asteroid" worksheets, with updates and corrections, are now bundled with the Electronic Edition of Topics in Astrodynamics.

    Regarding Herget's method, Project Pluto provides a Windows application called FIND_ORB that lets you "find" an orbit from the astrometric observations of a minor planet or comet. FIND_ORB uses Herget's method to solve for the preliminary orbit, but its capabilities do not stop there. It can fit a more precise orbit to the observations using the method of batch least squares differential correction. You can specify which planets are presumed to be "perturbers" of the comet's or asteroid's orbit as fitted to the observations. To download FIND_ORB, go to the Project Pluto website at

    Back | Top

    3. Orbit Propagation via State Space Analysis

    This worksheet uses the Earth escape trajectory of the Near-Earth Asteroid Rendezvous (NEAR) spacecraft, following its launch on 1996 February 17, to illustrate how state space analysis can be applied to the problem of propagating a trajectory through space and time.

  • More about Orbit Propagation via State Space Analysis

  • The latest version of the "Orbit Propagation" worksheet, with updates and corrections, is now bundled with the new Electronic Edition of Topics in Astrodynamics.

    Back | Top

    4. Effect of a Radial Impulse on a Circular Orbit

    This worksheet considers the effect of a 1.0 km/sec radial impulse on the orbit of an artificial Earth satellite at geostationary altitude (35786 km), both outward and inward.

  • More about Effect of a Radial Impulse

  • The latest version of the "Radial Impulse" worksheet, with updates and corrections, is now bundled with the new Electronic Edition of Topics in Astrodynamics.

    Back | Top

    5. Herget's Method with Cassini's Earth Flyby

    This set of worksheets applies Herget's method to determine the flyby path of the Cassini/Huygens space probe, launched 1997 October 15, following its return to Earth for a gravity assist on 1999 August 18. You will need to download, the "initiator" worksheet, and its associated input files observations.prn and sensors.prn, plus the "iterator" worksheet

    Once you have downloaded all four files, use Mathcad 8 or later to open the worksheet and click once on "Calculate Worksheet" in the Mathcad Math menu, which will cause to read the Cassini observations from observations.prn and the observers' coordinates from sensors.prn. Now open the worksheet and click once on "Calculate Worksheet", then scroll down to the RMS history area.

    Click three more times on "Calculate Worksheet", and watch as converges to the solution. Scroll around the worksheet to see all of what is involved: scroll upward from the RMS history to see the orbital elements solution; scroll downward to see the details of the residuals and convergence information. The worksheet is 15 pages long; be sure to scroll to the very end to see the references and summary analysis.

  • More about Herget's Method with Cassini's Earth Flyby

  • The latest versions of the "Herget's Method with Cassini's Earth Flyby" worksheets, with updates and corrections, are now bundled with the new Electronic Edition of Topics in Astrodynamics.

    Back | Top

    6. Sun Altitudes for Sextant Practice

    This worksheet provides an algorithm for predicting the sun's altitude and azimuth with respect to an Earth-fixed observer. Given your latitude and longitude, and an input file of local times of interest on a date of interest, will calculate for you the sun's altitude and azimuth at those times.

    Thus, if you take a sequence of sun altitude measurements with a sextant, the worksheet will compute for you what results you should get, as a check on your skill in making "sun shots".

    You need not have an interest in celestial navigation to find this worksheet useful. Want to know the sun's altitude, and its azimuth direction from your house, at the beginning of each daylight hour on some date of interest? The worksheet will allow you to calculate these data quickly and accurately. All you have to do is to (1) use, say, Windows Notepad to set up the local times of interest in the file times.prn, (2) change the year, month, and day in the worksheet itself, and (3) also change the latitude, longitude, height, and time zone offset in the worksheet. Then click on "Calculate Worksheet" in Mathcad's Math menu and scroll down to the page containing the worksheet's results.

  • More about Sun Altitudes for Sextant Practice

  • The latest version of the "Sun Altitudes" worksheet, with updates and corrections, is now bundled with the new Electronic Edition of Topics in Astrodynamics.

    Back | Top

    7. Sun-Sight Solutions Without Tables

    This worksheet presents a novel two-sight fix method and applies it to the angles.prn data in Richard R. Shiffman's "Sextant Noon-Day Sun Sightings" (Shiffman's worksheet can also be found in the Mathcad Library). It references the worksheet*, so you need to have also and its two input files, times.prn and angles.prn, in order for this worksheet, to work.

    (*If the reference in the very first page of the worksheet isn't working, go to Mathcad's Insert menu and click on "Reference", type in the name "", click on the "use relative reference" box, and then click "O.K". Now go to Mathcad's Math menu and click on "Calculate Worksheet". The reference to in the worksheet should now work.)

    Because references, it is able to take advantage of all of the procedural functions defined there: functions for calculating solar position in the ECI reference frame, referred to the mean equator and equinox of J2000.0; aberration (light-time correction); precession; nutation; ECI-to-topocentric transformation; refraction, etc. Thus the navigation solution is effected without reference to the Nautical Almanac, as well as without reference to sextant sight reduction tables.

  • More about Sun-Sight Solutions Without Tables

  • The latest version of the "Sun Sight Solutions" worksheet, with updates and corrections, is now bundled with the new Electronic Edition of Topics in Astrodynamics.

    Back | Top

    8. Rectilinear Two-Body Motion ("Earth Falls Into the Sun")

    If Earth were to stop in its orbit and fall into the sun, for how many days would it fall? Greg Neill, editor and publisher of The Orrery, posed and solved this problem in issue #43 (December 2001). Upon reading Greg's solution, I was reminded of a paper I presented at a conference of the American Institute of Aeronautics and Astronautics (AIAA) in 1986. This worksheet solves the problem of "how long would it take for Earth to fall into the sun" by means of the theory that I presented in that 1986 AIAA paper.

  • The latest version of the "Rectilinear Two-Body Motion" worksheet, with updates and corrections, is now bundled with the new Electronic Edition of Topics in Astrodynamics.

    Back | Top

    9. Gauss's Angles-Only Method with "Killer Asteroid"

    Asteroid 1997 XF11 received much attention in March 1998 when astronomers determined that it might hit Earth in October 2028. After the alarm was sounded, it was a matter of just one day before new information came in and it could be determined that a direct hit was not really very likely after all. In this worksheet I determine an orbital solution, via Gauss's angles-only method of orbit determination, using three observations taken in December 1997. In the worksheet I compare my solution with the definitive solution published by the Minor Planet Center in Cambridge, Massachusetts.

  • The latest version of the "Gauss's Angles-Only" worksheet, with updates and corrections, is now bundled with the new Electronic Edition of Topics in Astrodynamics.

    Back | Top

    10. Tracking Data Reduction for Galileo's Earth 1 Flyby

    The Galileo mission to Jupiter was one of the most exciting and successful endeavors in the history of robotic planetary exploration. Launched on 1989 October 18, the Galileo spacecraft orbited the Sun for more than six years, receiving one gravity assist from Venus and two from Earth, before achieving insertion into jovicentric orbit on 1995 December 7. After almost eight years of orbiting Jupiter and exploring the planet, its rings, and its moons, the Galileo spacecraft flew into Jupiter on command, and burned up in Jupiter's atmosphere, on 2003 September 21. Thus ended the Galileo mission.

    I was privileged to be able to assist in the near-realtime reduction of tracking data from Galileo's Earth 1 flyby on 1990 December 8. I published a related article in 1993 ["Algorithms for Reducing Radar Observations of a Hyperbolic Near-Earth Flyby," Journal of the Astronautical Sciences (April-June 1993), pp. 249-259]. The two worksheets in this Mathcad application, and, rigorously solve for the flyby trajectory of the Galileo spacecraft. They agree perfectly with the results I presented in my 1993 paper. (I did the calculations in the paper in 1992 by means of a Turbo Pascal 1.1 program that I developed and ran on an accelerated Apple Macintosh SE.)

  • The latest versions of the "Galileo's Earth 1 Flyby" worksheets and their associated data files, with updates and corrections, are now bundled with the new Electronic Edition of Topics in Astrodynamics.

    The tracking data are reduced by opening worksheet and clicking on "Calculate Worksheet". Then open and click on "Calculate Worksheet". Scroll down to page 14 of the worksheet to see the RMS error matrix. The RMS error for the first iteration of the differential correction (DC) should be 5.043 km. Now click a second time and watch the RMS error go to 4.584 km on the second iteration. Note that the DC has converged, on the second iteration, to the orbital solution that can be seen by scrolling down to pages 17 and 18. The results on page 18 show that the Galileo spacecraft flew to within about 961 km of Earth's surface at its Earth 1 flyby.

    Back | Top

    11. Calculating the Photoperiod in Plant Physiology

    What does calculating the photoperiod have to do with astronomy? The answer is this: the photoperiod is the duration of daylight each day, i.e., the time of sunset minus the time of sunrise. It is a function of both the geographic latitude (or the astronomical declination) and the season, expressed as a count of days since the beginning of the year.

    This contribution to the Mathcad Library consists of two worksheets, as follows.

    1. Calculating Photoperiod - derives simple formulas for the photoperiod, in minutes, and its diurnal rate, in minutes per day, as functions of declination and days since the beginning of the year.

    See also the chapter on photoperiodism in Plant Physiology, a textbook by Frank B. Salisbury and Cleon W. Ross (Wadsworth, 1992). The chapter on photoperiodism contains graphs prepared by Michael J. Salisbury circa 1983, using equations and data provided by Astroger back then. The "Calculating Photoperiod" worksheet generates these same graphs using new, analytical (vs. numerical) photoperiod and photoperiod rate formulas as derived in the worksheet.

    2. Photoperiod performs rigorous astronomical calculations to provide and plot the photoperiod, as a function of latitude and day count, for various latitudes. This worksheet was prepared in December 2003, in response to an e-mail request from Dr. Lilian B. P. Zaidan, plant physiologist at the Instituto de Botanica, Sao Paulo, Brazil.

  • The latest versions of the "Photoperiod" worksheets, with updates and corrections, are now bundled with the new Electronic Edition of Topics in Astrodynamics.

    Back | Top

    12. Modeling Blackbody Radiation

    2005 is the "World Year of Physics" and it marks the 100th anniversary of Albert Einstein's publication of papers on the photoelectric effect, Brownian motion, and the special theory of relativity. All four papers on these three topics appeared in the year 1905. But modern physics, i.e., modern quantum physics really began in 1901, when Max Planck propounded the notion that material bodies, but especially blackbodies*, emit and absorb thermal radiation in discrete quanta of energy, rather than continuously.

    Planck's hypothesis of quantized absorption and emission of radiation made it possible for him to derive a radiation law that applies to blackbody emission at all wavelengths and all frequencies, a universal law that succeeds in spectral regions where the prior radiation laws of Rayleigh, Jeans and Wien had failed. Planck received the Nobel Prize in physics in 1918 for his quantum theory of radiation.

    The Mathcad 12 worksheet, "Modeling Blackbody Radiation," revisits how Max Planck integrated the blackbody radiation curve for an arbitrary Kelvin temperature, T, over all possible wavelengths of thermal emission, to arrive at the Stefan-Boltzmann law. The Maple symbolic processing capability of Mathcad is invoked at key points of the derivation and Bernoulli numbers are used to evaluate the infinite series that is crucial to the derivation. Finally, Mathcad's X-Y Plot capability is used to plot the blackbody radiation curve for 2.725 degrees Kelvin.

  • More about Modeling Blackbody Radiation

  • The latest version of the "Blackbody Radiation" worksheet, with updates and corrections, is now bundled with the new Electronic Edition of Topics in Astrodynamics.

    (Note: there was a typographical error in the worksheet, in the value assigned to the speed of light, that has now been corrected. In case you downloaded the worksheet before the typo was corrected: the speed of light should have been assigned as c := 2.99792458 x 10^10 cm/sec in the worksheet. The second place to the right of the decimal should be a 9, not a 7.)

    *A blackbody is an ideal body in thermal equilibrium that absorbs all incident radiation and re-emits it as light energy distributed over the entire electromagnetic spectrum.

    Back | Top

    More about "Ephemeris of a Comet"

    The two "Ephemeris of a Comet" worksheets, (for Windows users) and (for Power Mac users) are Mathcad PLUS 6 worksheets that let you generate the Earth-relative sky coordinates of any celestial object, provided that you have orbital elements for the object. The worksheets first became available for downloading from MathSoft's website in September 1997.

    (Note: the worksheet is no longer available for download from the Mathcad Library. If you are a Power Mac user and have Mathcad PLUS 6 for Power Macintosh, e-mail me at the address at the very end of this web document and I'll send back to you as an e-mail attachment.)

    The worksheets use Comet Hale-Bopp (C/1995 O1) for their example, and they generate ephemeris points at 5-day intervals during the period 1997 March 17 - 1997 May 16. If you have Mathcad, PLUS 6 or some later Professional version (e.g., 7, 8 or 2000), you can use the live worksheet to generate an ephemeris for any comet that you have orbital elements for, for any time period of interest. (Actually, Mathcad Explorer will even let you modify the worksheets to run for your own test case, but you cannot save your modifications unless you have Mathcad PLUS 6 or later.)

    The UPM theory is the result of the author's research into universal variables methods of orbit propagation, for application to gravity-assist flybys (swingbys) of Earth by manmade spacecraft such as Galileo, NEAR, and Cassini. It is intended to be used with orbits having eccentricity 0.95 or greater. Such orbits are typically Earth escape and flyby or swingby trajectories.

    The full UPM theory is implemented in the tracking data reduction software of U.S. Space Command's Space Defense Operations Center (SPADOC) near Colorado Springs. It was used to reduce tracking data on the Galileo spacecraft's Earth 2 flyby trajectory (1992 December 8) and on the Near-Earth Asteroid Rendezvous (NEAR) spacecraft's Earth escape trajectory (1996 February 17). Other recent or upcoming applications of the UPM theory are the NEAR spacecraft's swingby of Earth on 1998 January 23 and the Cassini/Huygens spacecraft's gravity-assist flyby of Earth expected on 1999 August 18.

    For more on the full UPM theory, see "Algorithms for Reducing Radar Observations of a Hyperbolic Near-Earth Flyby," Journal of the Astronautical Sciences (April-June 1993).

    Back | Top

    More About "Herget's Method for an Asteroid"

    Paul Herget, director of the Cincinnati Observatory from 1943 to 1978, and first director of the IAU's Minor Planet Center (1947-1978), published in 1965 an ingenious method for determining a preliminary orbit for a comet or asteroid, from all available angles-only observations -- not from just three angles-only observations, as is typical of other preliminary orbit determination methods. Herget's method lends itself perfectly to the task of determining comet and asteroid orbits from observations taken via modern CCD ("charge-coupled device") astrometry.

    In Herget's time, images of comets and asteroids were obtained via photographic astrometry. Photographic astrometry employs a telescope with a camera that positions a glass photographic plate in the focal plane to record an image of a comet or asteroid against the star background. This method, which still finds use today for high-resolution work, is a tedious, time consuming, and expensive process. CCD astrometry facilitates the taking of multiple images of the same object (e.g., comet or asteroid) on multiple nights via CCD camera, telescope, and computer. The CCD camera replaces the photographic plate in the focal plane with a CCD array that does not have to be replaced after an image is recorded, nor does it require post-collection chemical development of the image, as does a photographic plate.

    How the CCD imaging process works is that the CCD camera uses the charge-coupled device array in the focal plane to record images that are stored on the computer as digital image files. Desktop computer programs such as Astrometrica and Charon permit the reduction of multiple CCD image files, collected over successive nights, to topocentric right ascension and declination measurements. Herget devised his method expressly to work with these topocentric right ascension and declination measurements, and as many of them as are available.

    Herget's method works by assuming that the direction lines -- from observer to asteroid -- for the first and last observations (obs) are exact, and by improving initial estimates of the observer-to-asteroid distances for the first and last obs. The two estimated distances are fitted to the remaining observations by treating their (the remaining observations') residuals as functions of the two estimated distances. The partials of the residuals with respect to the two estimated distances are calculated by numerical differencing and the fitting process is non-linear least squares.

    J.M.A. Danby gave Herget's method a particularly lucid exposition in the second edition of his textbook, Fundamentals of Celestial Mechanics (Willmann-Bell, 1988, pp. 235-238). This book is also highly recommended for its seminal exposition of Karl J. Stumpff's universal variables theory as based upon Stumpff's c-functions.

    Project Pluto has implemented Herget's method in the freely downloadable FIND_ORB (16-bit Windows, introduced in March 1997) and FIND_O32 (32-bit Windows, introduced in July 1998). FIND_ORB solves for a preliminary orbit by Herget's method (by taking "Herget steps") and optionally improves the preliminary orbit by a full, six-element differential correction, taking into account gravitational perturbations by the major planets and Earth's moon (each iteration of the differential correction is a "Full step").

    Don't limit your study of Herget's method to my worksheets. Get a copy of Danby's book from Willmann-Bell at and download FIND_ORB or FIND_O32 from the Project Pluto website.

    I should mention that I found out about FIND_ORB when it was described by Stuart J. Goldman in the June 1998 "Software Showcase" department of Sky and Telescope magazine. Note that while FIND_ORB is an easy-to-use program with a simple Windows interface, its "engine" is comprised of computer codes of great astronomical and mathematical sophistication.

    Project Pluto publishes the highly successful Guide and Charon programs for CCD astrometry, as advertised in Sky and Telescope (see the January 1999 issue, p. 139). Project Pluto is to be commended for making the FIND_ORB program available to the astronomical community at no charge via the Internet.

    Mathcad worksheets for "Herget's Method for an Asteroid", and, should be available for downloading from MathSoft's website sometime in January 1999. These two worksheets

  • Employ the UPM theory for path propagation, and
  • Use the UPM-extended Gaussian sector-to-triangle area ratio iteration.

    (The UPM-extended Gaussian iteration assumes that Gauss's hypergeometric X-function is a "quotient of c-functions", while the original Gaussian iteration assumes that Gauss's hypergeometric X-function is a "hypergeometric series".)

    The first worksheet,, is an "initiator", in that it sets up the test case input for Thus is run before, and its window should be kept open in the background when the window is opened. is the "iterator" in that each time you click on "Calculate Worksheet" in the Mathcad Math menu, hmc performs an iteration of Herget's method. To see this happen, scroll down to page 11, where the RMS matrix is accumulated, and watch the behavior of the RMS (root mean square of the residuals) as you iterate. For the test case supplied via hm1, hmc should converge on the eighth iteration to an RMS of 0.209 kilometers.

    I note in closing that although "Herget's Method for an Asteroid" takes up 23 printed pages, total, for the two worksheets, it illustrates only a portion of the total body of orbit computational code that Project Pluto has built into FIND_ORB.

    Back | Top

    Determining the Orbits of Comets and Asteroids

    Advances in CCD astronomy and in desktop computing technology have made it possible for amateurs to image comets and asteroids, and to determine their orbits. Why do so? I offer two principal reasons:

    (a) to help with the tracking and cataloging of asteroids and comets, which it is now known have impacted Earth in the past, and may do so again at some future date, with possibly disastrous consequences for life on Earth;

    (b) for the sheer intellectual challenge: you want to apply your knowledge of orbital mechanics, telescopes, computers, and CCDs, or you want to learn more about these areas of high technology in an "active" manner.

    In what follows I assume the reason for your interest in comet and asteroid orbit determination is (a) above; if it is (b), then you are at least interested in (a)!

    Where do I get observations of comets and asteroids?

    In CCD astrometry, as opposed to CCD astronomy, you (a) image comets and asteroids -- this much is in common with CCD astronomy; (b) reduce the images to topocentric right ascension (R.A.) and declination (DEC) measurements; (c) determine their orbits using the topocentric R.A. and DEC measurements, and/or (d) report your observations to the IAU's Minor Planet Center in Cambridge, Massachusetts.

    Thus in CCD astrometry your "make your own observations". You need a program such as the Guide program of Project Pluto to help you determine where to look, and you need a program such as Astrometrica or Project Pluto's Charon to reduce the CCD observations that you make. You can report your observations to the IAU's Minor Planet Center in Cambridge, Massachusetts, which uses them to determine and publish an improved orbit for the object. Recently now, thanks to Project Pluto, you can use their "observations to orbits" program, called FIND_ORB, to check the quality of your topocentric R.A. and DEC observations before you send them to the Minor Planet Center.

    The Minor Planet Center also provides a Computer Service at $6.00 per month ($72.00 per year) that is invaluable if you wish to do any serious work with comet and minor planet observations and orbits. I myself subscribe to this service and use the observations to construct test cases for my own research and teaching of theories of orbit determination and orbit improvement. The full web address is

    I should note here that the mathematics of comet and asteroid orbits is pretty much the same as that of interplanetary spacecraft when in coast phase. So if you are interested in the dynamics of interplanetary spacefaring, work with comet and minor planet orbits is not wasted.

    I'm Interested in CCD astronomy, how do I get started?

    You can find articles about CCD astronomy and the advertisements of vendors of telescopes and CCD cameras in the pages of Sky & Telescope. Richard Berry's book, Introduction to Astronomical Image Processing, published by Willmann-Bell is a good reference here.

    I'm already doing CCD astrometry -- what do you have to say to me?

    If you are already doing CCD astrometry, then you are already at least imaging comets and asteroids and reducing their images to topocentric right ascension and declination measurements. You may have read about the FIND_ORB program and how it can be used to determine an orbit from your observations (see Sky & Telescope, June 1998 Software Showcase, "Orbits from Observations").

    Since downloading FIND_ORB myself, I have been carrying out an ongoing e-mail dialogue with Project Pluto, the FIND_ORB program publisher. I have checked out and verified all aspects of the program, and have used the program to check my own work.

    FIND_ORB can be relied upon as an independent means of verifying the quality of your own CCD astrometric observations before sending them to the Minor Planet Center. But it is good for much more than that. If you generate simulated observations of an interplanetary trajectory, FIND_ORB can determine and propagate that trajectory. Given the high quality of implementation, FIND_ORB can be used to generate and/or check your own interplanetary orbital data. It is sure to become the standard for CCD observing experts who want to extend their expertise to include orbit determination from CCD observations.

    My own area of mathematical expertise encompasses the mathematics that is contained in FIND_ORB. I am doing research in this area, and I wish to help interested folks to understand how to use (and how not to misuse) the capabilities that this program provides. One project I have just completed in this area is my Mathcad worksheet on Herget's method. Herget's method is implemented as a preliminary orbit determination method in FIND_ORB; my Mathcad worksheets provide a detailed exposition of the method.

    Why did I choose to work with Mathcad? After many years of programming mathematical algorithms in Fortran, BASIC, Pascal, and C, and many years of supervising others doing the same, I have learned well that the codes themselves are not the best way to teach the math. Yet documenting the math via word processors was a task just as difficult as doing the original programming until Mathcad PLUS 6 came along (previous versions of Mathcad did not have the power of procedural function programming).

    While I do not claim that developing and/or reading Mathcad worksheets is a "piece of cake", I do say that Mathcad is the best way that I have found so far to specify, communicate, and validate orbital mechanics algorithms. So this is why I have adopted MathSoft's Mathcad as the medium by which I communicate my astronomical algorithms for orbit determination and orbit improvement in CCD astrometry.

    Of course I encourage all who are interested in orbit determination and orbit improvement to follow my example: get Mathcad and work with it. But I recognize that not all who are interested in CCD astrometry will want to make this investment of time. If you cannot justify getting and working with Mathcad, then at least download FIND_ORB from Project Pluto's website.

    Back | Top

    More About "Orbit Propagation via State Space Analysis"

    State space analysis is a way of modeling dynamical systems. It employs systems of first-order, ordinary differential equations in matrix form. It has been adopted with great success in the field of control systems engineering. I have written an article about this worksheet,, for publication in The Orrery, a bimonthly publication about "Models of Astronomical Systems" edited and published by Greg Neill / 4541 Anderson / Pierrefonds, Quebec / Canada H9A 2W6.

    The article shows how state space analysis is actually used to integrate the orbit of an artificial Earth satellite. It provides a QuickBASIC 4.5 program (Figure 2 in the article) and its output (Figure 1) for the NEAR Earth escape test case, as reference material for validation of the "orbit propagation via state space analysis" mathematics. The text of the Mathcad worksheet is also provided in the body in the article. The article shows that, by including calculations of the 6-by-1 vector of perturbative accelerations, P(X), where X is the state (cartesian position and velocity), the two-body state equations, Xdot = S(X) X, can be transformed into the perturbed state equations, Xdot = S(X) X + P(X).

    Subscriptions to The Orrery are presently $12.00 per year in the United States, and I believe that back issues are available. For further information write directly to Greg Neill at the address given, or e-mail him at

    (Note: Greg Neill ceased publication of The Orrery some time ago, and I have lost touch with him. But I leave in the references to him out of respect for his work. As a side note, Caxton C. Foster actually started The Orrery. There is a book by Caxton Foster, with that title, now available from Willmann-Bell. I was a contributor to Caxton Foster's book -- see Section 14.7, "Tracking the NEAR Launch.")

    Back | Top

    More about "Effect of a Radial Impulse"

    Suppose that you are in a spacecraft in a perfectly circular Earth orbit, and that you can maneuver by means of a powerful rocket motor that you can point in any direction you wish. Suppose also, for the sake of simplicity, that you can only operate the motor by firing it in one or more brief, but powerful pulses in any direction. This seems contrived, but in actual fact, orbital maneuvers are modeled and analyzed in this way. Realworld motor firings are controlled in such a way as to result in idealized impulsive velocity changes.

    Anyone who has studied the topic of orbital maneuvers knows that if you fire a pulse along your instantaneous velocity vector, in the direction of motion, then you will fall downward along an elliptical orbit whose apogee point is where you fired the pulse, and whose perigee point is 180 degrees away, and at a lower altitude than the fixed altitude of your original circular orbit. And if you fire a pulse of just the right magnitude again when you get to perigee, once more in the same direction as your instantaneous velocity vector, you will enter into a new circular orbit of lower altitude than your original orbit. The two-impulse sequence for maneuvering into the lower orbit is called a Hohmann transfer.

    Similarly, if you fire a pulse along your instantaneous velocity vector, but opposite to your direction of motion, you will climb upward along an elliptical orbit whose perigee point is where you fired the pulse, and whose apogee point is 180 degrees away, and at a higher altitude than the fixed altitude of your original circular orbit. If you now fire a pulse of just the right magnitude again when you get to apogee, once more in the direction opposite to your direction of travel, you will enter into a new circular orbit of higher altitude than your original orbit.

    The second of the two types of Hohmann transfer just described is quite familiar to orbital analysts, because it is used to maneuver Earth satellites into high-altitude orbits. The original orbit is called the parking orbit, the elliptical half-orbital segment is called the transfer orbit, and the final orbit is called the final orbit (of course!).

    Historically, impulses along the instantaneous velocity vector have proved to be the most useful kind of orbital maneuver. But what do you think would happen if you fire a pulse purely along the radius vector, and perpendicular to the instantaneous velocity vector (assuming again a circular orbit), either "straight up" or "straight down"? Will you then go straight up or straight down, or something else? The Mathcad worksheet answers this question.

    Back | Top

    More About "Herget's Method with Cassini's Earth Flyby"

    In the pair of worksheets associated with "Herget's Method for an Asteroid" we solved for the orbital elements of the asteroid (1035) Amata using Dan Burtz's CCD astrometric data, as collected by Dan at the U.S. Air Force Academy's Observatory while he was a graduate student at the University of Colorado, Colorado Springs.

    But in this pair of worksheets and we apply Herget's method to geocentric orbital motion, rather than to heliocentric orbital motion. The CCD astrometric observations were taken in this case by Gordon J. Garradd at Loomberah, Australia, and by Rob McNaught at Siding Spring Observatory, Australia. Bill J. Gray made the data available at the Project Pluto website, along with his own solution using his FIND_ORB program (click on to see Bill's analysis).

    What is most remarkable about the results is that the spacecraft Cassini was on the outbound hyperbolic arm of its Earth flyby trajectory, and at a distance that ranged from 50 - 350 Earth radii from the geocenter, when the 11 observations actually used in the Mathcad worksheet solution were taken. (The Herget's method solution uses 11 observations taken over a span of three days, and uses the two-body UPM theory as the orbit propagator.)

    Despite its use with only a relatively small number of CCD (angles-only) observations, and using only two-body mechanics (UPM), Herget's method converged to a solution yielding a time of perigee passage within about four minutes of JPL's definitive result, and a closest approach distance (1155 km) within about 16 km of JPL's final result (1171 km). Bill Gray's FIND_ORB solution, which accounted for perturbations by the sun and moon, did even better. Using all 15 CCD astrometric observations, Bill got to within better than a minute of JPL's time of perigee passage, and to within about 5 km of JPL's final closest approach distance estimate (1176 km vs. 1171 km).

    This pair of worksheets is the basis for a paper that I presented at the American Astronomical Society's Division on Dynamical Astronomy (AAS/DDA) conference at Yosemite National Park, April 9-12, 2000. The title of the paper is "Astrometry-Based Analysis of Cassini's Earth Flyby", and I presented it at 9:00 a.m. on Wednesday, April 12. To request a copy of the paper, e-mail me with your postal mailing address.

    Back | Top

    More About "Sun Altitudes for Sextant Practice"

    Recall that the sun's altitude is the angle, in degrees (0 degrees = on the horizon, 90 degrees = directly overhead), that it makes with the plane of the local horizon, while its azimuth is the angle that it makes with the direction "true north", as measured in the horizon plane clockwise from true north (0 degrees = due north, 90 degrees = due east, 180 degrees = due south, and 270 degrees = due west). Given that you know your latitude and longitude, how can you accurately calculate the sun's altitude and azimuth, i.e., its position in your local horizon reference frame, at times of interest during a day of interest? This is the problem that solves for you.

    There is a non-trivial modicum of dynamical astronomy in the calculations. Here is an outline of the steps for each time of interest:

    a. Calculate the sun's true Earth-centered, inertial (ECI) cartesian coordinates, referred to the mean equator and equinox of the J2000.0 epoch.

    b. Convert the sun's coordinates from true to astrometric by correcting for aberration (the light-time correction).

    c. Apply a precession matrix to refer the coordinates to the mean equator and equinox of date.

    d. Apply a nutation matrix to refer the coordinates to the true equator and equinox of date.

    e. Calculate the observer's ECI cartesian coordinates and subtract these from the sun's apparent coordinates in order to obtain the sun's cartesian coordinates in the topocentric, horizon-referenced reference frame.

    f. Transform the sun's topocentric, horizon-referenced cartesian coordinates to altitude, azimuth, and topocentric distance in the local horizon reference frame.

    g. Correct the sun's apparent altitude for atmospheric refraction.

    Back | Top

    More About "Sun-Sight Solutions Without Tables"

    My most recent spell of interest in celestial navigation, which has culminated in the worksheets and, was inspired by Dava Sobel's best-selling Longitude (Walker Publishing Company, 1995). Her book is about John Harrison, inventor the first handheld marine chronometer ("H-4", completed in 1759), and his struggles with Nevil Maskelyne, creator of The Nautical Almanac and Astronomical Ephemeris (the annual astronomical tables for navigation, first published in 1766 for use in the year 1767).

    Harrison was trying to win the Longitude Prize, created by the British Parliament in 1714 in order to stimulate a solution to one of the British Navy's (and every other contemporary navy's) most enduring problems: the lack of a way to compute a ship's longitude accurately when it is out of sight of land. In order to win the prize, Harrison invented a timepiece that could keep time accurately on a ship in storm-tossed seas, and through extremes of barometric pressure, temperature, and humidity.

    When one knows the longitude of the home port and the mean solar time difference between the home port and the location at sea, one can immediately calculate the local longitude. This is because Earth rotates 15 degrees of longitude per mean solar hour. One simply multiplies the time difference, in hours, by the conversion factor of 15 degrees per hour to arrive at the longitude difference. (If the "home port" were Greenwich, England, and one's marine chronometer kept Greenwich mean time, the longitude difference at each point of a voyage to the West Indies would be precisely the west longitude, since the longitude at Greenwich is zero degrees.)

    Thus Harrison was fully entitled to the prize no later than March 1762, following successful demonstration of his timekeeper's requisite accuracy during a sea trial consisting of a voyage from England to Jamaica and back. (John Harrison's son, William, had accompanied the Harrisons' fourth marine chronometer, H-4, on the voyage. H-4 lost only five seconds of time after 81 days at sea.)

    Yet Nevil Maskelyne, as a member of the Board of Longitude (the Board of Longitude constituted Parliament's duly appointed, independent panel of judges for the Longitude Prize), managed to keep the Harrisons from being awarded the Longitude Prize for another eight years, while he continued to work on his own astronomical solution to the longitude problem. We should note that Maskelyne's astronomical solution was the forerunner of the modern celestial approach: sight reduction tables, used together with the annual Nautical Almanac to reduce sextant sightings. (Today's Nautical Almanac, as issued jointly and annually by Great Britain and the United States, descends directly from Maskelyne's Nautical Almanac and Astronomical Ephemeris of 1767. So, while Dava Sobel's book turns most of us into partisans for Harrison and against Maskelyne, we should acknowledge that we are beneficiaries of both mens' genius.)

    Celestial navigation is, quite simply, the use of sextant sightings of celestial bodies to calculate one's latitude and longitude. The bodies useful for celestial navigation are the sun and moon, the brighter planets Venus, Mars, Jupiter, and Saturn, and the 57 navigation stars (the celestial positions of 57 bright stars, well distributed over the whole celestial sphere, are tabulated in the Nautical Almanac for use in celestial navigation).

    The marine navigator of today still uses a sextant to measure the apparent altitudes of the celestial bodies visible from his or her geographical location, and then uses sight reduction tables, together with the annual Nautical Almanac, to work up a position solution. Indeed, no finer instrument has emerged for this purpose than the modern sextant, virtually unchanged since World War II, which was independently invented in 1730 by an Englishman, John Hadley, and an American, Thomas Godfrey.

    In my worksheet, I apply a two-star fix method which I devised in 1982 (it was published in Navigation, Journal of the Institute of Navigation, at that time -- see for the reference). I test the method with "perfect" sun altitude measurements from the worksheet, then move on to the realworld, i.e., "imperfect" sun altitude measurements made by Richard R. Shiffman for his "Sextant Noon-Day Sun Sightings" worksheet. Using Mathcad's "regress" and "interp" functions, I smooth Shiffman's actual measurements and then show that pairwise solutions are possible by applying my two-star fix method to adjacent pairs of smoothed sun altitude measurements.

    Inspection of a plot of the smoothed measurements vs. the actual measurements suggests that Shiffman's measurement taking skill improved noticeably as he took more and more measurements. From this arises my conclusion that the last pair of adjacent, smoothed measurements yields the best solution available from the totality of data, i.e., from the 30 local time and sun altitude measurement pairs.

    Back | Top

    More about "Modeling Blackbody Radiation"

    Planck's radiation law is not just of historical interest. In 1964 Arno Penzias and Robert Wilson discovered radio noise emanating from all directions of the sky that is consistent with thermal emission from a blackbody at an equilibrium temperature of just a few degrees Kelvin. They deduced in 1965 that this radio noise is the cosmic microwave background (CMB). For this they were awarded a Nobel Prize in 1978 [1].

    By the 1960s there were two competing theories of the origin of the cosmos, the "steady state" theory and the "Big Bang" theory. Existence of the CMB was predicted by the Big Bang theory, but not by the steady state theory. So when the CMB was found by Penzias and Wilson, most physicists and astronomers came to accept the Big Bang theory and to reject the steady state theory.

    More recently, the Cosmic Background Explorer (COBE) spacecraft measured the CMB in all directions of space, from space (i.e., from Earth orbit). Its measurements of energy density vs. frequency fit almost perfectly on Planck's radiation curve for a blackbody at 2.725 degrees Kelvin. But small "ripples" in energy density were in fact found; these are believed to be evidence of variations in the early universe's energy density. Since these variations are thought to have seeded star and galaxy formation, it would have been a setback for the Big Bang theory had they not been found.


    [1] Mather, John C. and Boslough, John, The Very First Light, Basic Books, New York, 1996; pp. 49-50 and 64. John Cromwell Mather was the original proposer and project scientist for the COBE mission. The COBE satellite was launched on November 18, 1989.

    Back | Top

    Author's Background

    Astroger has pursued the following activities from 1967 to the present.

    Air Force Officer. Astroger began his space career in 1967 as a weather satellite orbital analyst for the Defense Meteorological Satellite Program. He taught mathematics at the U.S. Air Force Academy and, after serving a total of seven years in the Air Force, worked for 21 more years on Air Force space systems developmental projects. He earned a Bachelor of Science degree with high honors in chemistry from the University of Cincinnati and a Master of Arts degree in mathematics from the University of Nebraska at Omaha.

    Space Engineer. Astroger has held positions as (1) ballistic & orbital systems engineer, (2) 427M Space Computational Center (SCC) team leader for space applications software maintenance, (3) Space Defense Operations Center (SPADOC) Block 4B section supervisor for astrodynamic software development, and (4) SPADOC principal engineer for space surveillance applications. He held all of these positions while working for a single company that went by the names Philco-Ford, Aeronutronic Ford, Ford Aerospace & Communications Corporation, and Loral Aerospace Corporation during the years that he was employed there (1974-1995).

    Astroger was present in the NORAD Cheyenne Mountain Complex to assist with tracking data reduction for the Earth 1 and Earth 2 flybys of the Galileo spacecraft (December 1990 and December 1992, respectively), for the Mars Observer launch and Earth escape (September 1992), and for the NEAR launch and Earth escape (February 1996).

    Educator. In the course of his space career, Astroger has taught classes to cadets, to graduate space engineers, and to active duty Air Force operations personnel. He has presented papers at AAS/AIAA astronautics conferences and his work has been published in the Journal of the Astronautical Sciences.

    Educational Publisher. Astroger founded Astronomical Data Service (ADS) in 1976 to provide custom-prepared educational publications to science teachers. ADS distributes a mail order sales catalog annually in the fall. ADS accepts and fills orders by U.S. postal mail. Science teachers write to P.O. Box 26180, Colorado Springs, CO 80936, call (719) 597-4068 and leave a message, or click on the e-mail link at the bottom of this web page to request a catalog.

    Assistant Professor. Astroger held an appointment as Assistant Professor at CU-Colorado Springs during 1996-98. He taught astrodynamics and numerical methods to engineers working at the Lockheed Martin Astronautics Waterton Canyon facility, and then continued to advise his former students as they proceeded through the Master of Engineering in Space Operations program. Lockheed Martin engineers working at Waterton Canyon have been intimately involved with space missions such as Mars Pathfinder, Mars Global Surveyor, and Cassini; they build and operate the Titan and Atlas launchers.

    Book Author. Astroger's 378-page book, Topics in Astrodynamics, commenced publication on October 6, 2003.

    Amateur Radio Hobbyist. Astroger is a licensed amateur radio technician (KB0VJS) and operates his own ground equipment for direct image readout from U.S. and Russian polar-orbiting weather satellites. He has received his best images using a quadrifilar helix antenna built by his son, Jase (AB9BI), from plans in QST magazine.

    What is an "astroger"? See The Astroger Webpages.

    Back | Top

    (c) 1997-2011 by Astroger. Last updated 2011 April 3.

    Accesses:  myspace counters