James R. Scott

Glenn Research Center, Cleveland, Ohio

Taylor Series Trajectory Calculations Including

Oblateness Effects and Variable Atmospheric Density

NASA/TM—2011-217135

August 2011

AAS–11–612
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings that are preliminary or 

of specialized interest, e.g., quick release 

reports, working papers, and bibliographies that 

contain minimal annotation. Does not contain 

extensive analysis.

• CONTRACTOR REPORT. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to help@

sti.nasa.gov

• Fax your question to the NASA STI Help Desk 

at 443–757–5803

• Telephone the NASA STI Help Desk at

 443–757–5802

• Write to:

 NASA Center for AeroSpace Information (CASI)

 7115 Standard Drive

 Hanover, MD 21076–1320
James R. Scott

Glenn Research Center, Cleveland, Ohio

Taylor Series Trajectory Calculations Including

Oblateness Effects and Variable Atmospheric Density

NASA/TM—2011-217135

August 2011

AAS–11–612

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

2011 Astrodynamics Specialist Conference

cosponsored by the American Astronautical Society (AAS) and 

American Institute of Aeronautics and Astronautics (AIAA) 

Girdwood, Alaska, July 31 to August 4, 2011
Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

Available electronically at http://www.sti.nasa.gov

Trade names and trademarks are used in this report for identifi cation 

only. Their usage does not constitute an offi cial endorsement, 

either expressed or implied, by the National Aeronautics and 

Space Administration.

Level of Review: This material has been technically reviewed by technical management. 

This report contains preliminary fi ndings, 

subject to revision as analysis proceeds.
NASA/TM—2011-217135 1 

Taylor Series Trajectory Calculations Including 

Oblateness Effects and 

Variable Atmospheric Density 

James R. Scott 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

ABSTRACT 

Taylor series integration is implemented in NASA Glenn’s Spacecraft N-body Analysis 

Program, and compared head-to-head with the code’s existing 8th- order Runge-Kutta Fehlberg 

time integration scheme. This paper focuses on trajectory problems that include oblateness 

and/or variable atmospheric density. Taylor series is shown to be significantly faster and more 

accurate for oblateness problems up through a 4x4 field, with speedups ranging from a factor of 

2 to 13. For problems with variable atmospheric density, speedups average 24 for atmospheric 

density alone, and average 1.6 to 8.2 when density and oblateness are combined. 

INTRODUCTION 

We recently demonstrated the significant advantages of Taylor series integration in calculating 

spacecraft trajectories, i.e., order-of-magnitude speedup and simultaneous improvement in accuracy over a conventional 8th -order Runge-Kutta method.1

 Force models included central body 

gravitation (including Jଶ), other body gravitation, thrust, constant atmospheric drag and solar radiation pressure with constant illumination. 

In this paper we consider gravitational harmonics beyond Jଶ and drag due to variable atmospheric density. We show that Taylor series (TS) offers significant computational advantages in 

both cases, but also has limitations.

An alternative TS approach for handling oblateness and variable atmospheric density has been 

presented by Montenbruck.2

 The key differences here are the present implementation applies to 

spacecraft virtually anywhere in the solar system, can be used interchangeably with another integration method and employs the 1976 Standard Atmosphere3

. 

In the next section we summarize the TS methodology for handling oblateness effects and atmospheric density, and discuss its implementation in NASA Glenn’s Spacecraft N-body Analysis 

Program (SNAP)4

. Following that we present several example problems and discuss numerical 

results. Head-to-head comparisons are made with SNAP’s 8th-order Runge-Kutta Fehlberg 

(RKF) scheme.5
NASA/TM—2011-217135 2 

TAYLOR SERIES FORMULATION 

The governing equations of spacecraft motion are 

࢜  ൌ  ࢞

 (1)                                                                       ሻݐ ,଻ݔ ,଺ݔ ,ହݔ ,ସݔ ,ଷݔ ,ଶݔ ,ଵݔሺࢇ  ൌ   ࢜

