UFDC Home  |  Search all Groups  |  UF Institutional Repository  |  UF Institutional Repository  |  UF Theses & Dissertations |   Help

# A theory of the earth's precession relative to the invariable plane of the solar system

## Material Information

Title:
A theory of the earth's precession relative to the invariable plane of the solar system
Physical Description:
xi, 263 leaves : ill. ; 28 cm.
Language:
English
Creator:
Owen, William Mann
Publication Date:

## Subjects

Subjects / Keywords:
Precession   ( lcsh )
Astronomy thesis Ph. D
Dissertations, Academic -- Astronomy -- UF
Genre:
bibliography   ( marcgt )
non-fiction   ( marcgt )

## Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1990.
Bibliography:
Includes bibliographical references (leaves 259-262).
Statement of Responsibility:
by William Mann Owen.
General Note:
Typescript.
General Note:
Vita.

## Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 001659172
oclc - 24529617
notis - AHX0899
System ID:
AA00002116:00001

Full Text

A THEORY OF THE EARTH'S PRECESSION RELATIVE TO THE
INVARIABLE PLANE OF THE SOLAR SYSTEM

By

WILLIAM MANN OWEN, JR.

A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FTLORTDA IN PARTIAL FTTIFTT.T.MFNT

by

William Mann Owen, Jr.

DEDICATION

In the wee hours of March 13, 1960, a five-year-old boy was awakened by his grand-

father in

order to see a total eclipse of the moon.

Thirty years have passed;

has become a man, and his interest in astronomy, planted that night, has taken root and

blossomed. It is therefore to the memory of

R. Waldo Hambleton (1901-1981)

that this dissertation is lovingly dedicated.

ACKNOWLEDGMENTS

No enterprise of this nature can be accomplished without the support of many others.

My chairman, Dr. Heinrich K. Eichhorn, and my colleague at the Jet Propulsion Labo-

ratory, Dr. Jay

H. Lieske, reviewed early versions of the manuscript;

the quality of this

dissertation has been helped immeasurably by their critiques. The author benefited from

many conversations with them and with Drs. E. Myles Standish and X X Newhall of JPL

and with Dr. Jacques Laskar of the Bureau des Longitudes. Laskar's cooperation in pro-

viding his tabulation of the Earth's orbital parameters in advance of publication is greatly

appreciated.

The unwavering encouragement of my supervisors at JPL, Drs. Stephen P.

Synnott and James P. McDanell, was most helpful as well. Finally, the author is indebted

to Eva Eichhorn for translating the 1967 paper by Sharaf and Budnikova.

Financial support was provided at various times by the tuition assistance program of

the Jet Propulsion Laboratory, by a Graduate Council Fellowship from the University of

Florida, and by a Graduate Fellowship from the National Science Foundation.

Computer

support was provided by

The author is grateful to all these

institutions.

The NSF

requires in addition the following acknowledgment and disclaimer:

"This material is based on work supported under a National Science Foundation Gradu-

ate Fellowship.

Any opinions, findings, conclusions or recommendations expressed in this

nlhl'ir tinn taro thi~h ', n,-f th0 miih-,,hr c+n,1 Ar ncvf rnnr'ne..-n-.;l, .rfl,.-+ +. .;...,r "--.f +tr" MT"+n'n1' "

The TEX program (Knuth 1984) typeset the dissertation.

library (Pearson 1989) was used to create the figures. The hig

of this dissertation is due to these two programs; nevertheless,

The PGPLOT subroutine

h quality of the appearance

any defects in layout or

errors in substance are solely the responsibility of the author.

page

ACKNOWLEDGMENTS

LIST OF TABLES

S 5 0 5 5 5 5 0 5 5 5 5 S S S S S S S S S S S S v iii

LIST OF FIGURES

ABSTRACT

CHAPTERS

INTRODUCTION

S S S S S S S S *1

Background
Notation
Definitions

THE DETERMINATION OF THE INVARIABLE PLANE

Introduction

JPL Planetary Ephemerides .
The M04786 Planetary Ephemeris

. S . S
S S S S f

The Total Orbital Angular Momentum of the

olar System

The Uncertainty in the Total Orbital Angular Momentum
The Rotational Angular Momentum of the Solar System
The Adopted Orientation of the Invariable Plane .

THE SHORT-TERM THEORY

Introduction

Analytic Formulas for I, L,

and A

3.3.1
3.3.2
000\<-

Series Expansions for I, L, and A
The Expansion for I .
The Expansion for L ..
mi T ~l I' A

SS S S S. V
n- flV

THE LONG-TERM THEORY

. S. . 84

Overview

The Motion of the Ecliptic

4.3.1
4.3.2
4.3.3

The Equations of Motion for the Celestial Pole
The Vector Equation of Motion . .
The Equation of Motion in Component Form
Kinoshita's Expression for Luni-solar Precession

. .
.

The Numerical Integration of the Motion of the Celestial Pole
The Determination of the Precession Angles . .
The Chebyshev Representation of the Precession Angles
Sources of Error in the Long-Term Theory . .

CONCLUSIONS

The Short-Term and Long-Term Theories Compared . .
Comparing the Classical and Invariable Plane Precession Formulations
Summer .

REFERENCES

BIOGRAPHICAL SKETCH

LIST OF TABLES

page

Definitions of Symbols

. . . 5 9

Physical Constants from the M04786 Ephemeris

Vectors from the M04786 Ephemeris at 1989 August 25 04:00 ET

Asteroid Vectors at 1989 August 2

Relativistic Masses from the M047

5 04:00 ET

6 Ephemeris at 1989 August

25 04:00

Planetary Masses and Uncertainties

Covariances of the Planets'

Full-Precision

Orbital Angular Momenta

Values of Current Precessional Constants

* S S 5 0 5

* S S S S S S

Constants in the Short-Term Theory

Comparison of Rigorous Angles and Polynomial Approximations

Astronomical Constants for Kinoshita's Luni-solar Precession Model

Chebyshev Polynomial Coefficients for the Long-Term Theory .

Polynomial Coefficients from the Long-Term

Theory Near T

Chebyshev Coefficients for (A

9A, and

ZA from the Long-Term

Theory

NearT

S S S S S S S S S S S S S S S S S S S S S S S S

Polynomial Coefficients for CA

NearT

8A, and zA from the Long-Term Theory

* S S S S S S S S S S S S S S S 0 4 9 4 2 4 8

Polynomial Coefficients from the Long-Term and Short-Term Theories

. 245

LIST OF FIGURES

page

The Equatorial Coordinate System

. 4 7

The Ecliptic Coordinate System

The Effect of the Geocentric Lunar Orbit on the Invariable Plane

The Direction of the Total Angular Momentum (M04786 Ephemeris)

The Direction of the Total Angular Momentum (DE130 Ephemeris)

The Orientation of the Invariable Plane . . .

The Classical Precession Angles (, 0, and z

Precession Angles Using the Invariable Plane

The Equators and the Invariable Plane

. 6 6 . 6 S .

B S 6 0 .

Precession Between Two Arbitrary

Times

S S S 0 B J78

Spherical Coordinates for Q in the Eo System

S. S B 92

The Precession Angles for One Million

Years

S B 5 5 4 5 5 105

The Precession Angles CA, 0A, and ZA Near T

= +250

. S S 5 113

J2000 Equatorial Coordinates of the Ecliptic Pole and of the Pole of the
Invariable Plane . . .

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
A THEORY OF THE EARTH'S PRECESSION RELATIVE TO THE
INVARIABLE PLANE OF THE SOLAR SYSTEM

By

William Mann Owen, Jr.

December 1990

Chairman: Heinrich K

Eichhorn

Major Department: Astronomy

The most commonly used formulas for the rigorous application of general precession

employ three successive small rotations by which the equatorial coordinate system at epoch

is aligned with the equatorial system of date.

This dissertation presents an alternate method

for applying precession in which the rotation angles are intimately tied to the invariable

plane of the Solar System.

The work involved in constructing this theory divides neatly into three parts.

First,

the newest planetary ephemeris from the Jet Propulsion Laboratory, produced after the

Voyager

encounter with Neptune, allows one to infer the orientation of the invariable plane

with a standard error on the order of 0O'!04.

Second, the coefficients in the approximation

polynomials for the new angles are found in terms of their counterparts in the currently-

accepted IAU theory; this "short-term theory," like the IAU one, is valid for a few centuries

near the present time.

Third, the motion of the Earth's north pole is integrated numerically

over a million-year time span; values for both the standard precession angles and the new

ones are inferred at discrete times from the integration, and Chebyshev polynomials are fit

A comparison of the long-term and short-term theories reveals two possible improve-

ments to the IAU system of constants.

the rigid Earth, used in

The higher-order terms in Kinoshita's model for

the long-term theory, change the time derivative of Newcomb's

Precessional Constant from -0'00369 per Julian century to -0'00393/century at the stan-

dard epoch J2000.0. Laskar's investigation into the motion of the ecliptic, also used in the

long-term theory, yields -46'!8065/century for the rate of change of the obliquity versus

CHAPTER
INTRODUCTION

Background

The slow motion of the celestial pole (however one defines it) among the stars and

the concomitant slow westward drift of the location of the vernal equinox are the result of

the physical process known as astronomical precession.

The former effect results from the

gravitational interaction of the Sun and Moon with the oblate Earth; this produces a torque

which changes the direction of the Earth's rotational angular momentum vector. The short-

period portion of this phenomenon is known as "nutation" and will not be treated in depth

here;

the long-period

part is "luni-solar precession."

(As the longest nutation period is

18.6 years while the shortest period of luni-solar precession is about 25,000 years,

there is a clean division

between

the two manifestations of the same physical process.)

the gravity of the other planets in

the Solar System

perturbs the

Earth's

orbit, changing the orientation of its mean orbital plane (the ecliptic).

The influence of

this orientation change on the direction of the mean

vernal equinox is called

"planetary

precession."

Since the vernal equinox lies at the intersection of the Earth's mean equator

and the ecliptic, its location is affected by both motions.

The resultant of luni-solar and

planetary precession is accordingly called "general precession.

Luni-solar precession was discovered in the second century

B.C. by Hipparchus, who

2

motion of the ecliptic was not suspected until Newton applied his theory of gravitation to

the motions of the planets.

The theory of general precession is both dynamic and kinematic in nature.

The dy-

namics enters into the differential equations of motion for the ecliptic (or equivalently for

the pole of the ecliptic) and for the equator (or the celestial pole).

Once these equations are

integrated, the motion of the vernal equinox and of the equatorial and ecliptic coordinate

axes

becomes a problem in rotational kinematics.

Precession theory accordingly is devel-

oped in terms of angles by which one can rotate coordinate frames from their orientation

at one epoch to that at another.

The current theory of precession was developed by Lieske et al. (1977). Its time origin

is the beginning of the astronomical Julian

January 1 (Julian Ephemeris Date 2451545.0).

year 2000, or noon dynamical time on 2000

This zero epoch is denoted "J2000.0,

" and

time T in this theory is measured in units of Julian centuries of 36525 days of dynamical

time.

(The current distinction between

"terrestrial dynamical time"

and "barycentric dy-

namical time," a periodic difference of no more than two milliseconds, is of no consequence

in a theory whose rates do not exceed two degrees per century.) The paper of Lieske et al.

is based on Fricke's

comb's

(1971) determination of the speed of luni-solar precession and on New-

(1894) work on the motion of the ecliptic, the latter updated with modern values

of the masses of the planets.

While the rotation matrix approach to precession is fully rigorous, giving exact re-

sults for the effects of precession on equatorial coordinates,

the precession angles themselves are only approximate: Lies]

the formulas for evaluating

ke et at provide third-degree

3

reduction of modern precise astrometric observations. However, a blind application of the

polynomials over millenia will produce grotesquely wrong results.

Several authors have examined

precession

theory over much longer time intervals.

Here one generally expresses the orientation of the ecliptic as a sum of trigonometric terms

derived from long-period planetary perturbations; expressions for the accumulated general

precession in longitude follow.

This method was developed by Brouwer and van

Woerkom

(1950), corrected by Sharaf and Budnikova (1967), and extended to five million years by

Berger (1976). Berger's work was based on Bretagnon's

(1974) theory, which was complete

to second order in the planetary masses and to third degree in eccentricity and inclination.

The goal of this dissertation is to develop a theory of precession in which the invariable

plane of the Solar System is used as a fundamental reference plane.

Precession

angles

referred to the invariable plane can be used in both short-term and long-term applications.

With a suitable change in the origin of right ascensions, the new precession theory becomes

simpler computationally than the current theory.

The work described here divides neatly into three parts. In Chapter

of the invariable plane is found from

the orientation

the most recent planetary ephemeris produced

the Jet

Propulsion

Laboratory.

Chapter 3 presents formulas by which

the polynomial

coefficients of the new precession angles may be calculated from those of the current theory;

the results when the values in

the Lieske et al. (1977) paper are adopted may be called

"short-term theory."

Chapter 4 is devoted to the "long-term theory,"

developed by

numerical integration of the position of the celestial pole; the equations of motion are those

of Kinoshita (1977), while the "Numerical General Theory"

Notation

Table 1-1, found at the end of this chapter, presents a list of all the symbols used

throughout this dissertation with their definitions and the number of the section in which

the symbol first appears.

The table contains Latin symbols first, then Greek, then special

astronomical symbols, with the first two parts in alphabetical order.

Scalar quantities (including components of vectors and elements of matrices) are set in

italic type throughout this work.

Vectors are given in boldface roman type; individual rows

or columns of matrices are considered to be vectors for this purpose.

Matrices appear in

a boldface sans-serif font; the same font is used for both rotation matrices and covariance

matrices.

This scheme follows standard typographical practice for scalars and vectors; the

typography for matrices is that used by Goldstein (1980).

This paper follows commonly accepted notation for the precession angles.

The system

of notation for the coefficients of these angles pioneered by Lieske et al. (1977) is complete

and unambiguous;

at the same time the diacritical marks can be awkward.

Deviations

from their system are carefully noted both in Table 1-1 and at their first use in the text.

Equations are numbered by chapter in order of their first appearance.

When an equa-

tion appears more than once, all occurrences are labeled with the original equation number.

Definitions

Insofar

as possible, all the terminology used in this dissertation is currently in use by

the astronomical community.

An excellent glossary appears in section M of the annual

volumes of The Astronomical Almanac, published jointly by the United States Naval Ob-

servatory and the Royal Greenwich Observatory.

The reader is referred there for terms not

5

dissertation are right-handed: if the three axes be denoted (in order) by x, y, and z, and

if the x- and y-axes are placed on a piece of paper such that x increases to the right and y

increases upward, then the z-axis points out of the paper toward the observer.

The rectangular components of a vector are taken to form a column vector (a matrix

with only one column).

Vectors are also specified in terms of their spherical coordinates:

length (or magnitude) r, longitude angle A, and latitude angle ,3.

The transformation from

,A,/ } into rectangular components (x, y, z)T is

r cos A cos #
= rsinAcos fl
rsinf /

(1-1)

this is a vector equation.

The inverse transformation consists of the three scalar equations

r=V/

+y2+

(1-2)

A = plg(y,

(1-3)

f = tan

(1-4)

The notation "plg(y, xa)" in equation (1-3) was introduced by Eichhorn (1987/88)

as an al-

ternate to the more usual tan

-1(y/x).

The range of the plg function is 0

< plg(y, x)

< 360;

the quadrant depends on the signs of x and y individually.

language, A would be evaluated as ATAN2 (y, x). No such

In the Fortran programming

confusion exists for equations

(1-2) and (1-4); r is nonnegative, and /3 has the same range as the principal value of the

arctangent. Equation (1-4) is given in terms of the arctangent rather than the more usual

3 = sin

-(z/r) for numerical reasons:

the arcsine loses its precision

when its argument

approaches +1.

Finally, if r = 0, the equations for A and 0 are technically indeterminate:

a:
y
z

+y2).

the matrix is applied.

There are three elementary rotation matrices, denoted Ri(9),

producing a rotation by some angle 0 about one of the three coordinate

sign convention applies to the angles; for instance, a rotation about the

angle carries the x-axis toward the old y-axis.

axes. The standard

z-axis by a positive

Written out, the three elementary rotation

matrices are

0
COS0

Ri(9)

- sin 0

0
sin 0
cos 0

R2(0)

Ra(0)

cos9
0
sinG0

cos 0
-sin 9
0

- sin 0
0
cos 0

sin 0
cosO
COS0
0

Precession theory is intimately concerned with the orientation of the "mean equatorial

coordinate system.

" The

z-axis (denoted Q) of this system is directed to the mean Celestial

Ephemeris

Pole; the x-axis is directed toward the mean vernal equinox and is symbolized

by the

"ram's horns"

of Aries (T); the y-axis (yq

), directed toward right ascension 6h

and declination 0

(The word "mean" in this context

indicates that the effects of nutation are disregarded.) This system is also called the "mean

equator and equinox,

either "of epoch"

or "of date.

