UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations   Help 
Material Information
Subjects
Notes
Record Information

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 fiveyearold 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 (19011981) 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 SHORTTERM 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 LONGTERM 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 Lunisolar 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 LongTerm Theory . . CONCLUSIONS The ShortTerm and LongTerm 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' FullPrecision Orbital Angular Momenta Values of Current Precessional Constants * S S 5 0 5 * S S S S S S Constants in the ShortTerm Theory Comparison of Rigorous Angles and Polynomial Approximations Astronomical Constants for Kinoshita's Lunisolar Precession Model Chebyshev Polynomial Coefficients for the LongTerm Theory . Polynomial Coefficients from the LongTerm Theory Near T Chebyshev Coefficients for (A 9A, and ZA from the LongTerm 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 LongTerm 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 LongTerm and ShortTerm 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 "shortterm 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 millionyear 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 longterm and shortterm theories reveals two possible improve ments to the IAU system of constants. the rigid Earth, used in The higherorder terms in Kinoshita's model for the longterm 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 longterm theory, yields 46'!8065/century for the rate of change of the obliquity versus the currentlyadopted 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 longperiod part is "lunisolar precession." (As the longest nutation period is about 18.6 years while the shortest period of lunisolar 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 lunisolar and planetary precession is accordingly called "general precession. Lunisolar 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 lunisolar 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 thirddegree 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 longperiod 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 shortterm and longterm 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 "shortterm theory." Chapter 4 is devoted to the "longterm 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 11, 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 sansserif 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 11 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 righthanded: if the three axes be denoted (in order) by x, y, and z, and if the x and yaxes are placed on a piece of paper such that x increases to the right and y increases upward, then the zaxis 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 / (11) this is a vector equation. The inverse transformation consists of the three scalar equations r=V/ +y2+ (12) A = plg(y, (13) f = tan (14) The notation "plg(y, xa)" in equation (13) 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 (12) and (14); r is nonnegative, and /3 has the same range as the principal value of the arctangent. Equation (14) 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 xaxis toward the old yaxis. axes. The standard zaxis 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 zaxis (denoted Q) of this system is directed to the mean Celestial Ephemeris Pole; the xaxis is directed toward the mean vernal equinox and is symbolized by the "ram's horns" of Aries (T); the yaxis (yq ), directed toward right ascension 6h and declination 0 , completes the righthanded 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 11 shows the coordinate axes of this system, with circles marking the planes of the equator and the ecliptic. Ecliptic Equator Figure 11. Equator  The Equatorial Coordinate System   E I system; and again the yaxis completes the righthanded triad. The basis unit vectors are denoted T for the xaxis , YE for y, and E for Figure 12 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 ) singleletter 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 zaxis. (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 nonrotating system the B1950 (FK4) system above at epoch will be denoted "EME50" in accordance with JPL Table 11. 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 EarthMoon 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 11, 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 ;nnn +r 1 4, . 1 ..... k:s 1... L ri. T rTr n  1. Table 11, 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 lunisolar precession attributable to the Sun. Fundamental rate of lunisolar 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 lunisolar preces sion. Mass. The reciprocal mass of a planet (mo/m). The rotation matrix which transforms from mean pnnatnrial rnvr Table 11, 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 socalled 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 11, 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 zaxis. 1) The radius of a rotating body. = 1,2, or 3 for the , y, or 4.3.3 2) Kinoshita's rate of lunisolar 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 11, 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 lunisolar preces sion. 1) In cylindrical coordinates, the distance from a point to the zaxis. ) 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 longterm 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 yaxis 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 11 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 11, 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 11, 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 lunisolar precession; the rate at which the mean vernal equinox of date moves along the moving eclip tic of date. The speed of lunisolar 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 lunisolar 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 11, 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 computerreadable 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 1010 solar mass (Schubart 1974; Newhall et al. 1983) or 6 x 107 of about five less than this: about 6 x 108 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 deepspace 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 groundbased 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 xaxis 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 loworder 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 elevenbody 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), (21) 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) (22) Ti) the ak are the coefficients read from the file: and n is the decree of the highest nnlvnnminal dx dT dr akTkMr dT (23) 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 lowlevel subroutine. Users need only call a higherlevel 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 21. 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 1013 10"14 1014 1.6 x 1015 2.6 x 1015 0 km3/ AU3/day2 AU3/day2 AU3/day2 AU3/day2 AU3/day2 149597870.6094344 km 299792.458 81.300587 km/sec Table 21. task. The location of the Solar System barycenter is defined implicitly by equation (2) of Standish et al. (1976): (24) 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 (25) 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 (24) and (25) are interdependent; in practice, the Sun s position is not integrated at all but is inferred from equation (24). The orbital angular momentum of body i is then found by (Burkhardt 1982): (26) hi= *ri and the total orbital angular momentum is X Vi. (27) 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 (26) of the eleven bodies on the file are presented in Table relativistic masses, from equation ( those of the five asteroids appear in Table 25), are presented in Table 24. 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 (29) = 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 wellknown 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 periodan 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 23. 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 24. 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( 212) and is assigned a mass PB = Pe +[k. When the EarthMoon 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 21 shows the difference between the vector h computed by equation (2 7) and that computed using the EarthMoon 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.6year 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 21. 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 twelveyear synodic period of Jupiter and Neptune. Figure 22 displays the right ascension and declination of h as a function of time the entire sixteenyear span of the M04786 ephemeris. The twelveyear periodicity is quite obvious in the figure. (Corresponding plots using the DE130 ephemeris, shown to the same scale in Figure 23, 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 (27). 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, (221) 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 (222) 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 (224)  e2) Finally, the uncertainty in h must be rotated into EME50 equatorial coordinates, STyz MT SpqhM. giving (225) 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 26; 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 interplanet 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 (227) h +h/h2 2 9{hs,hyhz hxhz/h Table 26. 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 1013 = 2.5407 x 1014 2.5407 x 1014 4.4279 x 1015 radians2 (229) The resulting standard errors in the orientation of the total orbital angular momentum of the Solar System are oa cos 6S = 0'103865, (230) = '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 JupiterNeptune couple (Figure 22) 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. (232) In the typical case of solidbody rotation about the zaxis, transformation into cylindrical coordinates (s, , z) yields s2Q dm, (233) 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 (234) for constant p, = TrpRi5 MR2Q  IT. (235) An evaluation of equation ( 32), the general case will have the same functional form equation ( 235), 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 235), 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 rnxr1 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 (shortterm and longterm) 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 (236) When the vector h from equation (218) 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 24. 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 24 shows the angles Io and L0 as well as the spherical coordinates of h. These quantities are related by Lo = ao + 6h (239) 3t852572907 (240) lo = 90 (241) = 232008888075 = 2300'31'.99707. (242) = 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; (243) 66059'28'.003 0013; (244) 3051'09'.262 0"038; Lo = lo = (245) (246) 2300'31.997 + 0!013. These values will be used throughout Chapters 35. CHAPTER 3 THE SHORTTERM 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 currentlyaccepted shortterm 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 (31) 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 (32) 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) (33) and analogously for the other angles. Figure 31 shows the coordinate axes of both teams and the three "classical" angles. Reading the righthand side of equation (31) from right to left, the rotation Ra( ) places the initial yaxis, 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 yaxis 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). (34) Figure 32 shows the coordinate axes of epoch and of date from the same perspective as Figure 31; now the invariable plane and the five angles of equation (34) appear. magnified drawing of the region near the vernal equinoxes also appears as Figure 33. equation (34), the first (rightmost) rotation moves the xaxis to the intersection of the equator of J2000.0 and the invariable plane. The second rotation puts the yaxis (and consequently the entire xy 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 31. The Classical Precession Angles places the yaxis on the equator of date, and last a rotation about the vector Q (the final zaxis) positions the xaxis at the vernal equinox of date. The purpose of the shortterm 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 (34) 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 fourthdegree 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). (35) The desired angles may be isolated on the righthand 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). (36) One can expand both sides of equation (36), 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 (36) on the left by Ra( This rotation combines with R3( L) on the righthand side to yield R2(9) Ra3 (Lo + () I) Ra( A). (37) Next the two sides are expanded. The lefthand 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 (38) 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 righthand 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(39) 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]. (310) 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 (312) where plg(y, x) is defined in Section 13 to be the fourquadrant 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 (312) are therefore rigorously correct for all possible values of (, nations (310) 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 (311) 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 (310) through (312) above, their behavior with time depends on the behavior of the "classical" angles (, 0, and z that appear on the righthand 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. (313) These must be substituted into equations (310) through (312) 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). (314) 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); (316) 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). (317) 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 righthand side gives Ay= tan1 X0 + Ax dtan S(Ax) + 2! d2 tan (Ax)3 + 3! d3 tan S(Ax)4 + 4! d4 tan + 0 [(Aax)]. (318) Since xo = tan Yo, the first four derivatives can be expressed as dtan1 x dx d2tan1 x dx2 d3 tan1 x dx3 d4 tan1 x 1 + tan2 Yo 2tanYo (1 + tan2 yo)2 6 tan2 Yo (1 + tan2 yo)3 24(tan yo co2 = COS (319) (320) = 2 sinyocos3 yo; = 6 sin2 Yo0 COS4 Y0 2 cos6 yo; tan3 Yo) d4 /t i 9 \ (321) = 24(sin yo cos7 Yo sin3 yn cos5 /n ). (322)  f jsinLo)T2  ClC2 22   [(Cila + 1 = x1T +x2T2 + x3T3 + x4T4 + 0(T). (323) 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); (325) + O(T). (326) 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) (327) with the coefficients Yk given by = x2 ~2 sin Y cos3 ; xi1SInDoc0 CO yo; (329) = 3 cos 2x lX2 sin yo cos3 yo + x3(sin2 Y0o cos4 o i cos6 Yo); (330) = 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) (331) 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 = cos1 ([1 a 2T2  0102T3 (003 + 2*,2  4T4)] cosio + [01T + 02T2 + (0361 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)) (332) = 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 (333) 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) } (334) = cos 1 [cosIo + CT + c2T2 +c3T3 + cT4 + O(T5)], (335) where the Ck in equation (335) are defined by the various coefficients of Tk in equation (334). Now equation (335) has the form c= ciT+c2T2 + c3T3 + cT4 + 0(T5). (337) 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). (338) Since xo = cos lo, the leading term is simply Io. The first four derivatives become d cos x dx d2 cos1 x dx2 d3cos1x dx3 d4 cos1 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 (339) (340) (341) sin5 o 9 cos Io + 6 cos3 Io (342) sin7 Io When one substitutes into equation (338) the expression in equation (337) for c, along the values of the derivatives in equations (339) through (342), 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) (344) =Io+I T+I2T2 + 3T3 + 4T4 +0(Ts5). (345) 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 (334) into equation (344). After some tedious algebra, one obtains Cl sin Io = 01 sin Lo; c2 sin Io O? sin Lo sin Io sin Io (346) 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; (347) 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)]; (348) 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(302(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 + j902 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) (349) Equations (346) through (349) 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 (311): L = plg [cos 0 sin(Lo + () sin Io sin 0 cos Jo, cos(Lo + () sin 1o] + (311) = tan _1 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), (350) 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 (311) by n and the denominator by d. Then inserting the results of equations (314) through (317) 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) (351) 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). (352) If the leading terms in both n and d are factored out, the argument of the arctangent in equation (311) 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 (353) since sin Io $ 0. (The minus signs in the denominator cancel some of the minus signs in equation (352) 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; 354) n2 = (2 cot Lo  02 cot Io csc Lo; (355)  1/"3 641  12C1) cotLo ((1(2 + 102)  (03 l) cot Io CSC Lo; (356)  2 ,2 + 2 l2 Lo + ((4 + O(T5) y + O(TS) '  (02l 12) d = C(1 tan Lo; (358) tan Lo + 21; (359) A13)tanLo+ (2; (360) (361) The next step is to expand the quotient that forms the argument of the arctangent function in equation (311). 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) (362) + 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) (363) = 1 + qi T + q2T2 + q3T3 + q4T4 + O(Ts5). (364) This leaves equation (311) in the form l(qtan Lo) + (365) = tan1 { tan Lo + tan Lo[qiT + q2T2 + q3T3 + q4T4 +O(T5)]} +ziT+z2T2 + z3T3 + z4T4 + O(T5). (366) The expansion for the arctangent was given above in equations (327) through (331). 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 (367) 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. (370) It is easiest to leave the tan Lo with the qk. Now replacing nk and dk with their values from equations (354) through (361) 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 (371) = {[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; (372) 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 13l _ !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]; (373) 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 22) 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 jf 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]. (374) The final step is to substitute these values for the qk tan Lo into equations (367) through (370). For L1, we obtain Li = (qitanLo) Lo + zi (367) = (C1 Lo 091 cot Io sec Lo) cos Lo + zi = C1 + zl 01 cot Io COs Lo. (375) Next, we get for L2: L2 = (q2 tanLo)cos2 Lo (qi tan Lo)2 sin Lo cos3 Lo + z2 (368) = 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). (376) +( f 4 + 3(12( + cot Io Lo  0 6 65 The righthand side of equation (369) 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) (377) 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); (378) (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: (379)  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 (380) Finally, equation (370) for L4 contains four terms on its righthand 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]; (381)  [(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 (382) 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]; (383) 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)] (384) + 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 Lo2 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)]. (385) Equations (375), (376), (380), and (385) 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 timeconsuming 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 + ()] . (312) We will expand the righthand 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). (386) As in the previous subsection, denote the two arguments of the plg function in equation (312) 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 (314) through (317) 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); (387) 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). (388) 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 (312) in the form n1T+n2T2 diT  d2T2 + n3T3  d3T3 + n4T4 +O(T5) (389)  d4T4 + O(T5) with the coefficients nk and dk given by nl = 01 cos Lo csc lo; (390)  1 2 A =  2,1T +(03o 1^2  2(" 92  1011 (1  (lC(3 + 2C2 +3C1 S 0 l1)sminLo]jcsclo; (393) (394) d2 = 02 sin Lo cot Io + 09 (1 cos Lo cot Io + 1 9; d3 = (03 Of3 0i() sin Lo cot Io + (01 (395) (396) 91022 2li2n ~ 21 + (013 + 02(2 + 03(1 +(1+02 31(9103+"12 2 The next step,  0i(1i2)sinLocotIo 10 6l l  91 Ci) cos Lo cot Io lOf ). ~ 24 t/ ) (397) as before, is to expand the quotient that forms the argument of the arctangent function in (312). 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) (398) = 1+niT+(n2+ddn1)T2 + (n3 + din2 + d2n1 + dzni)T3 + (n4 + din3 + d2n2 + d3nl + d n2 + 2did2ni + dinl)T4 +O(T5) (399)  1+qiT+q2T 2 + q3T3 + q4T4 + O(T5). (3100) 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) (3101) = [02 cos Lo 0~i1 sin Lo + cot lo(0i sin Lo cos Lo)] csc lo;  .i 1. ,7Iy _L ,7Lyi r, . (3102) 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)] (3103) 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. (3104) This leaves equation (312) in the simple form tan 1 q (3105) = tan1 [qiT + q2T2 + q93T3 + q94T4 + O(T5)]. (3106) Once again we employ equations (327) through (331). This time, however. Xn is zero. = {(o3 + 3 csc lo } csclo} L0ocot I0o + )  01(1(2)cosLo  (81  Of3  (l2 A1 = ql (3107) = 1 cos Lo csc lo; (3108) A2 = q2 (3109) (3110) = [02 cos Lo 01(1 sin Lo + cot Io(9i sin Lo cos Lo)] csc lo; (3111) = {(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 (3112) A4 = 94 9l 92 (3113) = { (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). (386) 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 shortterm 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. (243) (244) 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, (245) and the inclination of the invariable plane to the mean equator of J2000.0 is Io= 90  60o = 23000'311!997. (246) 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 31. FullPrecision 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 31 2306'!21810 82828 +1'!09467786 21317 +0'e018202974269150 These values are given to sixteen digits, representing essentially the full precision of a Univac doubleprecision 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 31 and the above values for L0 and Io the coefficients Lk, k, and Ak have the values listed in Table 32. Table 32. Constants in the ShortTerm 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 32 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 lunisolar 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 correctat 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 33. Even after two centuries Table 33. 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 shortterm 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 (34) will now be denoted by In this way the initial and final time arguments of P are given explicitly, (3115) 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). (3116) 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 34. 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)], (311 which is geometrically evident from Figure 34. Equation (311 ) in fact has the same structure as equation (34), 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 (34) exactly. Also, if Ti and T2 be exchanged, equation (3118) 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 (3119) The last form of equation (3 119) is in the traditional twoargument 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 (3120) A more subtle numerical problem that occurs * 0 involves near cancellation in the offdiagonal 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 righthand side of equation (3118) 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 (3121) the result is a dramatic increase in the A(T )  A3T3 3  80 author believes that a straightforward mechanization of (3118), using (3120) to compute A(T2) A(Ti), is probably to be preferred. Summary and Discussion This chapter has developed the shortterm 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)], (3118) where the angles on the righthand side of the equation are computed by: L(T) = 351'09I.262  92'!7230T  1'!94824T2 + 0'!006539T3; (3122) I(T) = 23000'311!997  134'16685T + O'049754T2 + 0'!006173T3 (3123) A(T2) A(Tr) = (5116'.'1809 + 5'189432T1 0"'016908T2)t + (2'!92466  '0016908T1)t2  0'!005636t3. (3120) 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" fourthorder 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 fourthorder 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 longterm theory will have been presented. Equation (3118) 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 (3120) and (3123) 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 (3124) 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 lunisolar 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 lighttime 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 EarthMoon 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 EarthMoon 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 LONGTERM 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 seconddegree 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 polethat is, of the vector normal to the mean orbit of the EarthMoon barycenterhas 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 lunisolar precession is then used to integrate the motion of the celestial pole. The equations are close enough to linear that a simple fourthorder RungeKutta 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 longterm 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 lunisolar 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 EarthMoon barycenter, all reckoned with respect to 1 q are identically zero at J2000.0 (T 41 (43) (44) 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 (45) = plg(h, k); (46) sgnT sin 1 Vp2 q2 (47) = plg(psgnT ,qsgnT). In equations (46) and (48), 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 1001element arrays, which were interpolated as needed. The interpolating polynomial was selected to be the unique eleventhorder 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 lunisolar 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 socalled 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 lunisolar 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 (410) 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 