ݔ଻  ൌ   െ݉ሶ ሺݐሻ

where ሺݔଵ, ݔଶ, ݔଷሻ ൌ  ࢞ is the spacecraft position in Cartesian coordinates relative to an inertial 

frame centered at the central body, ሺݔସ, ݔହ, ݔ଺ሻ ൌ ࢜ is the velocity, ݔ ଻is the spacecraft mass, 

m is the mass flow rate and a is the acceleration, where 

 (2)                                                 ୭ୠ୭ࢇ ൅ ୭ୠୡࢇ ൅ ୱ୰୮ࢇ ൅ ୢࢇ ൅ ୦୲ࢇ ൅ ୭ୠࢇ ൅ ୡୠࢇ ൌ ࢇ

a is the sum of accelerations from the central body point mass, other body point masses, thrust, 

atmospheric drag, solar radiation pressure, oblateness effects of the central body, and oblateness 

effects of other bodies. 

 As described previously1

, the basic approach is to expand each component of the spacecraft 

state vector ࢄ ൌ ሺݔଵ, ݔଶ, ݔଷ, ݔସ, ݔହ, ݔ଺, ݔ଻ሻ in a local Taylor series which is then used to advance the 

solution to the next time level. Derivatives are obtained by reformulating Equations (1) into a 

desired canonical form and then making use of highly efficient differentiation arithmetic. A series expansion is made at each time level i, 

k

i

K

k

n n x (t) X (k) (t t )

0

   

 , (3) 

for n = 1,2, ..., N, where N is the number of equations in the reformulated system, and the Taylor 

coefficients ܺ௡ሺ݇ሻ are evaluated through simple recursive formulas. 

Oblateness Effects of the Central Body 

The gravitational potential due to a non-spherical shaped central body is expressed in the form 

of a truncated series whose coefficients are obtained from experimental observation. Nonspherical gravity accelerations are readily obtained as the gradient of the potential function. 

While simple in concept, the mathematical expressions are quite involved and will not be presented here. See Baker6

 (pp. 144-150), Battin7 (pp. 401-408) or Bate et al8 (pp. 419-423). 

Table 1. Number of Auxiliary Variables per Oblateness Case

 n 

 m 

0 2 3 4 

 0 2 33 52 59 

 1 59 74 90 

 2 68 90 113 

 3 99 129 

 4 138 
NASA/TM—2011-217135 3 

Due to the complexity of the oblateness equations, the number of series terms employed in the 

gravitational potential and the need to transform from earth-fixed coordinates to the equator and 

equinox-of-date framework, the number of auxiliary variables required to handle oblateness effects is significantly higher than that of other acceleration terms. Table 1 summarizes the number 

of auxiliary variables required up through a 4 x 4 field. Each auxiliary variable corresponds to 

an additional equation in the differential system. The number of auxiliary variables required for 

the central body point mass is shown as the n=0, m=0 case, where (n,m) = (degree, order) of the 

spherical harmonics. Note that the central body point mass requires only two additional equations, whereas Jଶ alone requires 33 additional equations. The (4,4) case requires 138 additional 

equations. Clearly, the increasing number of equations will make the TS method inefficient as 

oblateness resolution is increased beyond a certain point. 

Variable Atmospheric Density 

 The drag acceleration term ࢇୢ is 

       4 2 2 7

2

1

2

6

2

5 2 1

2

,1 1 4 2 2 a c x c x x c x x x c x / x d       (4) 

       5 2 1 7

2

1

2

6

2

5 2 1

2

,2 1 4 2 2 a c x c x x c x x x c x / x d       (5) 

     6 7

2

1

2

6

2

5 2 1

2

,3 1 4 2 2 a c x c x x c x x x / x d      (6) 

where c   CD  A   2

1

1 , ܣ is the spacecraft area, ߩ is the atmospheric density, CD is the 

