Implementation of
Gauss-Jackson Integration
for Orbit Propagation
1
Matthew M. Berry
2
and Liam M. Healy
3
Abstract
The Gauss-Jackson multi-step predictor-corrector method is widely used in numerical in-
tegration problems for astrodynamics and dynamical astronomy. The U.S. space surveil-
lance centers have used an eighth-order Gauss-Jackson algorithm since the 1960s. In this
paper, we explain the algorithm including a derivation from first principals and its relation
to other multi-step integration methods. We also study its applicability to satellite orbits in-
cluding its accuracy and stability.
Introduction
Determination and prediction of orbits requires an orbit propagator that finds the
phase space state of a satellite at one time based on its state at another time. This
function has traditionally been performed by an analytical computation such as
Hamiltonian normalization (general perturbations). While numerical integration
(special perturbations) can provide a wider range of force models, until recent years
the computation time required made it prohibitive for routine use in space surveil-
lance. With greatly increased computation speeds now widely available, the daily
processing of 10,000 or more satellites routinely tracked by space surveillance cen-
ters can conceivably be based on numerical integration, and a significant fraction
of this number are processed numerically now.
The Gauss-Jackson integrator [1] is a fixed-step predictor-corrector method. The
Gauss-Jackson software used in space surveillance originated in the 1960s as part 
1
Based on paper AAS 01-426 presented at the AAS/AIAA Astrodynamics Specialists Conference, Quebec
City, Canada, July 30?August 2, 2001.
2
Astrodynamics Engineer, Analytical Graphics Inc., 220 Valley Creek Blvd., Exton, PA 19341. E-mail:
mberry@agi.com. This work was completed while the first author was a Graduate Assistant in the De-
partment of Aerospace and Ocean Engineering at Virginia Tech and an employee of the Naval Research 
Laboratory.
3
Research Physicist, Naval Research Laboratory, Code 8233, Washington, DC 20375-5355. E-mail:
Liam.Healy@nrl.navy.mil.
The Journal of the Astronautical Sciences, Vol. 52, No. 3, July?September 2004, pp. 331?357
331
of the Advanced Orbit/Estimation Subsystem (AOES) [2]. Over the years, AOES
has given rise to a family of astrodynamics applications used in the various space
surveillance centers whose integrators are essentially unchanged, including a sys-
tem known as SpecialK [3]. SpecialK has been used by Naval Network and Space
Operations Command for special perturbations catalog maintenance for the last
six years. Another member of the family is Astrodynamics Support Workstation
(ASW) [4], used in Cheyenne Mountain to support the special perturbations cata-
log and in support of conjunction assessment for NASA. We have chosen to use
SpecialK as the implementation on which to focus, but our analysis has broad ap-
plicability to all Gauss-Jackson integrators, including the one in ASW. This paper
presents a study of the theory upon which the integrator is based, how it is imple-
mented, and its accuracy and stability.
After an overview of multi-step integrators, we start with an explanation of dif-
ference tables and operators that form the core of the derivation of all the multi-step
methods. We apply these to differentiation in the next section, followed by single 
integration where we derive the nonsummed Adams methods. The Adams method
has a summed form which is presented in the succeeding section. The next section,
on the St?rmer-Cowell integrator for second-order differential equations, makes
use of the Adams results. The summed form of the St?rmer-Cowell integrator,
called the Gauss-Jackson integrator, is presented next. Having avoided the issue of
how to start the integrator when backpoints are unavailable, we next address this
issue, in which the concept of a mid-corrector is introduced. To this point, the for-
mulation is in terms of the difference operators; for computational efficiency, the
ordinate form sums like terms but makes the result order-specific. The details are
included in the next section.
With the derivation of the integrators complete, we finish with two sections. The
first is a detailed description of implementation complete with pseudo-code that
should suffice for a reasonably skilled programmer to implement in any general-
purpose programming language. The second shows the effect of step size and order
of the Gauss-Jackson integrator on accuracy and stability.
Multi-Step Integrators
An ordinary differential equation (ODE) is an equation of the form
(1)
in which the derivative of a dependent variable y is a function f of that variable and
the independent variable t. In order to find a specific solution of a differential equa-
tion the initial conditions must be given. A given differential equation
along with initial conditions is known as an initial value problem. The solution to
an initial value problem, , can be found analytically for some differential equa-
tions. However, most differential equations must be solved numerically. Algorithms
that numerically solve initial value problems are known as numerical integrators.
Numerical integrators give an approximate solution at distinct values of t, known as
mesh points. The numerical solution at a mesh point is denoted , and due to
error in the numerical method will likely be different from the exact solution, .
The differential equation in (1) is called a first-order differential equation be-
cause the equation is for the first derivative. Differential equations may have a
higher order
yH20849t
n
H20850
y
n
t
n
yH20849tH20850
H20849t
0
, yH20849t
0
H20850H20850
dy
dt
H33527 fH20849t, yH20850
332 Berry and Healy
(2)
where n is the order of the differential equation. Initial value problems involving
higher-order ODEs require initial conditions for each derivative through .
Though higher-order ODEs may be rewritten as a system of first order equations
and solved with a standard numerical integrator, some numerical integrators are de-
signed to directly solve higher-order ODEs. Integrators that solve second-order
ODEs are called double-integration methods, while integrators that solve first-order
ODEs are called single-integration methods.
Multi-step integrators integrate forward from a particular time to the next mesh
point using function values at the current point as well as several previous mesh
points. The set of the previous points used as well as the current point is called the
set of backpoints. One disadvantage of multi-step integrators is that the initial set
of backpoints must be found through some startup procedure, which complicates
the method. The methods develop a Taylor series in the time separation between
mesh points, and this series must be truncated after some number of terms, which
thereby gives the order of the method. The order is one less than the number of
backpoints, and is not related to the order of the differential equation.
Multi-step integrators are known as predictor-corrector methods. The algorithms
first predict a y value for the next mesh point and the function is evaluated at
the predicted point. That predicted function value is added to the set of backpoints,
and a corrector formula is used with this revised set of backpoints to refine the pre-
dicted y value.
This procedure offers several variations on the implementation. First, when
adding the predicted value to the set of backpoints, the farthest backpoint from the
current point may or may not be dropped. If it is not dropped, the corrector is of
higher order than the predictor because the corrector uses one more backpoint than
the predictor. Otherwise, the predictor and corrector are of the same order. Also, a
second evaluation may or may not be performed after the corrector. A second eval-
uation improves the accuracy of the method, because improved function values are
used in the set of backpoints in the subsequent steps. However, the additional eval-
uation causes an increase in run time. Methods that use one evaluation per step are
known as Predict-Evaluate-Correct, or PEC, methods, and methods that perform a
second evaluation are called PECE methods. Some implementations perform addi-
tional iterations of the corrector to meet some tolerance, and are called or
methods.
The methods can be derived in a variety of ways, but the fixed-step methods can
be derived in a simplified form by using backward difference operators, which are
described in the next section. The derivations of the Adams, summed-Adams,
St?rmer-Cowell, and Gauss-Jackson methods follow. These derivations follow the
derivations presented by Maury and Segal [5], and the NORAD document known
as TP008 [6].
The function must be continuous and smooth through the set of back-
points. If there are any discontinuities in f, for example going through eclipse when
solar radiation pressure is considered, the integration must either be restarted, or
modified to handle the discontinuity [7].
In this paper we do not address variable step or error control, including 
s-integration, in numerical integration. The interested reader is referred to [8]
and [9] for information on this topic.
fH20849t, yH20850
PEH20849CEH20850
n
PH20849ECH20850
n
fH20849t, yH20850
y
H20849nH110021H20850
d
n
y
dt
n
H33527 fH20849t, y, yH11032,...,yH20850
H20849nH110021H20850
Implementation of Gauss-Jackson Integration for Orbit Propagation 333
Tables and Operators
Predictor-corrector integrators can be defined in terms of difference tables. For a
fixed-step method, assume that solution values are known on a discrete set of
equally-spaced mesh points, where h is
the step. As shorthand, these values are referred to as , so this set of
five values is , , , , . The corresponding function values are ,
where is the numerical solution at the mesh point (not the exact solution ).
The differential equations are solved by taking differences of the function values.
There are three kinds of differences, represented by different operators:
forward difference
backward difference
central difference
Note that first central differences cannot be calculated directly, as the half-step points
are not in the sequence of points, though the second central differences (see below for
definition of second differences) can be calculated. As previous authors ([10], [11])
have recognized, conversion between series in these operators is straightforward.
The backward difference derivation is used here, because when one is predicting or
correcting, i.e., any time other than startup, it is the most natural. Derivations using
other difference operators give algebraically equivalent methods, though differences
in their implementation may give slightly different results due to round-off error.
Not only can the differences of the function f be computed, the differences of the
differences, or second differences, can be computed, and so on. For instance,
the square of the backward difference operator is the operator applied to the first
difference
(3)
The relationships of the backward differences are illustrated in Table 1. The arrows
in the table point towards the difference; the upper component is always subtracted
from the lower. For example, .
In addition to the three difference operators, there is also a displacement opera-
tor E. The displacement operator is defined as
(4)
Powers of this operator act as expected; e.g., . The backward difference
operator can be defined in terms of the displacement operator
(5)
and alternatively the displacement operator can be defined in terms of the backward
difference operator
(6)
The summed integration formulas (see Summed Adams Method section), con-
tain negative powers of the operators; e.g., . The backward difference tableH11612
H110022
E H33527
1
1 H11002H11612
H11612 H33527 1 H11002 E
H110021
E
2
f
i
H33527 f
iH110012
Ef
i
H33527 f
iH110011
H11612
2
f
1
H33527 H11612f
1
H11002H11612f
0
 H33527 f
