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

MISSING IMAGE

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






























Copyright () 1990

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


s Navigation Section.


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.


















TABLE OF CONTENTS


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

the currently-adopted -46'.8150/century.
















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


about


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.)


In addition,


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"


of Laskar (1985) provided the









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


, completes the right-handed triad.


(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


This task is made


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


(1 Ceres) has about 6


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


x 10-7









of about five less than this:


about 6 x 10-8 radian, or 0'!01.


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.









task.


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


angular momentum is instead


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


plane about the apsis, about the latus rectum, and about about the orbit normal; and the


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-


rad).


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


radians2


(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


follow easily.


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


the four-quadrant arctangent.


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)


Additional expansions are required,


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


z is added,


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.


For the adopted invariable plane,


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


been obtained by Laskar


(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


(Laskar


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


detail: Laskar's


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)

Laskar's "Numerical General Theory"


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).


The standard orbital elements may be retrieved from these by


= 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


as the quadrant of


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