A Simple Method for Calibrating a Numerical Integrator
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.
Why I’m Not a Big Fan of COTS Software
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.
My Orbital Mechanics Recommended Reading List
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.)
The lost art of Fortran programming
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.
How I create computer programs
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.
Welcome to the Orbital and Celestial Mechanics blog
Hello,
I’ve decided to create a blogging capability in order to inform people about the computer programs I’m currently working on. This is also a place where I can read and respond to current and prospective users of the software I’ve developed.
I’ve been asked many times about why I offer so-called “freeware” computer programs. Actually, it’s a way that I can give back to an aerospace community which has been very good to me over a career spanning more than 30 years. The computer programming I do is not all about making money. It’s a way to keep my mind sharp (and to prevent it from totally turning to mush!) and it’s fun, too. I’m currently semi-retired, working part-time at Kennedy Space Center. After 30+ years in the business, it can be hard to just walk away. I love to solve problems, especially those that involve optimization of space flight mechanics problems.
The software that is offered for sale helps support my computer “habit”. The income helps pay for the Fortran compiler, graphics package and other software products that allow me to create computer programs (and the companion documentation) for solving practical problems. It also helps purchase hardware like a fast computer and a good printer. I’m currently using a “homemade” core I7 computer running Windows 7. Documentation is created with Word, MathType and Grapher.
I’m currently working on a Fortran program that will solve the integrated, n-body gravity-assist problem. The initial guess for this program is based on the patched-conic solution. The program then attempts to solve the problem by numerically integrating the full n-body equations of motion while constraining the B-plane coordinates at the flyby or gravity-assist planet and flying on to the destination or arrival body.
I will try to blog as often as possible. Your comments and questions are welcome. However, I will not do your homework nor will I provide software that you should create as part of your homework assignment.
David Eagle
www.cdeagle.com
Hello world!
Welcome to WordPress.com. After you read this, you should delete and write your own post, with a new title above. Or hit Add New on the left (of the admin dashboard) to start a fresh post.
Here are some suggestions for your first post.
- You can find new ideas for what to blog about by reading the Daily Post.
- Add PressThis to your browser. It creates a new blog post for you about any interesting page you read on the web.
- Make some changes to this page, and then hit preview on the right. You can always preview any post or edit it before you share it to the world.