drag coefficient and 2 c = rotation rate of Earth. For spacecraft orbits where ߩ varies, it is necessary to calculate the time derivatives of ܿଵ. The first derivative is just 

݀ܿଵ

݀ݐ ൌ െ 1 

ߩ݀ ܣ ஽ܥ  2

݀ݐ                                                                         ሺ7ሻ

where 

ߩ݀

ߩ߲   ൌ ݐ݀ߩ߲  ൅ ݐ߲

ଵݔ߲

ଵݔ݀

ߩ߲  ൅ ݐ݀

ଶݔ߲

ଶݔ݀

ߩ߲  ൅  ݐ݀

ଷݔ߲

ଷݔ݀

݀ݐ                                                 ሺ8ሻ

or 

ߩ݀

ߩ߲   ൌ ݐ݀

߲ݐ ൅  સρ ∙ ࢜                                                                        ሺ9ሻ

The derivatives of ߩ must be obtained from the atmospheric density model. We show how this 

can be done using the 1976 Standard Atmosphere. 

 Let the spacecraft altitude h be approximated by 

 (10) ௘ݎ െ| ࢞| ൌ ݄

where ݎ ௘is the radius of the earth. Then 

ߩ߲

ଵݔ߲

ߩ߲   ൌ

߲݄

߲݄

ଵݔ߲

ߩ  ߲ ൌ  

߲݄

ଵݔ

|࢞|

                                                                   ሺ11ሻ
NASA/TM—2011-217135 4 

ߩ߲

ଶݔ߲

ߩ߲   ൌ

߲݄

߲݄

ଶݔ߲

ߩ߲    ൌ  

߲݄

ଶݔ

|࢞|

                                                                     ሺ12ሻ

ߩ߲

ଷݔ߲

ߩ߲   ൌ

߲݄

߲݄

ଷݔ߲

ߩ߲    ൌ  

߲݄

ଷݔ

|࢞|

                                                                     ሺ13ሻ

so that 

ߩ߲   ൌ ߩસ

߲݄

1

|࢞|

  ࢞                                                                            ሺ14ሻ

Since డఘ

డ௧ ൌ 0 for the 1976 Standard Atmosphere, Equation (9) becomes 

ߩ݀

ߩ߲ ൌ ݐ݀

߲݄

1

|࢞|

ሻ15ሺ                                                                         ࢜ ∙ ࢞  

Equation (15) can be evaluated so long as డఘ

డ௛ can be determined for arbitrary h. Here డఘడ௛ can be 

approximated using finite differences. For h bracketed by ݄௜ and ݄௜ାଵ in the Standard Atmosphere discretization, ఘ೔శభି ఘ೔

௛೔శభି ௛೔

 approximates డఘ

డ௛ to at least first order, so ௗఘௗ௧

 is also known to at 

least first order. 

Equation (15) shows that ௗఘ

ௗ௧ is inversely proportional to the magnitude of the position vector. 

In addition, for a purely circular orbit, ࢞∙࢜≡0, so that ௗఘ

ௗ௧

 is zero. This is consistent with and 

follows from the fact that density is a function of altitude only in the 1976 Standard Atmosphere. 

Of course, there are no perfectly circular spacecraft orbits, but in general ࢞∙࢜ will be quite 

small for orbits with small eccentricity. The result is that ௗఘ

ௗ௧ will be small and higher time derivatives of density may be negligible in many instances. 

APPLICATION TO NEAR EARTH TRAJECTORIES 

We consider the following problems. (1) Earth-orbiting satellite with oblateness effects and 

zero atmospheric drag; (2) Spacecraft spiraling out of earth’s gravity well in a low thrust trajectory with oblateness effects and zero atmospheric drag; (3) Earth-orbiting satellite with atmospheric 

drag and zero oblateness; (4) Earth-orbiting satellite with variable atmospheric drag and oblateness effects; (5) Earth-orbiting satellite with constant atmospheric drag and oblateness effects. 

