Although I’m not a big fan of COTS software, I do use open-source or public domain __source code__ software from reputable organizations. I prefer source code because I can examine the underlying algorithms as well as the programming style.

The following is a brief list with Internet locations for several programs and software suites. There are many more. Surf the internet and I’m sure you’ll wash up on interesting sites.

**General Mission Analysis Tool (GMAT; NASA Goddard)**

gmat.gsfc.nasa.gov

**Solar Activity Site (NASA MSFC)**

sail.msfc.nasa.gov

**SPICE and MICE (NASA JPL)**

naif.jpl.nasa.gov/naif/

**Naval Observatory Vector Astrometry Software (NOVAS)**

www.usno.navy.mil/USNO/astronomical-applications/software-products/novas

**Standards of Fundamental Astronomy (SOFA; IAU)**

There are numerous resources for open-source and public domain numerical methods that are also useful for astrodynamics analysis.

This blog is a collection of resources that I have found useful for understanding, selecting and implementing numerical methods for aerospace trajectory optimization.

**Books**

*Optimal Trajectories for Space Navigation* by Derek F. Lawden, Butterworths, London, 1963.

*Spacecraft Trajectory Optimization*, Edited by Bruce A. Conway, Cambridge University Press, 2010.

*Advanced Design Problems in Aerospace Engineering, Volume 1: Advanced Aerospace Systems*, Edited by Angelo Miele and Aldo Frediani, Kluwer Academic Publishers, 2004.

*Practical Methods for Optimal Control and Estimation Using Nonlinear Programming* by John. T. Betts, SIAM, 2010.

**Survey Papers**

“A Survey of Methods Available for the Numerical Optimization of Continuous Dynamic Systems”, by Bruce A. Conway, Journal of Optimization Theory and Applications, 2011.

“Survey of Numerical Methods for Trajectory Optimization”, by John T. Betts, AIAA Journal of Guidance, Control and Dynamics, Vol. 21, No. 2, March-April 1998.

“A Survey of Numerical Methods for Optimal Control”, by Anil V. Rao, AAS 09-334, 2009.

**Software**

Applied Mathematical Analysis (Fortran object code libraries), www.appliedmathematicalanalysis.com

PROPT (MATLAB), tomdyn.com

DIDO (MATLAB), www.elissarglobal.com

GPOPS (MATLAB), www.gpops.org

SNOPT (MATLAB/DLL) scicomp.ucsd.edu/~peg/

One of the most common tasks in astrodynamics is the numerical solution of the orbital equations of motion. For motion which is not Keplerian or two-body, this is usually accomplished with some type of numerical method which propagates a spacecraft’s motion from one time epoch to another (earlier or later) epoch.

For elliptical orbits, it is often necessary to use a variable step size integrator. These techniques change the integration step size as a function of a spacecraft’s location in its orbit. Classical methods include Runge-Kutta algorithms that solve the first-order form of the differential equations of orbital motion. My favorite algorithm is a Runge-Kutta-Fehlberg (RKF) method because the integration coefficients have been “tailored” to solve the types of differential equations found in astrodynamics and celestial mechanics.

Many of these numerical methods allow the user to control how well the differential equations are solved via a truncation error tolerance or perhaps absolute and relative error controls. Setting these algorithm controls to correct values determines the performance of the numerical method.

If we formulate the differential equations of un-perturbed orbital motion in first-order form, it is possible to study the effect of algorithm controls on the orbit propagation using the following simple method.

The initial conditions for the x, y and z components of the position and velocity vectors are

r_x(0) = 1 – ecc; r_y(0) = 0; r_z(0) = 0;

v_x(0) = 0; v_y(0) = sqrt[ (1 + ecc) / (1 – ecc) ]; v_z(0) = 0

where ecc is the orbital eccentricity.

The system of six first-order differential equations (y_dot) is given by

y_dot(1) = v_x; y_dot(2) = v_y; y_dot(3) = v_z

y_dot(4) = –r_x / r^3; y_dot(5) = –r_y / r^3; y_dot(6) = –r_z / r^3