i
H11002 f
iH110021
H11002 H20849 f
iH110021
H11002 f
iH110022
H20850 H33527 f
i
H11002 2f
iH110021
H11001 f
iH110022
 H11612
2
f
i
H33527 H11612H11612f
i
H33527 H11612H20849 f
i
H11002 f
iH110021
H20850
 H9254f
i
H33527 f
iH110011/2
H11002 f
iH110021/2
 H11612f
i
H33527 f
i
H11002 f
iH110021
 H9004f
i
H33527 f
iH110011
H11002 f
i
yH20849t
n
H20850t
n
y
n
f
n
H33527 fH20849t
n
, y
n
H20850t
2
t
1
t
0
t
H110021
t
H110022
t
n
H33527 t
0
H11001 nh
...,t
0
H11002 2h, t
0
H11002 h, t
0
, t
0
H11001 h, t
0
H11001 2h,...,
y
n
334 Berry and Healy
(Table 1) can be extended to the left to get negative powers of the backward 
difference operator, as in Table 2. Extending the differencing to the left gives
. This relation can be changed into a recursion formula that
defines the inverse operator
(7)
Note that the initial ( ) term in the sum is arbitrary; the difference relation
holds no matter what this value is. This value is effectively an integration constant
. Thus, if the inverse operator is defined as a sum
(8)H11612
H110021
f
1
H33527
H20902
C
1
H11001 H20888
i
jH335271
f
j
,
C
1
H11002 H20888
0
jH33527iH110011
f
j
,
if i H11350 1
if i H11349H110021
C
1
i H33527 0
H11612
H110021
f
i
H33527 f
i
H11001H11612
H110021
f
iH110021
f
i
H33527 H11612
H110021
f
i
H11002H11612
H110021
f
iH110021
Implementation of Gauss-Jackson Integration for Orbit Propagation 335
TABLE 1. Backward Difference Table
0
1
2
3 H11612
3
f
3
H11612
2
f
3
H11612f
3
f
3
H11612
3
f
2
H11612
2
f
2
H11612f
2
f
2
H11612
3
f1H11612
2
f1H11612f1f1
H11612
3
f0H11612
2
f0H11612f0f0
H11612
3
fH110021H11612
2
fH110021H11612fH110021fH110021H110021
H11612
3
fH110022H11612
2
fH110022H11612fH110022fH110022H110022
H11612
3
fH110023H11612
2
fH110023H11612fH110023fH110023H110023
H11612
3
fiH11612
2
fiH11612fifii
TABLE 2. Inverse Backward Differences
0
1
2 f2 H11002 f1f2C1 H11001 f1 H11001 f2C2 H11001 2C1 H11001 2f1 H11001 f2
f1 H11002 f0f1C1 H11001 f1C2 H11001 C1 H11001 f1
f0 H11002 fH110021f0C1C2
fH110021 H11002 fH110022fH110021C1 H11002 f0C2 H11002 C1H110021
f
H110022 H11002 fH110023fH110022C1 H11002 f0 H11002 fH110021C2 H11002 2C1 H11001 f0H110022
H11612f
ifiH11612
H110021
fiH11612
H110022
fii
the difference relation holds. The constant is discussed further in the Imple-
mentation section.
Powers of the summation are understood to be multiple sums, analogous to
multiple differences. The presence of the operator in the Gauss-Jackson for-
mulas sometimes gives rise to the name second sum integration formula. The sec-
ond sum has an integration constant , analogous to .
Differentiation
Differentiation is the operator that turns a function into its deriva-
tive function. The differentiation operator D can be represented in terms of E by
noting that
(9)
for any real number p, using a Taylor series expansion. The exponential of an op-
erator is to be interpreted as its ordinary Taylor expansion
(10)
with powers of the operator well-defined. Thus we may identify the two operators
(11)
or taking the logarithm
(12)
This expression can be written in terms of the backward difference operator (6)
(13)
The integration methods can now be derived by taking the inverse of the differen-
tiation operator. The Adams method is derived first, and the other methods are
derived from it.
Adams Method
Indefinite integration is the operator inverse of differentiation
(14)
so its action on an arbitrary function is to give the antiderivative
(15)
The single-integration operator is computed as the inverse of the differentia-
tion operator (13)
(16)
With this operator Adams integration can be developed. The focus here is satellite
orbits, so the function f is replaced with acceleration r?.
D
H110021
H33527 H11002
h
logH208491 H11002H11612H20850
D
H110021
FH20849tH20850 H33527 D
H110021
fH20849tH20850
H20885
dt H33527 D
H110021
hD H33527 log E H33527 H11002logH208491 H11002H11612H20850
phD H33527 p log E
E
p
H33527 e
phD
e
t
H33527 1 H11001 t H11001
t
2
2!
H11001
t
3
3!
H11001 ...
 H33527 e
phD
fH20849t
0
H20850
 H33527 fH20849t
0
H20850 H11001 phDfH20849t
0
H20850 H11001
H20849 phH20850
2
2!
D
2
fH20849t
0
H20850 H11001
H20849 phH20850
3
3!
D
3
fH20849t
0
H20850 H11001 ...
 E
p
fH20849t
0
H20850 H33527 fH20849t
0
H11001 phH20850
D H33527 dH20862dt
C
1
C
2
H11612
H110022
H11612
H110021
C
1
336 Berry and Healy
The definite integral is computed by using the displacement operator on the
indefinite integral
(17)
For one step, p is 1. This integration corresponds to a predictor in a multi-step
numerical integration, because the equation finds a value of r? at a future time. If
, the predictor operator J can be written in terms of the backward differ-
ence operator as
(18)
To integrate from to , shift backward using the displacement operator (6)
(19)
This operation corresponds to the corrector in a multi-step numerical integration,
because it assumes that the velocity is already known and uses the correspon-
ding acceleration to recalculate a new, better, velocity. It is on this expression that
we focus initially because the operator has the simplest expansion.
The coefficients of the operators may be developed using a recursion relation
[12]. For convenience in developing expansions, we define an operator L and its
Taylor expansion with coefficients as
(20)
We shall develop the expansion of L recursively through the standard Taylor series
of the logarithm
(21)
which can be substituted in the definition of L
(22)
and then expanded in the unknowns 
(23)
Expanding and grouping by powers of H11612
(24) H11001
H20873
1
4
c
0
H11001
1
3
c
1
H11001
1
2
c
2
H11001 c
3
H20874
H11612
3
H11001 ...H33527 1
 c
0
H11001
H20873
1
2
c
0
H11001 c
1
H20874
H11612H11001
H20873
1
3
c
0
H11001
1
2
c
1
H11001 c
2
H20874
H11612
2
H20849c
0
H11001 c
1
H11612H11001c
2
H11612
2
H11001 c
3
H11612
3
H11001 ...H20850
H20873
1 H11001
1
2
H11612H11001
1
3
H11612
2
H11001
1
4
H11612
3
H11001 ...
H20874
H33527 1
c
i
L H33527 H11002
H11612
logH208491 H11002H11612H20850
H33527
1
1 H11001
1
2
H11612H11001
1
3
H11612
2
H11001
1
4
H11612
3
H11001 ...
H33527 H20888
H11009
iH335270
c
i
H11612
i
logH208491 H11002H11612H20850 H33527 H11002H11612 H11002
1
2
H11612
2
H11002
1
3
H11612
3
H11002
1
4
H11612
4
...
L H33527 E
H110021
J H33527 H11002
H11612
logH208491 H11002H11612H20850
H33527 c
0
H11001 c
1
H11612H11001c
2
H11612
2
H11001 c
3
H11612
3
H11001 ...
c
i
n
th
E
H110021
J H33527 H208491 H11002H11612H20850J H33527 H11002
H11612
logH208491 H11002H11612H20850
t
n
t
n
H11002 h
 H33527 H11002
H11612
H208491 H11002H11612H20850 logH208491 H11002H11612H20850
 J H33527