For all problems the spacecraft mass is 10,000 kg and the moon and sun are included as other bodies. The drag coefficient is 0.7 and the spacecraft area is 240 square meters. In Problems (1) 

and (3) - (5), the earth-orbiting satellite is at an inclination of 28.45 degrees with initial altitude 

415 km and the trajectory is propagated for 10 days. For Problem 2 the calculation ends when the 

semimajor axis equals 40,000 km. 

All calculations were run on a Dell Poweredge 6950, which has 8 AMD 2.8 GHz processers 

and 16 GB of RAM. TS calculations used a series with 20 terms for all problems. Both RKF 

and TS used a local error tolerance defined by ߬௜ ൌ   Min௡ ሾ |ݔ |௡߬ + τ] where ݔ ௡is the ݊௧௛ variable in the state vector/differential system, and τ denotes the error tolerance parameter. 
NASA/TM—2011-217135 5 

Problems with Oblateness 

We consider Problems 1 and 2, with spherical harmonics having degree and order (n,m) = 

(2,0), (2,2), (3,3) and (4,4). The tables below show (a) the head-to-head ratio of RKF/TS CPU 

time, (b) relative ܮଶ convergence of each method to its τ = 1.e-14 solution and (c) difference between TS and RKF solutions at τ = 1.e-14. Relative ܮଶ convergence is defined by ห࢞ െ  ࢌ࢞ห ̸หࢌ࢞ห

where ࢌ࢞ is the final spacecraft position at ߬ ൌ 1.e-14 and x is the final spacecraft position at 

larger ߬, as listed in the tables. Relative ܮଶ difference between RKF and TS is defined by 

 .|୊୏ୖ࢞| / |୊୏ୖ࢞ െ ୗ୘࢞|

Table 2.a CPU ratios, RKF/TS, for Problem 1 

 (n,m) (2,0) (2,2) (3,3) (4,4) 

 τ

1.e-10 8.61 3.07 2.29 1.61 

1.e-11 9.87 3.58 2.57 1.89 

1.e-12 12.94 4.12 3.10 2.22 

1.e-13 15.14 4.88 3.59 2.68 

1.e-14 17.41 5.83 4.31 3.17 

 Average: 12.80 4.29 3.17 2.31 

Table 2.b Relative ࡸ ૛convergence levels for Problem 1 

 RKF TS 

 (n,m) (2,0) (2,2) (3,3) (4,4) (2,0) (2,2) (3,3) (4,4) 

 τ

1.e-10 0.42e-06 0.42e-06 0.42e-06 0.42e-06 0.13e-07 0.36e-10 0.13e-09 0.24e-10 

1.e-11 0.34e-07 0.34e-07 0.34e-07 0.34e-07 0.56e-08 0.25e-10 0.71e-10 0.54e-10 

1.e-12 0.27e-08 0.27e-08 0.26e-08 0.27e-08 0.31e-09 0.79e-10 0.10e-09 0.94e-10 

1.e-13 0.22e-09 0.21e-09 0.18e-09 0.20e-09 0.35e-09 0.14e-09 0.47e-10 0.85e-10 

Table 2.c Difference between RKF and TS solution for Problem 1 

 (n,m) Relative ܮଶ Difference Absolute Difference (km) 

(2,0) 0.148e-09 0.100e-05 

(2,2) 0.322e-10 0.218e-06 

(3,3) 0.819e-10 0.554e-06 

(4,4) 0.146e-09 0.990e-06 

Table 3.a CPU ratios, RKF/TS, for Problem 2 

 (n,m) (2,0) (2,2) (3,3) (4,4) 

 τ

1.e-10 9.14 3.23 2.36 1.74 

1.e-11 10.60 3.74 2.76 2.03 

1.e-12 12.55 4.45 3.27 2.39 

1.e-13 15.32 5.23 3.86 2.83 

1.e-14 17.72 6.16 4.53 3.34 

 Average: 13.07 4.56 3.35 2.47 
NASA/TM—2011-217135 6 