For this formulation, the Keplerian or un-perturbed orbital period is exactly equal to 2 pi. In the last three equations r is the scalar magnitude of the position vector; r = sqrt(r_x^2 + r_y^2 + r_z^2).

To test an integration method, select an eccentricity value, say 0.75, set typical values for the algorithm controls, and run the propagator for an integer number of 2 pi “orbits”, printing the position and velocity vector components at “time” intervals of 2 pi. An examination of the state vector (position and velocity vector components) at the end of each complete orbit will indicate how well the integrator matches the initial conditions defined above as the propagation time increases. The goal is to select the “best” algorithm control for acceptable state vector differences for a user-defined propagation duration. A highly elliptical orbit exercises the step size features of a numerical method.

This method can also be used with dimensional values (perhaps kilometers for position and kilometers per second for velocity) for the initial state vector and Keplerian orbital period (seconds). This form also requires the value of the Earth’s gravitational constant (mu) in the three differential equations expressed as

y_dot(4) = –mu * r_x/r^3; y_dot(5) = –mu * r_y/r^3;

y_dot(6) = –mu * r_z /r^3

We have used this technique to calibrate an RKF method for (1) propagating orbital debris in the geosynchronous orbit region for one hundred years, (2) studying the lifetime of Earth satellites, (3) propagating lunar and interplanetary trajectories, (4) predicting impact locations for re-entering orbital debris, etc.

My first job after graduate school involved the independent verification and validation (IV&V) of the guidance software for the Inertial Upper Stage (IUS). The most important thing that task taught me was not to trust analysis performed by others until you perform that same analysis yourself using as many analytic and software techniques as possible.

A typical example was a problem we found in the state transition matrix subroutine. It was based on Goodyear’s method but was not working correctly. We proved this by creating a version using numerical derivatives and two versions based on other semi-analytic methods.

There are several good commercial-off-the-shelf (COTS) software products available for aerospace mission analysis. They usually work well for exactly the analysis they were designed to do.

However, some of the reasons I’m not fond of COTS software are as follows;

(1) inadequate information about the underlying algorithms, fundamental constants and computational assumptions

(2) incomplete or minimal testing of the software

(3) never quite does everything you need

(4) often difficult to add your own features

(5) bug fixes and updates are not always timely

Therefore, I tend to create the analysis capability I need and not rely on COTS software.

I’m often asked to recommend orbital mechanics books and technical reports. The following is a short list of the resources that I use frequently. Most of these books can be purchased from Amazon and several are available as eBooks.

“**Orbital Mechanics**”, Vladimir Chobotov, published by the American Institute of Aeronautics and Astronautics (AIAA), www.aiaa.org. This book includes clear explanations about fundamental orbital mechanics. It also contains a advanced chapter about low-thrust orbit transfer. Many of the chapters were written by staff members of the Astrodynamics Department at The Aerospace Corporation (which we fondly called Circle A). The book includes a CD ROM with several orbital mechanics applications.

“**Fundamentals of Astrodynamics and Applications**”, David Vallado, published by Microcosm Press. Companion software for free in MATLAB, C and Fortran is also available. This is a excellent text written by someone with lots of practical experience.

“**Orbital Mechanics for Engineering Students**”, Prof. Howard Curtis, published by Butterworth-Heinemann. Although written for college students, the practicing aerospace engineer will also find this book quite useful. MATLAB code for many of the examples and homework problems is available.

“**An Introduction to the Mathematics and Methods of Astrodynamics**”, Prof. Richard Battin, published by the AIAA. This is one of the most comprehensive textbooks about fundamental orbital mechanics.

“**Methods of Orbit Determination**”, P. R. Escobal, published by Krieger Publishing Company. This is a classic text with algorithms for predicting rise and set, shadow conditions, etc.

“**Adventures in Celestial Mechanics**”, Prof. Victor Szebehely, published by Wiley. Easy to read but quite informative.

There are also several technical documents in the Technical Reports section of my website (www.cdeagle.com). Among the resources in this section are the following;

“Orbit Theory and Applications”, E. Taylor

“Handbook of Orbital Perturbations”, Prof. Leon Blitzer