1
h
H20849E H11002 1H20850D
H110021
H33527
1
h
H20851H208491 H11002H11612H20850
H110021
H11002 1H20852D
H110021
H33527
1
h
H11612
H208491 H11002H11612H20850
D
H110021
p H33527 1
H20849E
p
H11002 1H20850D
H110021
r?
n
H33527 H20849E
p
H11002 1H20850r?
n
H33527
H20885
tnH11001ph
tn
r?dt
Implementation of Gauss-Jackson Integration for Orbit Propagation 337
allows for the development of a recursion relation for the coefficients
(25a)
(25b)
(25c)
or, in general
(26)
for with . The first few coefficients are
(27)
With the coefficients known, integration from to may be written in terms
of L as
(28)
A formula for the corrected value can be found by performing the integration and
using the coefficients 
(29)
This formula is the Adams-Moulton corrector formula in difference form. This
form is called the difference form because it uses the differences . An alternate
form, called the ordinate form, in which the differences are rewritten in terms of the
function values, is described in the Ordinate Forms section.
Returning to the predictor (18), the expansion of J can be computed by noting that
(30)
Each coefficient can be expressed as the sum of the coefficients ,
(31)
so they may be computed directly from (27)
(32)
i
H9253
i
0
1
1
1
2
2
5
12
3
3
8
4
251
720
5
95
288
6
19087
60480
7
5257
17280
8
1070017
3628800
H9253
i
H33527 H20888
i
kH335270
c
k
k H11349 ic
k
H9253
i
J H33527 H208491 H11002H11612H20850
H110021
L H33527 H208491 H11001H11612H11001H11612
2
H11001 ...H20850H20849c
0
H11001 c
1
H11612H11001c
2
H11612
2
H11001 ...H20850 H33527 H20888
H11009
iH335270
H9253
i
H11612
i
H11612
i
r?
n
H33527 r?
nH110021
H11001 h
H20873
1 H11002
1
2
H11612H11002
1
12
H11612
2
H11002
1
24
H11612
3
H11002
19
720
H11612
4
H11002
3
160
H11612
5
H11001 ...
H20874
r?
n
c
i
H20885
tn
t
n
H11002h
r?dt H33527 hE
H110021
Jr?
n
H33527 hLr?
n
t
n
t
n
H11002 hc
i
i
c
i
0
1
1
H11002
1
2
2
H11002
1
12
3
H11002
1
24
4
H11002
19
720
5
H11002
3
160
6
H11002
863
60480
7
H11002
275
24192
8
H11002
33953
3628800
c
0
H33527 1n H11350 1
c
n
H33527 H11002H20888
nH110021
iH335270
c
i
n H11001 1 H11002 i
 0 H33527
1
3
c
0
H11001
1
2
c
1
H11001 c
2
 0 H33527
1
2
c
0
H11001 c
1
 1 H33527 c
0
338 Berry and Healy
and the integral is then
(33)
or
(34)
This formula is the Adams-Bashforth predictor formula in difference form, so called
because when applied to the velocity and acceleration at it finds the velocity at the
next mesh point .
In practice, given m known accelerations
(35)
the backward difference table (Table 1) for the derivatives can be computed through
, a predicted value for is made based on (34), then that value is corrected
with (29).
Although our purpose in deriving the nonsummed Adams-Bashforth and Adams-
Moulton coefficients is to use these results in the Gauss-Jackson method, these
formulae stand on their own as a widely-used numerical integration method of
first-order differential equations. For greater accuracy, summed methods may be
preferable, and we treat the summed form of the Adams method next.
Summed Adams Method
An alternative form of single integration, called the summed form, can be devel-
oped using the summation operator . The Adams-Moulton corrector formula (29)
may be written so both velocity terms are on the left side
(36)
The difference of the two velocity terms may be written using the backward differ-
ence operator
(37)
The corrected value of the velocity can be found by applying to both sides
of (37), [5]
(38)
This expression is the summed Adams corrector formula. The summed form uses
the same coefficients as the Adams-Moulton corrector, but they have been shifted
by one place. The formula is called the summed form because of the presence of
the summation operator.
r?
n
H33527 h
H20873
H11612
H110021
H11002
1
2
H11002
1
12
H11612H11002
1
24
H11612
2
H11002
19
720
H11612
3
H11002
3
160
H11612
4
H11001 ...
H20874
r?
n
H11612
H110021
H11612r?
n
H33527 h
H20873
1 H11002
1
2
H11612H11002
1
12
H11612
2
H11002
1
24
H11612
3
H11002
19
720
H11612
4
H11002
3
160
H11612
5
H11001 ...
H20874
r?
n
r?
n
H11002 r?
nH110021
H33527 h
H20873
1 H11002
1
2
H11612H11002
1
12
H11612
2
H11002
1
24
H11612
3
H11002
19
720
H11612
4
H11002
3
160
H11612
5
H11001 ...
H20874
r?
n
H11612
H110021
r?
nH110011
H11612
mH110021
r?
nH11002mH110011
,...,r?
n
t
nH110011
t
n
r?
nH110011
H33527 r?
n
H11001 h
H20873
1 H11001
1
2
H11612H11001
5
12
H11612
2
H11001
3
8
H11612
3
H11001
251
720
H11612
4
H11001
95
288
H11612
5
H11001 ...
H20874
r?
n
H20885
tnH11001h
t
n
r?dt H33527 hJr?
n
Implementation of Gauss-Jackson Integration for Orbit Propagation 339
A similar process can be performed on the Adams-Bashforth predictor for-
mula (34)
(39)
giving the summed Adams predictor. The summed form gives the predicted or cor-
rected value by integrating directly from epoch, while the nonsummed form inte-
grates from point to point. According to Henrici ([13], p. 327) the summed form is
preferred for reducing the propagation of roundoff error.
St?rmer-Cowell
To find a corrector formula for position, the corrector may be applied to the
Adams-Moulton corrector (28)
(40)
Application of the corrector to both velocity terms, e.g., , gives an
expression for position
(41)
To find the coefficients in the expansion of 
(42)
in terms of the Adams coefficients c (26), note that (42) is the square of (20)
(43)
so the first nine coefficients are
(44)
These are the Cowell corrector coefficients. The corrector formula is given by using
these coefficients in (41)
(45)
This formula is the Cowell corrector formula [5]. It is a double-integration formula,
since it computes the position given acceleration.
As with single integration, the predictor formula is found by shifting the step at
which the operator is applied, in other words applying the displacement operator E
to both sides of (41)
(46)
Writing the displacement operator E in terms of H11612, (6), and applying to 
(47)
EL
2
H33527 H208491 H11002H11612H20850
H110021
L
2
H33527 H208491 H11001H11612H11001H11612
2
H11001 ...H20850H20849q
0
H11001 q
1
H11612H11001q
2
H11612
2
H11001 ...H20850 H33527 H20888
H11009
iH335271
H9261
i
H11612
i
L
2
r
nH110011
H33527 2r
n
H11002 r
nH110021
H11001 h
2
EL
2
r?
n
r
n
H33527 2r
nH110021
H11002 r
nH110022
H11001 h
2
H20873
1 H11002H11612H11001
1
12
H11612
2
H11002
1
240
H11612
4
H11002
1
240
H11612
5
H11001 ...
H20874
r?
n
i
q
i
0
1
1
H110021
2
1
12
3
0
4
H11002
1
240
5
H11002
1
240
6
H11002
221
60480
7
H11002
19
6048
8
H11002
9829
3628800
q
i
H33527 H20888
i
kH335270
c
k
c
iH11002k
L
2
H33527 q
0
H11001 q
1
H11612H11001q
2
H11612
2
H11001 ...
L
2
r
n
H33527 2r
nH110021
H11002 r
nH110022
H11001 h
2
L
2
r?
n
hLr?
n
H33527 r
n
H11002 r
nH110021
hL
H20885
tn
t
n
H11002h
r?dt H33527 hLH20849r?
n
H11002 r?
nH110021
H20850 H33527 h
2
L
2
r?
n
r?
nH110011
H33527 h
H20873
H11612
H110021
H11001
1
2
H11001
5
12
H11612H11001
3
8
H11612
2
H11001
251
720
H11612
3
H11001
95
288
H11612
4
H11001 ...
H20874
r?
n
340 Berry and Healy
shows that the predictor coefficients can be written as a sum of the corrector
coefficients, just as with single integration. The predictor coefficients can be com-
puted from (44)
(48)
Using these coefficients in (46)
(49)
gives the St?rmer predictor formula [5].
Gauss-Jackson
The Gauss-Jackson method is the summed form of the St?rmer-Cowell method.
According to Herrick ([10], p. 12), the method was named Gauss-Jackson because
of its appearance in a 1924 paper by Jackson [1]. The Gauss-Jackson method is also
referred to as the second-sum method. Sometimes, the name is associated with a
central difference derivation, because that is how Jackson showed it. However, as
we mentioned, the derivation by any of the difference operators are equivalent;
what is unique about Gauss-Jackson is the presence of the second sum.
The Gauss-Jackson corrector can be derived from the Cowell corrector by not-
ing that the position terms in (45) may be combined using a second backward
difference
(50)
Applying the second sum to both sides gives an equation for position [5]
(51)
Note that , so the first two terms can be simplified to
(52)
This expression is the Gauss-Jackson corrector formula.
A similar process on the St?rmer predictor (49) gives a summed predictor
(53)
which is the Gauss-Jackson predictor formula. Note that for the predictor the sec-
ond sum operators acts on , while for the corrector it acts on .
Startup Formulas
In order to calculate the differences , the difference table (Table 1) must be cal-
culated, which depends on having values of acceleration at the backpoints. To
calculate the difference points must be known. The initial value problemN H11001 1N
th
H11612
i
r?
nH110021
r?
n
r
nH110011
H33527 h
2
H20875
H11612
H110022
r?
n
H11001
H20873
1
12
H11001
1
12
H11612H11001
19
240
H11612
2
H11001
3
40
H11612
3
H11001 ...
H20874
r?
n
H20876
r
n
H33527 h
2
H20875
H11612
H110022
r?
nH110021
H11001
H20873
1
12
H11002
1
240
H11612
2
H11002
1
240
H11612
3
H11002
221
60480
H11612
4
H11001 ...
H20874
r?
n
H20876
H208491 H11002H11612H20850r?
n
H33527 E
H110021
r?
n
H33527 r?
nH110021
r
n
H33527 h
2
H11612
H110022
H20873
1 H11002H11612H11001
1
12
H11612
2
H11002
1
240
H11612
4
H11002
1
240
H11612
5
H11002
221
60480
H11612
6
H11001 ...
H20874
r?
n
H11612
H110022
H11612
2
r
n
H33527 h
2
H20873
1 H11002H11612H11001
1
12
H11612
2
H11002
1
240
H11612
4
H11002
1
240
H11612
5
H11002
221
60480
H11612
6
H11001 ...
H20874
r?
n
r
nH110011
H33527 2r
n
H11002 r
nH110021
H11001 h
2
H20873
1 H11001
1
12
H11612
2
H11001
1
12
H11612
3
H11001
19
240
H11612
4
H11001
3
40
H11612
5
H11001 ...
H20874
r?
n
i
H9261
i
0
1
1
0
2
1
12
3
1
12
4
19
240
5
3
40
6
863
12096
7
275
4032
8
33953
518400
H9261
i
Implementation of Gauss-Jackson Integration for Orbit Propagation 341
only gives the initial conditions at epoch, so a startup procedure is required to
find position and velocity, and then acceleration, at the other backpoints. One
possible method is to use a single-step integrator, such as Runge-Kutta, to find
the values at the initial set of backpoints. Alternatively, one may use the two-body
motion solution, which may be computed analytically. In the latter case particu-
larly, we use an iterative method that takes an initial guess of the backpoints and
refines them with corrector formulas. These corrector formulas, which correct a
value using not only previous points, but also future known points, may be called
mid-corrector formulas.
Where two-body motion is close to the actual motion modeled, the initial set of
guess values can be found by using the analytic two-body solution. In other cases,
a good analytic approximation of the full motion would be necessary. For a method
that uses the difference, total points are needed, so N points in addition
to epoch must be found. To reduce the error that comes from the two-body solution,
instead of propagating N steps forward, NH208622 steps are taken both forward and back-
ward. These points are numbered using a symmetric index, from H11002NH208622 to NH208622,
with epoch numbered 0. This system requires that N be even, as it is for the inte-
grators in this study. If N were odd an additional point would be needed either be-
fore or after epoch.
After the set of guess values is found, each value other than epoch is corrected
with a mid-corrector formula, except for the last value, which uses the corrector
formula. The formulas are numbered using the letter j as an index with symmetric
index numbers. So the mid-correctors are numbered to ,
and the corrector is numbered . Similarly the predictor can be considered
the formula.
The mid-correctors are found from the corrector by applying the inverse dis-
placement operator repeatedly. Using the corrector operator L as an
example, the relation of operators to coefficient sets can be seen on the diagram for
an eighth-order method ( ):
mid-correctors
corrector L
predictor
The coefficients for each mid-corrector may be computed from the next one by 
application of the operator . For an arbitrary power series in H11612 with coeffi-
cients d, multiply the series
(54)
In other words, the coefficients for a particular value of j are just the differences of
coefficients of adjacent values of powers of H11612 for the next higher value of j. As
shown in the derivation of the Adams and St?rmer-Cowell predictors from their
H20849d
0
H11001 d
1
H11612H11001d
2
H11612
2
H11001 ...H20850H208491 H11002H11612H20850 H33527 d
0
H11001 H20849d
1
H11002 d
0
H20850H11612H11001H20849d
2
H11002 d
1
H20850H11612
2
H11001 ...
1 H11002H11612
L
1 H11002H11612
j H33527 5
j H33527 4
H208491 H11002H11612H20850Lj H33527 3
bapply H208491 H11002H11612H20850
H110021
H208491 H11002H11612H20850
2
Lj H33527 2
aapply 1 H11002H11612H11080
H11080
H11080
H11080
H11080
H11080
H208491 H11002H11612H20850
7
Lj H33527 H110023
H208491 H11002H11612H20850
8
Lj H33527 H110024
N H33527 8
E
H110021
H33527 1 H11002H11612
j H33527 NH208622 H11001 1
j H33527 NH208622
j H33527 NH208622 H11002 1j H33527 H11002NH208622
N H11001 1N
th
342 Berry and Healy
?
?
?
?
?
?
?
correctors, the coefficients are also the sums of all the coefficients of the next lower
j, through the particular power of H11612
(55)
All the coefficients, mid-correctors, corrector, and predictor, taken together may
be viewed as an array, with the rows indexed by j, and the columns by the expo-
nent i of H11612. For an -order integrator, there are terms in the series expan-
sion of the operator for each j, and there are values of j; therefore, the full
array of coefficients has elements. For example, the eighth-order
integrator has coefficients.
The summed Adams corrector formula (38) can be written using the coefficient
array 
(56)
From this, the last mid-corrector formula is obtained from the coeffi-
cient differences (54) of the corrector
(57)
The general mid-corrector formulas can be written as
(58)
where . The coefficients for each value of j are obtained by
taking differences of the coefficients for . The mid-corrector coefficients are
thus computed recursively
(59)
for .
The predictor formula (39) may be rewritten in terms of the coefficient array 
(60)
The eighth-order coefficients are presented in Table 3.H9252
ji
r?
nH110011
H33527 h
H20873
H11612
H110021
r?
n
H11001 H20888
N
iH335270
H9252
N
2
H110011, i
H11612
i
r?
n
H20874
H9252
H11002
N
2
H11349 j H11349
N
2
H11002 1
H9252
ji
H33527
H20877
H9252
jH110011, i
H11002 H9252
jH110011, iH110021
H9252
jH110011, i
if i H11022 0
if i H33527 0
j H11001 1
H9252
ji
H11002
N
2
H11349 n H11349
N
2
H11002 1
r?
n
H33527 h
H20873
H11612
H110021
r?
n
H11001 H20888
N
iH335270
H9252
ni
H11612
i
r?
N
2
H20874
 H33527 h