Table 3.b Relative ࡸ ૛convergence levels for Problem 2 

 RKF TS 

 (n,m) (2,0) (2,2) (3,3) (4,4) (2,0) (2,2) (3,3) (4,4) 

 τ

1.e-10 0.13e-04 0.13e-04 0.13e-04 0.13e-04 0.32e-05 0.21e-07 0.15e-08 0.10e-08

1.e-11 0.11e-05 0.11e-05 0.11e-05 0.11e-05 0.87e-07 0.18e-08 0.32e-09 0.86e-09

1.e-12 0.83e-07 0.83e-07 0.83e-07 0.83e-07 0.30e-07 0.77e-09 0.25e-09 0.12e-08

1.e-13 0.59e-08 0.58e-08 0.63e-08 0.61e-08 0.58e-08 0.13e-09 0.85e-09 0.16e-09

Table 3.c Difference between RKF and TS solution for Problem 2 

 (n,m) Relative ܮଶ Difference Absolute Difference (km) 

 (2,0) 0.526e-07 0.210e-02 

 (2,2) 0.662e-07 0.265e-02 

 (3,3) 0.753e-07 0.301e-02 

 (4,4) 0.769e-07 0.308e-02 

Average TS speedups range from 2.31 to 12.8 for Problem 1, and from 2.47 to 13.07 for Problem 2. Tables 2.b and 3.b show that TS is more converged than RKF in all but one case. CPU 

times versus convergence levels are plotted in Figures 1 and 2. 

The small differences in RKF and TS solutions shown in Tables 2.c and 3.c are due to the fact 

that TS calculates the motion of the sun and moon, whereas RKF obtains their motion from 

ephemeris files. The differences are bigger in Problem 2 because the trajectory is propagated 

over 112 days versus 10 days for Problem 1. 

Figure 1.a CPU time vs. valid digits for Problem 1, (n,m) = (2,0), (2,2) 

0

2

4

6

8

10

12

14

5 6 7 8 9 10 11 12

CPU Time (sec)

- Log (Convergence Level)

RKF Prob. 1 - (2,0)

TS Prob. 1 - (2,0)

RKF Prob. 1 - (2,2)

TS Prob. 1 - (2,2)
NASA/TM—2011-217135 7 

Figure 1.b CPU time vs. valid digits for Problem 1, (n,m) = (3,3), (4,4) 

Figure 2.a CPU time vs. valid digits for Problem 2, (n,m) = (2,0), (2,2) 

0

2

4

6

8

10

12

14

5 6 7 8 9 10 11 12

CPU Time (sec)

- Log (Convergence Level)

RKF Prob. 1 - (3,3)

TS Prob. 1 - (3,3)

RKF Prob. 1 - (4,4)

TS Prob. 1 - (4,4)

0

10

20

30

40

50

60

70

4 5 6 7 8 9 10 11

CPU Time (sec)

- Log (Convergence Level)

RKF Prob. 2 - (2,0)

TS Prob. 2 - (2,0)

RKF Prob. 2 - (2,2)

TS Prob. 2 - (2,2)
NASA/TM—2011-217135 8 

Figure 2.b CPU time vs. valid digits for Problem 2, (n,m) = (3,3), (4,4) 

Problems with Atmospheric Drag 

Table 4. Effect of Time Derivative Density Term on Results 

 Difference in RKF and TS spacecraft position (km) 

 Without Variable With Variable 

 Density Model Density Model 

 τ

 1.e-10 0.0306 0.00058 

 1.e-11 0.0292 0.00189 

 1.e-12 0.0256 0.00153 

 1.e-13 0.0230 0.00148 

 1.e-14 0.0205 0.00124 

 We present results for an earth-orbiting satellite with atmospheric drag but no oblateness 

(Problem 3). We first consider the effect of the time derivative density term ௗ௖భ

ௗ௧

. Table 4 shows 