“Astrodynamics”, G. S. Gedeon

*“If I have seen further it is by standing on ye sholders of Giants.”*

**Isaac Newton, February 1676**

(Note that *sholders* is the original spelling, not a typo.)

My first and only Fortran programming class was at a small college with a big computer. The computer was actually part of the business school and I was one of only a few people who seemed to use it only a regular basis. The computer was fast which allowed me to make programming mistakes fast and fix bugs even faster.

The source code was punched on small cards by a machine that resembled the old teletype devices. It almost walked around the room while punching the cards. Now I have totally dated myself!

I was excited about programming in Formula Translation. I could finally solve some fairly sophisticated math and engineering problems again and again by simply changing the input to the software.

Engineering and math students graduating from college today seem to be quite proficient in MATLAB, C and C++. I personally can barely spell C! Unless we are developing deliverable software, people tend to program in a language with a high comfort level. Many young engineers have heard of Fortran but few have actually seen Fortran source code. I am not aware of Fortran classes offered to engineering students.

In aerospace engineering we do lots of things in MATLAB. But we also do many more types of analysis in Fortran. The so-called “industry standard” aerospace programs like POST (Program to Optimize Simulated Trajectories) and OTIS (Optimal Trajectories by Implicit Simulation) are both written in Fortran. I don’t know about other engineering disciplines, but in aerospace the majority of heritage software is written in Fortran. Also, I don’t know many people who are interested in porting thousands of lines of source code to another high level language.

Fortran was written for number-crunching. It has evolved over the years to take advantage of such features as modern programming constructs, dynamic memory allocation and modules. It’s easy to read and easy to use. Fortran is a compiled language that executes orders of magnitude faster than MATLAB which is an interpretive language. If you know MATLAB you should be able to learn Fortran quickly.

**We need engineers who are not intimidated by Fortran and large scientific simulations.**

For example, where I work we have a task to “hook” the Fortran version of the GRAM 2010 atmosphere model into a trajectory simulation. Since both GRAM 2010 and the trajectory simulation are written in Fortran, the integration and checkout should be straight forward. The problem is the lack of engineers who are proficient with Fortran.

In this blog I will talk about how I create different computer programs using MATLAB, Fortran and Numerit Pro.

Most of these computer programs rely on a suite of subroutines (Fortran) or functions (MATLAB and Numerit Pro). These are typically short utility routines that perform a single type of calculation (coordinate transformation, evaluate a system of differential equations, etc.). The collection of Fortran routines I call the “Orbital Mechanics with Fortran” toolbox. Likewise, for MATLAB the suite is called “Orbital Mechanics with MATLAB”.

These support routines have been used for many years and are for all practical purposes, “bullet proof”. This means that they will always work correctly provided I give them the right information in the correct units. This is very important when developing new applications. If the application or main program crashes, it’s likely something in the call to these routines and not the routines themselves. This approach makes creating complex scientific simulations much easier.

An excellent example is the computer program that solves the n-body, integrated Earth-to-Mars flight mechanics problem. This is a fairly complicated trajectory optimization problem.

The solution is achieved using the following major computational steps;

(1) solve the two-body, patched-conic interplanetary Lambert problem for the energy (C3), declination (DLA) and asymptote (RLA) of the outgoing hyperbola

(2) compute the orbital elements of the geocentric launch hyperbola and the components of the interplanetary injection delta-v vector

(3) perform geocentric orbit propagation from perigee of the geocentric launch hyperbola to the Earth’s sphere-of-influence (SOI)

(4) perform an n-body heliocentric orbit propagation from the Earth’s SOI to the B-plane at Mars encounter

(5) target to the user-defined B-plane coordinates by minimizing a heliocentric trajectory correction maneuver (TCM) at the Earth’s sphere-of-influence

The analysis results at the conclusion of each step can be displayed and checked for accuracy before going on to the next step. Of course, it is the responsibility of the programmer/analyst to determine is the results from each step are correct.

You might think of this programming approach as “divide and conquer”.

Finally, remember that numbers computed by computers are not always the right numbers. They simply do what you tell them.