The "epoch" is the initial time for

precession, usually J2000.0, and the

"date"

refers to an arbitrary final time for precession.

Figure 1-1 shows the coordinate

axes

of this system,

with circles marking the planes of the

equator and the ecliptic.

Ecliptic

Equator

Figure 1-1.

Equator
-

The Equatorial Coordinate System

- -

E

I

system; and again the y-axis completes the right-handed triad.

The basis unit vectors are

denoted T for the

x-axis

, YE for y, and E for

Figure 1-2 depicts the ecliptic coordinate

system.

The vectors E and Q appear in

both figures; the angle between them is

obliquity of the ecliptic.

For the sake of brevity, Eichhorn

's (1987

) single-letter notation for coordinate

teams will often be used.

The mean equatorial system of date is called the "Q system," and

the ecliptic coordinate system of date is likewise called the

"E system.

" Their counter-

parts at epoch carry a subscript zero: Qo and E0,

respectively.

The orientation of the Qo

system is defined implicitly by the right ascensions and declinations at J2000.0 of the 1535

fundamental stars in the FK5 star catalog (Fricke et al. 1988).

Two other equatorial systems are still in frequent use.

The "B1950 (FK4) system

is defined by the positions and proper motions of the FK4 star catalog (Fricke and Kopff

1963) at epoch B1950.0 (the beginning of Besselian year 1950, or Julian Ephemeris Date

2433282

.42345905).

Observations subsequent to those used in the FK4 have revealed a small

but significant systematic error in the FK4 proper motions in right ascension.

Consequently

this system, although intended to be inertial, is in slow rotation about its

z-axis.

(No doubt

future observations will reveal a similar flaw

, albeit of smaller magnitude,

in the FK5.)

By contrast

the planetary ephemerides developed at the Jet Propulsion Laboratory

(JPL) are integrated numerically in a coordinate system that is inertial by construction.

The latest available ephemeris, produced for the

Voyager project after the Neptune en-

counter, has its coordinate

B1950.0.

axes

aligned

This non-rotating system

the B1950 (FK4) system above at epoch

will be denoted

"EME50" in accordance with JPL

Table 1-1.

Definitions of Symbols

Symbol

Section

Definition

The rotation matrix which transforms from the mean ecliptic and
equinox of J2000.0 into the mean equator and equinox of date.

4.2.3

1) The smallest principal moment of inertia of the Earth.
2) Subscript denoting the accumulated precession angle rather than
its rate.

The semimajor axis of an orbit, usually with a subscript to indicate
the orbiting body.

Coefficients of Chebyshev polynomials obtained from JPL planetary
ephemerides.

Prefix denoting a Besselian year.

B1950

B1950.0

C

4.2.3
2.4

1.3

1.3

4.2.3

1) The intermediate principal moment of inertia of the Earth.
2) Subscript denoting the Earth-Moon barycenter.

The mean equatorial coordinate system at epoch B1950.0; used with
a following (FK4) to denote the rotating system of the FK4
star catalog.

The instant of time corresponding to the start of the Besselian year
1950; Julian Ephemeris Date 2433282.42345905.

The largest (polar) principal moment of inertia of the Earth.

1) The speed of light in vacuo.
2) The product sin irA cos 11A

3.3.1

The coefficient of Tk in the expression for the cosine of the angle I.

The Jacobian

of the transformation from rectangular coordinates

into right ascension and declination.

The polynomial forming the denominator in the quotient n/d.

The coefficient of Tk in the denominator polynomial d.

Table 1-1, continued.

Symbol

Section

Definition

The coordinate system defined by the ecliptic and mean equinox of
date.

The coordinate system defined by the ecliptic and mean equinox of
epoch.

The orbital eccentricity, often used with a subscript to denote the
orbiting body.

EME50

G

The nonrotating coordinate system aligned with the B1950 (FK4)
system at epoch B1950.0.

The universal constant of gravitation.

The Jacobian of the transformation from Set III

coordinates into

rectangular coordinates in the orbital system.

4.3.3

2.4

2.5

2.6
4.2

3.3.1

2.4
9 _;

The ratio of the moments of inertia of the Earth.

An angular momentum vector.

The unit vector in the direction of the orbital angular momentum
of a planet.

1) The magnitude of angular momentum.
2) One of the two rectangular components giving the eccentricity
and longitude of perihelion of the Earth's orbit.
3) The stepsize in a numerical integration.

1) The rotational moment of inertia of a body.
2) The inclination of the invariable plane to the mean equator of
date.

The inclination of the invariable plane to the mean equator of epoch.

The coefficients of Tk in the polynomial expansion of I.

1) Running index over all bodies on a JPL planetary ephemeris file.
9'\ Tb0i ;nn-n +r 1 4, .- 1 ..... k:s- 1... L ri. T rTr n --- -1.

Table 1-1, continued.

Symbol

Section

Definition

Prefix denoting a Julian year.

J2000

The mean equatorial coordinate system at epoch J2000.0, as realized
by the FK5 star catalog.

J2000.0

The standard epoch for precession theory:

Julian Ephemeris Date

2451545.0, or 2000 January

1 12:00 dynamical time.

1) The Gaussian gravitational constant.
2) A running index, or the power of time T.
3) One of the two rectangular components giving the eccentricity
and longitude of perihelion of the Earth's orbit.

4.3.3

4.3.3

Fundamental rate of luni-solar precession attributable to the Sun.

Fundamental rate of luni-solar precession attributable to the Moon.

The angle, measured along the equator of date, from the mean vernal
equinox of date to the intersection of the invariable plane and
the equator of date; the right ascension of this point.

Same as above, but using the equator and mean vernal equinox of
epoch.

3.3.2

2.5

2.5

4.3.3

The coefficients of T in the polynomial expansion of L.

The rotation matrix which transforms from the EME50 system into
the orbital system (p, q, h) of a planet.

The mean anomaly of a planet.

The lunar coefficients in Kinoshita's expression for luni-solar preces-
sion.

Mass.

The reciprocal mass of a planet (mo/m).

The rotation matrix which transforms from mean pnnatnrial rnvr-

Table 1-1, continued.

Symbol

Section

Definition

The vector directed

toward the

ascending node of the ecliptic of

J2000 on the ecliptic of date (for T
scending node if T < 0.

> 0) or toward the de-

vector directed

toward the

ascending node of the invariable

4.3.3

3.3

plane on the mean equator of J2000.

1) The degree of the highest Chebyshev polynomial whose coefficient
is included in the data records of a JPL planetary ephemeris
file.
2) The polynomial forming the numerator in the quotient n/d.
3) The mean motion of a planet.

The coefficient of Tk in the numerator polynomial n.

4.3.3

The rate of regression of the nodes of the Moon
(a negative quantity).

Order of terms which are omitted from an equ

s orbit on the ecliptic

ation, as in O(T5).

The rotation matrix which transforms from the mean equator and
equinox of epoch to the mean equator and equinox of date; the
precessionn matrix."

The analog to P above with the intersection of the invariable plane
and equator replacing the vernal equinox.

4.3.1

Newcomb's

"Precessional

Constant."

This is not the same

as the

so-called

4.3.1

"constant of precession

p below.

The coefficient of Tk in the polynomial for Newcomb's

"Precessional

4.3.1

Constant."

The unit vector directed toward perihelion of a planet.

1) One of the two rectangular components giving the inclination and
node of the ecliptic of date on the ecliptic of J2000.
2) The speed of general precession in longitude, also known as the
"constant of precession."

Table 1-1, continued.

Symbol

Section

Definition

The speed of general precession in longitude at J2000.

The rotation matrix which

transforms from the mean equator and

invariable plane of J2000 to the mean equator and vernal equi-
nox of date.

The vector directed toward the Celestial Ephemeris Pole.

The coordinate system defined by the mean equator and vernal equi-
nox of date.

The coordinate system defined by the mean equator and vernal equi-
nox of epoch.

The unit vector directed

toward 90

true anomaly in

the orbital

plane of a planet.

1) The quotient of two polynomials, q

= n/d.

2) One of the two rectangular components giving the inclination and
node of the ecliptic of date on the ecliptic of J2000.

The coefficient of Tk in the quotient of two polynomials.

Ri(9)

The elementary rotation matrix which rotates the coordinate system

by angle 0 about axis i, where i
z-axis.

1) The radius of a rotating body.

= 1,2,

or 3 for the

, y-, or

4.3.3

2) Kinoshita's

rate of luni-solar precession.

The position of a body relative to the barycenter of the Solar System.

1) The magnitude of an arbitrary vector.
2) The magnitude of the difference of the barycentric position vectors
of two bodies.

The covariance matrix of the Set III elements and mass of a planet.

Table 1-1, continued.

Symbol

Section

Definition

The covariance matrix of the angular momentum of a planet referred
to B1950.0 equatorial coordinates.

The covariance matrix of the orbital angular momentum of the Solar
System in terms of the right ascension and declination of the
angular momentum vector.

4.3.3

The solar coefficients in Kinoshita's expression for luni-solar preces-

sion.

1) In cylindrical coordinates, the distance from a point to the

z-axis.

) The quantity sin 7rA sin HA*

The rotation matrix which transforms from the EME50 system (as
realized by the DE130 planetary ephemeris) into the J2000
system (as realized by the DE202 planetary ephemeris).

The transpose of a matrix or column vector.

Time, measured in Julian centuries past J2000.

The central time in a time interval of the table of long-term Cheby-
shev coefficients.

1) The earliest time covered by one record of a JPL planetary ephem-
eris file.

) The initial

time when one desires to apply precession

between

two arbitrary times.

1) The latest time covered by one record of a JPL planetary ephem-
eris file.
2) The final time when one desires to apply precession between two
arbitrary times.

Tk(x)

The Chebyshev polynomial of the first kind of degree k.

The time interval T2 T1

when one desires to apply precession be-

tween two arbitrary times.

The velocity of a body relative to the barvcenter of the Solar Svstem.

Table 1

-1, continued.

Symbol

Section

Definition

The first Cartesian coordinate axis, or the projection of a vector in
that direction.

The unit vector in the y direction of the E system.

The unit vector in the y direction of the Q system.

The second Cartesian coordinate axis, or the projection of a vector
in that direction.

1) The third Cartesian coordinate axis, or the projection of a vector
in that direction.
2) The accumulated angle, measured along the equator of date, from
the intersection of the equator of date and equator of epoch

to the y-axis of the mean equatorial system of date

denoted

this is

as ZA by Lieske et al. (1977)

Another notation for ZA

above.

The coefficients of Tk in the polynomial expansion for

ficients

the coef-

z3 are denoted z4 and z' respectively by Lieske

et al. (1977).

1) Right ascension, in particular of the angular momentum vector
of the Solar System.
2) Any one of the precession angles.

The right ascension of the angular momentum vector of the Solar

System in the Qo

system.

The coefficient of Tk (a polynomial coefficient) for any one of the
precession angles.

The coefficient of Tk(

/3 1.3

7 2.6

r) (a Chebyshev coefficient) for any one of the

precession angles.

The generalized latitude angle of a vector.

The coefficient of moment of inertia of a body.

Table 1-1

continued.

Symbol

Section

Definition

3.3.3

The coefficient of Tk in the polynomial expansion for A.

The difference between the inclination of the invariable plane to the
true equator of date and the inclination of the invariable plane
to the mean equator of date.

The Set III parameter specifying a rotation of an orbital plane about
its line of apsides.

The Set III parameter

specifying a rotation of an orbital plane about

the latus rectum.

Nutation in obliquity.

Nutation in longitude.

Declination

, in particular of the angular momentum vector of the

Solar System.

The declination of the angular momentum vector of the Solar
in the Qo system.

1) The difference between the angle A (q.v.)

System

as computed from the

rigorous equation and from the polynomial approximation.
) The difference between the angle A measured to the true equator
of date and the analogous angle measured to the mean equator
of date.

The difference

between

the angle

(q.v.) as computed

from

rigorous equation and from the polynomial approximation.

The difference

between

the angle

L (q.v.) as computed from

rigorous equation and from the polynomial approximation.

The obliquity of the ecliptic; the inclination of the ecliptic of date
to the equator of date.

4.3.1

The obliquity of the ecliptic at J2000.

The accumulated angle, measured along the equator of J

2000,.

from

/

Table 1-1, continued.

Symbol

Section

Definition