the difference in RKF and TS spacecraft positions after the satellite trajectory has been propagated for 10 days, where the TS trajectory is propagated with and without ௗ௖భ

ௗ௧ . The results clearly 

show the improvement in accuracy due to inclusion of the time derivative density term. 

Table 5 summarizes RKF and TS performance on Problem 3. Constant density results are included for comparison. The speedups are quite large here due to the fact that TS can handle the 

drag equations with only four auxiliary variables. Figure 3 plots the CPU times versus convergence levels. 

0

10

20

30

40

50

60

70

4 5 6 7 8 9 10 11

CPU Time (sec)

- Log (Convergence Level)

RKF Prob. 2 - (3,3)

TS Prob. 2 - (3,3)

RKF Prob. 2 - (4,4)

TS Prob. 2 - (4,4)
NASA/TM—2011-217135 9 

Table 5.a CPU ratios, RKF/TS, for Problem 3 

 Constant Density Variable Density 

 τ

 1.e-10 23.38 17.39 

 1.e-11 28.82 19.85 

 1.e-12 34.84 23.30 

 1.e-13 40.77 28.22 

 1.e-14 47.90 33.41 

 Average: 35.14 24.43 

Table 5.b Relative ࡸ ૛convergence levels for Problem 3 

 RKF TS 

 Constant Density Variable Density Constant Density Variable Density 

 τ

1.e-10 0.43e-06 0.44e-06 0.68e-08 0.17e-06 

1.e-11 0.35e-07 0.35e-07 0.10e-09 0.13e-06 

1.e-12 0.27e-08 0.28e-08 0.37e-09 0.46e-07 

1.e-13 0.20e-09 0.20e-09 0.57e-10 0.36e-07 

 

Table 5.c Difference between RKF and TS solution for Problem 3 

 Constant Density Variable Density 

Relative ܮଶ Difference Absolute Difference (km) Relative ܮଶ Difference Absolute Difference (km)

 0.117e-10 0.794e-07 0.182e-06 0.124e-02 

Figure 3. CPU time vs. valid digits for Problem 3 

0

2

4

6

8

10

12

14

4 5 6 7 8 9 10 11

CPU Time (sec)

- Log (Convergence Level)

RKF Prob. 3 - Variable Density

TS Prob. 3 - Variable Density

RKF Prob. 3 - Constant Density

TS Prob. 3 - Constant Density
NASA/TM—2011-217135 10 

Problems with Atmospheric Drag and Oblateness

We consider Problems 4 and 5. RKF and TS performance are summarized below and results 

are plotted in the figures which follow. Average TS speedups range from 1.6 to 8.2 for Problem 

4 and from 2.2 to 11.9 for Problem 5. The largest absolute difference between RKF and TS 

shown in Table 6.c is 81.8 meters. It is also worth noting that RKF and TS both predict a reduction in the orbit semimajor axis of 2.709 km (for the (4,4) case). This compares with the dragonly case (Problem 3) in which both methods predict a reduction in semimajor axis of 1.906 km. 

Atmospheric drag combined with earth’s oblateness thus leads to a significantly greater loss in 

spacecraft altitude. 

Table 6.a CPU ratios, RKF/TS, for Problem 4

 (n,m) (2,0) (2,2) (3,3) (4,4) 

 τ

 1.e-10 5.89 2.12 1.55 1.13 

 1.e-11 6.84 2.47 1.81 1.32 

 1.e-12 7.96 2.88 2.09 1.53 

 1.e-13 9.19 3.41 2.45 1.84 

 1.e-14 11.03 3.98 2.88 2.17 

 Average: 8.18 2.97 2.16 1.60 

 

Table 6.b Relative ࡸ ૛convergence levels for Problem 4

 RKF TS 

 (n,m) (2,0) (2,2) (3,3) (4,4) (2,0) (2,2) (3,3) (4,4) 

 τ

1.e-10 0.42e-06 0.43e-06 0.43e-06 0.43e-06 0.40e-05 0.20e-04 0.14e-04 0.12e-04 