H20875
H11612
H110021
r?
N
2
H110021
H11001
H20873
H11002
1
2
H11001
5
12
H11612H11001
1
24
H11612
2
H11001
11
720
H11612
3
H11001
11
1440
H11612
4
H11001 ...
H20874
r?
N
2
H20876
 H33527 E
H110021
h
H20875
H11612
H110021
r?
N
2
H11001
H20873
H11002
1
2
H11002
1
12
H11612H11002
1
24
H11612
2
H11002
19
720
H11612
3
H11002
3
160
H11612
4
H11001 ...
H20874
r?
N
2
H20876
 r?
N
2
H110021
H33527 E
H110021
r?
N
2
H20849 j H33527
N
2
H11002 1H20850
r?
n
H33527 h
H20873
H11612
H110021
r?
n
H11001 H20888
N
iH335270
H9252N
2
, i
H11612
i
r?
n
H20874
H9252
9 H11003 10 H33527 90
H20849N H11001 1H20850H20849N H11001 2H20850
N H11001 2
N H11001 1N
th
H33527 H20888
k
H20873
H20888
k
iH335270
m
i
H20874
H11612
k
H33527 H20849m
0
H11001 m
1
H11612H11001m
2
H11612
2
H11001 ...H20850H208491 H11001H11612H11001H11612
2
H11001H11612
3
H11001 ...H20850
H20849m
0
H11001 m
1
H11612H11001m
2
H11612
2
H11001 ...H20850
1
1 H11002H11612
Implementation of Gauss-Jackson Integration for Orbit Propagation 343
Analogously, the Gauss-Jackson corrector (52) can be written using a coefficient
array 
(61)
The mid-correctors are then written
(62)
where . The coefficients of the mid-correctors are computed as
the difference (54) of the corrector coefficients; the entire set of coefficients are
thus computed recursively
(63)
for . Finally the predictor is included in the array
(64)
The eighth-order coefficients are presented in Table 4.
Note that the sum introduced by the term in the mid-correctors or corrector
is accumulated through the point before the point being corrected; this is given by
the term in (52).r?
nH110021
H11612
H110022
H9251
ji
r
nH110011
H33527 h
2
H20873
H11612
H110022
r?
n
H11001 H20888
N
iH335270
H9251
N
2
H110011, i
H11612
i
r?
n
H20874
H11002
N
2
H11349 j H11349
N
2
H11002 1
H9251
ji
H33527
H20877
H9251
jH110011, i
H11002 H9251
jH110011, iH110021
H9251
jH110011, i
if i H11022 0
if i H33527 0
H11002
N
2
H11349 n H11349
N
2
H11002 1
r
n
H33527 h
2
H20873
H11612
H110022
r?
nH110021
H11001 H20888
N
iH335270
H9251
ni
H11612
i
r?
N
2
H20874
r
n
H33527 h
2
H20873
H11612
H110022
r?
nH110021
H11001 H20888
N
iH335270
H9251
N
2
, i
H11612
i
r?
n
H20874
H9251
344 Berry and Healy
TABLE 3. Eighth-Order Summed-Adams Difference Coefficients, H9252
ji
i
01 2 3 4 5 6 7 8
25713
89600
1070017
3628800
5257
17280
19087
60480
95
288
251
720
3
8
5
12
1
2
5
H11002
8183
1036800
H11002
33953
3628800
H11002
275
24192
H11002
863
60480
H11002
3
160
H11002
19
720
H11002
1
24
H11002
1
12
H11002
1
2
4
425
290304
7297
3628800
13
4480
271
60480
11
1440
11
720
1
24
5
12
H11002
1
2
3
H11002
7
12800
H11002
3233
3628800
H11002
191
120960
H11002
191
60480
H11002
11
1440
H11002
19
720
H11002
3
8
11
12
H11002
1
2
2
2497
7257600
2497
3628800
191
120960
271
60480
3
160
251
720
H11002
31
24
17
12
H11002
1
2
1
H11002
2497
7257600
H11002
3233
3628800
H11002
13
4480
H11002
863
60480
H11002
95
288
1181
720
H11002
65
24
23
12
H11002
1
2
0
7
12800
7297
3628800
275
24192
19087
60480
H11002
2837
1440
3131
720
H11002
37
8
29
12
H11002
1
2
H110021
H11002
425
290304
H11002
33953
3628800
H11002
5257
17280
138241
60480
H11002
1011
160
6461
720
H11002
169
24
35
12
H11002
1
2
H110022
8183
1036800
1070017
3628800
H11002
11603
4480
520399
60480
H11002
22021
1440
11531
720
H11002
239
24
41
12
H11002
1
2
H110023
H11002
25713
89600
10468447
3628800
H11002
1354079
120960
1445281
60480
H11002
45083
1440
18701
720
H11002
107
8
47
12
H11002
1
2
H110024
j
Ordinate Forms
The formulas given above all involve multiple differences and sums. In order to
compute these formulas to a particular order, all the elements of the appropriate
difference tables need to be computed. Calculating the differences can be time-
consuming and is unnecessary because the formulas can be re-expressed in terms
of the function values themselves.
As shown above in equation (3), powers of the difference operators can be re-
duced to linear combinations of the original function values. The backward differ-
ence operator raised to any power can be expressed as a linear combination of the
function value at the mesh points using binomial coefficients
(65)
for . A sum over arbitrary coefficients of the backward difference operator
can thus be written in terms of the values of the function
(66)
 H33527 H20888