Another notation for (A above.

The coefficients of Tk in the polynomial expansion for (A; the coef-
ficients (2 and (3 are denoted ({ and (' respectively by Lieske
et al. (1977).

1) The longitude angle in cylindrical coordinates.
2) The dihedral angle between the equator of J2000 and the equator
of date; the angle between the Celestial Ephemeris Poles of
J2000 and of date; denoted OA by Lieske et al. (1977).

Another notation for OA above.

The coefficients of Tk in the polynomial expansion for 6A

the coef-

ficients 02 and 03 are denoted 09 and 0(' respectively by Lieske
et al. (1977).

The angle, measured in the ecliptic of date, from the mean vernal
equinox of date to the intersection of the ecliptic of date and
the ecliptic of J2000; denoted AA by Lieske et al. (1977).

The generalized longitude angle of a vector.

The rest mass of a body multiplied by G.

The relativistic mass of a body multiplied by G.

The angle, measured in the ecliptic of J2000, from the mean vernal
equinox of J2000.0 to the intersection of the ecliptic of date
and the ecliptic of J2000.0; the J2000 ecliptic longitude of the
ascending node of the ecliptic of date; denoted HA by Lieske
et al. (1977).

The angle between

the ecliptic of date and the ecliptic of J2000,

taken to be negative for T

< 0; denoted

A by

Lieske et al.

(1977).

The longitude of the Earth's perihelion in the E0 system.
mil 1 P r i l j

Table 1-1, continued.

Symbol

Section

Definition

4.3.2

4.3.1

4.3.1

The colatitude angle of the Celestial Ephemeris Pole referred to the
ecliptic and mean equinox of J2000.0.

The instantaneous speed of luni-solar precession; the rate at which
the mean vernal equinox of date moves along the moving eclip-
tic of date.

The speed of luni-solar precession at J2000.

The angle, measured along the ecliptic of J2000, from the mean ver-
nal equinox of J2000 to the intersection of the ecliptic of J2000
and the equator of date; the accumulated luni-solar precession;
denoted 'A by Lieske et al. (1977).

4.3.1

4.3.1

4.6

The instantaneous speed of planetary precession; the rate at which
the mean vernal equinox of date moves along the moving mean
equator of date.

The speed of planetary precession at J2000.

The angle, measured along the equator of date, from the vernal equi-
nox of date to the intersection of the mean equator of date and
the ecliptic of J2000; the accumulated planetary precession;

denoted XA

4.3.1

by Lieske et al. (1977).

1) The rotation rate of a body.
2) The rotation rate of the Earth.

The argument of perihelion

of a planetary orbit,

referred

to the

EME50 coordinate system.

The inclination of the ecliptic of J2000 to the mean equator of date;
denoted WA by Lieske et al. (1977).

The Sun.

The planet Mercury.

The planet Venus.

Table 1-1, continued.

Symbol

Section

Definition

The planet Jupiter.

The planet Saturn.

The planet Uranus.

The planet Neptune.

The planet Pluto.

The Earth's Moon.

The unit vector directed toward the mean vernal equinox of date.

The longitude of the ascending node.

CHAPTER

THE DETERMINATION OF THE INVARIABLE PLANE

2.1. Introduction

The invariable

plane of the Solar System is rigorously defined

as that plane which

contains the center of mass of the Solar System and is perpendicular to the total angular

momentum of the Solar System.

This definition holds for both classical and relativistic

physics.

The "total angular momentum" must rigorously include rotational angular momentum

as well as orbital angular momentum, and it must account for all the mass in the Solar

System: satellites, asteroids, and comets in addition to the Sun and planets. Nevertheless,

the major contributor to the total angular momentum of the Solar System is the orbital

angular momentum of the planets, in particular of Jupiter.

The goal of this chapter is to determine the orientation of the invariable

equivalently to estimate the direction of the total angular momentum.

plane, or

much easier thanks to the availability of computer-readable planetary ephemeris files, and it

can provide reliable results now that the masses of the outer planets have been determined

from spacecraft encounters.

The angular moment of the largest asteroids must be included.

The largest asteroid

x 10-10 solar mass (Schubart 1974; Newhall et al. 1983) or 6

x 10-7

of about five less than this:

This is comparable to the

uncertainty found in Section 2.5 and therefore must be included.

Smaller asteroids individually will have an effect correspondingly less.

Taken as an

ensemble, their angular momentum should lie very close to the normal to the invariable

plane; thus while the magnitude of the total angular momentum of the Solar System would

be increased, the direction of the total angular momentum would be virtually unchanged.

So small asteroids, and by implication comets and smaller bodies, can be safely ignored at

this time.

Rotational angular momentum will also be ignored, but for an entirely different reason:

the large uncertainty in the Sun's

rotational angular momentum.

This topic is explored in

more depth in Section 2.6.

The approach adopted here will be to compute the orbital angular momentum of each

planetary system (comprised of the planet and its satellites) under the assumption that

each planetary system is a point mass. Accordingly, this chapter will first describe the most

recent planetary ephemeris file produced at the Jet Propulsion Laboratory.

The relativistic

equations giving the orbital angular momentum of the Sun and planets are presented next.

Then the total angular momentum will be found and its uncertainty estimated.

For most planets, the product of the universal gravitational constant G and the body's

mass is known much more precisely than the mass itself.

This arises because only their

product (typically denoted by t) enters into the equations of motion.

The value of G must

be determined in the laboratory by measuring the gravitational force between two bodies

of known mass; modern measurements (e.g., Luther and Towler 1982) give its value to only

22

Throughout this remainder of this chapter, the word "mass" will therefore be taken

to mean

- Gmi rather than m, itself.

"Mass"

accordingly is measured in

units of

length3 timee.

passing,

the square of the Gaussian constant of gravitation k gives

the solar mass, exactly

1 AU3/day2

if p is specified in km3/sec2

, the length of the AU

follows.)

JPL Planetary Ephemerides

The planetary ephemeris files produced by the Jet Propulsion Laboratory provide the

foundation for all NASA's deep-space missions.

In addition to providing the position and

velocity

of the Sun,

Moon, and nine

planets,

these files also include

the values of the

physical constants that were used in their generation.

These constants include the speed of

light, the length of the astronomical unit (expressed in kilometers), and the masses of the

planets.

The constants are read by the programs that integrate the equations of motion of

natural satellites or of a spacecraft; this insures compatibility between files.

Planetary ephemeris files are produced by the "Solar System Data Processing System"

(SSDPS), which was originally designed by Charles L. Lawson (1981) and is now maintained

E. Myles Standish, Jr., and X

X Newhall.

The SSDPS has three main functions:

reduces astrometric observations (both ground-based and from spacecraft) to determine a

set of initial conditions for the Solar System bodies; it integrates the equations of motion;

and it transforms the integrator output into an easily interpolated ephemeris file.

reader is referred to Newhall et al. (1983), Standish (1990a), or Newhall (1989) for a more

complete description of the SSDPS; an overview of the system for our purposes here is

presented below.

pole at a standard

epoch,

either

B1950.0 (as realized

by the

FK4 catalog) or

J2000.0

(FK5).

The x-axis

is directed to the corresponding (FK4 or FK5) "catalog equinox"

zero point for measuring right ascensions) at that epoch.

Since the equations of motion

integrated by the SSDPS are expressed in an inertial

coordinate system (i.e.,

there are

no Coriolis accelerations in the differential equations),

construction, an inertial coordinate system. When the

the resulting ephemeris realizes,

B 1950.0 epoch is used, the resulting

inertial coordinate system is termed EME50 (an abbreviation for "Earth mean equator of

1950"

) in order to distinguish it from the slowly rotating coordinate system of the FK4

catalog

The first

step in

the production of a new planetary

ephemeris

is the estimation of

the masses and initial positions and velocities of the principal bodies in the Solar System.

The a priori model is an earlier ephemeris.

The data (Standish 1990a) include ground-

based optical observations,

radar ranging data, and spacecraft ranging data.

The latter two

classes provide accurate synodic periods independent of transit circle observations; this in

turn gives the planets' inertial mean motions. Since transit observations are made relative

to fundamental stars, the offset and drift of the fundamental catalog (FK4 or FK5) vernal

equinox are therefore visible in the residuals

, and these quantities can be estimated.

The second step in the preparation of an ephemeris is the numerical integration.

equations of motion (Newhall et al. 1983) are fully relativistic in

their treatment of the

gravitational forces. Each planet other than the Earth is taken to be a point whose mass

the sum of the masses that comprise that planetary system. Each point mass therefore de-

fines the barycenter of its planetary system.

(This

scheme is consistent with the other JPL

24
Several asteroids are included as perturbing bodies; their orbits are given analytically, but

they themselves are not integrated.

The Earth and Moon, however, are included as separate bodies.

Here the low-order

gravity harmonics of both Earth and Moon are modeled, as are the Earth tides and lu-

nar librations.

These accelerations are vital for the production of an accurate geocentric

ephemeris for the Moon.

Therefore the SSDPS integration is effectively of an eleven-body

system: the nine planets, the Sun, and the Moon.

The integrated positions and velocities

are written to a file for further processing.

The third step produces the planetary ephemeris file itself.

To save storage space, and

to make interpolating easier, a planetary ephemeris file contains a Chebyshev representation

of the motions of the planets (Newhall

1989).

The final ephemeris file contains many

records,

with each record covering a fixed

time interval.

Each record in

turn contains

coefficients of Chebyshev polynomials used in the interpolation of position and velocity. If

T1 < T < T2, where T is the desired interpolation time and T1 and T2 are the start and

end times of the record bracketing T

, then the x component of position is found by

x=Z

akTk(T),

(2-1)

where the Tk are the Chebyshev polynomials of the first kind (see, for instance, Rivlin

1974); r

, in the range |r|

< 1, is a dimensionless measure for time defined by

- (Ti + T2)

(2-2)

-Ti)

the ak are the coefficients read from the file: and n is the decree of the highest nnlvnnminal

dx
dT

dr
akTkMr dT

(2-3)

The y and z components are found in the same manner. In practice, the values of Tk(r) and

Tk(r) are found recursively in a low-level subroutine.

Users need only call a higher-level

subroutine, passing the time and desired planets, and retrieving the position and velocity

vectors.

The Chebyshev coefficients are fit to positions and velocities output from the numeri-

cal integrator; the fitting process, which is constrained to match both position and velocity

across a record boundary, is described by Newhall (1989).

The Chebyshev representation

has been found to match the integrated trajectory to within a tolerance of 0.5 mm in posi-

tion; this is many orders of magnitude below the accuracy to which the planets'

positions

are currently known.

The M04786 Planetary Ephemeris

The particular ephemeris file used in

this work, known internally

as M04786 after

the identification number of the run which created it, has a slightly different history from

the scheme outlined above. It was prow

Determination Program (ODP; Moyer

duced by Jacobson et al. (1990) using JPL

1971) in

's Orbit

the final steps of postflight analysis of

the 1989 Voyager encounter with Neptune (Stone and Miner 1989).

Its precursor, known

within JPL as DE130, was created by Standish (1987) to serve as an a priori ephemeris

for the Neptune encounter. In accordance with Voyager project requirements, DE130 (and

therefore M04786) is on the EME50 system; a J2000 equivalent is known as DE202.

These

ephemerides supersede DE200, which is still used in the production of the annual volumes

These five asteroids were found

to have the largest effect on

Viking lander ranging

observations.

As Voyager encountered Neptune, the data it collected (both radiometric and optical)

helped

to render Neptune's ephemeris more accurate.

The planetary ephemeris file used

during encounter operations was accordingly updated from

time to time.

Due to time

constraints,

the full SSDPS was not executed.

Rather, Neptune's orbit and mass were

merely linearly corrected based on parameters estimated by the ODP.

The masses and other physical constants used by the M04786 ephemeris are presented

in Table 2-1. As mentioned above, the "masses" all include a factor of G, and the satellites

of the superior planets are included with their primaries.

Physical Constants from the M04786 Ephemeris

Mass of Mercury, p/
Mass of Venus, p#
Mass of Earth, Fe
Mass of Mars system, Pc
Mass of Jupiter system, #p

Mass
Mass

of Saturn system, ph
of Uranus system, Pi6

Mass of Neptune system, p,
Mass of Pluto plus Charon, PB

Mass
Mass

of the Sun, pe
of the Moon, p(

22032.08045038011 km3/

324857.4782

760278 km3/

398600.4420277941 km3/

4282

8.28654533971 km3/sec2

126712700.9827376 km3/sec2
37940448.53389772 km3/sec2
5794559.117031542 km3 /sec2
6836534.716231512 km3/sec2
1020.864919671097 km3/sec2
132712439800.9096 km3/sec2
4902.799066232991 km3/sec2

Mass of Earth plus Moon, pB
Mass of 1 Ceres

Mass of

Pallas

Mass of 4 Vesta

Mass of 7 Iris
Mass of 324 Bamberga
Length of astronomical unit
Speed of light in vacuo, c
Earth/Moon mass ratio, Pe/Ht

403503.241094027

1.54674588
3.8448884762

7.2858525

10-13
10-"14
10-14

1.6 x 10-15
2.6 x 10-15

0 km3/

AU3/day2
AU3/day2
AU3/day2
AU3/day2
AU3/day2

149597870.6094344 km

299792.458
81.300587

km/sec

Table 2-1.

The location of the Solar System barycenter is defined implicitly by equation (2) of

Standish et al. (1976):

(2-4)

where r2 is the position of body i relative to the Solar System barycenter and p

relativistic mass (corrected for both special and general relativity effects).

is its

One calculates

p* using equation (3) of the same report:

2_ 1
- -
2c2 2c2

rir

(2-5)

In this equation, pi is the rest mass of body

Vi =

vii = i'rJ is the speed of body i relative

to the Solar System barycenter; c is the speed of light; and ri = |r,-rj[ is strictly speaking

the magnitude of the difference of the barycentric position vectors of body i and another

body j. Newhall et al. (1983) point out that equations (2-4) and (2-5) are interdependent;

in practice, the Sun

s position is not integrated at all but is inferred from equation (2-4).

The orbital angular momentum of body i is then found by (Burkhardt 1982):

(2-6)

hi= *ri

and the total orbital angular momentum is

X Vi.

(2-7)

The sixteen bodies are the nine planets, the Sun, and the Moon, all interpolated from the

planetary ephemeris file; and the five asteroids,

whose positions and velocities are found

16
,ir.1
i=1

pzri = 0

this will presumably give the most accurate position and velocity for Neptune.

ing positions,

The result-

velocities, and angular moment from equation (2-6) of the eleven bodies on

the file are presented in Table

relativistic masses, from equation (

those of the five

asteroids appear in Table

2-5), are presented in Table 2-4.

In Table

-3; and the

the solar

position and velocity are found from equation (

The total of these angular moment

5.5138791989915901
-8.1599811430560712
1.9241747810631041

x 1016
x 1017
x 1018

km5/

The spherical coordinates of h are

.0907754056836565

x 1018

km5/sec3

(2-9)

= 2732865725677

660972373413

= 6658'20'54429.

The right ascension and declination here are still expressed in the EME50 coordinate system

of the ephemeris; the uncertainty in these numbers is discussed in the next section.

A straightforward application of equation (

7) will not yield an unvarying result.

nodes on the ecliptic of the Moon's geocentric orbit regress with a well-known period of

18.6 years; this regression is the result of the gravitational influence of both the Earth's

oblateness and the Sun.

Since the Moon

s geocentric orbital plane thus changes,

the con-

tribution of the Moon to h also varies.

(The direction of the Earth's

rotational angular

momentum varies due to nutation with the same period-an equal and opposite effect-

thus ensuring that the total angular momentum is conserved.) Therefore it is preferable to

= 18h15m27.774162:

Table

Vectors from the M04786 Ephemeris at 1989 August 25 04:00 ET

Body

Mercury

Venus

-16158199.07503421
-60540281.38189575
-30792661.60242493

x -48044757.13184071
y -89820030.22659636
z -37464034.81962667

37.68574922489486
-6.219471575576255
-7.199241280123081

31.19917668472052
-13.50936340740929
-8.058115422625290

5.383098408817692
-2.812992786916346
5.248044464575123

7.071007368738365
-5.054773396094979
1.121201755436503

x 1012

x 1013

x 1013

Earth

132534047.2682675
-66391831.78962909
-28790114.33862247

13.80094834015177
23.89128363392751
10.35901588646045

3.138869950535301
-7.056237004025107
1.627357555590895

x 1010
x 1014

Mars

Jupiter

Saturn

-244898357.5855488
35884297.47698200
23023213.81212496

71436178.22611272
699498249.8393588
298310659.1621094

303660932.6890458
-1353840709.129601
-573039586.0391354

-3.157305918724009
-19.87024113963734
-9.036182675380297

-13.17134485303879
1.568629228634501
0.994759615211981

8.938087519218495
1.928822496211554
0.411069166116023

5.705573328457442
-9.788994312814919
2.132629013701246

2.887703775435171
-5.068779699170997
1.181645297602056

2.082055996257906
-1.990622886004587
4.813297847774068

x 1013
x 1014

x 1016
x 1017

X 1016
x 017
x 1017

Uranus

176076406.2304637
-2646457517.863089
-1162055888.894550

6.742392857724443
0.124628946369494
-0.040540240190316

1.460888460184428
-4.535922436215454
1.035221192520042

Neptune

Pluto

x 843360977.1913863
y -4101540837.662588
z -1701376011.969666

x -3070732451.816337
y -3201106574.346603
z -80576284.87092758

5.305123897076706
1.016154761387040
0.282701827614558

4.079629959404747
-3.842341036034878
-2.450284528360149

3.892369040033341
-6.333659888390561
1.546162148446206

7.691217076519461
-8.016739555119871
2.537679451419032

x 1015
x 1016
x 017

x 1012
x 1012

x -206331.5970958987

46449.36015359861
17559.05332732530

0.009329389609560
-0.002139033914718
-0.001087745178433

-1.720703970228159
-8.045154208580101
1.062509844429565

Moon

132604706.2412162

12.79303203667523

8.598343384697650

x 1015

x 1012
x 1012

Table 2-3.

Asteroid Vectors at 1989 August 25 04:00 ET

Body

Ceres

186281992.1788968
343470578.6841611
123850658.9230570

436755840.8648572
38833040.99610688
-39639922.64470550

Pallas

-16.35952815373246
5.116105802899675
5.747722186944968

-5.507766068499227
15.64886644512920
-2.692689866417276

92992435509.65316
-214825510453.6537
455899280135.5158

8893557968.250119
24044304801.65627
121544809850.9764

Vesta

125506706.2166484
-275465663.8456852
-126367091.5783639

19.46247210077522
7.392776730344995
0.401434801390846

26912693424.32996
-82010229009.18285
205502281954.1065

x -409332580.4006283
y -24051697.50313012
z -51893151.65421642

-1.124700364629995
-15.12047499555895
-6.362926495829524

-453228983.1614052
-1827087511.594522
4421892639.311001

Bamberga

x -524855815.2412224
y -53995531.06175173
z -99176075.74679006

2.253161649691973
-10.70612252024783
-6.668306875225594

-818264398.8558691
-4341672429.418912
6694172275.600386

Table 2-4.

Relativistic Masses from the M04786 Ephemeris at 1989 August 25 04:00 ET

Body

Mercury
Venus

(km3/sec2)

22032.08040253906
324857.4782712153

Earth
Mars system
Jupiter system
Saturn system
Uranus system
Neptune system
Pluto plus Charon
Sun
Moon
Earth plus Moon

Ceres
Pallas
1 7n/, 4-.-

39860
42828
12671
37940
57945
68365

0.4418123855
.28653388460
2700.9849606
448.53290854
59.117018051
34.716225980

1020.864919713455
132712439800.7594
4902.799065514813
403503.2408779003

69.36936487995916
17.24378096264524
00 C/?fl'7/? oflflA cn7

Pere + P(r(

2-12)

and is assigned a mass

PB = Pe +[k.

When the Earth-Moon barycenter is treated

as a single fictitious

body,

the total orbital

5.513

8792089473392

-8.1599811342261851

1.924174

7793873432

x 1016
x 1017
x 1018

km5/

with corresponding spherical coordinates

= 2.0907754037994348

= 273?865725688

= 66?972373417

X 1018 km

= 18hl5m27.774165;

= 6658'20'.'54430.

Figure 2-1 shows the difference between the vector h computed by equation (2-

7) and

that computed using the Earth-Moon barycenter instead of the Earth and Moon individu-

This plot was actually generated from the longer DE130 ephemeris in order to show

the effect more clearly; the results from the M04786 ephemeris are nearly identical.

18.6-year oscillation, with an amplitude of 17 parcsec,

dominates. An annual perturbation

in the lunar orbit is also apparent.

Neither equation (

) nor equation (

-14) gives the desired result

. If the M04786

ephemeris is interpolated at various other times, h is seen to be a function of time, with

^ +#

-?2 WLW.

were integrated using the DE130 value of Neptune's

mass, 6828879.09654 km3/

(Stan-

dish 1987).

The Voyager data cause the estimate of Neptune's mass to increase by 0.112%,

to the value shown in

Table 2-1.

When

the estimate of Neptune's

mass increased

gravitational attraction by

Neptune on Jupiter ought also to have increased.

However,

the ODP corrected only Neptune's

ephemeris; the operational software cannot change the

orbits of the other planets when the mass of one planet changes.

Therefore the M04786

ephemeris contains, in effect, a pair of forces (Jupiter on Neptune and vice versa) that are

not equal and opposite. The total angular momentum of the system is therefore not strictly

conserved

The magnitude of the imbalance varies with the twelve-year synodic period of

Jupiter and Neptune.

Figure

2-2 displays the right

ascension and declination of h

as a function of time

the entire sixteen-year span of the M04786 ephemeris.

The twelve-year periodicity is quite

obvious in the figure. (Corresponding plots using the DE130 ephemeris, shown to the same

scale in Figure 2-3, are quite stable by contrast.) The average value of these coordinates

will be taken as the direction of the total orbital angular momentum of the Solar System;

the magnitude is given by equation (2-

Consequently we have

5.5138790907107613
8.1599810902131458
1.9241747812877163

1016
1017
o107
io'8/

km5/sec

.0907754037994348

= 2732865725626

= 662972373550

x 1018 km5

= 18h15m277774150;

= 6658'20'.54478.

18hl5m27 743

27?742

277741

18h 15m27.7740

6658'20':.'5455

0

O
20t'5450

0
*j

o 20"5445
")

AA0S>R'20"SzthO4.

18hl5m27.7816

27.7815

27.7814

18h15m27.7813

6658'20'3115

o
o 20'3 110

0
4.*)

o 2023105
U)

R0FtR'9fl"R i f

-- -. ..T- I | -I J I I I I I I I I

I

The Uncertainty in the Total Orbital Angular Momentum

The uncertainty in

the total orbital angular momentum h of the Solar System, as

found in the previous section, may be estimated through applying the standard formulas

of error propagation (Bevington 1969) to equation (2-7).

The SSDPS solves not for initial

position and velocity, but rather for changes to the "Set III" orbital elements of Brouwer and

Clemence (1961): Aa/a, Ae, Ap, Aq, eAw, and

Aw + AM.

These parameters represent

the fractional change in semimajor axis; change in eccentricity;

rotations of the orbital

change in the mean argument of latitude.

These are referred to each planet's osculating

orbit at a standard epoch; the transformation from Set III elements into changes in the

position and velocity vectors is well defined.

The orbital angular momentum for a planet is given by

h = /ji/(/p +u)a(1 -e2)h,

(2-21)

Where h, the unit vector in the direction of the positive orbit normal, defines the third as
where h, the unit vector in the direction of the positive orbit normal, defines the third axis

of the orbital coordinate system.

The first and second unit vectors that define the orbital

coordinate system are respectively the normalized "Laplace vector" p, directed toward the

perihelion, and 4 h x p.

The rotation matrix from EME50 equatorial coordinates into

the orbital system (p, q, h) is

M = R3(w) RI(i

(2-22)

Let the

covariance matrix from

the SSDPS for a planet

be denoted by SIm;

37
where the matrix H consists of partial derivatives of h with respect to the seven parameters:

9{hp, hq, hh}

O{p,Aa/a,Ae,

Ap,Aq

,Aw+AM}

0
= 0
h/p+ h/2(g

(2-24)

- e2)

Finally,

the uncertainty in h must be rotated into EME50 equatorial coordinates,

STyz

MT SpqhM.

giving

(2-25)

The orbital elements here can

be taken to be the osculating elements at any epoch

(which are easily accessible from converting position and velocity vectors) without introduc-

ing undue error in the calculation; after all, the goal is merely to determine the uncertainty

in h.

The major contribution to Suz will be from the uncertainties in the planetary masses.

These uncertainties are usually expressed in units of inverse solar masses, so that a planet

inverse mass

_= P/pi.

Accordingly

- aim

Values for

omj, from Standish (1990b),

are presented in Table

along with the

resulting a,1 from equation (

Except for Pluto, the mass uncertainties are typically

five or

six orders of magnitude less than the masses themselves. By contrast

the uncertainty

is Ap or

Aq is typically O'.'01, at most 0"'05 for

Pluto, or O(10-

Although

Table

Planetary Masses and Uncertainties

Body

Mercury
Venus
Earth plus

Moon

6023600
408 523.7
328 900.55

0.025

Mars system
Jupiter system
Saturn system
Uranus system
Neptune system
Pluto plus Charon

The covariance matrix Sxyz

3098 708

1 047.349
3497.09
22902.94
19412.24

135000000

0.001

7000000

for each planet is given in

Table 2-6; as expected,

major contributor is the uncertainty in the masses.

The uncertainty in the total angular

momentum of the Solar System is simply the sum of these matrices, given at the end of

the table, as one can ignore correlations between parameters for one planet and those of

the others.

The inter-planet correlation coefficients for the outer planets (which carry the

bulk of the angular momentum) are rarely above one percent. By contrast, the parameters

for the inner planets are very highly correlated, but the inner planets' masses are so low

and their contribution to h so small that these correlations can safely be ignored.

The final step is to transform the uncertainty in h into uncertainties in the adopted

right ascension a and declination 6 of the normal to the invariable plane.

Since tana =

hy/hx and sin 6 = hz/h, the Jacobian of the transformation is

9{a,6}

-h,/(h,+hh,)

'Vh+ h2

h+/(h + h)
_ z 2 1h242

(2-27)

h +h/h2 2

9{hs,hyhz

-hxhz/h

Table 2-6.

Covariances of the Planets' Orbital Angular Momenta

Body

x (km10/s6)

y (kml/s6)

z (kml/s6)

Mercury

Venus

5.0008
-2.6084
4.8705

1.2601
-8.5806
1.8988

x 1016
x 1017
x 1017

x 1015
x 10is
x 1016

-2.6084
1.3606
-2.5405

-8.5806
6.1258
-1.3555

X 1017
x 1018
x 1018

x 1015
x 1016
x 1017

4.8705
-2.5405
4.7436

1.8988
-1.3555
3.0092

x 1017
x 1018
X 10ls

x 1016
x 017
x 1017

Earth

1.6179
2.8183
1.2981

3.3970
-5.8144
1.2667

Mars

x 1013
x 1012
x 1012

x 1014

x 1016

2.8183
2.0857
-4.1848

-5.8144
9.9629
-2.1703

x 1012
x 1015
x 1015

x 1015
x 1016
x 1017

1.2981
-4.1848
9.9230

1.2667
-2.1703
4.7282

x 1012
x 1015
x 1015

x 1016
x 1017
x 1017

Jupiter

Saturn

Uranus

Neptune

Pluto

Total-

3.4311
-1.4203
3.0757

1.5633
-1.3513
3.2732

6.2578
-1.8624
4.6589

2.5811
-2.3328
5.7654

1.7137
-1.7859
5.6481

1.9075
-3.3045
9.2912

x 1021
x 1022
x 1022

x 1022
x 1023
x 1023

x 1019
x 1020
x 1020

x 1020
x 1021
x 1021

x 1023
x 1023
xl1023

x 1023
x 1023
x 1023

-1.4203
2.3870
-5.4560

-1.3513
1.2945
-3.1272

-1.8624
6.3122
-1.4266

-2.3328
3.8529
-9.3596

-1.7859
1.8612
-5.8863

-3.3045
1.7642
-4.3693

x 1022
x 1023
x 1023

x 1023
x 1024
x 1024

x 1020
x 1021
x 1022

x 1021
x 1022
x 1022

x 1023
x 1023
x 1023"

x 1023
x 1024
x 1024

3.0757
-5.4560
1.2766

3.2732
-3.1272
7.5630

4.6589
-1.4266
3.2619

5.7654
-9.3596
2.2868

5.6481
-5.8863
1.8616

9.2912
-4.3693
1.0963

x 1022
x 1023
x 1024

x 1023
x 1024
x 1024

x 1020
x 1022
x 1022

x 1021
x 1022
x 1023

x 1023
x 1023
x 1024

x 1023
x 1024
x 1025

2.2940 x 10-13
= 2.5407 x 10-14

2.5407 x 10-14
4.4279 x 10-15

(2-29)

The resulting standard errors in the orientation of the total orbital angular momentum of

the Solar System are

oa cos 6S

= 0'103865,

(2-30)

= '0t01373.

The correlation

between a and 6 is +0.7972,

indicating that the

"error ellipse"

is quite

elongated; the semimajor and semiminor

axes

of the error ellipse are 0 .04166 and 0"!01343,

and the major axis has a position angle of 73?53. It is apparent that the effects of the Moon

nodal regression (Figure

-1) and of the unbalanced Jupiter-Neptune couple (Figure 2-2)

are both quite small in comparison to these errors.

The major obstacle to shrinking the error even further is the remaining uncertainty

in Pluto's mass and to a lesser extent in its orbit.

Hubble Space Telescope observations

should resolve Pluto and Charon, thereby giving a reliable semimajor axis of their relative

orbit and improving the estimate of the sum of their masses. Increases in the precision of

Pluto

s orbit estimate can come only after decades of additional observation. Nevertheless,

before the

Voyager

encounter with Neptune, the mass of Neptune was believed known

percent; this alone would have produced an uncertainty of about

25" in the results.

The improvement in the knowledge of the orientation of the invariable plane realized by

the Voyager mission is evident.

The Rotational Angular Momentum of the Solar System

)dm.

(2-32)

In the typical case of solid-body rotation about the

z-axis,

transformation into cylindrical

coordinates (s,

, z) yields

s2Q dm,

(2-33)

where

- 0 is the rotation rate.

For spherically symmetric bodies of radius R,

whose

density p is a function only of the internal radius r

the integral takes the form

7rpr Sdr

(2-34)

for constant p,

= -TrpRi5

MR2Q

- IT.

(2-35)

An evaluation of equation (

-32),

the general

case

will have the same functional form

equation (

2-35),

except that the constant

will be replaced by a coefficient 7 whose value

depends on the body's density and angular speed as a function of position; R and 0 can

be taken to be a reference radius and rotation rate, respectively.

Within the Solar System, the Sun has the largest rotational angular momentum; its

enormous mass (over a thousand times Jupiter's) and radius (nearly ten times Jupiter'

more than compensate for its slow rotation rate.

An approximate evaluation of equation

2-35),

using the polar rotation rate and setting 7

= 0.2 to account in part for central con-

densation, gives h

= 3.5 X 1016 km5/sec3 for the solar rotational angular momentum.

is about one percent of the total orbital angular momentum of the planets and obviously

chn,-,1 B,0 ;nrllsiio in tho r-=lriiltitnn nC +tb0 n rnnl- tfvflnn rn+ Cha foa1 rnxr-1 fl;l Kxi]-n rnl

(xy

model

(cf. Iben

1967),

but due

to uncertainties in

the opacity,

energy generation, and

convection models the results are probably good to not much more than

three or four

digits.

Furthermore, the photosphere of the Sun does not rotate uniformly; the sidereal rota-

tion period varies from 26.2 days at the solar equator to 36.6 days at the poles (Howard and

Harvey

1970).

These rates can only be obtained through measurements of Doppler shifts

at the solar limb; the more precise techniques applied to the planets (analysis of periodic

radio emissions or landmark tracking) cannot be used.

Although recent developments in

helioseismology can in principle determine the rotation characteristics of the solar interior,

results to date are in at best qualitative agreement (for example, Duvall et al. 1984, Brown

1985, and Duvall et al. 1986).

For both

these reasons, the rotational angular momentum of the Sun is not known

to the precision necessary to include it in a precise determination of the invariable plane.

The prudent course of action is therefore to ignore all rotation and to define the "working

model"

which

of the invariable plane (as opposed to the "conceptual"

definition) as that plane

contains the Solar System barycenter and is normal to the total orbital

angular

momentum of the Sun, planets, and largest asteroids.

This decision is also justified on the grounds that rotation (of the planets

as well

of the Sun) does not enter into the equations of motion of the planetary barycenters as

currently formulated in the SSDPS (Newhall et al. 1983).

The Sun is assumed spherically

symmetric; there is no detectable perturbation, even on Mercury's perihelion, due to solar

oblateness.

(Consequently there is also no mechanism

whereby the Sun's spin axis can

>

model"

for h is expected to remain constant; and indeed, as seen in Section 2.4, this is

reasonably true of a planetary ephemeris file.

The two theories (short-term and long-term) that are developed in the next two chap-

ters of this work do not depend on the special physical nature of the invariable plane. Any

plane that does not change its orientation relative to an inertial coordinate system will

suffice. So there is no reason to insist on the inclusion of rotational angular momentum.

The Adopted Orientation of the Invariable Plane

Since rotational angular momentum is not included in the working definition of the

invariable plane,

the total orbital angular momentum found above becomes the normal

vector to the invariable plane.

However,

the vector

h found above was specified in EME50 coordinates.

In order

to transform to J2000, one must not use the standard procedure (Aoki

et al.

1983) for

transforming star coordinates from the B1950 (FK4) system to the J2000 (FK5) system:

because the M04786 ephemeris is already inertial, there is no need for the equinox drift

term. Rather

, Standish (1987) published a 3 x 3 rotation matrix T that transforms positions

and velocities from the DE130 ephemeris (the predecessor of M04786) to the corresponding

positions and velocities in DE202:

+0.9999
= +0.0111
+0.0048

256795509812
814782944231
590037723526

-0.0111814782756923

+0.99993748494

86071

-0.0000271702869124

-0.004859

-0.000027
+0.999988

0038154553
1625775175 .
1946023742

(2-36)

When the vector h from equation (2-18) is multiplied on the left by

vector will be expressed in J2000 coordinates.

T above, the resulting

This is the result we seek:

* '0

Invar.
Plane

atOt

I

Figure 2-4.

The Orientation of the Invariable Plane

The invariable

plane can also be described

by giving the right ascension

L0 of its

ascending node on the Earth's mean equator of J2000 (which direction is defined by n =

where Qo = (0,0,1)T in J2000 equatorial coordinates), and the inclination Io of

the invariable plane to the Earth's mean equator of J2000. Figure 2-4 shows the angles Io

and L0 as well as the spherical coordinates of h.

These quantities are related by

Lo = ao + 6h

(2-39)

3t852572907

(2-40)

lo = 90

(2-41)

= 232008888075 = 2300'31'.99707. (2-42)

= Ohl5m24 .617498;

45
uncertainties permit; accordingly the above results for these angles will be rounded to the

nearest 0'!001, and the rounded values will be adopted for the orientation of the invariable

plane:

= 27351'09t!262 0'!038;

(2-43)

66059'28'.003 0013;

(2-44)

3051'09'.262 0"038;

Lo =
lo =

(2-45)

(2-46)

2300'31.997 + 0!013.

These values will be used throughout Chapters 3-5.

CHAPTER 3
THE SHORT-TERM THEORY

Introduction

The precession matrix P, which transforms from the Qo system of epoch into the Q

system of date, can be built up from elementary rotation matrices in several ways. Lieske

et al. (1977), who present the currently-accepted short-term theory, construct P most easily

in terms of the three angles

and z (Newcomb 1906; Andoyer 1911) by

P = R3(-z) R2(0) R3(-).

These three angles are approximated in that paper by cubic polynomials,

and final times, but in the initial time T

(3-1)

not in the initial

and the difference t between the initial and final

times.

When

the initial time is the standard epoch J2000.0, Lieske et al. (1977) denote

these angles by (A

OA, and

respectively.

The coefficients of the various powers of time

are denoted there as

CA-

CT3

(3-2)

and similarly for the other two angles.

In this chapter,

intermediate epoch

the tilde and subscript

does not enter until

A will

the end

be suppressed for clarity:

of the chapter,

since an

there is no ambiguity

. 1 1 1 1 P /** I

C1T -cl2

47

measured in Julian centuries past J2000. Furthermore, the system of primes and subscripts

used in the paper by Lieske et al. is replaced by a

power of T

system of subscripts only indicating the

Thus within this chapter we will use the notation

+ (4T4

+ O(T5)

(3-3)

and analogously for the other angles.

Figure 3-1 shows the coordinate

axes

of both

teams and the three

"classical"

angles.

Reading the right-hand side of equation (3-1) from right to left,

the rotation Ra(

) places

the initial y-axis, marked yoq0 in the figure, at the intersection of the two equators.

second rotation moves the pole from Qo to Q, and the third rotation puts the y-axis at its

final location, along yq

The ecliptics of epoch and of date are suppressed in the figure for

clarity.

The precession matrix P can also be expressed by the following sequence of rotations:

L) RI(-

A) Ri(Io) R3(Lo).

(3-4)

Figure

3-2 shows the coordinate

axes

of epoch and of date from the same perspective

as Figure 3-1; now the invariable plane and the five angles of equation (3-4) appear.

magnified drawing of the region near the vernal equinoxes also appears as Figure 3-3.

equation (3-4),

the first (rightmost) rotation moves the x-axis to the intersection of the

equator of J2000.0 and

the invariable

plane.

The second rotation puts the y-axis (and

consequently the entire

x-y plane) into the invariable plane.

These two rotations depend

on the orientation of the invariable

plane at the standard

epoch.

The last three

=-(T IT+

= R3(-

I) Ra(-

Eq. of

epoch

T T0 Eq. of date

Figure 3-1.

The Classical Precession Angles

places the y-axis on the equator of date,

and last a rotation about the vector Q (the final

z-axis) positions the x-axis at the vernal equinox of date.

The purpose of the short-term theory is to provide analytical and numerical expressions

for the coefficients of the angles I

, L, and A.

These are obtained by equating the expressions

for P in equations (3

-1) and (3-4) above; solving for I, L, and A in terms of 0o, L0,

and z; and expanding the solutions in powers of time.

The theory is developed to T4 even

though Lieske et al. (1977) go only

as far

asT3

The extra term will give an indication of

the sufficiency of the new theory to model precession with only cubic polynomials:

coefficients of T4 were large even in the absence of fourth-degree coefficients for (,

if the

0, and

z, then the theory would be inadequate.

-

Invar.
Plane

of epoch

0 Eq. of date

Figure 3-

Precession Angles Using the Invariable Plane

Equator of epoch

Equator of date

Analytic Formulas for I

L, and A

Equating the two expressions for the precession matrix above gives

R2(9) R3(

-A) Ri(Io) R3(Lo).

(3-5)

The desired angles may be isolated on the right-hand side of equation (3-

5) by multiplying

both sides on the right by R3(

Lo) Ri(-

lo), producing

-z) R2(0) Ra(

- Lo) RI(

-I) Ra(

-A).

(3-6)

One can expand both sides of equation (3-6), obtaining the matrix elements sin A sin I

cos A sin I, sin L sin I

, and cos L sin I in terms of (, 0,

and Lo.

The angles themselves

The algebra is simplified considerably, however, if one first multiplies both

sides of equation (3-6) on the left by Ra(

This rotation combines with R3(-

L) on the

right-hand side to yield

R2(9) Ra3

(Lo + ()

-I) Ra(

-A).

(3-7)

Next the two sides are expanded.

The left-hand side of equation (3-

7) becomes

R2(9) R3

(Lo + )

K cos 0 cos(Lo + ()

- cos

0 sin(Lo +

() cos o

- sin 0 sin lo0

cos 0 sin(Lo +

) sin Io

sin(Lo + C)

cos(Lo +

cos Io

- cos(Lo +

) sin Io

(3-8)

sin 0 cos(Lo + C)

- sin 0 sin(Lo + ()
+ cos 0 sin Io

cos lo

sin 0 sin(Lo + ) sin Io
+ cos 0 cos Io

the right-hand side is

- sin 0 cos lo

= R3[

= R3(

L) RI(

-I) Ra(

= R3(

L) RI(

I)Ra(-A)

cos(L z) cos A
- sin(L z) cos I sin A

cos(L
- sin(L -

- z) sin A
z) cos I cos A

sin(L

- z)sinm

si
+ cos

S- z)cos A
- z) cos sin A

sin(L
+ cos(L -

- z)sin A

z) cos I

cos A

- cos(L

- z) sin

S(3-9)

sin I sin A

sin I cos A

cos I

The angle I is then found from equating the (3,3) components:

= cos-' [cos 0 cos Io + sin 0 sin(Lo + () sin lo].

(3-10)

Then

, assuming that I

0 so that the factors sin I can be cancelled

, the (1,3) and (

components yield L by

= pig [cos 0 sin(Lo + () sin Io

- sin 0 cos I0o

cos(Lo + () sin I1

and the (3,1) and (3,

) components give A by

= plg [sin 0 cos(Lo + (),

cos 0 sin Io

- sin 0 sin(Lo + () cos Io

(3-12)

where plg(y,

x) is defined

in Section

1-3 to be

There is

no sign ambiguity in either case because I

is restricted to the first or second quadrants;

therefore sin I

> 0 always.

(The extreme case of I

= 0 never occurs in practice, as the

Earth's

through

equator is always inclined at least 19 to the invariable plane.) Equ

(3-12) are therefore rigorously correct for all possible values of (,

nations (3-10)

and all

physically reasonable values for To and Lo.

-z)]Ri(

(, 0, and z, guaranteeing that A < 0. Similarly, for small posil

(corresponding to T > 0), both arguments are positive, and A

Equation (3-11) for L reduces to first order to L0 + C

the time origin. Finally, if T

tive values of (, 0, and z

will be positive as well.

therefore L is continuous near

= 0, one recovers I = Io, L = L0, and A = 0.

Series Expansions for I

L, and A

Given that the angles I, L, and A are found by equations (3-10) through (3-12) above,

their behavior with time depends on the behavior of the "classical" angles (, 0, and

z that

appear on the right-hand sides. (The angles Io and L0 are constant.) If the classical angles

are modeled as polynomials in time, then the new angles can also be so expressed.

Let the classical precession angles be approximated by the polynomials

=E

where T

(kT

9kTk

z=E

ZkTk

denotes time in Julian centuries from the standard J2000.0 epoch.

(3-13)

These must be

substituted into equations (3-10) through (3-12) and the trigonometric functions approx-

imated by polynomials.

For the sine,

sin 0 = 0

+0(0')

= (OT

S+02T2

+03T3

+ 4T4) (1

+ 02T2)3

+ O(T5)

= O1T+02T2

+ (93 -

+(04 02T4

+ O(Ts5).

(3-14)

Similarly for the cosine,

Using these two expansions, the expansion for the functions of (Lo + ) follows:

sin(Lo + () = sin

cos Lo + cos ( sin Lo

= sinLo+(1CicosLo)T+(

cosLo

+ [(C3

- 6f ) cos Lo

sin Lo]T

+ [(C4 C2)cOS

- (C

14 ) sin Lo]T4

+ O(T5);

(3-16)

cos(Lo + () = cos C cos Lo sin ( sin Lo

= cos

- (Ci

sin Lo)T

- ( f2 cosLo +

sin Lo)T2

cos Lo + (C3

1- )sin Lo]T3
-- 6% 8111)

- 24C)cosLo+ (C4(4- 21

) sin Lo]T4

+ 0(T5).

(3-17)

The last series expansion that will be used more than once is for the arctangent func-

tion. Let yo + Ay = tan

-(xo + Ax). A Taylor expansion of the right-hand side gives

Ay=

tan-1 X0 + Ax

dtan

S(Ax)
+ 2!

d2 tan

(Ax)3
+ 3!

d3 tan

S(Ax)4
+ 4!

d4 tan

+ 0 [(Aax)].

(3-18)

Since xo = tan Yo, the first four derivatives can be expressed as

dtan1 x
dx

d2tan-1 x
dx2

d3 tan-1 x
dx3

d4 tan-1 x

1 + tan2 Yo

-2tanYo

(1 + tan2 yo)2

6 tan2 Yo

(1 + tan2 yo)3

24(tan yo

co2
= COS

(3-19)

(3-20)

= -2 sinyocos3 yo;

= 6 sin2 Yo0 COS4 Y0 2 cos6 yo;

-tan3 Yo)

--d4 /t i 9 \

(3-21)

= 24(sin yo cos7 Yo sin3 yn cos5 /n ).

(3-22)

- f jsinLo)T2

- ClC2

22 -

- [(Cila + 1

= x1T

+x2T2

+ x3T3

+ x4T4

+ 0(T).

(3-23)

The next three powers of Ax,

complete to T4

are then

(Ax)

= xrT2

+ 2x1x2T3

ix3)T4

+O(T);

(Ax)3

(Ax)4

= xiT3

+3x2

= ixT4

+O(T);

(3-25)

+ O(T).

(3-26)

When these expressions for (Ax

Taylor expansion,

k and the values of the derivatives are inserted into the

the resulting equation is of the form

= yIT+ y2T2

+ y3T3

+ y4 T4

+O(T5)

(3-27)

with the coefficients Yk given by

= x2

-~2 sin Y cos3 ;
-xi1SInDoc0 CO yo;

(3-29)

= -3

cos

-2x

lX2 sin yo cos3 yo + x3(sin2 Y0o cos4 o i cos6 Yo);

(3-30)

= X4 COS

x4(sin yo cos7

-xlx3) sin Yo co
- sin3 o cos5
-sin yo cos 1

Y0 + 3xx2(sin2 yo cos4 Y0 cos6 Y0)

(3-31)

be developed below

3.3.1

but because they will be used but once, they will

as the need arises.

. The Expansion for I

The inclination I of the invariable plane to the equator of date is given by equation
/'*5_1 (\ ^k ^rn

= a1

x2T4

-(Z+

55
First the expressions for cos 0, sin 0, and sin(Lo + () must be substituted into the argument

for the arccosine:

I = cos-1 ([1 a 2T2

- 0102T3 (003 + 2*,2

- 4T4)] cosio

+ [01T + 02T2 + (03-61 T3 + (4 1 2)T4]
x (sin Lo + (Ci( cos Lo)T + ((2 cos Lo (2 sin Lo)T2

+[(C3
+ [((4

- 6() cosLo (1(iC2 sin Lo]T3
- 2f(2)cosLo ((iC3 + _J )sinLoT}si In o

+O0(TS))

(3-32)

= cos

-(1 -(l !022

- 02T3 (0103 + 22 4 )T cosIo

+ {(01 sin Lo)T + (02 sin Lo + 0i(i cosLo)T2
+ [(03 1 11C2) sin Lo + (01(2 + 02(1)cosLo]T3

+ [(04 -2 12
+ (091( 3+02

12 11i(2) sin Lo

+" 3(1 103( 1iC13) cosLojT4 sin Io

+ O(Ts5))

= cos

(3-33)

-1 { coslo + (01sin Lo sin Io)T

+ (02 sin Lo sin lo + 01(1 cos Lo sin Jo 02 cos Io)T2
+ [(3 ? -- 2 C) sin Lo sin Io + (01 (2 2C1) cos Lo sin Io
O102 cos lo]T3
+ [(04 }02 02(2 01(12)sin Losin Io
+ (9C3 +92 2 +3(C 1,a 01,i() cosLosinlo

+ (204 0103 22 )cos0o]T4 + O(T5) }

(3-34)

=- cos

-1 [cosIo + CT + c2T2 +c3T3 + cT4 + O(T5)],

(3-35)

where the Ck in equation (3-35) are defined by the various coefficients of Tk in equation

(3-34).

Now equation (3-35) has the form

c= ciT+c2T2

+ c3T3

+ cT4

+ 0(T5).

(3-37)

We therefore perform a Taylor expansion of the arccosine function, as was done above for

the arctangent:

I = cos

-1 Xo + c

dcos

c2
+2!

c3
+3!

d3 cos

c4
+ .

d4 cos

+ 0(c5).

(3-38)

Since xo = cos lo, the leading term is simply Io.

The first four derivatives become

d cos- x
dx

d2 cos-1 x
dx2

d3cos-1x
dx3

d4 cos-1 x

(1 -2)1/2

X
x
(1 X2)3/2

1+2a:2
(1 2)5/2
(1 x)/

9x + 6x3

- x2)

1
sinT0o

cos I0
sin3 To

1 + 2 cos2 0

(3-39)

(3-40)

(3-41)

sin5 o

9 cos Io + 6 cos3 Io

(3-42)

sin7 Io

When one substitutes into equation (3-38) the expression in equation (3-37) for c, along

the values

of the derivatives in

equations (3-39) through

(3-42), one obtains an

equation giving the coefficients Ik of the various powers of T in terms of the ck:

I= Io (ciT

+ c2T2

+ c3T3

smTo)

- j(cT+c c2T2

+c3T3

+ c4T ) s3 0
Vsm I'o/

d2 cos

=Io- CT
sm Io

c2
- (
\sin Io

C COS 10
2sin3 I0 /

c C
- sn +
\sin Io

CiC2 cos To

c?(1 +

2 cos2

[o) T

sin3 1o

(c) +

2clC3)cos Io

2 sin3 Io

c (9 cos o + 6 cos3 Jo0) 4
24 sin7 1o

+ 0(T)

(3-44)

=-Io+I T+I2T2

+ 3T3

+ 4T4

+0(Ts5).

(3-45)

This last equation also shows, as expected

that I --* Io as T

-+0.

The last step is of course to obtain the coefficients Ik in terms of the coefficients Ok and

in the approximation polynomials for 0 and

by substituting the various coefficients in

equation (3-34) into equation (3-44). After some tedious algebra, one obtains

Cl
sin Io

= 01 sin Lo;

c2
sin Io

O? sin Lo sin Io

sin Io

(3-46)

c cos I0
2 sin3 1o

02 sin Lo sin Io + O i1 cos Lo sin Io 12 cos Io
_ _ _2 1 c sI

sin Io
(0l sin Lo sin Jo)2 cos Jo
2 sin3 Io

= 02 sin Lo 01(i cos Lo + j.of

Lo cot Io;

(3-47)

C3sin
sin Io

C C2 COS l0

cl(1 +

sin3 Jo

cos2 Io)

sin5 1o

(03 6 I1iC2) sinLo sin lo + (01G2 + 21) cos Lo sin 0o 012 cos lo
T

58
= (-3 + 0 + 0(3) sin Lo (01(2 + 2i) cos Lo
+ cot lo(0102 cos2 Lo 02C1 sin Lo cos Lo)
+ cot2 o[3(1 sin Lo sin3 Lo)] csc2 o(0 sin3 Lo)
= (-03 + 011i2) sin Lo (01(2 + 02(1)cosLo + 9 cOs2 Lo sin Lo
+ cot lo(0102 cos2 Lo 0Ci(1 sin Lo cos Lo)
+ cot2 lo[f( sin Lo cos2 Lo)]; (3-48)

Sc4_ (cj +2cic3)coslo
TA -- -- --
sin -o 2 sin3 o
c c2(1 + 2 cos2 Io) c4(9 cosIo+ 6 cos3 Io)
2 sin5 1o 24 sin7 Io
= [(4 212 22 1 1 2) sin osino
+ (01(3 + 02(2 + 03(1 C1 1i13) cos Lo sinlo
+ (gOf 0j03 420) cos o / sin Jo
[(02 sin Lo sin o + 0iCi (cos Lo sin 1o j09 cos o)2
+ 2(0i sin Lo sin Io)[(93 lo3 ^ 2) sin Lo sin Io
+ (02(1 + 01i2)cos Lo sin Io 0192 cos Jo]] coso/(2 sin3 o)
[(01 sin Lo sin lo)2 (02 sin Lo sin Io + 01 (1 cos Lo sin Io cos lo)
x (1 + 2 cos2 o)]/(2 sins Io)
[(0i sinLosin o)4(9cos o10 + 6cos3 Io)]/(24sin7 Io)
= [- 04 + 01(1(2 + 02(01 + (12)] sin Lo
+[-01(3-02(2- 03(1 11+ 6l(i 1 + C2)] cosLo
+ cot o [ + ( + (2) sin Lo + (0103 2t0 12) cos2 Lo
cotlto[ +~i ( +9 9 lsnL (1322111

(20102(1 + 022)sin Locos Lo]
+ cot2" Io [01202 (2 sin Lo sin3 Lo) + 03 (1( cos Lo sin2 Lo cos Lo)]
+ cot3 Io [0(9f(- + sin2 Lo 1 sin4 Lo)]
+ csc2 o [ 101 sin2 Lo(02 sin Lo + 01(1 cos Lo)]
+ cot Io csc2 Io jje sin2 o(2 3 sin2 Lo)
= (-94 1(2 1C2) sin Lo + j-902 sin Lo cos2 Lo 91i sin2 Lo cos Lo
I F .- It 11 n r 1 n .,- /^ n\2 7-

+ cot3 I0 [(-

S+ 3 sin2 Lo -

sin4 Lo)

(3-49)

Equations (3-46) through (3-49) are the desired results for the coefficients of I

3.3.2.

The Expansion for L

The right ascension of the ascending node of the invariable plane on the equator of

date, denoted by

L, is given by equation (3-11):

L = plg [cos 0 sin(Lo + () sin Io sin 0 cos Jo, cos(Lo + () sin 1o] +

(3-11)

= tan

_-1 cos cos(

Lo

) sin I0 sin 0 cos Io\
+ () sin o

In order to derive the coefficients Lk such that

L= Lo + L1T

+L2T2

+L3T3

+ L4T4

+O(T5),

(3-50)

one must develop both the numerator and denominator of the argument of the arctangent

as approximation polynomials, obtain their quotient, and expand the arctangent itself.

Denote the numerator of the argument of the arctangent in equation (3-11) by n and

the denominator by d.

Then inserting the results of equations (3-14) through (3-17) gives

the following:

n = cos 0 sin(Lo + () sin Io sin 0 cos Io

- 0102 T3

- (0103 + i922

x {sin Lo + ( cos Lo)T

cos Lo

- 141 sin Lo)T

+[((3
+ [(14

- R?) cos Lo

- (1(2 sin Lo]T

cosL

241) sin Lo]T } sin 1o

l- 2T2
~- 2VllA

- N 4)1

- 1i)

- (1(3 +

-(03

- ia) coslo]T3

+ {[4

- 0102l 1

- )] + 0 (2] cosLo sin o

+ [-0103

S12 +.4" (0 +-4)]sinLosinmo

-(04

- 2 1 2f) COS to} 4

+O(T5)

(3-51)

d = cos(Lo + C) sin lo

= { cos Lo

- ((~1 sin Lo)T

- (lcos

Lo+

sin Lo)T2

- [(1(2 cos Lo + (Ca

j fil) sin Lo]T3

- [(ia+3

(14) cosLo + (44

- ((2) sin Lo]T4} sinlo + 0(T5)

= cos Lo sin Io ((1 sin Lo sin lo )T

- (1('2cosLo +

sin Lo) sin loT2

cosLo + ((3

- 1() sin Lo] sin IoT3

--[((1(Ca

-24 4lcos

) sin Lo]sin IoT4

+0(T5).

(3-52)

If the leading terms in both n and d are factored out, the argument of the arctangent in

equation (3-11) takes the form

sinLo sin Iofl1 + nT

cos Lo sin Io[1

-diT

+n2T2
- d2 T2

+ n3T3
-a 3T
- da T

+n4T4
- d4T4

+ O(T)]
+ O(T5)]

1 + n1T + n2T2
= tan Lo -dT -d2T2
V1 -diT -d2/

n3T3
-d3T3

+n4T4
-d4T4

(3-53)

since sin Io \$ 0.

(The minus signs in the denominator cancel some of the minus signs in

equation (3-52) and also facilitate expansion of the quotient later.)

Now the coefficients

nk and dk are given by

= (1 cot Lo 9 cot Io csc Lo;

3-54)

n2 = (2 cot Lo

- 02 cot Io csc Lo;

(3-55)

- 1/"3
641

- -12C1) cotLo ((1(2 + 102)

- (03 l) cot Io CSC Lo;

(3-56)

-- 2 ,2 +

-2 l2

Lo + ((4

+ O(T5) y
+ O(TS) '

- (02l 12)

d = C(1 tan Lo;

(3-58)

tan Lo + 21;

(3-59)

A13)tanLo+ (2;

(3-60)

(3-61)

The next step is to expand

the quotient that forms the argument of the arctangent

function in equation (3-11). Denote the quotient by q; one obtains by simple long division

1+nlT+n2T2
q 1 diT d2T2

+ n3T3
- d3T3

+ n4T4
- d4T4

+O(T5)
+ O(TP)

(3-62)

+ d1)T +(n2

+d2 +dinl +d#)T2

+ (rn3 + d3 + 2did2

+ d2n1

+ din2

+ dinl + d3)T3

+ (n4 + d4 + 2did3 + d2n2

+ d + 2did2n1 + 3d(d2

+ din3 + d1n2 + dIn + d4)T4

+ O(T5)

(3-63)

= 1 + qi T + q2T2

+ q3T3

+ q4T4

+ O(Ts5).

(3-64)

This leaves equation (3-11) in the form

-l(qtan Lo) +

(3-65)

= tan-1 { tan Lo + tan Lo[qiT + q2T2

+ q3T3

+ q4T4

+O(T5)]}

+ziT+z2T2

+ z3T3

+ z4T4

+ O(T5).

(3-66)

The expansion for the arctangent was given above in equations (3-27) through (3-31).

Here L0 plays the r6le of x0, and qk tanLo are to be substituted for Xk.

When

the expansion becomes:

L1 = (qitan Lo)

Lo + z1

(3-67)

L = tan

d2 = (2

d3 =(

d4 = ((4 1_2

)tan Lo + (1 3 +K 2C2 14l

= 1 + (nl

,

b *

L4 = (q4 tan Lo) cos2 Lo [(q2 + 2qlq3) tan2 Lo] sin Lo cos3 Lo

+ 3q1q2 tan3 Lo(sin2 Lo cos4 Lo

-- 3 COS6

+ (qi tan Lo)4(sin Lo cos

- sin3 Lo cos

Lo) + z4.

(3-70)

It is easiest to leave the tan Lo with the qk.

Now replacing nk and dk with their values from equations (3-54) through (3-61) gives

ql tan Lo = (i + di) tan Lo

= [Ci(1 cot L0 01 cot Io csc Lo + Ci(1 tan Lo] tan Lo

= (C 1 cot Io secLo + C(1 tan

= (1 sec Lo 01 cot Io sec Lo;

q2 tanLo = (n2 + d2 + din1 + df)tanLo

(3-71)

= {[C2cot Lo

- ( + C(1)

- 02 cot Io csc

Lo] + ((2 tan Lo + 4l2)

+ ((1 cot Lo O8 cot Io csc Lo)((1 tan Lo) + ((1 tan Lo)2} tan Lo

= [C2

- 2(~1 + C)tanLo

- 02 cot Io

sec Lo

tan2 Lo + 2C? tan Lo)

+ ((1 tan Lo

- 0 (1 cot lo sin Lo tan Lo) + (2 tan" Lo

02 cot Io

sec Lo + (12 tan Lo

- 909 tan Lo0

- 01(1 cot Io tan Lo

sec Lo;

(3-72)

q3 tanLo = [In3 + n2d1 + ni(dt + d2) + d3 d3 + 2did2]tanLo

- 2IC(1)cotLo

+ [(2 cot Lo

- (w +
!(2
-- 21,

+0102)

2 cot Io

- 1) cot o cscLo
tan Lo)

csc Lo((1

+ ((1 cot Lo 01 cot Io csc Lo)[((1 tan Lo)2

+ [(Ca

- 6(3) tanLo + (12] + (Ci1 tan Lo)3

+ 2((1 tan Lo)(C2 tan Lo + zi2)} tan Lo

= [(3

- !-1
6 %1

- 2 ,11

- (( (2 + 192)tanLo

- (03

- 6
~6Gi)

cot Io

sec Lo]

+ [(1(2 tanLo

- i(0? + ( tan

- 02 ( cot In tan Lo sec

tan Lo + (1)]

Lo-

1-3l
_ !S1

-(o03

s--

63

= ((3 + 3 1 l f1) sec2 Lo + 2(1(2 tan Lo sec2 Lo 0102 tan Lo
+ (3 tan2 Lo sec2 Lo
1 3
+cot Io sec Lo[-03 + 1 ? 112 (2(1 + 1i(2) tan Lo
01f tan2 Lo]; (3-73)
q4 tanLo = (n4 + d4 + 2did3 + d2n2 + d2 + 2dld2n1 + 3d d2
+ din3 + d n2 +dni + dl)tanLo

= ({[4 2(2(1 (2) 01021] cot Lo
[9103 + (1C3 + C82 2-2) 1*2 *(0 + l4)]
-(04 -102 ) cot Io cscLo }
+ [((4 1' (2) tan Lo + (1(3 + 22 -- 41
+ 2((C tan Lo)[(C3 k3 tanLo + (1(2)
+ ((2 tan Lo + )1)[1 cot Lo 2 + (12) 02 cot Io csc Lo]
+((2tanLo+ 2 )2
+ 2(C1 tan Lo)((2 tan Lo + 1f)((1 cotLo 1 cot o sc Lo)
1/-
+ 3((1 tanLo)((2 tanLo + 2(1)
+ (1Ci tan Lo)[(C3 63( 10 (C) cot Lo
((1(2 + 0102)- (03- 13) cot o cscLo]
+ (C tan Lo)-(2 cot Lo -0 + (12) 02 cot Io csc Lo]
+ ((1 tan Lo)3((1 cotL0 01 cot Io cscLo) + ((C tan Lo)4) tanLo

= [(4 2(0 2) -0 021
-t019,(1l3+ 2121)-- 002-
[0103 + 3 + + 2) + (f)] tanLo
(04 0 02) cot Iosec Lo]
2+ [(14 22)tan2 Lo + ((1(3 + (22 f tanLo]
+(2(1iC3 tan3 Lo j-f tan3 Lo+2C(2 tan2 Lo)
+ [ tan Lo (2(10 + C12) tan2 Lo 02(2 cot o tanLosecLo
+ (12 t1C2)tanLo- 0212cotlo secLol
+ ((22 tan3 Lo + (122 tan2 Lo + 4?1 tanLo)
+ (2(1(2 tan2 Lo + (4 201 2 cot 10 tan2 Lo sec Lo 01i3 cot Jo sec Lo)
i /o/'2i +.__.4 r i 3/-4 ..... r \

+ ((1f tan3 Lo 01( cot Io tan3 Lo sec Lo) + (1 tans Lo
C

= ((4 + 1

--01021- 1422)sec2 Lo

+ 2(1( 3

- 012 12)tanLo

tan2 Lo

Lo + C14 tan3 Lo sec2 Lo

- ( 2 + 0103 14) tan Lo

sec Lo[(0120

2-- 4
2 lC3

+(1

- 212C 1CC2)
-02(2 -93a1)tanLo

- (02 ( + 201(1(2)tan2 Lo

- 01iC3 tan3 Lo].

(3-74)

The final step is

to substitute these values for the qk tan Lo into equations (3-67)

through (3-70). For L1, we obtain

Li = (qitanLo)

Lo + zi

(3-67)

= (C1

Lo 091 cot Io sec Lo) cos Lo + zi

= C1 + zl 01 cot Io COs Lo.

(3-75)

Next, we get for L2:

L2 = (q2 tanLo)cos2 Lo (qi tan Lo)2

sin Lo cos3 Lo + z2

(3-68)

= C2

- 02 cot Io sec Lo + (2 tan Lo sec2 Lo

- 2 1 tan Lo

- 011C cot Io tan Lo sec

Lo] cos2 Lo

- (C

01 cot Io sec Lo)2 sin Lo cos3 Lo +

- 02 cot Io cos Lo + (2 tan Lo

- 102 sin Lo cos Lo 11 cot Io sin Lo
--21l

- (C? tan Lo

z2--

- 201(1 cot Io sin Lo + 01 cot2 1o sin Lo cos Lo)

02 cot Io cos Lo + 01 (1 cot Io sin Lo

- 0f sin L0 cos Lo(cot2

I0+o).

(3-76)

+( f 4
+ 3(12(
+ cot Io

Lo-

- 0- 6

65
The right-hand side of equation (3-69) for L3 has three terms involving qk, plus z3.

Let us include z3 in the first term (which will contain (3); the three terms become, in order:

L3a = (q3 tanLo) cos2 Lo +

- 092l )

Lo +2((42 tanLo

Lo 6102tanLo

+ (Cf tan2 Lo

+ cot losecLo[-a03 +

-01(2 tan2 Lo]}

-(92C1 +1C2)tan Lo

Lo+z3

- 02 jl + 2(1C

- 011i2)

tan Lo
cosLo

- 019 2 sin Lo cos Lo + (Cf tan2 Lo

- (02( + 012) sinLo

- ~iC( sin Lo tan Lo];

2(qiq2 tan2 Lo) sinLo cos3 Lo

((1 sec2 Lo 0 cot lo secLo)

(3-77)

x (C2

Lo 02 COt Jo

sec Lo + 2 tan Lo sec Lo 92 tan Lo
--2 1tnL

0(1i cot JIo tanLo
x sin Lo cos3 Lo

sec Lo)

- 2C

tan Lo

- 2(C tan2 Lo

- 0(1 sin2 Lo

+ cot o [2(01(2 + 02(1)sin Lo + 40i1C sin Lo tan Lo 0 sin2 Lo cos Lo]

+ cot

1o(-20102 sin Lo cos Lo

- 20121C sin2 Lo);

(3-78)

(qi tan Lo)3(sin

Lo cOS4 Lo -

1 cos6 Lo)

= ((C sec2 Lo 0i cot Io sec Lo)3(sin2 Lo cos4 Lo cos6 Lo)

= ( tan2 Lo

+ cot

-1+cot
--341 -cotlI0(

2 o(3021 sin2

-301(i'C sin Lo tan Lo + 30i12 cos Lo)

313
Scot3 Io(-0c si2 Lo cos Lo + isf cos3 Lo).

The desired coefficient L3 is the sum of these three:

(3-79)

- 11

= (C+ R

z3 + 3 (l

+ cot lo[(-03 + 1-

L3b =

_ 12 2

+ cot

2 lo0 124(sin2 Lo

- cos Lo)

- 20102 sin Lo cos Lo

+ cot3 Io[(03( cos3 Lo

- sin Lo cos

(3-80)

Finally, equation (3-70) for L4 contains four terms on its right-hand side:

L4a = (q4 tan Lo) cos2 Lo + z4

= {((4 + ( 2

- O 02 41

+ 2((3

tan2 Lo

l12
-- 21l

-- (fl)tanLo
Lo + (4 tan3 Lo

- (? 02 + 003 )tanLo

+(1Of4

2-^4
2- 04l
- 02(

- 0?v(s
501 r3
-- 1 1~

- 02(2 01(1(2)
- 02(2 03(1)tanLo

- (02(2 +

1 1(2) tan2 Lo

- 1(1 tan3 Lo]}

Lo + z4

= C4

+ z4 + 1C(2

- 9102 (1

1 2
-- 2 l

- (22 + 0103

- 10f) sin Lo

+ 14
--(3,l +

- 2)9 tanLo + 3(2(2 tan

Lo + (f tan3 Lo

+ cot Io[((01202

-01 1 2)cos Lo

+(0Ci

- 01(3

- 0242

- 031)sinLo

01(1(2) sin Lo tan Lo 01(3 sin Lo tan2 Lo];

(3-81)

- [(22 + 2qiq3) tan2 Lo] sin Lo cos3 Lo

- 02 COt Io

sec Lo + (12 tan Lo

L 20 tan Lo

- 01 (1 cot I0o tan Lo sec Lo

- 01 cot Io

x {((a + 4

sec Lo)

- 1i ( l) sec2 Lo

+ 2(1(2 tan Lo

+ cot lo secLo[-03 + j0

Lo 0102 tan Lo + (3 tan2 Lo sec Lo

- 212
- 2tl,

- (02( + 01(2)tanLo

tan Lo + 022 cot2 Io sin Lo

- 01(2 tan2 Lo

cos Lo + (4 tan3

}) sin Lo cos3 Lo

Lo + 1 sin

+ 2 cot2
01(1 cot2

Io sin2 Lo tan Lo

- 202 (2 cot lo sin Lo + 2(t2

(2 tan2 Lo

- 0?2 sin2 Lo

+ ot I

+3(,(2
+ cot Io

sec Lo[( O1 8

cos Lo

L4b =

sec Lo

Lo cos Lo

Lo)].

- 02121

5- /-3
-- 6i1i

- (02Cl2 +

= (

+ C( tan3 Lo

+cotlo[(-0a3l + tCi

- (0?13 -C- 1 2)sin Lo

- (011 +j01(12)sinLotanLo -

01( sinLo tan2 Lo

- cot Io[(01(3 + 031

- 0 02 sin2

Lo cosLo + 01(?13 tan

sin Lo + 201 (1

sin Lo tan Lo

Lo sin Lo

- cot2 o[(-0103 + I

- -Of 2) sin Lo cos Lo

- ( i02(1 + 2)sin2

- Of sin2

Lo tan Lo] }

+ 201 02 1i) sin

- 14- sin3 L0 cos Lo

-r4 9(2) tanLo 9f2Sin2
--3"1 + lS / a ^ + 01l 1 sln2

Lo tan Lo 6(C tan2 Lo

- 3c4 tan3 Lo

+ cot Io sin Lo[201(3 + 202 2

+ 2039

- 391(1 + O 31C

- 30f02 sin Lo cos Lo

- O1i1 sin2 Lo + (801i(1C2 + 402 (f) tan Lo + 601i13 tan2 Lo

+ cot2 Io sin Lo[(-O2 + 34l

- (40102C1 + 202

- 20103

- 0 ff1) cos Lo

)sinLo 30 (2 sin Lotan Lo];

L4c = 3q1q2 tan3 Lo(sin2 Lo cos4 Lo cos6 Lo)

= 3(1Ci sec2 Lo 01 cot Io sec Lo)2

(3-82)

Lo 02 cot Io

01i1 cot Io tan Lo
X (sin2 Lo cos4 Lo 1

sec Lo + C2 tan Lo sec2 Lo tan L

sec Lo)

cos6 Lo)

01 cot Io cos Lo)

- 02 COt I0 COS

Lo + ( tan Lo j0 sin Lo cos Lo 01(1 cot Io sin Lo)

x (3 tan2 Lo 1)

- iO22 sin Lo cos Lo + (1 tan Lo
--2 isioo1ian

+ cot Io(-02(C cos Lo

- 01 ?13 sin Lo

- 201(1(2 cos Lo + 03 sin Lo

- 20163 sin Lo)

+ cot

Io(20i02 C

Lo +

0~2Q sin Lo co

s2 Lo 201(3 Sin Lo)

+cot3 Io(-00 2 COS3 L0o

- 013i sin Lo co

s2 Lo)](3tan2 Lo -1)

A

AL i 9,. r T* r ifl/ A-/ I r r

- (0()

=(

+ (-

- 2ClC3

x ((2

= ((1 -

x ((2

= [C1

^- -< -fc

+ cot2 Io[(-20102(1

- O(C2)cos2 Lo+ i- sin Lo cOS3 Lo 30(2C sin Lo cosLo

+ (601021

+ 3022) sin2 Lo

- 1f sin3 Lo cos

Lo + 909fC sin2

Lo tan Lo]

+ cot3 Io[0(02 cos3 Lo + 06(i sin L0o cos2 Lo

- 302 sin Lo cos Lo

- 303C sin3 Lo];

(3-83)

L4d = (qi tan Lo)4(sin Lo cos7

Lo sin3 Lo cos5 Lo)

Lo 01 cot Io sec Lo)4 (sin Lo cos7 Lo sin3 Lo cos5 Lo)

= ((Ci 1 cot Io cos Lo) (tanLo tan3 Lo)

= (1 (tan Lo tan3 Lo)

+ cot Io[401(13(sin Lo tan2 Lo
+ cot2 Io[60 (f(sin Lo cos Lo

+ cot3 Io[40(3 (sin3 Lo

- sin Lo)]
- sin2 Lo tan Lo)]

- sin Lo cos2 Lo)]

(3-84)

+ cot4 Io[01 (sin Locos3 Lo sin3 Lo cos Lo)].

Their sum gives L4:

-12 I2
- -~l

- 0102 (1

1 2
- 2

- sin2 Lo) + 0~

- 0l03+ 9f) sin Lo cos Lo

sin2 Lo

+ cot Io{[-04 + 01(1(2 + 0212C + 222 (

1 sin3 Lo cos Lo
- 3 sin2 Lo)] cos Lo

+ [(a + 02(2 + 0a3l

-&(OC + 0C)+0C9(sin2

Lo-2

Lo)] sin Lo}

+ cot2 Io[(20C12

- 02 +

- 201 03) sin Lo cos Lo

- (2012 + 12C2)(cos2 Lo sin2 Lo)

+2a sin Lo cos Lo(cos2 Lo
+ cot3 o10[002 cos Lo(cos2 Lo

+ cot4 2o[1 sin Lo cos Lo(cos2 Lo

- 3 sin2 Lo)]

- 3 sin2 Lo)

- 0(i sin Lo(3 cos2 Lo sin2 Lo)]

- sin2 Lo)].

(3-85)

Equations (3-75), (3-76),

(3-80), and (3-85) therefore contain the final results for the first

four coefficients of the approximation polynomial for L. No attempt has been made here to
J J A i'r 1 1 "/ ) r 9 1 \

= ((i

12 2
-L- (f1 1

L4 = (4 + z4

69

relatively time-consuming on a computer, and there will be no significant cancellation in

the subtraction.

3.3.3.

The Expansion for A

The angle A is measured along the invariable plane; its origin is the ascending node

of the invariable plane on the mean equator of J2000.0, and its endpoint is the ascending

node of the invariable plane on the mean equator of date.

The exact expression for A was

found at the beginning of this section:

A = plg [cos(Lo + () sin 0, sin Io cos 0 cos Io sin 0 sin(Lo + ()] .

(3-12)

We will expand the right-hand side of this equation in a manner similar to the expansion

of L in the previous subsection, resulting in an expression of the form

A = A1T + A2T2

+Aa3T3

+ A4T4

+ O(T5).

(3-86)

As in the previous subsection, denote the two arguments of the plg function in equation

(3-12) by n and d respectively, as these are the numerator and denominator of the argument

of the implied arctangent.

(Although the notation is the same, the values of the symbols

are obviously quite different in this subsection.)

Substituting the expansions for the sine

and cosine of 0 and (Lo + () from equations (3-14) through (3-17) gives:

n = cos(Lo + () sin 0

= [cosLo

- K1

-(C sin Lo)T
(2cos Lo + ((3

sin Lo)T2

- 3C?)sin Lo]T3

- [((1iG+ -

- C-)fcosLo +(4

'

) sin Lo]T4

O(T5)]

-(J^zcos

Lo+

-1_

^ Jj

,

- O0 l 0141

) cos Lo

- (0I(3 + 02(2 + 03(1

1a3
--6 1l<1~

- -01(43) sin Lo]T4

+ O(T5);

(3-87)

d = sin Io cos 0 cos Io sin 0 sin(Lo +

= sin Io [1

- 0102 T3

_ 11)T4

-(0193 + 2I2

+ O(T5)]

1_

- cosIo[0iT + 02T2

+(o4

- 12)T4

+O(Ts5)]

x { sinLo + ((i

cosLo)T+ ((2cosLo 2( sinLo)T2

+[(C3
+ [((4

- i COs Lo

- (1(2 sin Lo]T3

) cos Lo

- ((Cl3 +

-24 sin LojT4

+O(T5)}

= sin Io

- (01 sinLocos Io)T

(02 sin Lo cos o + 01(1 cos Lo cos To + j sin Io)T

- [(03
- [(04

6 1
- 2!
l2 "I

+(01(3+02

- -i0(j)sinLocoslo + (O1

- o0(1
2 + 3C l

+ 02 (1) cosLo cos o + 1i 02 sin Io]T3

- l ,1 2)sin Lo cos Io

- 01(1
l03/-
-- lql

cosLo sin lo

+ (0103 + 2!2

- 01)sinIo]T4

+O(T5).

(3-88)

The absence of a constant term in the numerator is consistent with the statement made

earlier that

0OatT

= 0. The algebra is simplified, however, if we divide both numer-

ator and denominator by the the constant term in the denominator, namely sin Io; this is

permissible because (as stated before) o0

This leaves the argument of the arctangent

in equation (3-12) in the form

n1T+n2T2

-diT

- d2T2

+ n3T3

- d3T3

+ n4T4

+O(T5)

(3-89)

- d4T4 + O(T5)

with the coefficients nk and dk given by

nl = 01 cos Lo csc lo;

(3-90)

- -1 2

A =

- 2,1T

+(03o

-1^2
-- 2("

92

-- 1011 (1

- (lC(3 + 2C2 +3C1

S 0- -l1)sminLo]jcsclo;

(3-93)

(3-94)

d2 = 02 sin Lo cot Io + 09 (1 cos Lo cot Io + 1 9;

d3 = (03 Of3

0i() sin Lo cot Io + (01

(3-95)

(3-96)

91022
--2li2n

~ 2-1

+ (013 + 02(2 + 03(1

+(1+02
31-(9103+"1-2 2

The next step,

- 0i(1i2)sinLocotIo

10
--6l l

- 91 Ci) cos Lo cot Io

-lOf ).
~ 24 t/ )

(3-97)

as before, is to expand the quotient that forms the argument of the

arctangent function in (3-12).

Once again we shall denote the quotient by q.

By an even

simpler long division than before,

n1T+n2 T2

-diT

-d2T2

+ n3T3

+

-d3T3

n4T4 +0(T6)
-d4T4 +O(T)

(3-98)

= 1+niT+(n2+ddn1)T2

+ (n3 + din2 + d2n1 + dzni)T3

+ (n4 + din3 + d2n2

+ d3nl + d n2 + 2did2ni + dinl)T4

+O(T5)

(3-99)

- 1+qiT+q2T

2 + q3T3

+ q4T4

+ O(T5).

(3-100)

The expressions for the qk are given by:

9qi = n1

= 01 cosLo csc Io;

q2 = n2 + dilnt

= [(02 cosLo 01(i sin Lo) csc o] + (Oi sin Lo cot Io)(9i cos Lo csc Jo)

(3-101)

= [02 cos Lo 0~i1 sin Lo + cot lo(0i sin Lo cos Lo)] csc lo;

-- .i 1. ,7-Iy _L ,7-Lyi- r, .

(3-102)

di = O1 sin Lo cot lo

+ 02 ( 1) COS Lo cot Io + 102;

I)-

q=

d4 4 (04

A

+ (0i sin Lo cot lo)2 (01 cos Lo csc Io)

-_ 01(1) cos

- (01(2 + 02(i)sinLo

+ cot 1o [20192 sin Lo cos Lo + 9 (1( (
+ cot2 Io(9f sin2 Lo cos Lo) } csc Io;

- sin2 Lo)]

(3-103)

q4 = n4 + din3 + d2n2 + d3n1 + d n2 + 2dld2nl + d1ni

= [(4

- 12

- 02(2

3 +02

- 011iCi(2) cos Lo

+031

- -1

--61

l) sin Lo

+ (01 sin Lo cot Io){[(03

- 101(1) cos0o

+ (02 sin Lo cot Io + 01(i cos Lo cot Io + j0 )[(62

i0f3
- 61

- (91(2 + 02(1)sinLo

cos Lo

- 01 (1 sin Lo)

- 011'i) sin Lo cot Io + (01(2 + 02(1) cos Lo cot Io + 0102]

x (0i cosLo csc lo)

+ (0x sin Lo cot o)2 [(02 cos

- 011 sin Lo)csc Io]

+ 2(01 sin Lo cot lo )(02 sin Lo cot Io + 9 i( cos

x (01 cos Lo csc lo)
+ (01 sin Locot lo)3(09 cos Locsc lo)

= {((4 + 022

- (0143 + 0242 0331 + 3 C1o

- 1i() sin Lo

+ cot Io[(20103

-2t(2

+022 + 3 1)sinLocos

+ (20102(1 + (42)(cos2 Lo

- sin2 Lo)]

+ cot2 Io[3102 sin2 Lo cos Lo + 0f (1i sin Lo(2 cos2 Lo sin2 Lo)]

+ cot3 Io(Oe sin3 Lo cos Lo)}

csc Io.

(3-104)

This leaves equation (3-12) in the simple form

tan 1 q

(3-105)

= tan-1 [qiT + q2T2

+ q93T3

+ q94T4

+ O(T5)].

(3-106)

Once again we employ equations (3-27) through (3-31).

This time, however.

Xn is zero.

= {(o3 + 3

csc lo }
csclo}

L0ocot I0o + )

- 01(1(2)cosLo

- (81

- Of3

- (l2

A1 = ql

(3-107)

= 1 cos Lo csc lo;

(3-108)

A2 = q2

(3-109)

(3-110)

= [02 cos Lo 01(1 sin Lo + cot Io(9i sin Lo cos Lo)] csc lo;

(3-111)

= {(03 + 361

- 1i()cos

+ 02 1i) sin Lo

+ cot o [20102 sin Lo co

s Lo + C1(cos2

- sin2 Lo)]

+ cot2 Jo(93 sin2 Lo cos Lo)} csc Io
- -6" cos3 Lo csc3 Io

= (3

- i(l+- 3 1 COs2

Lo) cos Lo

+ 02 (1) sin Lo

+ cot 10o [20102 sin Lo cos Lo + OPCi(cos

-sin Lo)]

+ cot

2 Io[O3(sin2 Lo cos Lo

- cos Lo)]}

csc 1o

(3-112)

A4 = 94 9l 92

(3-113)

= { (04 + 12 02

- 91(1

) cos Lo

- (013 + 02C2 +C O3(1 + 31

+ cot lo[(20103

- 619i1) sin Lo

- 20112 ( 22 1) sinLocosLo

+(20102(1 122)(cos2 Lo

+ cot

- sin

2 Io[3002 sin2 Lo cos Lo + 913 sin Lo(2 cos2 Lo -sin2 Lo)]

+ cot3 Io(0O sin3 Lo cos Lo)} csclo
- (0(02 cos3 Lo 1Ca cos2 Lo sin Lo + 94 sin Lo cos3 Lo cot lo) csc3 Io

= {(04

- (01(3 + 022 3(1 31

- 1il) sin Lo + 01i cos2 Lo sin Lo

+ cot Io [(201i 03

- 201( + 022 +

0 ) sin Lo cos Lo

+ (2091021 122)(cos2 Lo

- sin2 Lo)

- 01 sin Lo cos3 Lo]

+ cot2 lo[02 02 cos

Lo(3s2Lo

- cOS2

+ 0(1 sin Lo(3cos2 Lo -

sin2 Lo)]

I _- f E I n [ / a P T- 1 .

+ 0202 sin2 Lo) cos Lo

A3 = q3 -

- (01

- (01

- 2 12

- o2(42

-i 01 1

A = AzIT+A2T2

+ A3T3

+ A4T4

+ O(T5).

(3-86)

These expressions,

together with

the corresponding equations for

Ik and

Lk that were

developed

in the

previous

two subsections,

complete the analytical

development of the

short-term theory of general precession relative to the invariable plane of the Solar System.

Numerical Results and Verification

The right ascension &o and declination So of the angular momentum vector of the Solar

System were found in Chapter

to be

= 2735'109262;
= 6659'2811003.

(2-43)
(2-44)

From these angles, the right ascension of the ascending node of the invariable plane on the

mean equator of J2000.0 is

Lo = ao 270

= 351'091262,

(2-45)

and the inclination of the invariable plane to the mean equator of J2000.0 is

Io-= 90

- 60o = 23000'311!997.

(2-46)

Now in order to obtain numerical

values for the coefficients for the polynomials in

L, I, and

we also require numerical

values

for the coefficients of

0, and

IAU currently recommends the values of these quantities which were published by Lieske

et al. (1977) in

their Table 5.

However, these coefficients have been rounded, and using

'1 1 1 1* i 1 1 i i 1 1 1 I 1 r 1 i .

Table 3-1.

Full-Precision

Values of Current Precessional Constants

2306'1218108282

+0'.301

87986821595

+0'e017997590049701

Lieske

2004'13109489144
-0':42665200261276
-0'.041832591413886

s) from that printout are presented in Table 3-1

2306'!21810

82828

+1'!09467786

21317

+0'e018202974269150

These values are given to sixteen

digits, representing essentially the full precision of a Univac double-precision number.

The formulas for the coefficients in

the polynomial expansions of L,

and A

were

found in

the previous section in terms of the corresponding coefficients of the precession

angles (,

0, and

z and of the initial

angles

and lo.

Given

the numerical

values in

Table 3-1 and the above values for L0 and Io

the coefficients Lk,

k, and Ak have the

values listed in Table 3-2.

Table 3-2.

Constants in the Short-Term Theory

3051'09'!262
-96"!7230

-1'94824
0'006539
0'!0000881

2300'313997
-134'6685
0'!49754
O'1006173

-0'!0000188

01'000
5116"'1809

2'!92466

-0!005636
-0'!0000736

It must be noted here that the numbers in the last row of Table 3-2 were obtained

under the assumption that (4

=z4

as the theory developed by Lieske et al. (1977)

only goes to the third power of time.

There is no contribution from these coefficients until

the fourth power of time in any event.

The values for

L4, I4, and A4

arise solely from

products of coefficients with subscripts 1 to 3. Th4

be said to be deterministicc" in that they appear

e results in the last row of Table 3

as a consequence of the transformation

to the invariable plane; their true values cannot be found without first having values for

76

While L0 and I0 are intended to represent the orientation of the invariable plane in

the Qo system, there is no such restriction in the equations themselves.

Any nonrotating

plane will do. In particular,

if Lo is set to zero and I is set to Co,

the obliquity at epoch,

the invariable plane is replaced by the ecliptic of J2000.0.

In this instance, the resulting

L(T) is none other than the negative of the accumulated planetary precession )(A(T) (in

Lieske's

of J2000

notation)

. denoted

; I(T) is the angle between

)A(T) by Lieske; and A(T)

the mean equator of date and the ecliptic

is the accumulated luni-solar precession,

CA(T).

When this test

case

was run, the output coefficients duplicated Lieske's

results for

these angles to better than a microarcsecond.

This increases one's

confidence that both

the theory itself and the computer program that calculates the coefficients are correct-at

least in those terms independent of sin L0.

A further verification is possible by comparing the angles obtained by evaluating the

rigorous formulas for

and A in equations (3

-10) through (3

-12) with

the results

obtained from evaluating the time polynomials.

(Since the

"rigorous"

equations depend

on ,

0, and

, their absolute accuracy depends on the accuracy with which these angles

themselves are modeled.)

If these comparisons are made at regular time intervals on either

side of the origin, the differences between the rigorous angles and the polynomial approx-

imations should not show any systematic trend except proportional to the fifth or higher

powers of time.

This is indeed the

case

for various random values of L0 and Io in addition

to the two

cases

mentioned above.

the differences SL

, 61, and 6A between the rigorous

angles and polynomial approximations are presented in Table 3-3. Even after two centuries

Table 3-3.

T

Comparison of Rigorous Angles and Polynomial Approximations

-2.0
-1.5
-1.0
-0.5

0"0000078626

0'!0000001

8587

0'f0000002438
0'0000000076
O'fOoooooo000000000ooo
-010000000075
-0!0000002398
-00000018127
-0'0000076039

-0'0000034404

-0!000000

-0t0000001086
-0't0000000034
0O0o00000000o
0'0000000034
00000001107
0'!0000008443
0'10000035746

'0f0000103361
0'0000024600
0"!0000003249
0'0000000102
0ft000000000
-00000000102
-0"!0000003290
-0'0000025064
-0'!0000105968

Precession Between Two Arbitrary Times

The short-term theory as it has been presented thus far provides only for processing

from the standard epoch (J2000.0) to another date. In this section there will be developed

expressions for precession from one arbitrary time to another.

For the sake of clarity, the matrix P from equation (3-4) will now be denoted by

In this way the initial and final time arguments of P are given explicitly,

(3-115)

as is the dependence

, I, and A on the final time T.

The precession matrix from one arbitrary time T1

to a second time T2 can be thought

as first processing from T1

back to J2000.0, then from J2000.0 to T2.

Therefore

P(T1 -* T2)= P(0 -* T2) P(T1 0).

(3-116)

Because the precession matrix is orthogonal, this is equivalent to

- I -

P(0 -+ T) = R3[-L(T)] R [-I(T)] R3[-A(T)] R (Io) R3(Lo).

. of J2000.0

T(T,)

L(T1)

Eq. of epoch

T(T2-)_ L(T2)

of date

Figure 3-4.

Precession Between Two Arbitrary Times

P(Ti -- T2)

L(T2) R1

(T2)] R3[-A(T2)] Ri(Io) R3(Lo)}

x{R3

= {R3[-

L(T1)] R1

L(T2)

I(Ti)] R3 [-

I(T2)] R3

AT1)] Ri(Io) R3(Lo)}T

A(T2)] Ri(Io) R3(Lo)}

x {Ra(

= Ra3[-

Lo)RI(-

L(T2j)] R1

lo)R3

I(T2)] R3

(Ti)] R[I(T1)

R3[L(Ti)]}

- A(T2)] R1 [1(T1)]Ra[L(T1)],

(3-11

which is geometrically evident from Figure 3-4.

Equation (3-11

) in fact has the same

structure

as equation (3-4), except that Lo0

and Io

are replaced

by their values at Ti,

and A is replaced by the accumulated angle from T1

to T2. IfTi is zero,

equation (3-

reduces to equation (3-4) exactly. Also, if Ti and T2 be exchanged, equation (3-118) shows

that resulting matrix P(T2

-- T1) is the transpose of P(Ti1

as in fact it must be.

-* T),

= {R3

- A(T)

= (AiT2 +

T2 A3T 23+A424)

- (A Ti AT + A3T1 + A4T4 )

= [AI(T1 i+ t)+

A2(T + t)

2 A3(T1 t)3

+ A4(T+ t)4

- (A T1 + A2T2

A4T4)

= Ait+ A2(
+ A4(4T:

+ 6Tit2 +

A3(3T12 + 3Tit2

4Tit3

+t3)

+t4)

= (A

+2A2

Ti + 3A3T2

+4A4T13)

+ 3A3T1 6A4T12)t2

+(A2
+ (A3

+ 4A4T1)t3

(3-119)

The last form of equation (3

-119) is in the traditional two-argument form used by Lieske

et al.

(1977),

with T1 replacing their T inside the parentheses.

Numerically, with the expansion restricted to the third power of time,

A(T2)

A(Ti)

= (5116'.1809 +

'189432T1

- 0'.016908T2)t

+ (2'.92466

- 0"'!016908Tr1)t2

- 0'!005636t3

(3-120)

A more subtle numerical problem that occurs

-* 0 involves near cancellation in

the off-diagonal elements of P

. Clearly, in the limit as t

- 0, P must approach the identity

matrix

However, the individual terms in any one element will not vanish; rather, they will

tend to cancel.

Certainly the right-hand side of equation (3-118) can be multiplied out to

give the individual terms of the P

by judicious applications of the identity

cos x

the gross cancellations can be avoided.

sin2(lx-

However

(3-121)

the result is a dramatic increase in the

A(T )

- A3T3 3 -

80
author believes that a straightforward mechanization of (3-118), using (3-120) to compute

A(T2) A(Ti), is probably to be preferred.

Summary and Discussion

This chapter has developed the short-term aspect of precession theory using the in-

variable plane of the Solar System as a reference plane.

The classical precession matrix

that rotates from the Earth mean equator and equinox of epoch (T1 centuries past J2000)

to the Earth mean equator and equinox of date (T2 T1 + t centuries past J2000) is the

product of five elementary rotations:

P = R3[-L(T2)] RI[-I(T2)] R3[A(Ti) A(T2)

R1 [I(Ti)] R3 [L(Ti)],

(3-118)

where the angles on the right-hand side of the equation are computed by:

L(T)

= 351'09I.262

- 92'!7230T

- 1'!94824T2

+ 0'!006539T3;

(3-122)

I(T)

= 23000'311!997

- 134'16685T + O'049754T2

+ 0'!006173T3

(3-123)

A(T2) A(Tr)

= (5116'.'1809 + 5'189432T1 0"'016908T2)t

+ (2'!92466

- '0016908T1)t2

- 0'!005636t3.

(3-120)

This version of the matrix P agrees with the classical matrix (Lieske 1979) to better than

0"0001 for IT I

< 1 century,

with virtually all the difference due to the neglected "deter-

ministic" fourth-order terms.

The accuracy with which either theory tracks the actual processional motion of the

Earth depends on the time span over which the model for the precession angles (A, 0A, and

ZA is valid. Laskar (1986) found a fourth-order term in the accumulated general precession

only 2.5 centuries on either side of J2000.0.

This topic will be addressed in more detail in

Chapter 5, after the long-term theory will have been presented.

Equation (3-118) for the matrix P is unnecessarily

complicated, since any rotation

matrix can be expressed as the product of three, not five, elementary rotations.

The best

way to eliminate the extra two rotations would be to refer the origin for the right ascension

system not to the vernal equinox, as has been done by

European astronomers since the

Renaissance, but instead to the intersection of the Earth mean equator and the invariable

The precession matrix would then be written as

P' = Ri[-I(T2)] R3[A(T1) A(T2)] R1 i[I(Ti)];

now only equations (3-120) and (3-123) are required to evaluate the matrix.

This change would cause all right ascensions to be decreased by

variations in right ascension to be decreased by dL/dT

(3-124)

L(T) and all annual

. (There would be no change to the

declination system.) In particular, at the standard epoch J2000.0 one would subtract L0 =

15m24s.6175 from all right ascensions, and add

= 6118153/century to all centennial

variations min right ascension.

Such a change would involve a radical redefinition of the standard coordinate systems

that have been

used in astronomy for centuries.

The r6le of the vernal equinox as the

origin for right ascensions would be lost; in effect the ecliptic itself would be replaced by

the invariable plane.

The origin of right ascension would still be defined by the intersection

of two planes; however, only one of them (the equator) would be moving.

The resulting

theory,

when

put into

practice, is

simpler,

because only two different

polynomials

plane.

are

82

This is not to imply that the concept of the ecliptic has completely outworn its use-

fulness. Any analytical theory of precession must use the ecliptic.

over time, both the Sun and Moon are found in the ecliptic. The s

When suitably averaged

;ecular torque that causes

luni-solar precession will therefore be a function of the location of the ecliptic and of the

obliquity.

One of the major objections to using the invariable plane in this fashion will be: "But

we can't observe the invariable plane, and we can observe the ecliptic directly!" Both halves

of this claim cannot simultaneously be true if identical assumptions are used to evaluate

their truth.

Can we observe the invariable plane? There is no flashing beacon in the sky that marks

the plane, true, but that is also true of the ecliptic.

We do observe the apparent locations

of the planets; centuries of observations have furnished reliable mean motions and mean

orbits.

Spacecraft encounters have given

us both reliable initial conditions and dramatic

improvements in our estimates of the planetary masses.

So while we do not observe the

invariable plane per se, we certainly do have enough information at hand, all gleaned from

observation, to enable us to determine its orientation.

Can we actually observe the ecliptic? The Sun itself does not lie in the ecliptic.

we observe is the apparent direction to the Sun.

What

This is affected by the light-time correction;

by annual, monthly, and diurnal aberration; and by diurnal and monthly parallax.

(The

monthly effects are caused by the motion of the Earth about the Earth-Moon barycenter.)

A glance at an ephemeris of the Sun as expressed in ecliptic coordinates will confirm this:

the Sun's

geocentric latitude can reach almost one arcsecond in magnitude. All these effects

83

will determine an osculating orbit, to be sure; but this is not the ecliptic, as the ecliptic

is defined as the mean orbit (whatever that is) of the Earth-Moon barycenter.

therefore account somehow for planetary perturbations in

One must

deriving the orientation

motion of the ecliptic from observations.

Therefore both planes require for their determination the entire system of planetary

masses and orbits as well as models for reducing apparent positions to true ones. No object

(as a general rule) lies in either plane.

but the conclusion is unavoidable: if

The details of the reduction are of course different,

the ecliptic can be said to be observable, then so is

the invariable plane; and if the invariable plane is said to be unobservable, then so is the

ecliptic.

The only reason therefore to prefer one to the other is convenience.

simplicity of the precession

Is the relative

theory using the invariable plane outweighed by centuries of

observations (not to mention tradition) referred to the vernal equinox? The answer to this

question must be given by the astronomical community as a whole.

CHAPTER 4
THE LONG-TERM THEORY

Overview

As pointed out by Simon Newcomb (1906) and reinforced by Lieske et at (1977), it is

not possible to develop rigorous analytic formulas for the precession angles that are valid

for all time.

The reasoning behind this statement is that the equations of motion of both

the ecliptic and the mean equator cannot be integrated analytically in their most general

form. Lieske et al. took the velocity of the ecliptic pole vector to be parabolic in time; they

obtained it by fitting a second-degree polynomial through the values of the perturbation

function at three epochs.

This provided an approximate expression, as a cubic in

time,

for the ecliptic pole.

The cubic polynomial for the ecliptic pole was next substituted into

Newcomb's equation of motion for the mean celestial

pole;

another analytic integration

supplied an approximation for the location of the mean celestial pole; and the locations of

both poles yielded the desired precession angles.

Accurate precession angles over long periods of time must, however, be obtained by

numerical integration.

The orientation of the ecliptic pole-that is, of the vector normal

to the mean orbit of the Earth-Moon barycenter-has already

(1990) over an interval of +500,000 years from

the standard epoch J2000.0,

as part of

a revision of his

"Numerical

General

Theory"

of the Solar System

1985,

1986,

85

Given the ecliptic pole and the eccentricity of the Earth's orbit, the angular speed

of the celestial

pole can

be found.

This work uses the equations of motion in the form

developed by Kinoshita (1975, 1977), a more complete accounting of the torques on the

Earth than Newcomb's model.

Kinoshita's result for the speed of luni-solar precession is

then used to integrate the motion of the celestial pole.

The equations are close enough to

linear that a simple fourth-order Runge-Kutta integrator (Abramowitz and Stegun 1964,

Section 25.5.10),

with a constant step size, sufficed to give numerical accuracies on

order of a nanoarcsecond after 500,000 years.

Coordinates of both poles were written to a

file at regular intervals.

Given the poles of the ecliptic and equator, it is a straightforward matter to obtain the

desired precession angles.

These are found, not by the historical developments of spherical

trigonometry, but by the simpler methods of vector algebra. All of the classical precession

angles, as well as the angles L, I, and A using the invariable plane,

were determined at

each output point.

Since the resulting output file is quite large (one output point per century for a million

years), the results were condensed

by fitting Chebyshev polynomials through

the values

for each angle, in a manner similar to that employed in the production of JPL planetary

ephemerides

(Newhall

1989).

These Chebyshev

coefficients form

the final result of the

long-term theory.

The individual sections of this chapter describe the various facets of this procedure in

theory for the ecliptic; the development of the equations of motion of the

celestial pole, including Kinoshita's equation for the speed of luni-solar precession; their

The Motion of the Ecliptic

The problem of finding the orientation of the ecliptic, and indeed of the orbits of all

the planets, is an old one in the field of celestial mechanics. The work of Jacques Laskar

builds on the strong French tradition started by Laplace, Le Verrier, and Lagrange and

continuing today with Duriez (1977, 1979)

and Bretagnon (1982).

(NGT) uses a combination of analytical and

numerical methods to evaluate the mean orbital elements of the first eight planets (Laskar

1985). First the disturbing function was developed to second order in the planetary masses

and fifth degree in eccentricity and inclination.

This process yielded some 153,824 monomial

terms in the differential equations for the eccentricity and inclination of the planetary orbits.

These equations were then integrated numerically.

The initial orbital elements in Laskar's

work were based on

Bretagnon's (1982)

VSOP

theory;

the planetary masses are those

currently recommended by the International Astronomical Union (1976).

A subsequent paper by Laskar (1988) isolated the most prominent frequencies in the

eccentricity and inclination variables for the planets. In this technique, Laskar took a Fast

Fourier Transform of each component under study, found the frequency with the strongest

amplitude, subtracted that term from the original, and repeated the process to find the

next most significant term.

for the outer planets. Such

His work stopped after 80 terms for the inner planets and 50

Sa formulation would have been very useful for this work, but

for the fact that Laskar's Fourier representation does not preserve the initial conditions.

I therefore requested from Laskar a file of the eccentricity and inclination

variables

themselves. He very graciously sent me such a file, not from the original NGT but from a

= e sin w

= e cos w.

p = sin 1TIA sin IAA,

q = sin 1f7TACOSIHA,

where e is the eccentricity, w the longitude of perihelion,

longitude of the ascending node (of the descending node,

the Earth-Moon barycenter, all reckoned with respect to 1

q are identically zero at J2000.0 (T

4-1

(4-3)

(4-4)

rAI the inclination, and 11A the

the E

< 0) of the mean orbit of

System. Therefore p and

=0).

= v/h2 + k2

(4-5)

= plg(h, k);

(4-6)

sgnT sin

-1 Vp2 q2

(4-7)

= plg(psgnT

,qsgnT).

In equations (4-6) and (4-8), the notation follows Eichhorn (1987/88

w and of HA is sensitive to the signs of the numerator and denominator; this is the ATAN2

function of Fortran.

The factors sgn T

in equations (4-

7) and (4

) provide continuity for

TA and HA

as T passes through zero,

as p and q change signs then. (Otherwise the derivative

of rA would change sign instantaneously, and HA itself would suffer a discontinuity of 180

Equation (4-

) cannot be evaluated at T

= 0; then HA assumes its limiting value.

The equations of motion will turn out to be evaluated most easily if the unit vector

directed to the ecliptic pole is expressed in terms of its rectangular coordinates.

Relative

88

In order to save computer time during the integration, the p and q values for each

record were transformed at the start into s and c. These were then stored along with h and

k in 1001-element arrays, which were interpolated as needed. The interpolating polynomial

was selected to be the unique eleventh-order polynomial passing through twelve consecutive

points bracketing the interpolation time. (To improve the numerical stability of the process,

the Chebyshev representation was chosen in lieu of the simple polynomial representation.)

Except at the very beginning or end of the file, six points out of the twelve preceded the

interpolation time and six followed it.

This scheme guarantees that the interpolated result

will be the tabular value at each point in

the table;

the discontinuity

between regimes

enters only at the first derivative, and this should be small due to the high degree of the fit.

Subroutines DPFIT and DCPVAL, of JPL

s MATH77 library (JPL Applied Mathematics

Group 1987), performed the polynomial construction and evaluation, respectively.

The Equations of Motion for the Celestial Pole

The development of the equations of motion for the mean celestial pole (the vector

normal to the mean equator of date) proceeds in several logical steps.

equation of motion in vector form.

One begins with the

From this relatively simple equation the derivatives of

the components of the pole vector are found, since it is the components that are actually

integrated.

Finally

Kinoshita's expression for the rate of luni-solar precession is adopted,

and the appropriate constants of his theory are found.

4.3.1.

The Vector Equation of Motion

The principal source of the motion of the celestial pole is a torque produced by the

89

precession may be obtained by considering the Sun and Moon to be smeared into a ring

around the ecliptic, and the oblateness of the Earth likewise to be concentrated into a ring

around the equ

the ecliptic, e.

ator.

These two rings have a mutual inclination equal to the obliquity of

The resulting torque on the equatorial ring is in the direction of the line

of intersection of the two planes (the vernal equinox), and its magnitude is proportional

to sin 2e. A detailed derivation is given in chapter 8 of Woolard and Clemence (1966) and

need not be repeated here.

The second contribution to the motion of the celestial pole is the so-called

precession."

"geodesic

Although de Sitter (1916) first showed that general relativity predicts that

spinning bodies in a gravitational field will undergo a forced precession, only two decades

later did

de Sitter and

Brouwer (1938) include its effect

in luni-solar precession.

effect of geodesic precession is also a torque in the direction of the vernal equinox, but its

magnitude is proportional to sin e instead of sin 2e.

Since sin2e = 2sin ecos e, the magnitudes of the two torques may be combined into

an expression of the form

= (P cos

pg) sin

(4-10)

where Q is the celestial pole vector, a unit vector; P is Newcomb's

"Precessional Constant"

(Newcomb 1906, p.

228); and pg

is the rate of geodesic precession.

Neither P nor pg is quite

constant, as both depend slightly on the eccentricity of the Earth's

orbit, which varies over

time.

Now the dot product of two unit vectors gives the cosine of the angle between them,

and the magnitude of the cross product of two unit vectors is the (positive) sine of the an1le