1.e-11 0.35e-07 0.36e-07 0.34e-07 0.37e-07 0.12e-05 0.13e-04 0.86e-05 0.73e-05 

1.e-12 0.10e-08 0.36e-08 0.28e-08 0.48e-08 0.94e-05 0.92e-05 0.44e-05 0.44e-05 

1.e-13 0.92e-09 0.22e-08 0.73e-09 0.48e-09 0.77e-05 0.36e-05 0.20e-05 0.19e-05 

Table 6.c Difference between RKF and TS solution for Problem 4

 (n,m) Relative ܮଶ Difference Absolute Difference (km) 

 (2,0) 0.881e-06 0.597e-02 

 (2,2) 0.121e-04 0.818e-01 

 (3,3) 0.856e-05 0.579e-01 

 (4,4) 0.732e-05 0.495e-01 

Table 7.a CPU ratios, RKF/TS, for Problem 5 

 (n,m) (2,0) (2,2) (3,3) (4,4) 

 τ

 1.e-10 8.14 2.98 2.12 1.62 

 1.e-11 9.44 3.43 2.52 1.84 

 1.e-12 11.87 4.03 2.92 2.18 

 1.e-13 13.70 4.72 3.44 2.53 

 1.e-14 16.42 5.71 4.07 3.05 

 Average: 11.91 4.17 3.02 2.24 
NASA/TM—2011-217135 11 

Table 7.b Relative ࡸ ૛convergence levels for Problem 5

 RKF TS 

 (n,m) (2,0) (2,2) (3,3) (4,4) (2,0) (2,2) (3,3) (4,4) 

 τ

1.e-10 0.42e-06 0.42e-06 0.42e-06 0.42e-06 0.29e-07 0.13e-09 0.69e-10 0.88e-10 

1.e-11 0.34e-07 0.34e-07 0.34e-07 0.34e-07 0.39e-08 0.45e-10 0.67e-10 0.39e-10 

1.e-12 0.27e-08 0.27e-08 0.27e-08 0.26e-08 0.34e-09 0.48e-10 0.39e-10 0.75e-10 

1.e-13 0.20e-09 0.15e-09 0.21e-09 0.21e-09 0.51e-10 0.14e-09 0.71e-10 0.45e-10 

Table 7.c Difference between RKF and TS solution for Problem 5 

 (n,m) Relative ܮଶ Difference Absolute Difference (km) 

 (2,0) 0.843e-10 0.570e-06 

 (2,2) 0.282e-10 0.191e-06 

 (3,3) 0.700e-10 0.474e-06 

 (4,4) 0.563e-10 0.381e-06

 

Figure 4.a CPU time vs. valid digits for Problem 4, (n,m) = (2,0), (2,2) 

0

2

4

6

8

10

12

14

3 4 5 6 7 8 9 10

CPU Time (sec)

- Log (Convergence Level)

RKF Prob. 4 - (2,0)

TS Prob. 4 - (2,0)

RKF Prob. 4 - (2,2)

TS Prob. 4 - (2,2)
NASA/TM—2011-217135 12 

Figure 4.b CPU time vs. valid digits for Problem 4, (n,m) = (3,3), (4,4) 

Figure 5.a CPU time vs. valid digits for Problem 5, (n,m) = (2,0), (2,2) 

0

2

4

6

8

10

12

14

3 4 5 6 7 8 9 10

CPU Time (sec)

- Log (Convergence Level)

RKF Prob. 4 - (3,3)

TS Prob. 4 - (3,3)

RKF Prob. 4 - (4,4)

TS Prob. 4 - (4,4)

0

2

4

6

8

10

12

14

5 6 7 8 9 10 11

CPU Time (sec)

- Log (Convergence Level)

RKF Prob. 5 - (2,0)

TS Prob. 5 - (2,0)

RKF Prob. 5 - (2,2)

TS Prob. 5 - (2,2)
NASA/TM—2011-217135 13 