N
mH335270
z
Nm
r?
nH11002m
 H33527 H20888
N
mH335270
H20849H110021H20850
m
r?
nH11002m H20888
N
iH33527m
H9256
i
H20873
i
m
H20874
 H20888
N
iH335270
H9256
i
H11612
i
r?
n
H33527 H20888
N
iH335270
H9256
i H20888
i
mH335270
H20873
i
m
H20874
H20849H110021H20850
m
r?
nH11002m
H9256
i
i H11350 0
H11612
i
r?
n
H33527 H20888
i
mH335270
H20873
i
m
H20874
H20849H110021H20850
m
r?
nH11002m
Implementation of Gauss-Jackson Integration for Orbit Propagation 345
TABLE 4. Eighth-Order Gauss-Jackson Difference Coefficients, H9251
ji
i
012 3 4 5 6 7 8
3250433
53222400
8183
129600
33953
518400
275
4032
863
12096
3
40
19
240
1
12
1
12
5
H11002
330157
159667200
H11002
407
172800
H11002
9829
3628800
H11002
19
6048
H11002
221
60480
H11002
1
240
H11002
1
240
0
1
12
4
45911
159667200
641
1814400
1571
3628800
31
60480
31
60480
0H11002
1
240
H11002
1
12
1
12
3
H11002
3499
53222400
H11002
289
3628800
H11002
289
3628800
0
31
60480
1
240
19
240
H11002
1
6
1
12
2
317
22809600
0H11002
289
3628800
H11002
31
60480
H11002
221
60480
H11002
3
40
59
240
H11002
1
4
1
12
1
317
22809600
289
3628800
1571
3628800
19
6048
863
12096
H11002
77
240
119
240
H11002
1
3
1
12
0
H11002
3499
53222400
H11002
641
1814400
H11002
9829
3628800
H11002
275
4032
23719
60480
H11002
49
60
199
240
H11002
5
12
1
12
H110021
45911
159667200
407
172800
33953
518400
H11002
6961
15120
73111
60480
H11002
79
48
299
240
H11002
1
2
1
12
H110022
H11002
330157
159667200
H11002
8183
129600
1908311
3628800
H11002
20191
12096
172651
60480
H11002
347
120
419
240
H11002
7
12
1
12
H110023
3250433
53222400
H11002
427487
725760
7965611
3628800
H11002
45601
10080
347539
60480
H11002
371
80
559
240
H11002
2
3
1
12
H110024
j
with the ordinate coefficients defined as
(67)
Notice that these coefficients, unlike the difference coefficients , depend on N,the
order of the integration method. Thus a variation of order to control errors while in-
tegrating is feasible in difference form, but is more difficult in ordinate form. Here
the index m numbers the backpoints so the current point is and each previ-
ous point has a higher positive value of m.
As an example, consider the summed Adams corrector (38). For a fourth-order
method, , the difference coefficients and the ordinate coefficients are
(68)
An alternative formulation includes the term with the sum rather than with
the coefficients. The example summed Adams-Moulton corrector (38) is now written
(69)
With this formulation the coefficients are different from (68)
(70)
All the single-integration corrector and mid-corrector coefficients are computed
this way; however, the term is included in the coefficients for the predictor. The
eighth-order summed Adams coefficients in ordinate form are shown in Table 5,
and the Gauss-Jackson coefficients are shown in Table 6. For Gauss-Jackson, the
term , is included with the coefficients and not the sum . Although this
alternate formulation is not necessary, it makes the programming easier, as described
in the next section. Both the Gauss-Jackson and the summed Adams coefficients
may be compared with [6], Table III (pages 25?26) and [12] Appendix C.3. Ordi-
nate forms for the Adams, summed Adams, St?rmer-Cowell, and Gauss-Jackson
predictors and correctors for orders between 1 and 14 may also be found in [5].
In this section, the are generic difference coefficients such as and , and the
z are the corresponding ordinate coefficients, which have an extra subscript to in-
dicate the order. The actual eighth-order ordinate coefficients a or b do not have an
order subscript even though they are dependent on the order; the two subscripts are
the point of integration j and the point of reference k.
H9252H9251H9256
H11612
H110022
1
12
r?
n
H11612
0
a
jk
b
jk
H11612
0
m
H9256
i
z
4m
0
0
H11002
49
288
1
H11002
1
12
77
240
2
H11002
1
24
H11002
7
30
3
H11002
19
720
73
720
4
H11002
3
160
H11002
3
160
m H33527 0
r?
n
H33527 h
H20875
H11612
H110021
r?
n
H11002
1
2
r?
n
H11001
H20873
H11002
1
12
H11612H11002
1
24
H11612
2
H11002
19
720
H11612
3
H11002
3
160
H11612
4
H11001 ...
H20874
r?
n
H20876
H11002
1
2
r?
n
H11612
0
m
H9256
i
z
4m
0
H11002
1
2
H11002
193
288
1
H11002
1
12
77
240
2
H11002
1
24
H11002
7
30
3
H11002
19
720
73
720
4
H11002
3
160
H11002
3
160
z
4m
H9256
i
N H33527 4
m H33527 0
H9256
z
Nm
H33527 H20849H110021H20850
m
H20888
N
iH33527m
H9256
i
H20873
i
m
H20874
346 Berry and Healy
Implementation of Gauss-Jackson Integration for Orbit Propagation 347
T
ABLE 5.
Eighth-Order Summed-Adams Coeff
icients in Ordinate F
o
rm,
b
jk
k
j
01
2
3
4
3288521 1036800
H11002
40987771
3628800
10219841
403200
H11002
135352319
3628800
167287
4536
H11002
9839609
403200
5393233
518400
H11002
9401029 3628800
25713 89600
5
H11002
19087 89600
427487 725760
H11002
3498217 3628800
500327 403200
H11002
6467 5670
2616161 3628800
H11002
24019 80640
263077
3628800
H11002
8183
1036800
4
H11002
8183
1036800
H11002
57251
403200
1106377 3628800
H11002
218483 725760
69
280
H11002
530177
3628800
210359
3628800
H11002
5533
403200
425
290304
3
425
290304
H11002
76453
3628800
H11002
5143
57600
660127
3628800
H11002
661
5670
4997
80640
H11002
83927
3628800
19109
3628800
H11002
7
12800
2
H11002
7
12800
23173
3628800
H11002
29579
725760
H11002
2497
57600
2563
22680
H11002
172993
3628800
6463
403200
H11002
2497
725760
2497
7257600
1
2497
7257600
H11002
1469
403200
68119
3628800
H11002
252769
3628800
0
252769
3628800
H11002
68119
3628800
1469
403200
H11002
2497
7257600
0
H11002
2497
7257600
2497
725760
H11002
6463
403200
172993
3628800
H11002
2563
22680
2497
57600
29579
725760
H11002
23173
3628800
7
12800
H11002
1
7
12800
H11002
19109
3628800
83927
3628800
H11002
4997
80640
661
5670
H11002
660127
3628800
5143
57600
76453
3628800
H11002
425
290304
H11002
2
H11002
425
290304
5533
403200
H11002
210359
3628800
530177
3628800
H11002
69
280
218483 725760
H11002
1106377 3628800
57251
403200
8183
1036800
H11002
3
8183
1036800
H11002
263077
3628800
24019 80640
H11002
2616161 3628800
6467 5670
H11002
500327 403200
3498217 3628800
H11002
427487 725760
19087 89600
H11002
4
H11002
1
H11002
2
H11002
3
H11002
4
348 Berry and Healy
T
ABLE 6.
Eighth-Order Gauss-J
ackson Coeff
icients in Ordinate F
o
rm,
a
jk
k
j
01
2
3
4
103798439 159667200
H11002
24115843
9979200
18071351
3326400
H11002
159314453
19958400
25162927
3193344
H11002
8660609 1663200
6322573 2851200
H11002
11011481 19958400
3250433
53222400
5
3250433
53222400
572741
5702400
H11002
8701681
39916800
4026311
13305600
H11002
917039
3193344
7370669
39916800
H11002
1025779
13305600
754331
39916800
H11002
330157
159667200
4
H11002
330157
159667200
530113
6652800
518887
19958400
H11002
27631
623700
44773
1064448
H11002
531521
19958400
109343
9979200
H11002
1261
475200
45911
159667200
3
45911
159667200
H11002
185839
39916800
171137
1900800
73643
39916800
H11002
25775
3193344
77597
13305600
H11002
98911
39916800
24173
39916800
H11002
3499
53222400
2
H11002
3499
53222400
4387
4989600
H11002
35039
4989600
90817
950400
H11002
20561
3193344
2117
9979200
2059
6652800
H11002
317
2851200
317
22809600
1
317
22809600
H11002
2539
13305600
55067
39916800
H11002
326911
39916800
14797
152064
H11002
326911
39916800
55067
39916800
H11002
2539
13305600
317
22809600
0
317
22809600
H11002
317
2851200
2059
6652800
2117
9979200
H11002
20561
3193344
90817
950400
H11002
35039
4989600
4387
4989600
H11002
3499
53222400
H11002
1
H11002
3499
53222400
24173
39916800
H11002
98911
39916800
77597
13305600
H11002
25775
3193344
73643
39916800
171137
1900800
H11002
185839
39916800
45911
159667200
H11002
2
45911
159667200
H11002
1261
475200
109343
9979200
H11002
531521
19958400
44773
1064448
H11002
27631
623700
518887
19958400
530113
6652800
H11002
330157
159667200
H11002
3
H11002
330157
159667200
754331
39916800
H11002
1025779
13305600
7370669
39916800
H11002
917039
3193344
4026311
13305600
H11002
8701681
39916800
572741
5702400
3250433
53222400
H11002
4
H11002
1
H11002
2
H11002
3
H11002
4
Implementation
To find both position and velocity given a function for acceleration, either 
a single-integration method must be used once for velocity and the method 
used again for position, or a single-integration method must be used along with
a double-integration method. Assuming that summed and nonsummed forms are
not used together, the four integrators derived above offer four implementations:
two Adams, two summed Adams, Adams and St?rmer-Cowell, or summed
Adams and Gauss-Jackson integration. For convenience, Gauss-Jackson integra-
tion is understood to mean Gauss-Jackson and summed Adams integration, and
St?rmer-Cowell integration is understood to mean St?rmer-Cowell and Adams
integration. Because Gauss-Jackson integration involves sum terms, its imple-
mentation is somewhat complicated, and so is described in detail in this section.
Overview
In order to use the predictor-corrector formulas to integrate forward in time with
an order method, points must already be known. To preserve the order
of the method, these points must have also been found using an order
method. Since initially only one point, epoch, is known, a startup procedure is nec-
essary. If initial estimated points are given, the mid-corrector and corrector
formulas can be used to refine them. By applying the mid-correctors iteratively
until the points converge, the resulting points are accurate through order.
For the eighth-order methods used in this study, a Taylor expansion of the two-
body solution to fifth order ( f and g series, [14] equation (4?68)) is used for the
initial estimate, though a single-step integrator, such as Runge-Kutta could be used.
The analytic solution generates eight positions and velocities in addition to epoch:
four before and four after epoch. The accelerations at all nine points are then found
from the positions and velocities with full perturbations, and corrector and mid-
corrector formulas are applied to the eight non-epoch accelerations, which gener-
ates corrected positions and velocities. These corrected positions and velocities are
used to find corrected accelerations, and the cycle repeats, until the accelerations
converge, i.e., on two successive iterations, the absolute difference is less than a
specified tolerance. This process is denoted SECECE...CE, where ?S? is startup
estimate, ?E? is evaluate, and ?C? is correct. This process is done for the eight ini-
tial points other than epoch.
Once the startup points have been corrected to satisfactory accuracy, the regular
predictor-corrector process can start. First the predictor equation is used to find the
position and velocity of the first point after startup, and the force model is evaluated
to find the acceleration of that point. This predicted acceleration is then used in the
corrector formula to find a corrected position and velocity. The corrected position
and velocity may be compared to the predicted values, and if they do not match to
some tolerance, another acceleration is evaluated, which is then used in the correc-
tor. Note that the procedure is different from startup, where convergence of acceler-
ation is checked. This cycle would then be repeated until the position and velocity
converge, or until some specified maximum number of iterations is reached. This
process gives a implementation, with a maximum value of n specified.
Single Integration
The first step in startup of single integration is to use the mid-corrector formu-
las to correct the eight estimated velocities around epoch. In difference form, these
PH20849ECH20850
n
N
th
N H11001 1
N
th
N H11001 1
N H11001 1N
th
Implementation of Gauss-Jackson Integration for Orbit Propagation 349
formulas are (58); using the ordinate coefficients b, the formulas are
(71)
with . Note the term is excluded from the coefficient sum in an-
ticipation of its inclusion in the acceleration term, as in (69). The integration con-
stant (see Tables and Operators) may be determined by making sure that the
initial conditions ( )
(72)
are satisfied; solving for 
(73)
The acceleration term may be included with the sum to make software imple-
mentation easier. Define the term , which replaces the first two 
terms in (71), and define a new integration constant Then (71) may 
be rewritten
(74)
where n ranges from H110024 to 4, with
(75)
Table 7 shows the relationship between the summed accelerations and near epoch.
For a given value of n, notice the symmetry of the formulas in the last column;
changing the sign of n merely changes the sign of the term added to . This for-
mulation makes programming simpler and is the reason that the term is
moved from the coefficients to the sum. The are computed successively in order
to compute the startup velocities. While the formulas are correct for , no cor-
rection is performed on epoch.
n H33527 0
s
n
H11002r?
0
H208622
CH11032
1
H11032
s
n
s
n
H33527
s
nH110021
H11001
r?
nH110021
H11001 r?
n
2
CH11032
1
H11032
s
nH110011
H11002
r?
nH110011
H11001 r?
n
2
if n H11022 0
if n H33527 0
if n H11021 0
r?
n
H33527 h
H20873
s
n
H11001 H20888
4
kH33527H110024
b
nk
r?
k
H20874
CH11032
1
H11032
H33527 C
1
H11002
r?
0
2
.
s
n
H33527 H11612
H110021
r?
n
H11002
1
2
r?
n
C
1
H33527
r?
0
h
H11002 H20888
4
kH33527H110024
b
0k
r?
k
H11001
r?
0
2
C
1
H33527 H11612
H110021
r?
0
r?
0
H33527 h
H20873
H11612
H110021
r?
0
H11002
r?
0
2
H11001 H20888
4
kH33527H110024
b
0k
r?
k
H20874
n H33527 0
C
1
H11002
r?
n
2
H110024 H11349 n H11349 4
r?
n
H33527 h
H20873
H11612
H110021
r?
n
H11002
r?
n
2
H11001 H20888
4
kH33527H110024
b
nk
r?
k
H20874
350 Berry and Healy
?
?
?
?
?
?
?
TABLE 7. Relationship Between Summation Term and s
n
0
1
2 CH110321 H11001
1
2
r?0 H11001 r?1 H11001
1
2
r?2C1 H11001 r?1 H11001 r?2r? 2
CH110321 H11001
1
2
r?0 H11001
1
2
r?1C1 H11001 r?1r? 1
CH110321C1r? 0
CH110321 H11002
1
2
r?0 H11002
1
2
r?H110021C1 H11002 r?0r?H110021H110021
CH11032
1 H11002
1
2
r?0 H11002 r?H110021 H11002
1
2
r?H110022C1 H11002 r?0 H11002 r?H110021r?H110022H110022
s
n H33527 H11612
H110021
r?n H11002
1
2
r?nH11612
H110021
r?nr?nn
Once the startup points have been corrected, integration may proceed forward
with application of the predictor, followed by one or more applications of the cor-
rector. The predictor finds the velocity from the coefficients in row 5 of the
b array, and the nine most recent accelerations . The formula also includes
a term which is the sum of all the accelerations since epoch and 
(76)
with . This expression differs from the difference form given in (60). For
programming, the predictor can be written using 
(77)
The Adams-Moulton corrector formula in difference form (69), for any point
after startup, can be rewritten using the eighth-order ordinate coefficients b.
These coefficients are applied to the eight previous accelerations and the current
predicted acceleration
(78)
with . Note that (78) is the same as (71) for . The corrector can be writ-
ten using 
(79)
Depending on the implementation, this equation may be applied through several
evaluate and correct (EC) iterations, but only the last acceleration, , changes,
so the summation for the first eight need only be done once, then saved through
the iterations.
Note that when applying the predictor, the first formula in (75) uses one half
of the corrected acceleration and one half of the predicted acceleration in the
sum. If the acceleration term is not shifted, then the entire acceleration comes
from the predicted value. Therefore, the shift has a presumed advantage of using
half the corrected acceleration and half the predicted, rather than all predicted
acceleration. A possible further improvement would be to use instead of
in (78) with b redefined appropriately. Then this summation would have
only corrected values.
Double Integration
The mid-correctors, predictor, and corrector for the second integration ordinate
form are respectively
(80)
(81)
(82)n H11350 4 r
n
H33527 h
2
H20873
H11612
H110022
r?
nH110021
H11001 H20888
4
kH33527H110024
a
4k
r?
nH11001kH110024
H20874
n H11350 4 r
nH110011
H33527 h
2
H20873
H11612
H110022
r?
n
H11001 H20888
4
kH33527H110024
a
5k
r?
nH11001kH110024
H20874
H110024 H11349 n H11349 4 r
n
H33527 h
2
H20873
H11612
H110022
r?
nH110021
H11001 H20888
4
kH33527H110024
a
nk
r?
k
H20874
H11612
H110021
r?
n
H11612
H110021
r?
nH110021
k H33527 4
r?
n
H33527 h
H20873
s
n
H11001 H20888
4
kH33527H110024
b
4k
r?
nH11001kH110024
H20874
s
n
n H33527 4n H11350 4
r?
n
H33527 h
H20873
H11612
H110021
r?
n
H11002
r?
n
2
H11001 H20888
4
kH33527H110024
b
4k
r?
nH11001kH110024
H20874
n H11022 4
r?
nH110011
H33527 h
H20873
s
n
H11001
r?
n
2
H11001 H20888
4
kH33527H110024
b
5k
r?
nH11001kH110024
H20874
s
n
n H11350 4
r?
nH110011
H33527 h
H20873
H11612
H110021
r?
n
H11001 H20888
4
kH33527H110024
b
5k
r?
nH11001kH110024
H20874
C
1
H11612
H110021
r?
nH110028
...r?
n
r?
nH110011
Implementation of Gauss-Jackson Integration for Orbit Propagation 351
corresponding to the difference form (62), (64), (61). To find the integration con-
stant for double integration, , use in (80)
(83)
Table 2 shows that the second sum on the point immediately before epoch is the dif-
ference of the two constants of integration
(84)
which than allows us to solve for , because is known (73)
(85)
For software implementation, a term is defined recursively
(86)
The mid-correctors, predictor, and corrector can now be written in terms of 
(87)
(88)
(89)
As with first integration, the summation in the corrector may be split, because the
first eight accelerations do not change during the EC iteration.
Procedure
A step-by-step procedure for the operation of the integrator can now be written:
Startup
11. Use f and g series to calculate eight positions and velocities surrounding epoch
labelled with epoch.
12. Evaluate nine accelerations from these positions and velocities, and those of
epoch labelled 
13. While the accelerations have not converged:
(a) Calculate (75) and (86).
(b) For each point , :
i. Calculate (75) and (86).
ii. Calculate using Table 5.
iii. Calculate using Table 6.H20888
4
kH33527H110024
a
nk
r?
k
H20888
4
kH33527H110024
b
nk
r?
k
S
n
s
n
n HS33527 0n H33527 H110024...4
S
0
s
0
r?
H110024
...r?
4
.
r
0
r
H110024
...r
4
n H11350 4 r
n
H33527 h
2
H20873
S
n
H11001 H20888
4
kH33527H110024
a
4k
r?
nH11001kH110024
H20874
n H11350 4 r
nH110011
H33527 h
2
H20873
S
nH110011
H11001 H20888
4
kH33527H110024
a
5k
r?
nH11001kH110024
H20874
H110024 H11349 n H11349 4 r
n
H33527 h
2
H20873
S
n
H11001 H20888
4
kH33527H110024
a
nk
r?
k
H20874
S
n
S
n
H33527
S
nH110021
H11001 s
nH110021
H11001
r?
nH110021
2
C
2
H11002 C
1
S
nH110011
H11002 s
nH110011
H11001
r?
nH110011
2
if n H11022 0
if n H33527 0
if n H11021 0
S
n
H33527 H11612
H110022
r?
nH110021
C
2
H33527
r
0
h
2
H11002 H20888
4
kH33527H110024
a
0k
r?
k
H11001 C
1
C
1
C
2
H11612
H110022
r?
H110021
H33527 C
2
H11002 C
1
r
0
H33527 h
2
H20873
H11612
H110022
r?
H110021
H11001 H20888
4
kH33527H110024
a
0k
r?
k
H20874
n H33527 0C
2
352 Berry and Healy
?
?
?
?
?
?
?
iv. Calculate (74) and (87).
v. Evaluate the updated acceleration using the appropriate force model.
(c) Test convergence of the accelerations.
Predict
14. Calculate (86).
15. Calculate .
16. Calculate .
17. Calculate (77) and (88).
Evaluate?Correct
18. Evaluate the acceleration .
19. Increment n.
10. While and have not converged and the maximum number of corrector it-
erations is not exceeded:
(a) Calculate (75).
(b) If first iteration:
i. Calculate .
ii. Calculate .
(c) Calculate and .
(d) Sum the results of 10(a), 10(b), and 10(c) to obtain (79) and (89).
(e) Test convergence of and ; if not converged, and this is not the last
allowable corrector iteration, evaluate .
11. Go to Step 4 unless this is the final time desired.
Similar procedures are used for implementing two Adams, two summed Adams,
or Adams with St?rmer-Cowell.
Effects of Step Size and Order on Accuracy and Stability
In general, the accuracy of a numerical integrator improves as the step size is
reduced and the order of the method is increased. One exception to this rule is
that round-off error can occur if the step size becomes too small. Another ex-
ception occurs when the integrator becomes unstable. Increasing the order may
cause instability because higher-order methods are susceptible to instability at
large step sizes [15] (pp. 139, 146).
To assess the effect of step size and order on accuracy, we compare ephemeris
generated using step sizes of 30, 60, 120, and 240 seconds, with 6th, 8th, 10th, 12th,
and 14th order Gauss-Jackson integrators. The integrators follow the procedure
given in the Procedure section, modified for the order of the integrator. The tests
compare ephemeris generated at one minute intervals over 72 hours. For step sizes
greater than one minute, a fifth-order Hermitian interpolator ([16], p. 28) is used to
find the intermediate points. The perturbation forces include WGS-84
geopotential [17], Jacchia 70 [18] drag, and lunar and solar forces. Note that all of
these perturbations are continuous forces. To simplify this study, discontinuous
forces such as solar radiation pressure are not considered. The two test cases used
are the International Space Station (ISS, catalog number 25544), with a period
24 H11003 24
r?
n
r
n
r?
n
r
n
r?
n
a
4k
r?
n
b
4k
r?
n
H20888
3
kH33527H110024
a
4k
r?
nH11001kH110024
H20888
3
kH33527H110024
b
4k
r?
nH11001kH110024
s
n
r
n
r?
n
r?
nH110011
r
nH110011
r?
nH110011
H20888
4
kH33527H110024
a
5k
r?
nH11001kH110024
H20888
4
kH33527H110024
b
5k
r?
nH11001kH110024
S
nH110011
r?
n
r
n
r?
n
Implementation of Gauss-Jackson Integration for Orbit Propagation 353
of 92.05 min. and eccentricity of 0.001, and CRRES (20712), with a period of
607.28 min. and eccentricity of 0.716.
In the tests, a metric for integration accuracy is defined using an error ratio 
defined in terms of the RMS error of the integration [19]. First define position
errors as
(90)
The RMS position error can be calculated
(91)
The RMS position error is normalized by the apogee distance and the number of
orbits to find the position error ratio
(92)
The error ratios for ISS are given in Table 8, and the error ratios for CRRES are
given in Table 9. In the comparison the ephemeris generated by the 14th order
method with 30 second step size is considered to be the reference. The tables show
that step size affects accuracy more than the order of the method does. The trend is
similar for both the circular and the eccentric objects. In Table 8 the asterisk de-
notes that the integration has become unstable. It is obvious when instability has
occurred because the eccentricity of the integrated ephemeris starts to grow to the
point of the orbit becoming hyperbolic. The instability at a large step size for high
orders shows that if computation time is a higher priority than accuracy, a lower
order method is preferable, because a larger step size can be used without the threat
of instability.
Herrick [10] claims that the corrector is not usually required in Gauss-Jackson
integration. This saves computation time, because only one evaluation is required
(PE) per step. A predictor-only form of Gauss-Jackson can be created by skip-
ping steps 10(b)?10(e) of the procedure above. Tables 10 and 11 list error ratios for
predictor-only results for ISS and CRRES, respectively. Once again, the 14th order,
30 second step size results, with the corrector, are considered to be the reference.
Comparing Table 10 to Table 8, we see that for the circular satellite the results
with and without the corrector are similar for smaller step sizes and lower orders. 
H9267
r
H33527
H9004r
RMS
r
A
N
orbits
H9004r
RMS
H33527
H20881
1
N
H20888
N
iH335271
H20849H9004r
i
H20850
2
H9004r H33527 H20841r
computed
H11002 r
reference
H20841
354 Berry and Healy
TABLE 8. Error Ratios for ISS (e H11549 0.001)
Step Size
Order
(sec) 6 8 10 12 14 
?
**1.2 H11003 10
H110024
1.3 H11003 10
H110024
1.1 H11003 10
H110024
240
1.1 H11003 10
H110027
8.8 H11003 10
H110028
1.1 H11003 10
H110027
1.1 H11003 10
H110027
9.7 H11003 10
H110028
120
9.0 H11003 10
H1100210
9.7 H11003 10
H1100210
1.1 H11003 10
H110029
1.5 H11003 10
H110029
2.4 H11003 10
H110029
60
1.5 H11003 10
H1100213
3.8 H11003 10
H1100213
1.5 H11003 10
H1100212
1.0 H11003 10
H1100211
30
However, when the corrector is not used instability occurs at lower orders and
smaller step sizes than when the corrector is used. For the eccentric satellite, shown
in Table 11, there is a somewhat more noticeable loss in accuracy when the correc-
tor is not applied. However, the integration of the eccentric satellite remains stable
in more cases than the circular case.
In general, when the corrector is applied at least two evaluations are performed
(PECE). This means that a predictor-only implementation (PE) has roughly half the
computation time as a predictor-corrector implementation, because evaluations ac-
count for most of the computation time when using a realistic force model. Compu-
tation time is also related to the step size; when the step size is doubled the
computation time is halved. Comparing Table 8 to Table 10, and Table 9 to Table 11,
we see that 30 second predictor-only results are more accurate than 60 second
predictor-corrector results, though these have roughly the same computation time.
This justifies Herrick?s claim that the corrector is not necessary; however, if it is not
used, care must be taken to avoid instability. The inclusion of the corrector (PEC)
with no final evaluation gives an improvement in accuracy over the predictor-only
method, with little cost in computation time, since there are no extra evaluations.
However this method has the same stability problems as the predictor-only method.
An alternate method performs a second evaluation in which only the two-body
force is re-evaluated; the perturbation force from the first evaluation is re-used. Be-
cause the two-body evaluation is a simple calculation, this method has essentially
the same run-time as performing only one evaluation per step. However, in most
Implementation of Gauss-Jackson Integration for Orbit Propagation 355
TABLE 9. Error Ratios for CRRES (e H11549 0.716)
Step Size
Order
(sec) 6 8 10 12 14 
?
9.0 H11003 10
H110025
2.0 H11003 10
H110024
4.0 H11003 10
H110024
1.9 H11003 10
H110025
9.3 H11003 10
H110024
240
8.7 H11003 10
H110028
2.6 H11003 10
H110028
2.2 H11003 10
H110027
7.6 H11003 10
H110027
1.1 H11003 10
H110027
120
1.6 H11003 10
H1100210
1.6 H11003 10
H1100210
1.2 H11003 10
H1100210
3.9 H11003 10
H1100211
5.3 H11003 10
H110029
60
1.1 H11003 10
H1100214
9.4 H11003 10
H1100214
2.5 H11003 10
H1100213
1.4 H11003 10
H1100211
30
TABLE 10. Error Ratios for ISS, Predictor only (e H11549 0.001)
Step Size
Order
(sec) 6 8 10 12 14 
*
**
***
***4.9 H11003 10
H110024
240
1.2 H11003 10
H110027
1.3 H11003 10
H110026
120
9.2 H11003 10
H1100210
1.6 H11003 10
H110029
4.7 H11003 10
H110029
60
1.5 H11003 10
H1100212
3.5 H11003 10
H1100213
1.9 H11003 10
H1100212
1.4 H11003 10
H1100211
30
cases this method has the same stability advantage as performing two full evalua-
tions, though the partial evaluation is not as accurate as a full evaluation. Results
using a partial evaluation are available in [9].
Summary
This paper derives the Gauss-Jackson integration method widely used in special
perturbations software for space surveillance. Both single and double integration
are derived from difference tables. The ordinate form simplifies the calculations be-
cause it avoids the need to calculate tables of differences; however, it requires new
sets of coefficients. We show how the theory is actually implemented in the soft-
ware. Finally, there is a discussion of stability and accuracy which shows that step
size has more of an effect on accuracy than the order of the method does, and that
higher orders are subject to instability at large step sizes.
Coefficients may be computed to any order with code that is available at the
permanent link http://hdl.handle.net/1903/2202.
Acknowledgments
We thank Keith Atkins, Shannon Coffey, Felix Hoots, and Alan Segerman for helpful com-
ments on the original paper. Fred Lutze provided comments on this version of the paper, and we
dedicate it to his memory.
References
[1] JACKSON, J. ?Note on the Numerical Integration of ,? Monthly Notes,
Volume 84, Royal Astronomy Society, 1924, pp. 602?606.
[2] KOSKELA, P. E. ?Astrodynamic Analysis for the Advanced Orbit/Ephemeris Subsystem,?
Technical Report U-4180, Philco-Ford Corporation, Aeronutronic Division, Newport Beach,
California, September 1967.
[3] NEAL, H. L., COFFEY, S. L., and KNOWLES, S. ?Maintaining the Space Object Catalog with
Special Perturbations,? in F. Hoots, B. Kaufman, P. Cefola, and D. Spencer, editors, Astro-
dynamics 1997 Part II, Volume 97 of Advances in the Astronautical Sciences, pp. 1349?1360.
[4] CASALI, S. ?Conjunction Analysis Approach for the Astrodynamics Support Workstation
(ASW),? presented at the Core Technologies for Space Systems Conference, Colorado Springs,
Colorado, November 19?21, 2002.
[5] MAURY, J. L. and SEGAL, G. P. ?Cowell Type Numerical Integration as Applied to Satellite
Orbit Computation,? Technical Report X-553-69-46, NASA, 1969, NTIS #N6926703.
d
2
xH20862dt
2
H33527 f H20849x, tH20850
356 Berry and Healy
TABLE 11. Error Ratios for CRRES, Predictor only (e H11549 0.716)
Step Size
Order
(sec) 6 8 10 12 14 
*
**1.4 H11003 10
H110023
1.0 H11003 10
H110022
1.3 H11003 10
H110022
240
4.3 H11003 10
H110026
2.0 H11003 10
H110026
2.3 H11003 10
H110025
4.2 H11003 10
H110025
120
5.8 H11003 10
H110027
1.4 H11003 10
H110029
1.9 H11003 10
H110029
2.0 H11003 10
H110029
8.8 H11003 10
H110028
60
5.4 H11003 10
H1100213
2.1 H11003 10
H1100213
1.8 H11003 10
H1100213
6.6 H11003 10
H1100212
3.5 H11003 10
H1100210
30
[6] NORAD ?Mathematical Foundation for SCC Astrodynamic Theory,? Technical Report TP
SCC 008, Headquarters North American Aerospace Defense Command, 1982. Obtainable
from Defense Technical Information Center, Cameron Station, Alexandria, VA 22304-6145,
as AD #B081394.
[7] WOODBURN, J. ?Mitigation of the Effects of Eclipse Boundary Crossings on the Numerical
Integration of Orbit Trajectories Using an Encke Type Correction Algorithm,? presented as
paper AAS 01-223 at the AAS/AIAA Space Flight Mechanics Meeting, Santa Barbara, CA,
February 11?14, 2001.
[8] BERRY, M. and HEALY, L. ?The Generalized Sundman Transformation for Propagation of
High-Eccentricity Elliptical Orbits,? Advances in Astronautics, San Diego, CA, February 2002.
[9] BERRY, M. M. A Variable-Step Double-Integration Multi-Step Integrator, Ph.D. thesis,
Virginia Polytechnic Institute and State University, Blacksburg, Virginia, April 2004. Available
at permanent link http://scholar.lib.vt.edu/theses/available/etd-04282004-071227/.
[10] HERRICK, S. Astrodynamics: Orbit Correction, Perturbation Theory, Integration, Volume 2,
Van Nostrand Reinhold Company, New York, 1972.
[11] FRANKENA, J. F. ?St?rmer-Cowell: Straight, Summed and Split. An Overview,? Journal of
Computational and Applied Mathematics, 62: 129?154, 1995.
[12] WALLNER, R. N. ?Earth Gravitational Error Budget Initial Report,? Technical report, Kaman
Sciences Corporation, 1500 Garden of the Gods Road, Colorado Springs, CO 80907, November
1994, System Memorandum KSPACE 95-18.
[13] HENRICI, P. Discrete Variable Methods in Ordinary Differential Equations, John Wiley and
Sons, New York, 1962.
[14] VALLADO, D. A. Fundamentals of Astrodynamics and Applications, Volume 12 of The Space
Technology Library, Microcosm Press and Kluwer Academic Publishers, El Segundo, Califor-
nia and Dordrecht, The Netherlands, second edition, 2001.
[15] MONTENBRUCK, O. and GILL, E. Satellite Orbits, Springer, New York, 2000.
[16] DAVIS, P. J. Interpolation and Approximation, Dover, New York, 1975.
[17] National Imaging and Mapping Agency ?Department of Defense World Geodetic System 1984,
its Definition and Relationships with Local Geodetic Systems,? Technical report, National Imag-
ing and Mapping Agency, 1997. Available at permanent link http://earth-info.nga.mil/GandG/
tr8350_2.html. This publication covers changes made to the standard since 1984 which are not
incorporated into this software.
[18] JACCHIA, L. G. ?New Static Models of the Thermosphere and Exosphere with Empirical Tem-
perature Models,? Technical Report 313, Smithsonian Astrophysical Observatory, 1970.
[19] MERSON, R. H. ?Numerical Integration of the Differential Equations of Celestial Mechanics,?
Technical Report TR 74184, Royal Aircraft Establishment, Farnborough, Hants, UK, January
1975, Defense Technical Information Center number AD B004645.
Implementation of Gauss-Jackson Integration for Orbit Propagation 357