Figure 5.b CPU time vs. valid digits for Problem 5, (n,m) = (3,3), (4,4) 

CONCLUSION 

Taylor series integration is implemented in a high fidelity trajectory propagation code and applied to near earth trajectory problems. Head-to-head comparisons are made with 8th-order 

Runge-Kutta Fehlberg integration. Results show that TS is faster by a factor of 2-13 for oblateness problems up through a 4x4 field, while simultaneously improving accuracy. For problems 

including both oblateness and variable atmospheric density, TS is faster by a factor of 1.6 to 8.2. 

We conclude that TS is a highly competitive integration method for realistic trajectory problems 

involving oblateness and atmospheric density. 

REFERENCES 

1 Scott, J.R. and Martini, M.C., “High Speed Solution of Spacecraft Trajectory Problems Using Taylor Series Integration,” Journal of Spacecraft and Rockets, Vol. 47, No. 1, Jan.-Feb. 2010, pp. 199-202; also presented as paper no. 

AIAA 2008-6957 at the AIAA/AAS Astrodynamics Specialist Conference and Exhibit, 18-21 Aug. 2008, Honolulu, 

Hawaii.

2

Montenbruck, O., “Numerical Integration of Orbital Motion Using Taylor Series,” American Astronautical Society/AIAA Spaceflight Mechanics Meeting, American Astronautical Society Paper 92-195, Feb. 1992.

3 “1976 U.S. Standard Atmosphere”, NOAA, NASA, USAF, U.S. Government Print Office, Washington, D.C., October 1976. 

4

Martini, M.C., “S.N.A.P. 2.3 User’s Manual, Spacecraft N-Body Analysis Program,” Analex Corporation, Glenn Research Center, Cleveland, OH, 2005. 

0

2

4

6

8

10

12

14

5 6 7 8 9 10 11 12

CPU Time (sec)

- Log (Convergence Level)

RKF Prob. 5 - (3,3)

TS Prob. 5 - (3,3)

RKF Prob. 5 - (4,4)

TS Prob. 5 - (4,4)
NASA/TM—2011-217135 14 

5

Fehlberg, E., “Classical Fifth-, Sixth-, Seventh-, and Eighth-Order Runge-Kutta Formulas with Stepsize Control,” 

NASA TR R-287, Oct. 1968. 

6

Baker, R., Astrodynamics: Applications and Advanced Topics, Academic Press, New York and London, 1976. 

7

Battin, R., An Introduction to the Mathematics and Methods of Astrodynamics, Revised Edition, AIAA Education 

Series, J. Przemieniecki Editor-in-Chief, American Institute of Aeronautics and Astronautics, Inc., Reston, VA, 1999. 

8

Bate, R., Mueller, D., and White, J., Fundamentals of Astrodynamics, Dover Publications, Inc., New York, 1971. 
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 

control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

01-08-2011 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Taylor Series Trajectory Calculations Including Oblateness Effects and Variable Atmospheric 

Density 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Scott, James, R. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 136905.02.09.03.04 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-17871 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2011-217135 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Categories: 13 and 64 

Available electronically at http://www.sti.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

Taylor series integration is implemented in NASA Glenn’s Spacecraft N-body Analysis Program, and compared head-to-head with the

code’s existing 8th- order Runge-Kutta Fehlberg time integration scheme. This paper focuses on trajectory problems that include oblateness 

and/or variable atmospheric density. Taylor series is shown to be significantly faster and more accurate for oblateness problems up through 

a 4 by 4 field, with speedups ranging from a factor of 2 to 13. For problems with variable atmospheric density, speedups average 24 for 

atmospheric density alone, and average 1.6 to 8.2 when density and oblateness are combined.

15. SUBJECT TERMS

Taylor series; Spacecraft trajectories; Runge-Kutta method; Trajectory analysis; Spacecraft N-body Analysis Program (SNAP); 

Numberical integration; High order 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

20 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18