THE EARTH
GRAVITATIONAL
FROM
SATELLITE
POTENTIAL
DERIVED
MOTION
YOSHIHIDE KOZAI Tokyo Astronomical Observatory, Mitaka, Tokyo, Japan
(Received April 2, 1966) Abstract. Recent works on the determination of the earth gravitational potential by its dynamical effects to the motion of artificial satellite are reviewed. Future programs to improve the results and to detect from observations effects due to earth tide and any other time variations of the potential to the satellite motions are discussed. 1. Introduction
Even before the first artificial satellite was launched, geodesists and astronomers believed that tracking of satellites would be very useful for geodesy, although many of them thought at that time that by tracking two or three satellites the shape of the earth would be known accurately. However, it turned out that the potential of the earth, that is, the shape of the geoid was not yet known completely, even though many satellites have been tracked for more than eight years. The first unexpected discovery by the tracking was made by O'K~ErE, ECI~FLS and SQUW,Es (1959), who concluded that long-periodic variation with a eighty-day period in the perigee height for 1958 fi2, Vanguard 1, was caused by north-south asymmetry in the potential, that is, zonal spherical harmonics of odd orders. It was also found even in 1958 that the coefficient of the zonal spherical harmonics of the second order, J2, determined from satellite motions was considerably less than the value derived by classical geodesy. The flattening of the international ellipsoid being 1/297.0 and the value of J2 for this being 1.092 x 10 .3 ( J = 1.638 x 10 -3) (JEFFREVS, 1959), whereas that derived by satellite motions was less than 1.0827 x 10 -3, from which the flattening has been computed as 1/298.25. However, the values of J2 determined by various authors were scattered between 1.0822 x 10 .3 and 1.0833 x 10 -3, three years ago, although this discrepancy has been reduced recently. The discrepancies in the determination of Jz as well as in other coefficients were due mainly to the assumption for deriving these values, that the values of Jj decrease rapidly with increasing j. If this assumption were correct, the tracking of two or three satellites would be enough to determine all the values of Jj which are not negligibly small. The fact that Jj converges very slowly with increasing j is also one of important discoveries by satellites. By using existing data higher order terms which cannot be neglected have to be disregarded to determine coefficients in the potential. Therefore, we need more data for more satellites with various inclinations for satellite geodesy. Space Science Reviews 5 (1966) 818-879; 9 D. Reidel Publishing Company, Dordrecht-Holland
THE EARTH GRAVITATIONAL POTENTIAL
819
The first determination of the ellipticity of the equator by satellite motions was made by IzsA~ (1960). However, for this case also, it was found that coefficients of non-zonal harmonics do not converge rapidly, and that the shape of a geoid is much more complicated than a tri-axial ellipsoid. The determination of coefficients of nonzonal terms is much more difficult than that of zonal terms because of differences in the periods of perturbations, and in order to determine them it is necessary to improve coordinates of tracking stations from observations simultaneously. Even if we can derive the potential very accurately, the question whether this is constant or not still remains. We know that the earth is deformed by tides. And NEWTON (1965) succeeded in detecting effects of the tides by observing the satellites. The geocentric gravitation constant, GM, has also been determined accurately by tracking Ranger satellites, as well as the lunar gravitation constant which gave the mass of the moon. Theories of motion of satellites will not be discussed here, since Coo~ (1963) already reviewed them in Space Science Reviews. And for the methods of tracking there was an article by VEIS (1963). Readers should also refer to papers by SHAPIRO (1963) and KAULA (1962) and to books by M~J~LLER (1964) and BERROTH and HOFMANN (1960). 2. Motion of Satellites
When the earth is assumed to be a point-mass, that is, to have density distribution of spherical symmetry, the equations of motion for an artificial satellite moving around the earth have rigorous solutions indicating that the orbit for this problem is an ellipse fixed in space (Figures 1 and 2). The position of the orbital plane is determined by the inclination, i, to the equator of the earth and the longitude of the ascending node, ~2, measured from the equinox. The direction of the major axis of the ellipse is determined by the argument of perigee, co, the angle measured from the ascending node to the perigee. The semi-major axis, a, and the eccentricity, e, fix the size and the shape of the ellipse, respectively. orbital plane
Fig. 1. Orbital plane and the equator.
820
YOSH1HIDE KOZAI
O
Fig. 2.
Ellipse.
The mean motion, n, that is, the mean angular velocity of the satellite, is related to the semi-major axis as,
(1)
n2a 3 = GM,
where G is the gravitation constant and Mis the mass of the earth. The mean anomaly, l, is defined as, l=_l ndt+~, (2) where ~c is constant. The six constants a, e, i, ~c, co, and O are called Keplerian orbital elements. Once they are given, positions of the satellite can be computed as functions of time. The actual procedure for computing geocentric equatorial coordinates of the satellite is as follows: (i) Compute the mean motion and the mean anomaly by Equations (1) and (2). (ii) The eccentric anomaly, u, and the true anomaly, f, corresponding to the value of l are derived by the equations, u - e sin u = l, tan f 2
|
/l + e tan u /
~/1 - e
(3)
2
(iii) The geocentric distance, r, is computed as, r-
P -a 1 + ecosf
(1 - e cos u),
(4)
where p is the parameter of the ellipse and is defined as p = a (1 - e2).
(5)
821
THE EARTH GRAVITATIONAL POTENTIAL
(iv) Then the geocentric rectangular coordinates of the satellite referred to the equinox and the equator can be derived as,
x/r y/r z/r
i = cos (L + (2) + 2 sin 2 ~ sin L sin ~ , i = sin (L + O) - 2 sin 2 - sin L cos f2, 2 = sin i sin L ,
(6)
where L is called the a r g u m e n t of latitude, that is, L -- f + co.
(7)
(v) The geocentric coordinates are transformed to the topocentric coordinates expressed by the topocentric distance, p, the right-ascension, c~, and the declination, 6, as x - X = p cos 6 cos ~, / y-Y=pcos6sine, } (8) z - Z = p sin 6,
/
where X, Y, and Z the earth is rotating ordinates (X, Y, Z ) coordinates of the equator as
are the geocentric equatorial a r o u n d the z-axis, X and Y are referred to the equinox can observer (4, t/, ~) referred to
(i) (V =
si /~c
coordinates of the observer. Since functions of time. However, the cobe c o m p u t e d from the geocentric the Greenwich meridian and the
cos/~G 0
i)(1) ,
(9)
where /~G is Greenwich sidereal time. Z ( North Pole )
x
x
lit~tY or
/ y
GreenwichMeridian Fig. 3. Moving coordinates (X, Y) and fixed coordinates (x, y).
822
YOSHIHIDE KOZAI
When the equations of motion for the satellite are solved by using the actual potential of the earth, Keplerian orbital elements are no longer constant, but are functions of time. However, formulae (2) through (8) can be used with the instantaneous values of the orbital elements to compute the coordinates of the satellite, although instead of (1) modified Kepler's third law of the form of (26) given later should be used to compute the semi-major axis. The quantities e, 6, and p in the expressions (8) and dp/dt can be measured by tracking satellites. Most of the tracking data used in satellite geodesy come from (i) optical measurements of angular coordinates (e, 6), (ii) range (p) and/or range-rate (/5, Doppler shift) measurements by radio waves, and (iii) by radio interferometry. One of the typical cameras for the optical tracking is the Baker-Nunn camera of the Smithsonian Astrophysical Observatory (VEIS, 1963). (See Figure 4.) Accuracies within 2 seconds of arc in position (ct cos b and c5) and 1 millisecond in timing can be obtained by using this type of cameras. In other words, accuracy within 10 meters in the geocentric coordinates, that is, on an average 10 .6 in (x/r, y/r, z/r) can be obtained. Radio methods are not so precise, being limited at present to a precision of 30 meters (NEWTON, 1964). However, the radio methods have a compensating advantage, since the measurements can be made at any time of day and in any weather. On the
Fig. 4. Baker-Nunn Camera.
THE EARTH GRAVITATIONAL POTENTIAL
823
other hand, the optical measurements can be made only when the satellite is visible, that is, when it is night at the observing station, the satellite is sun-lit, and the sky is clear. Therefore, radio data are sometimes more useful than optical data in satellite geodesy. In order to derive as much information as possible from observations of satellites we need theories of motion which can compute positions of the satellite with the same accuracy as that of observations, that is, accuracy of 10 .6 . The orbit of the satellite is disturbed by non-spherical terms in the potential of the earth, the lunar and the solar gravitational forces, the solar radiation pressure, the atmospheric drag, and so on. (i) The gravitational potential of the earth is expressed by
1/= --
1-
J,
P~ (sin (p)
Y n=2
+
p~ (sin q0 (C,,,, cos m2 + S,,~ sin miC),
(10)
n=2m=l
where the potential is evaluated at a point whose geocentric spherical coordinates are, r, the geocentric distance, qg, the geocentric latitude, and 2, the longitude measured eastwards. In the expression (10) R is the equatorial radius of the earth, P, is the Legendre polynomial of the order n, and p~ is the normalized associated Legendre polynomial defined as, (n + m)! (1 - x) "/2 dmp" (x) ( 2 n + 1) ~ - - rn)! dx m Of numerical coefficients J,, C,,m and S,,m in the expression (10) J2 takes the largest value of the order of 10 .3 . Other coefficients take values of order of 10 .6 or less. (ii) The disturbing forces due to the sun and the m o o n have factors (nQ/n) 2 and (n~/n)Zm~, respectively, where no, n~, and m> are, respectively, the mean motions of the sun and the m o o n and the mass of the m o o n in units of the mass of the earth. When the mean motion of the satellite, n, is expressed in degrees per day, we have
(no~n) z = 0.9715/n 2 , (n~/n) z m~ = 2.1354/n 2 . For usual satellites n is between 3000 ~ and 5000 ~ per day, Therefore, n -2, and the disturbing factors are of the order of 1 0 - 7 . (iii) The acceleration due to the solar radiation pressure to the satellite close to the earth depends on the area-to-mass ratio, A I M s, of the satellite. When A/Ms is given in cgs units, the acceleration is 4.6 x 10 -8 A / M s times as large as the total acceleration due to the earth. For satellites used for geodesy A / M s is of the order of 10- i.
YOSHIHKOZAI IDE
824
3. Reduction of Observations
When approximate orbital elements are known accurate values of the orbital elements are determined by a differential orbit improvement program, that is, by comparing observed coordinates and/or range-rate of a satellite with computed values based on the approximate values of the orbital elements (VEIs and MOORE, 1960; MERSON, 1963 ; SOCHILINA, 1963; GAPOSCHKIN, 1964; ESCOBAL, 1965). The (O-C), de, dr, and dp, in the topocentric right-ascension, the declination and the range are expressed by corrections in the geocentric rectangular coordinates of the satellite approximately as, cos 6 de d6
~
dp/pj
where B is a rotation matrix such as, - sin e sin 6 cos e \ cos 6 cos e
B = ~-
(dx)
1 =
B P
dy
(11)
,
dz
cos -sin6sine cos 6 sin e
0)
cos sin
.
Coefficients in formulas expressing the corrections dx, dy, and dz approximately as linear functions of the corrections in the orbital elements are derived by differentiating formulas (6) which express x, y, and z as functions of the orbital elements with respect to these orbital elements. As the inclination, i, and the longitude of the ascending node, f2, appear explicitly and the argument of perigee, o9, appears only through L = f + e ) , the differentiation with respect to these elements can be made easily. The other elements, that is, the semi-major axis, a, the eccentricity, e, and the mean anomaly, l, appear in (6) through the radius, r, and the true anomaly, f. And the corrections in r and f are expressed approximately as linear functions of the corrections in a, e, and I as follows:
r
dr=-ada+~f~,
12ae
_e~Sinf'dl-ac~
/ (12)
where p = a(1 - eZ). The corrections in ! and a are transformed to those in the mean anomaly at the epoch, K, and the mean motion, n, as dl=d•+tdn, da/ a = - ~ dn/n .
THE EARTH GRAVITATIONAL POTENTIAL
(costa
825
Then the (O-C)'s in the observations including range-rate are expressed by dE, a vector whose components are the corrections in the orbital elements as,
de5
dp/p
= AdE,
(13)
df)/p / with a known matrix A which is a function of the orbital elements and time. As the correction vector dE consists of six components, the equation (13) has six unknowns and, therefore, it is necessary to have at least six observations to solve it. Usually, dE is solved by the method of least squares by using a sufficient number of observations covering a certain interval of time, during which dE can be assumed to be constant. The orbital elements of the satellite are changing very rapidly owing to perturbations. Any one of the orbital elements, E, is expressed as
E = E ~ + bE s + (~El + 8Ep,
(14)
where E ~ is a constant term and 6E s, (SE1, and 6Ep indicate secular, long-periodic, and short-periodic terms, respectively. Of short-periodic perturbations produced by many sources only those due to d2-term in the gravitational potential of the earth are very important and can be computed with sufficient accuracy by analytical or numerical methods with an approximate value of J2 if approximate values of the orbital elements are known. If necessary, short-periodic perturbations due to d~-, J 3 - , and J4-terms and non-zonal terms in the potential should be computed and included in 6Ep. Then, from the short-periodic perturbations in the orbital elements, 6Ep, evaluated for each epoch of observations deviations in the coordinates of the satellites due to 6Ep are computed and subtracted from the observed positions and/or range-rate. By doing so the observed positions of the satellites can be transformed to positions of a fictitious satellite moving on a mean orbit which does not change rapidly owing to the short-periodic perturbations. Orbital elements of this mean orbit, /~= E - f E p , are still changing slowly because of secular and long-periodic perturbations, and can be expressed as polynomials of time during a couple of days after some lunar and non-zonal long-periodic terms, which can be computed with sufficient accuracy, have been subtracted from the observations if necessary. In the expressions of the mean orbital elements, F,, k
= Y~ Ej. t j,
(15)
j=0
it sometimes happens that not only the constant term Eo but also EI, Ez,... should be corrected. For this case the correction vector dE in (13) should be written as, k
dE = ~ dEj. t J', j=0
(16)
826
YOSHIHIDE KOZAI
where dE now consists of five components, do, dr2, di, de, and dl instead of six, as by the following relations d~c and dn are expressed by d/j; d~c=d/o,
dn=dla,
2dnl =dl2,....
The constant term Eo in (15) is the sum o f E ~ and (6E~ + 6Ep) at t = 0 . Now the equation (13) is written in the following form; cos ~ dc~ / d6 dp/p
k -- 2 (A~'tJ)'dEi,
(17)
j =o
d~/p/ where Aj takes an identical form for any value o f j except for the coefficients of dlj, for which contributions from dn through da should be added except for j = 0 . The values of k on the right-hand side of (17) may be different for different elements. On the left-hand side of (17) (O-C)'s based on approximate orbital elements, E, which depend on time, are given. From the equation (17) not only the corrections for the mean orbital elements at t = 0 , dE o, but also those for the secular motions including the mean motion, dE a, and for the acceleration, dE 2, and so on can be solved, if a sufficient number of observations are well distributed for a certain period of time. The secular motions thus derived are not always accurate enough, except for the mean motion. However, they were sometimes used for determining coefficients of zonal spherical harmonics of even orders in the gravitational potential of the earth. Thus the improved mean orbital elements corresponding to the mean epoch, t = 0, are derived from the observations, say, for four days. When a sufficient number of observations for any one satellite is distributed for a long interval of time, say, for more than three months, sets of the mean orbital elements, say, every two days are derived, and values of secular motions and amplitudes of long-periodic perturbations can be determined accurately by using these data. To evaluate these, known long-periodic perturbations such as due to the luni-solar gravitational attractions and to the solar radiation pressure, should be computed for each epoch and for each element, and should be subtracted from the mean elements, E. Also, long-periodic perturbations with arguments 2on, 3co, ..., due to zonal spherical harmonics can be computed with sufficient accuracy, and can be eliminated, since their amplitudes are small. Thus we have the mean orbital elements of the second stage, E, for each epoch. Time variations of E thus derived are produced mainly by secular perturbations due to zonal spherical harmonics of even orders, due to luni-solar forces, due to air-drag and by long-periodic perturbations with argument co due to zonal harmonics of odd orders. However, the variations in the mean motion and the semi-major axis are produced only by the air-drag, unless long-periodic variations are produced by resonance with non-zonal terms in the potential of the earth. Therefore, from these variations the air-drag effects can be estimated by comparing values of the mean motion and the
827
THE EARTH GRAVITATIONAL POTENTIAL
semi-major axis at each epoch with those, n o and ao, at the mean epoch t = 0 as
Aa
=, a - aoo[)3
(is)
The drag-effects in the mean anomaly, l, can also be estimated by comparing observed values of I with the computed values based on the mean motion n o and the mean anomaly, ~c, at t = 0 as, A1 = I - (not + to), (19) if long-periodic perturbations with argument co in the mean anomaly can be eliminated approximately. Between Al and An the following relation should hold: t
Al = j a n
dt.
,J 0
The mean height of the perigee, q = a ( 1 - e ) , and the mean inclination are also changing secularly owing to drag. However, they vary more slowly than the semimajor axis, the mean motion and the eccentricity do, and they can be expressed approximately as linear functions of An as Aq = b "Aa = Ai = c . A n / n ,
z3 b" a "An/n ,]
(20)
with constant coefficients b and c, which can be computed theoretically if parameters concerning the atmosphere are known. The value of b is very small if the perigee is high enough and the eccentricity is not very small. The variation of the inclination, Ai, is caused by the motion of the atmosphere following the rotation of the earth, and is therefore small (STERNE, 1959; COOK and PHMMER, 1960). The rate of the secular variation due to the drag for the eccentricity is derived by the first relation in (20) as, Ae = (1 - e - b).Aa/a = -
2(1 - e -
b) An~n,
(21)
The secular variations Aa, Ae, and and Ai cause secondary secular variations in the argument of perigee and the longitude of the ascending node as, 1col ACO=3 n
7- e
---{lYe +
1Ol{~-e - - AY2=3 n +e
+
8eb 1-e 2
8eb - - - 3 c 1-e 2
15csin2i) 4~ } t
a
n
i
(22) AI,
where col and f21 are secular motions for co and f2, respectively. When the density and the scale-height of the atmosphere at the perigee are constant, that is, independent of time, Aa and An, and, therefore, Aq and Ae, are proportional to t, and Al, Aco, and AY2 are proportional to t 2.
828
YOSHIHIDE KOZAI
Of course, as formulae (20) through (22) are approximate linear expressions, they cannot hold for large values of An and At. After subtracting the air-drag effects from the mean orbital elements of the second stage, E, the mean orbital elements of the third stage are derived. These orbital elements should be expressed as, a = ao~
/'/~
t/o,
i = i o + A ~ sin co, e = eo + Ae sin co,
(23)
co = coo + COlt + A~ cos CO, f2 = f2o + O~t + A~ cos co. The expression for the mean anomaly is not given here, since it cannot be used for determining coefficients J, in the potential of the earth because of large drag influences. When the eccentricity takes a very small value, the amplitude of cos co in the expression for co in (23), A~,, becomes very large because of a small divisor l/e, and the expressions Ae and A~oin (23) do not hold in this form. For this case the eccentricity and the argument of perigee are better expressed as e sin co = eo sin (coo + tort) + Ae ,~ e cos co e0 cos (coo +cott)J
(24)
The final mean elements, io, e0, coo, and f2o at the mean epoch t = 0 , the secular motions, cos and f2,, and the amplitudes of the long-periodic terms with argument co, Ai, Ae, A~ and A9 can be determined by the method of least squares by the equations (23). However, when the eccentricity and the argument of perigee are expressed by (24), approximate values of A~ and co~ must be known to determine values of coo, e3l, eo, and A~. They can be determined by the following equations. By using the equations e sin co = eo sin coo (cos COl"t)+ e o cosco o (sin (~l't) ] e cos co
+ (e cos co) t dcoI + A~, t) I eo cos COo(cos Nl"t) - eo sin ~oo (sin c~i - (e sin CO- -4e) t dcoa ,
(25)
where c~I and Ae are approximate values of co, and A~, four unknowns eo sin coo, eo cos COo,Ae, and dco1, the correction for c9s are determined. Thus we have the mean orbital elements free from any kind of perturbations. However, some of these values might be different when different theories of perturbations due to Jz-term in the potential are applied in the reduction of observations even if the same observations are used. For example, when M~RSON'S theory (1961) is adopted to compute periodic perturbations, smoothed elements which are modified orbital elements corresponding to L = 0 , that is, at the ascending node, are derived instead of the mean elements.
829
THE EARTH GRAVITATIONAL POTENTIAL
Merson's smoothed inclination is the inclination in a modified sense when the satellite passes through the ascending node. This value is variable and is a periodic function with argument co. Therefore, the orbital elements based on Merson's theory cannot be reduced by the same way as explained here. The inclination in K~NG-HELE'S theory (1958) is defined as the highest latitude attained by the satellite and is constant. However, it is different from the mean inclination. When BROUWER'S (1959) or KOzAFs (1959b) theory is applied, mean orbital elements are derived by the method described here. However, both theories do not give identical values for some of the orbital elements. The mean orbital elements of the first stage which are free from the short-periodic perturbations and which are derived directly from the differential orbit improvement program differ, as outputs, from each other by long-periodic terms, having J2 as a factor with argument 2e). This means that BROUWER and KozAI gave different expressions for the long-periodic perturbations with argument 2~o, and the differences are absorbed in the expressions of short-periodic perturbations. Therefore, after the periodic terms with argument 2m given in their papers have been subtracted, namely, when the mean orbital elements of the second stage have been derived, both theories will give identical values for the orbital elements, except for the semi-major axis. That is, the same values for the constants appearing in the formulae (23) and (24) except for ao are derived from the same set of observations. Brouwer's mean semi-major axis, aB, is the mean value of the expression for the semi-major axis with respect to the mean anomaly. On the other hand Kozai's mean semi-major axis, arc, is nearly equal to the mean value of the geocentric distance of the satellite. This has been a conventional definition of the semi-major axis adopted for natural satellites. Between a B and arc the following relation holds: J2 (1 - 3 sin 2 i)~/1- ~ } arc = ar~ 1 - ~3 -p~ For both Brouwer's and Kozai's theories Kepler's third law does not hold in the form (1), but, for example, for Kozai's theory Kepler's third law is written as,
23
noa o = G M
{13- ~J( 21 - ~ - s i n
2 i)x/1-e
.
(26)
Comparison between BROUWER'S elements (1959) and M[USEN's elements (1959) based on Hansen's method and adopted at the U.S. National Aeronautics and Space Agency has been made by BAILIEand BRYANT (1960). Also BATRAKOV'Stheory (1963) and Zhongolovitch's theory have been compared by Brouwer's theory (ZHONGOLOVITCH, 1960; ZHONGOLOVITCHand PELLINEN,1962). Of course, differences between any sets of elements are not larger than a quantity of the order of J2.
830
YOSHIHIDE KOZAI
4. E x p r e s s i o n s for the S e c u l a r M o t i o n s
Of the quantities determined with equations (23) and (25) f r o m observations, the secular motions, co1 and O1, are expressed theoretically by the following forms:
(-~
2
=
4
J2jG2j/P 2j + J2 G2, 2/Po +
A~LS,
j=1
f~,I(-noOo) =
Y', j=l
(27)
JzjHzj/p 2j + ,12112, 2/p4o + Af2LS,
where the semi-major axis is expressed in units of equatorial radius of the earth and Po = ao (1 - eoz),~ 0 o = cos i o . Here
Gzj and
Hz: are functions of e 2 and 02 and are given here for j = 1, 2 . . . . , 7 ;
3 U2 = ~, //4 =
15 ~(3-702
G4 = H6 -
3 G2=--~(1--502), ) (2+3e2),
15 128 {4 (3 - 3602 + 4904 ) + 9e 2(1 - 1402 + 2104)},
105
(5 - 3002 + 3304 ) (8 + 40e 2 + 15e4),
1024 105
G6 = 02H6 - ~ /48 -
315 32768
G8 = 02H8
(5 - 10502 -~- 3150" - 23106).(8 + 20e 2 + 5e4),
(35 - 38502 + 100104 - 71506)'(16 + 168e 2 + 210e 4 + 35e6), 315 262144
(35 - 126002 + 69300'* - 1201206 + 643508) 9
9(64 + 336e 2 + 280e 4 + 35e6), 3465 Hi~ - 4194304 (63 - 109202 + 491404 - 795606 + 419908) (128 + 2304e 2 + 6048e 4 + 3360e 6 + 315e8),
Glo = 02Hlo
3465 8388608(63 - 346502 + 300300'* - 9009006 + 10939508
- 4618901~
+ 1152e z + 2016e 4 + 840e 6 + 63e8),
9009 HIE -- 67108864(231 -- 577502 + 3927004 -- 10659006 + 12435508 5200301o)'(256 + 7040e z + 31680e 4 + 36960e 6 + 11550e 8 + 693e1~
THE EARTH GRAVITATIONAL POTENTIAL
831
9009
G,2 = 02H12 268435456(231 - 1801802 + 22522504 - 102102006 + 20785050 s - 19399380 ~~ + 6760390 I2) (512 + 7040e 2 + 21120e 4 + 18480e 6 + 4620e s + 231ei~ 45 045 Hi4 = 2147483648(429 - 1458602. + 1 3 8 5 6 7 0 * - 55426806 + 10623470 s - 9657700 ~~ + 334305012) (1024 +
39936ez + 274560e 4 + 549120e 6
+ 360360e 8 + 72072e t~ + 3003e12), 45 045
GI4 = 02H14 4294967296( 429 - 4504502 + 76576504 - 484984506 + 145495350 s - 223092870 I~ + 16900975012 - 5014575014) (1024 + 19968e z + 91520e 4 + 137280e 6 + 72072e s + 120t2e I~ + 429e12). Here and in the following the suffix 0 to be attached to e and 0 has been omitted for simplicity. Secutar motions due to the lunar and solar attractions are denoted by AcoLs and AfaLS in (27) and are expressed as (KozAL 1959a), AmES = 89 (-- 13+2 502 + AOLs ~ ( l + ~ - e ) ,
e2)'}
(28)
where 3
1
4 n Z x / 1 - e~ In 2 (1 - ~}sinZ~) +
n~m} {1 - 88sinZJ (t + cos%) - ~ sin2e cos / a}].
In this expression n o and n~ are mean motions of the sun and the moon, respectively, m~ is the mass of the moon in units of the mass of the earth, ~ is the obliquity, that is, the inclination between the ecliptic and the equator, and J is the inclination of the lunar orbit with reference to the ecliptic and is about 5 ~ When n is given in degrees per day, ~, the dimensionless quantity, takes the following value: t.762
!g
nZ ~ / 1 - e2.
(29)
Functions G2, 2 and//2, 2 in (27) also depend on ez and 02. As different definitions are taken for io, eo, and ao, forms for G2, 2 and H2, 2 by different authors are different from each other. Therefore, definitions of the orbital elements should be found to compute their values. Expressions for G2, 2 a n d / / 2 , 2 by some authors are given as follows:
832
YOSHIHIDE
H2,2 KING-HELE
(1958)
BROUWER (1959) MERSON
(1961)
3
(7-1902)+0(e2), 3
-
KOZAI
32
{4(1-1002)-(9-502)eZ+12(1-302)x/1-e~},
3 32 {4 (3 - 1402) + (25 - 5302) e2},
G2, 2
3 {2(1 -502) (5-t-4302)+(25 - 12602 128 + 1450'*) e 2 - 2 4 ( 1 - 5 0 2 ) (1-302) x/1-e2},
BROUWER (1959)
-
MERSON (1961)
-
3 128
{200- 72402 + 5850'* + (61 - 44602 + 5850'*) e2}.
KozAI'S expressions (1959b) for G2, 2 and/-/2, 2 are derived by changing the numerical coefficients of , f ( 1 - e 2) in the corresponding expressions by Brouwer from + 12 to - 2 4 and from - 2 4 to +48, respectively. Values of o92 and f2t in degrees per day for n = 14 revolutions per day (a= 1.14)
degrees per day 15"
I0~
5~
-~~
- I 0 ~ 0I o
J---~,
' ~ t5
2e
~~
i 6o ~
i 75 ~
'
Inclination
Fig. 5.
Daily secular motions for the argument of perigee and the longitude of the ascending node as functions of the inclination for n = 14 revolutions per day and e = 0.1.
THE EARTH GRAVITATIONAL POTENTIAL
833
and e = 0.1 as functions of the inclination are shown in Figure 5. If the motion of the satellite is retrograde, the secular motions also change their signs. The value of cot vanishes for any plausible value of n near i = 63 ~ which is called the critical inclination. For i < 6 3 ~ col is positive, whereas for i > 6 3 ~ cot is negative. And for i = 0 ~ and e=0.0, that is, for the equatorial circular satellites col is twice as large as the value of - f21Approximate values of cot and ~21 for other values of the mean motion are derived by multiplying a factor (n/14) 7/3 by the values given in Figure 5. The last two terms in the expressions (27), that is, the luni-solar terms and the J~-terms, take small values for usual satellites, and can be computed with sufficient accuracy with the approximate value of J2. And they can be subtracted from the observed values of the secular motions. The values for col and ~21 thus corrected can be expressed by linear forms of Jzj with known coefficients, and by using these equations values of Jzi can be determined if a sufficient number of equations are available. Of Jzj, J2 takes a value of the order of 10 -3, and other values are of the order of 10 -6 or less. Therefore, it is necessary to compute the coefficients of .12 in (27) with six significant figures and those of other Jzj with three figures to compute values of 0) 1 and s with six figures. This can be done only when the values of eo, i0 and ao are determined with more than six significant figures. When many observations are distributed for a very long interval of time for satellites with large perigee height, that is, with small air-drag effects, values of cot and s a can be determined from the observations with higher accuracy than the accuracy for eo and io. For this case theoretical values of col and ~2t cannot be computed as accurately as the observed values. This means that even if observations for a very long interval of time are available the accuracy for determining J2~ cannot be increased so much. To eliminate correlations between various coefficients in (23) and (25) for determining cot, f2t and the amplitudes of the periodic terms with argument co, observations covering one revolution period of the argument of perigee should be used. It need not be necessary to use observations for more than this period of time. When we try to express values of col and ~21 theoretically with more than seven significant figures, many more coupling terms such as J~, JzJ4, JZ/J2, JZ/J z should be added to the right-hand sides of (27) besides J2Z-terms (KozM, 1962) and more accurate values of ao, io, and eo should be known. In order to derive the mean orbital elements with such a high accuracy short-periodic perturbations with factors j2, J3 and J4 should be computed and eliminated in the orbit improvement program. Otherwise, additional terms would not be necessary in the formulae (27). When we know all the forces acting on the satellites, and they can be evaluated accurately, one set of data, col, ~?t, ao, i0, and eo, are enough for one satellite. However, when a couple of sets of data are derived from observations during different periods, it usually happens that they are not consistent with each other. Therefore, as many sets of data as possible are needed even for one satellite to eliminate accidental or possibly systematic errors in the determinations of the secular motions and the mean elements.
834
YOSHIHIDE KOZAI
If there are errors Ae, Ai, and An in eo, it, and n o, the computed value of (21, for example, has errors expressed as,
As
4eo.Ae
7 An
tan i. A i + - - 3 n
-
f21
1 - e0z
(30)
Besides these errors, if there is an error in the value of GM, the geocentric gravitation constant, the value of the semi-major axis which is computed from the mean motion by Kepler's third law such as (26) has the following error:
Aa a
1 A(GM)
3 GM '
and, therefore, - ~A (GM)/GM should be added to the right-hand side of (30). When satellites with low perigee height must be used, that is, when air-drag effects are rather large, procedures explained in the previous section cannot be followed for determining the values of 031 and f21. However, if approximate values of dj are known, mean orbital elements at each epoch can be compared with computed values derived by integrating disturbing functions containing only secular and long-periodic terms which are functions of the mean orbital elements at each epoch. Then corrections de) 1 and dr21 based on the assumed values of dj can be derived, and the)' are expressed as linear functions of corrections to Jj as:
d03t = ~ P2j'dd:j, ) j=l
(31)
dr21 = ~ Q2J'dJ2j" i j=l / As values of dJ2~ are estimated to be in the order of l0 -6, it is not necessary to compute P2i and Q2j very accurately even'for j = 1, that is, it is not necessary to put very- accurate values of ao, io and e o in (31). Actually, a o and eo vary rather rapidly when the air-drag effects are large, and to compute values for P2j and Q2) in the formulae (31) approximate mean values can be used. MERSON'S theory (1961) does not express secular motions and long-periodic variations for the argument of perigee and the longitude of the ascending node in the integrated form (23), but gives instantaneous rates of changes for the smoothed elements as functions of other smoothed elements which are not constant. The instantaneous rates of changes for 03 and f2 can be derived directly by the orbit improvement program, although their accuracies are not so high as those derived by use of equations (23). As CooK's theory (1962) of the luni-solar perturbations gives also instantaneous rates of changes of elements, by applying his theory corrections due to the luni-solar effects can be made. By inserting approximate values for Jj with odd values of j, equations of condition to determine values for Jzj are derived. When k equations of condition are given, values of J2j up to j = k can be derived theoretically. Although, to avoid ill-conditioning in solving these equations data for various satellites with various sets of orbital elements should be used, ranges for
THE EARTH GRAVITATIONAL POTENTIAL
835
values of inclinations, eccentricities a n d semi-major axes for the satellites to be used for this purpose are limited. F o r the satellites for which very accurate orbital elements are available the eccentricities are n o t larger t h a n 0.25, the inclinations are higher t h a n 28 ~ a n d the semi-major axes are between 1.15 a n d 1.7 earth radii. Also there are gaps in the i n c l i n a t i o n between 50 ~ a n d 67 ~ a n d between 70 ~ a n d 85 ~. Of course, observations for very low satellites are available. However, as air-drag effects are very large for such satellites, their data are n o t adequate for d e t e r m i n i n g values of J2j. Even if perigee height is n o t low, satellites, which have large area-to-mass ratios (A/Ms), such as Echo satellites, are inadequate, since for these satellites effects of the solar r a d i a t i o n pressure are too large to be evaluated theoretically with sufficient accuracy, a n d the air-drag effects are very large. As positions of perigee are very difficult to determine from observations for satellites with nearly circular orbit, accuracies for the secular m o t i o n s of the arguments of perigee are very poor for such satellites. Also, effects by solar radiation pressure on the a r g u m e n t of perigee c a n n o t be c o m p u t e d as accurately as the a r g u m e n t of perigee is determined from observations, especially, for such nearly circular satellites because of a small divisor e appearing in the expression of the perturbations.
5. Z o n a l Harmonics of Even Orders I n Table I coefficients derived recently by four authors for zonal spherical h a r m o n i c s of even orders are given. Satellites used by various authors for d e t e r m i n i n g coefficients of n o t only zonal h a r m o n i c s of even orders b u t also those of odd orders a n d n o n - z o n a l harmonics are listed with approximate orbital elements in Table II.
TABLE I Coefficients of Zonal Harmonics of Even Orders (in units of 10-6) J2 COOK (1965)
J4
1082.65
1.61
4- 1 0
I
KOZAI (1964)
4- 1 5
J8
1.52 3
0.57 4-7
0.44 4-11
2
1082.68
-- 1.61
0.71
0.13
1
1082.63 4-1
-- 1.63 4-2
•
0.59 3
-- 0.15 4-3
1082.64 4-1
-- 1.65 4-2
0.65 • 3
-- 0.27 • 5
1.70
0.73 4-40
-- 0.46 •
1082.64 • 8
•
YlZ
J14
0.73 40
--
•
J10
4-
1082.64 4-2
1
2
SMrrH (1965)
JG
0.09
-- 0.31
--
0.16 4-5
-- 0.29 4-5
--
0.05 4-5
-- 0.36 • 4
0.18 4- 6
0.17 • 29
-- 0.22 • 10
0.19 4-11
836
YOSHIHIDEKOZAI TABLE ff Satellites used for Geodesy Zonal harmonics of even orders
61 v 64-5A 61 at/1 59 a 1 59 r/ 58fl2 61 fi 1 62 a e 60 t 2 62fly 60~r2 62fl/z 1 64-15A 62o 1 63-24A 61 o 1 61o2
Explorer 11 Saturn5 Transit 4B Vanguard2 Vanguard 3 Vanguard 1 Explorer 9 Telstar Echo 1 rocket Relay 1 Tiros2 Anna 1B Ariel2 Ariel 1 Tiros7 Transit 4A Greb3Injun 1 Alouette Transit 5B
62flc~ 1 63 38C 63 49B 61 c~fi1 Midas4
i
e
a
28~ 31.5 32.4 32.8 33.4 34.2 38.8 44.8
0.086 0.032 0.010 0.165 0.189 0.190 0.105 0.242
1.18 1.08 1.16 1.30 1.33 1.36 1.25 1.52
47.2 47.5 48.5 50.1 51.6 53.8 58.2 66.8
0.011 0.284 0.008 0.007 0.073 0.057 0.002 0.008
1.25 1.69 1.11 1.18 1.13 1.13 1.10 1.15
66.8 80.5 89.9 90.0 95.9
0.008 0.003 0.004 0.004 0.012
1.15 1.16 1.17 1.17 1.57
C
KC
S
x
Those of odd orders
KO
GN KCS S
x
x
Non-zonal harmonics GN
I
KA
X
X
X
X
x x x
x
X x
• • X x
• x
•
x x
x • •
• •
•
X
x
x
x
x
X
•
• • •
•
• • •
•
•
•
x x
• •
•
•
x
•
x
C: COOK,KC: KING-HELEand CooK, S: SMITH,KO: KozAI, GN: GUIERand NEWTON, KCS: KING-HELE,COOKand SCOTT,~: [ZSAK, KA: KAULA
The values given in Table I were derived by the secular m o t i o n s for the longitude of the ascending node only except for values by KOZAI (1964). Table I shows that the values derived by the various authors are consistent with each other except for values for Js. COOK (1965) used four rather high satellites with inclinations 34~ 44~ 660.8 a n d 950.9. By using only high satellites, for which effects of higher order harmonics are negligible, Cook intended to separate the effects of lower order harmonics to those of high orders. Therefore, he derived three coefficients, J2 t h r o u g h J6, from four data. Seven satellites with inclinations 280.8, 320.8, 38~ 47~ 530.8, 66~ a n d 950.9 were employed by KrNG-HELE a n d COOK (1965) in their determination. The orbital elements for these satellites were derived from Baker-Numa a n d M i n i t r a c k observations. Values of f2~, io, eo, a n d ao were n o t determined by the m e t h o d of least squares, b u t were derived from a few data by assuming that the atmospheric density at the perigee is constant. For, from values of ~2 corresponding to ~o= 90 ~ a n d 270 ~ where the longperiodic term with a r g u m e n t co vanishes, the value of f2~ was derived, a n d from the values of n, e, a n d i corresponding to co = 180 ~ n o, %, a n d io were determined.
THE EARTH GRAVITATIONAL POTENTIAL
837
The average relative accuracy of Q1 adopted by King-Hele and Cook is of order of 10 - 4 , t h a t is, the average error in ~ i is 3~ 10 -4 per day. They divided each equation of condition by -noOo/p ~ so that the coefficient of J2 is always 1.5. They selected four equations out of the seven equations and solved them with four unknowns, J2, J4, J6, and Js, on the assumption that J10, J12, ... were zero. There are twenty-four possible sets of four equations. However, some sets of equations had to be excluded, because the inclinations of the satellites used in them are not widely distributed or two of the four equations have large standard deviations. Of twenty-four sets of solutions eight solutions had standard deviations less than 0.10x 10 -6, 0.20 x 10 - 6 , and 0.40x 10 - 6 , for "]2, J4, J6, and Js respectively. As King-Hele and Cook found that the scatter of these eight solutions was extremely small, the solution 1 by them in Table I was derived by taking the mean values of the eight solutions. They thought that the solution 1 was the most satisfactory representation of the potential. The solution 2 with six unknowns was derived by solving the equations by minimizing the largest residuals. This method is different from the method of least squares. The maximum residual based on these solutions is 0.68 times as large as the standard deviation assigned to the corresponding equation of condition. The solution derived by the method of least squares with six unknowns was, J2 = 1082.61 x 10 - 6 , _ 12
"]4 = -
1.46 x 10 - 6 , -4- 30
J6 =
0.44 • 10 .6 , + 50
J8 =
0.53 • 10 -6 , + 80
J10 =
- 0.23 x 10 - 6 , _+ 65
J12 =
0.06 • 10 - 6 . _+ 66
(32)
The sum of the squares of the residuals for the solution (32) is 0.81, whereas that for the solution 2 is 1.30. However, the maximum residual based on the solution (32) is 0.82 times as large as the standard deviation of the corresponding data. King-Hele and Cook preferred the solution 2 to the solution derived by the method of least squares. KozAI (1964) used secular motions of both the argument of perigee and the longitude of the ascending node for nine satellites (1961 ol and 02 being regarded as one satellite) with inclinations 28~ 32~ 330.4, 44~ 47~ 470.5, 500.1,660.8 and 95~ The inclinations for these satellites are not distributed so uniformly as those for the satellites used by the other authors. Observations used were precisely reduced BakerN u n n observations, except for 1961v with inclination 280.8, for which Baker-Nunn observations reduced at each camera station were used. In order to eliminate accidental and systematic errors in the data, as many sets of data as possible were derived even for one satellite at different epochs, and ( O - C ) based on approximate vMues of Jz~ was computed for each value of the secular motion. It was found that the scattering of the ( O - C ) was sometimes much larger
838
YOSHIHIDE KOZAI
than their assigned standard deviations, and by taking their mean value for each satellite Kozai thought that he derived reliable data. Since 1961 ol and 1961o2 with inclination 66~ have almost identical orbital elements, these two satellites were regarded as one, and produced only one set of data. The standard deviations finally assigned to the secular motions for the ascending node were on an average, 10-5 degrees per day, that is, of order of 10-6. The secular motions for the argument of perigee could not be determined accurately for satellites whose eccentricities are very small. For nearly circular orbits perturbations due to the solar radiation pressure in the argument of perigee cannot be evaluated as accurately as the argument of perigee can be determined from the most precise observations. Therefore, the standard deviations assigned to the arguments of perigee take various values, depending on the eccentricities of the satellites. For the nine satellites the standard deviations for the perigee take values between 2 ~ x 10 -5 and 1~ x 10 -3 per day. Kozai constructed eighteen equations of condition in the forms (31) which express (O - C) based on the approximate values of Jz~ as linear functions of the corrections, dJzi, to J2j. Then he solved the equations with six and seven unknowns by the method of least squares after assigning a weight reciprocally proportional to the square of the standard deviation to each equation. The weights thus assigned to the equations of condition are therefore widely spread over the region between 2 and 5 x 10 -4, and actually one third of the equations were neglected in the determination. The results were given in Table I. Kozai's weighting system is a conventional one, when the standard deviation assigned to each equation expresses the error in the equation. However, as values of J16, J18 . . . . had to be neglected to solve the equations, errors caused by neglecting these coefficients should be added to the standard deviations computed from those estimated for con, f2~, eo, io, and n 0. Since these errors due to neglecting higher order terms could not be estimated in any way and were not taken into account by Kozai, his standard deviations cannot be regarded as the real errors in the equations. This might cause systematic errors in the determination of Jzj, since the inclinations of the satellites Kozai used were not uniformly distributed and, therefore, effects due to neglecting J~6, J18 . . . . could not be cancelled out. SMITH (1965) used seven satellites with inclinations 31~ 38O.8, 440.8, 51~ 58~ 66~ and 80O.5 and solved seven simultaneous linear algebraic equations with seven coefficients as unknowns. Standard deviations assigned to the observed values of the secular motions are 2~ 10 -5 per day, or of order of 10 -5 on average. 6. Z o n a l H a r m o n i c s o f O d d Orders
Coefficients of zonal spherical harmonics of odd orders in the potential of the earth can be evaluated by using the amplitudes of the long-periodic perturbations with argument co, namely, Ae, Ai, A~o and Ae determined by Equations (23) or (25) from observed orbital elements. The amplitudes are expressed approximately as linear
THE EARTH GRAVITATIONAL POTENTIAL
functions
839
J2j+ 1 a s ,
of
oo
sin i Ae = - J2a (1 - -5O 2)
J2j+ 1 B2j+ 1 P2j+ 1/P2(j- 1 ) , j=0
e0 At
V/1 - e i s i n i e0
J2j+l { - sin2 i (1 - 502)'C2~+1
An = J2P sin i(1 - 502) 2
(33)
j=0
+ (9 - 502 ) B2j+I} P2j+~/P2~ sin i Ao, = - OA, - eJ2p(1 ~ 502)
J 2 j + l B z j + I Qaj+I/P 2(j-l) , j=0
where the semi-major axis is expressed in units of the equatorial radius of the earth and the suffix 0 to be assigned to 0, i, e and p have been omitted. Coefficients B2j+ 1 and C2j+ t in (33) are polynomials of 02 and are not divided by 0 and sin i, and P2j+I and Qzj+ 1 are polynomials of e z not divided by e 2. A l t h o u g h B3 and C3 have factors ( 1 - 5 0 2 ) and ( 1 - 5 0 2 ) 2, respectively, for other values of j B z j + I a n d C2j+I do n o t have such factors. Therefore, the amplitudes become very large w h e n 1 - 502 vanishes, that is, at the critical inclination. These linear expressions (33) do n o t hold near the critical inclination, neither for very small eccentricity cases and the nearly equatorial orbit because of divisor e and sin i. The coefficients in (33) take the following form" B 3
=
] -- 502 ,
P3
= Q3
Bs
=1-
P5
5 = 32 (4 + 3e2),
C 3 =
10 (1--502)
2 ,
1
2 1402+2104 ,
c , = 28 (1 - 302),
5 Qs = 3 2 ( 4 + 2 5 e
2 + 6 e 4)
B 7 = 5 - 13502 + 49504 - 42906 ,
C7 = 2 (135 - 9900 z + 128704), P7
-
35 2O48 35
Q7 --
- -
2048
(8 + 20e 2 + 5e4),
(8 + 124e 2 + 145e 4 + 20e6),
B9 = 7 - 30802 + 200204 - 400406 + 243108 ,
840
YOSHIHIDE K O Z A I
88 (7 -- 9102 + 27304 -- 22106),
C 9
=
e9
105 ~ - --(64 65536
Q9
-
105
+ 336e 2 + 280e 4 + 35e6),
(64 + 1776e 2 + 4760e 4 + 2485e 6 + 210eS),
65536 B~I = 21 - 136502 + 1365004 - 4641006 + 6298508 - 2 9 3 9 3 0 l~ C l l = 130 (21 - 42002 + 214204 - 387606 + 2 2 6 1 0 8 ) , 1155 P1, - 4 1 9 4 3 0 4 (128 + 1152e 2 + 2016e 4 + 840e 6 + 63e8), 1155 Qll
-
B13 =
(128 + 5504e 2 + 2 6 2 0 8 e 4 + 3 0 0 7 2 e 6 + 8967e s + 504e1~
4194304
33 - 297002 + 4207504 - 213 18006 + 47965508 -
4 9 0 3 1 4 0 *~ + 185725012 ,
C13 = 60 (99 - 280502 + 2131804 - 6395406 + 8171908 - 3 7 1 4 5 0 1 ~ 3003 P13
--
QI3
--
(512 + 7040e 2 + 2 1 1 2 0 e 4 + 18480e 6 + 4620e 8 + 231e1~
67108864 3003
(512 + 3 1 3 6 0 e 2 + 2 3 2 3 2 0 e 4 + 4 6 7 2 8 0 e 6 + 300300e 8
67108 864 + 5 7 9 8 2 e 1~ + 2 3 1 0 f 2 ) .
0.002,
,
,
,
O.O01
O.O00
d
-
0.0010
o
15o
I
30 ~
I
45~
f
60 ~
l
75 ~
90 ~
lnclinotion
Fig. 6. The amplitude, Ae, for a = 1.14 and e = 0.1 as functions of the inclination.
THEEARTHGRAVITATIONALPOTENTIAL
841
The amplitudes in (33) are of the order of 10 -3, as J2j+l are all of the order of 10 . 6 and the divisor Jz in these expressions is of the order of 10-3, T h e amplitudes Ae, A~o, and An based on SMITH'S value o f J2j+l (1963) for n-----14 revolutions per day ( a = 1.14) and e = 0 . 1 are s h o w n as functions of the inclination in Figures 6 and 7. The amplitude At can be c o m p u t e d by multiplying - e 0 / ~ / i l - e 2) sin i to As. As was m e n t i o n e d earlier these curves do not represent correct amplitudes near the critical inclination, although it is true that the amplitudes become very large near the critical
A~u
Af~
r
J~o
c,~5
f2
~
/
o~o -<5
o?os
~oo
i
9o~o
5
o~lo
- i~o 0~
15"
30 ~
45 ~
60 ~
7~, ~
"90 ~
Inclinotion
Fig. 7. The amplitudes, Ao~ and An, for a --~ 1.14 and e = 0.1 as functions of the inclination. Scales given on the left side correspond to A~o and those on the right side to An.
inclination. The amplitudes, As~ and - A ~ , take also very large values near i = 0 ~ because of a divisor sin i, although An + Ao is very small. It should be noticed that Ao,, A~, a n d A ~ have factors l/e, e0/sin i, and e0, respectively. The equations o f condition to determine the values of J2j + 1 are expressed as linear functions, and it is sufficient to c o m p u t e the coefficients of J2j + 1 with three significant figures. T h e values of Jzj+ 1 recently obtained by various authors are given in Table III. GOIER and NEWTON (1965) used four satellites with inclinations 320.4, 500.1, 66~ a n d 900.0. These are Transit and Geodetic satellites, and their eccentricities are all less t h a n 0.0l. T h e orbital elements were derived f r o m D o p p l e r observations by using
842
YOSHIHIDEKOZAI
TABLE III Coefficients of Zonal Harmonics of Odd Orders (in units of 10-6) J3 GUIER and NEWTON
--
(1965) KING-HELE, COOK
J9
2.68 •
-- 0.03 •
-- 0.59 •
0.18 • 2
--
2.56 • 9
-- 0.15 -- 16
-- 0.44 • 30
0.12 • 28
2
-- 2.73
0.17
-- 0.94
0.53
1
--
2.56 •
-- 0.19 •
-- 0.38 +2
-- 2.55 •
-- 0.21 •
-- 0.33 •
2.44 •
-- 0.18 ~ 3
-- 0.30 • 3
2 SMITH (1963)
J7
1
and SCOTT(1965)
KOZAI (1964)
J5
--
--
0.68
0.04 • -- 0.05 •
Jla
Jll
•
0.30 4 0.30 • 4
-
0.11 •
numerical integrations to compute perturbations. As the eccentricities are very small for these satellites, An and A~ are also very small because of a factor, e, and, therefore, were not used for the determination o f J2j'+ 1. They adopted formulas (24) to determine the value of Ae. Therefore, GUIER and NEWTON'S values of A e, A / B in their notation, were derived by using both the eccentricity and argument of perigee. In determining Ae they simultaneously also derived the factor of the solar radiation pressure effects f r o m the observations. The average standard deviation of their values of A e is 10 -6. By using four data for the four satellites four unknowns J3 t h r o u g h 3"9 were solved, and their standard deviations were computed f r o m those of A~. KIN6-HELE, C o o n and SCOTT (1965) used the amplitude Ae for six satellites with inclinations 28~ 32~ 47~ 51%, 66~ and 95~ They derived values of eo and A e by the equation e = eo + Ae sin o9, even for satellites with very small eccentricities. The smallest value of the standard deviations assigned to Ae was 2 x 10 -6, and the largest one was 11 x 10 -6. They divided each equation of condition by sin i/2Jaa so that the coefficient o f Ja is always - 1. To derive values of J2i+ 1 from these equations of condition they followed similar procedures as those by KING-H~LE and C o o k (1965) to determine the values of J2j. They solved the six equations in fours for J3, -/5, -/7, and J9- Out of fifteen possible sets o f solutions thus derived nine sets have standard deviations less than 0.15 x 10-6, 0.35 x 10 -6, 0.40 x 10 -6, and 0.40 x 10 -6, for J3, Js, JT, and J9, respectively. By taking the mean of the nine values the following solution similar to the solution 1 by KING-
THE EARTH GRAVITATIONAL POTENTIAL
843
HELE and COOK (1965) for J2i was obtained: J3 = -2.61
x 10 .6
J5 =
-0.19
+_ 11 J7 = -0.56
x 10 -6
+22 x 10 -6 ,
J 9 -=
• 27
(34)
0.10 x 10 .6 .
• 27
However, King-Hele, Cook and Scott did not recommend this solution for this case, since it does not represent the value of Ae for 1961 ol with inclination 66~ Instead of the solution (34) they recommended the solution I given in Table III as four-coefficient representation. The solution 1 was obtained by the method of least squares by assigning equal weight to all the equations of condition that had been divided by sin i/2Jaa. For five-coefficient representation King-Hele, Cook, and Scott recommended solution 2 in Table III. This was obtained by minimizing the largest residual. The solution by the method of least squares with five unknowns was, J3 = - 2.78 x 10 .6 , +_ 18
Js = 0.24 • 33
J7 = -
J9 =
1.08 • 55
x 10 .6 ,
0.67 + 48
x 10 - 6 , "
x 10 -6 ,
(35)
J l l = - 0.80 • 10 .6 . _ 61 The sum of the squares of the residuals for this solution (35) is 1.46, whereas that for the solution 2 is 1.87. And the largest residual for the solution (35) is 0.80a, and that for the solution 2 is 0.62er, where a is the standard deviation of the equations of condition and is assumed to be 0.1 x 10 -6. KozM (1964) used the same satellites for this analysis as those for determining J2~', namely, the nine satellites with inclinations 28~ 32~.8, 33~ 440.8, 47~ 50~ 660.8, and 95~ For satellites with small eccentricity Ae was derived with (25) and Ai and Ara with (23). And for other satellites the four amplitudes including Ao were derived with (23). In Figures 8 through 12 observed orbital elements which are changing owing to zonal harmonics of odd orders are plotted as functions of the argument of perigee and also of time for 1959~1 and 1960z2. The eccentricities are 0.1645 and 0.0115 and the inclinations are 320.88 and 47~ for 1959el and 196012, respectively. These orbital elements were determined from precisely reduced Baker-Nunn observations. The other perturbations including the secular perturbations have already been subtracted from the observed values. As these diagrams show, Ai, Ao,, and Ae can be determined, besides Ae, with two or three significant figures for 1959el, and can be used for the determination of J2j+ IKOZAI (1964) derived many sets of values for Ae, A~,, Ai, and Aa even for one satellite by using observations during different periods of time in order to eliminate
844
YOSHIHIDE
KOZAI
e
L0
270~ 0.1650
0~
I
90 ~
I
180~
u
lille
270 ~
I
0~
i
a
t
i
e o
e 9
9
0.1640
o
e9
9
ooe
I
I
Oot.O.O
oe
I
Nov. 0.0
t
Dec.O.O
d on.O.O
1962
Fig. 8.
1963
Observed variations of the eccentricity due to zonal harmonics of odd orders for 1959ct 1.
i
32~
LO
270 ~
0~
90~
180~ 270~
0~
90 ~
l9 OO
Ool 9
I oe t
l9
32 .8 8 o
oO
Oi I
o0
i
320.870
Fig. 9.
, Nov. 0.0 Oct.O.O 1962 I
I Dec.O.O
I
dan.O.O
196.3
Observed variations of the inclination due to zonal harmonics of odd orders for 1959al.
THE EARTH
GRAVITATIONAL
845
POTENTIAL
al 270
q~
0~
~
I
90 ~
I
180 ~
270
!
I
I
OO
~
9o ~
I
9O 9
9
OOO
0~ 0 DJ
9
9o
oo
o~
eOO 9
oo
-
9149
9
0~
O
- 0.010
J
'
Oct.QO
'
Nov.O.O
a
Dec.O.O
Jan.O.O 1963
1962
Fig. 10.
Observed variations of the longitude of the ascending node due to zonal harmonics of odd orders for 1959ctl.
1.0 270 o o.~
Oo
i
n e.
90 ~ ,
27o ~
~so ~
oo
9o ~
. loegeo
~e
~n 0
~o
o
~ ~ ~9
_ 0~
, O c t 0.0 1962
Fig. 11.
, NoV. 0.0
I
De6.O.O
J an.O.O 196.5
Observed variations of the argument of perigee due to zonal harmonics of odd orders for 1959c~1.
846
YOSHIHIDE KOZAI w
o~
270~
e
90 ~
oo
270 ~
Jso ~
J
0.0120
0.0115
0,0110 i 00
9
~176149149
Oiler m
0,0i05'
|
I
I
I
Apr,O.O
Nar, O.O
MoyO.O
I
I
dulyO.O
duneO.O
Aug 0.0
1962
W 4o
270 ~
2 7 0~
180 ~
90 ~
0~
01~ 9
o 9149
QO D
O
0 0
d
2~
0
o~
9
_2 ~
9
9o ~
- ~a
I
Mar 0.0
I
Apr.O,O
I
NoyO.O
I
June 0.O
I
July O.0
I
Aug.O.O
1962
Fig. 12.
Observed variations of the eccentricity and the argument of perigee due to zonal harmonics of odd orders for 1960z2.
847
THE EARTH GRAVITATIONAL POTENTIAL
accidental and systematic errors in the amplitudes. By taking their mean value one equation of condition was obtained. Thus Kozai used thirty-three equations of condition for the nine satellites. After assigning a weight reciprocally proportional to the square o f the standard deviation to each equation o f condition, the two sets o f solutions listed in Table I I I were obtained by the m e t h o d of least squares. The standard deviations of the values for Ae were between 2 x 10 . 6 and 2 x 10- s, those for Ai were between l ~ 10 . 4 and 9 ~ 10 -4, those for A~, are 5 ~ 10 . 4 and 1~ x 10 -2, and those for A~ were 2 ~ x 10 -4 and 3 ~ x 10 -3. Kozai's m e t h o d has the same disadvantages as mentioned for the case of determination o f even order harmonics, although he used m a n y and accurate data. SMITH (1963) used Ae for three satellites with inclinations 280.8, 48~ and 660.8. The standard deviations for A e were between 8 x 10 -6 and 16 x 10 -6. By using three equations Smith determined three coefficients given in Table III. 7. Discussion of the Results The values o f coefficients o f zonal spherical harmonics derived by different authors are not completely consistent with each other, although the discrepancies have recently been greatly reduced. In order to compare these solutions with each other for J2~, deviations of the values in the two solutions by KING-HELE and COOK (1965), the two solutions by KOZAI (1964), and the solution by SMITH (1965) from the mean value o f five are c o m p u t e d for each J2j and are given in Table IV in units of 10 -9. To compute these values more figures were taken than those given in Table I if they were given in the original papers. Next, daily secular motions for the argument o f perigee and the longitude of the ascending node based on the values of 32j given in Table IV are computed for n = 14 revolutions per day ( a = 1.14) and e = 0 . 1 , and plotted as functions of the inclination in Figures 13 and 14. However, in Figure 14, instead o f O 1 itself f ~ sin i is given, where sin i is regarded as a kind of weight, namely when the inclination
TABLE IV Deviations from the mean values for J2~ (in units of 10-9) J2
J4
J6
J8
,/lo
Jl~
J14
KC 1 2
33 -- 7
12 101
61 --80
192 501
147 58
-- 73 236
--73 --74
KO1 2
--17 -- 2
-- 6 --28
--57 -- 4
-- 87 --208
-- 97 4
-- 58 --121
--74 105
S
-- 7
-- 79
80
-- 398
-- 112
16
116
KC: KING-HELEand COOK(1965), KO: KozAI (1964), S: SMrrH (1965)
848
YOSHIHIDE KOZAI
is very small, the longitude of the ascending node and its secular motion cannot be determined accurately from observations. In the upper parts of Figures 13 and 14 values of inclinations of the satellites used are marked for each author. For 1961 c~61, for which i=95~ the mark is plotted at i=84~.1. In Figures 15 and 16 coefficients of ,/2 through J14 to compute co1 and f21 sin i, respectively, are shown for the same values of n and e, namely for n = 14 and e = 0. I. Similar diagrams for other values of the mean motion, n, are roughly obtained by
doo4
~oo2
o~00o
_ o~
-~004
0~
L5~
30 ~
450
60*
75 ~
90 ~
Fig. 13. Daily secular motions of the argument of perigee, col, based on various values of J2j given in Table IV for n = 14 revolutions per day and e = 0.1 as functions of the inclination. The marks in the upper part show the inclinations of the satellites used for the determination of J2j.
multiplying a factor, (1 +4j/3)-th power of (n/14), by the coefficient for J2j. These factors are given in Table V for n = 16, 15 . . . . , 8 revolutions per day. Figure 13 shows that large discrepancies in col exist near i = 0 ~ 35 ~ 50 ~ 70 ~ and 90 ~ and Figure 14 shows that near i = 10 ~ 40 ~ 60 ~ and 80 ~ discrepancies in f21 sin i become very large. Discrepancies in cos for small inclinations are particularly large, although co1 cannot be accurately determined for these cases. However, the secular motion of the longitude of the perigee, col + f2~, can be determined rather accurately when the eccentricity is not so small, even if the inclination is very small. The values of col + f2t based on the values of J2j given in Table IV for
THE EARTH GRAVITATIONAL POTENTIAL
849
n = 14 and e = 0.1 as functions o f the inclination less than 25 ~ are shown in Figure 17, and the coefficients o f Jzj for c o m p u t i n g t h e m are shown in Figure 18. Discrepancies in co~ + (2~ c o m p u t e d by values of Jzj determined by various authors are very large for equatorial satellites.
~ c l
-_,~, I
f
oooo2oII It~o~ I 2-S.....I I
I
II
.....
f/
I
'
I
I
/
I
I
1
~176176176 / ~176176176176 / ,/
OOOCO5 I/
/~,
/
~L
/ 9
doooc
,,
",,
/
~i',
>--, ', //TW~
~:Y/
"-/?><~_\i' I
,"
X\, //,,,\ /,\~ ,,'\x. ~' ',kJ~t
( /
,:
"~"
",
/
-EO00~
",...;
,,:
,,"
/' /'~ k
/
/
lit kk
/
/
-o~
,I 0~
4'
15 ~
I 50 ~
f
I
I
45 ~
60 o
75 ~
90 ~
Fig. 14. Daily secular motions of the longitude of the ascending node times sin i, Ot sin i, based on various values of J~; in Table IV for n ~ 14 revolutions per day and e = 0.1 as functions of the inclination. The marks in the upper part show the inclinations of the satellites used for the determination of J2;.
Figure 18 shows that the coefficients o f J2j in col +f21 increase for n = 1 4 revolutions a s j increases near i = 0 ~ Therefore, higher order terms cannot be neglected to c o m p u t e co t + f2 t accurately when the m e a n m o t i o n is larger than twelve revolutions per day. At least two nearly equatorial satellites, one high and one low, will be necessary to determine J2j accurately f r o m the secular m o t i o n s o f the longitude o f the
850
YOSHIHIDE KOZAI
~oooo~ I \a~
J.
J,4
J~
10000 !
-20000
0~
15~
20000
30~
45~
60~
75~
90~
Jio
J1z
i
12
_,oooo
/ /~
/t
-? 0 0 0 0
o
Fig. 15.
J"
Jl~
,~~
~o ~
4%~
6~~
~~
90 ~
Coefficients o f J2j to c o m p u t e 091 for n = 14 revolutions and e = 0.1 as functions o f the inclination.
THE EARTH
GRAVITATIONAL
851
POTENTIAL
3000
2000
I000
- I000
- 2000
-3000
0~
i
I
i ~
i
30 ~
ooo -
i
;
45 ~
60 ~
,
r
,
5~
90 ~
J6
_,ooW' J6 - 3000~-, 0~
I
I
I
15 ~
30 ~
45 ~
I
60 ~
I
75 ~
90 ~
Fig. 16. Coefficients of J~j to compute I2~ sin i for n = 14 revolutions and e = 0.1 as functions of the inclination.
perigee. I f such satellites with e = 0.1 ~ 0.2 are available in future, tracking o f these satellites as well as satellites with inclinations 10 ~ 40 ~ 60 ~ and 80 ~ will become very useful for the determination of J2j, although, unfortunately, existing tracking stations near the equator are very few. As only K o z M (1964) used perigee data in the determination o f Jzj, the discrepancies in 091 are usually much larger than those in (21 and become sometimes as large as 0~ even for the satellites for which the argument o f perigee have been determined more accurately than +0~.0001. Although coefficients of J2s for corn-
852
KOZAI
YOSHIHIDE
TABLE V Factors to be multiplied by coefficients of J ~ in Figures 15, 16 and 18 as functions of the mean motion n
J~
J4
J6
J8
dlO
Jl~
Jl~
16 15 14 13 12 11 10 9 8
1.37 1.17 1.00 0,84 0.70 0.57 0.46 0.36 0.27
1.63 1.29 1.00 0.76 0.57 0,41 0.29 0.20 0.13
1.95 1.41 1.00 0.69 0.46 0.30 0.19 0,110 0,061
2.33 1.55 1.00 0.63 0.38 0.22 0.119 0.061 0.029
2.78 1.70 1.00 0.57 0.31 0.16 0.076 0,034 0,014
3.33 1,86 1.00 0.51 0.25 0.114 0.048 0.019 0.0065
2.97 2.04 1.00 0,46 0.20 0.083 0.031 0.0104 0.0031
o~
o?oos
",
/ /
ooooo
/
/ }
_ o.~oo5
/
/ / /
-~olo
/
/
/ KG
/ KO
- o~,o~5
I
--"--~--~
2
/
I
.......
2
/
$
............
/ / /
- # 0 2 0 oo
Fig. 17.
&o
i'oo
t'5.
i
ao~
asO
D a i l y secular motions o f the longitude o f the perigee, col -F ~2~, based o n values o f J2j in T a M e I V f o r n = 14 revolutions and e = 0.1 as functions o f the inc|ination.
853
THE EARTH GRAVITATIONAL POTENTIAL
p u t i n g ~ot d e c r e a s e m o r e s l o w l y t h a n t h o s e f o r c o m p u t i n g f2 I, t h e d a t a f o r t h e a r g u ment of perigee are useful for this kind of work. V a l u e s o f J2j b y d i f f e r e n t a u t h o r s d e v i a t e n o t o n l y b e c a u s e o f d i f f e r e n t m e t h o d s a n d o f d i f f e r e n t satellites u s e d b u t a l s o b e c a u s e o f s y s t e m a t i c d i f f e r e n c e s i n t h e a d o p t e d
0~
Fig. 18.
5~
1 0~
15"
20 ~
25 ~
Coefficients of J2j to compute coi + -Q1 for n 14 revolutions and e = 0.1 as functions of the inclination. =
TABLE V]' Deviations from the mean values for J2j,1 (in units of 10 -9 )
J~
J5
J7
GN
-- 91
69
-- 96
KCS 1 2
25 --145
-- 53 267
57 --443
26 39
-- 88 --113 --
KO1 2 S
146
82
G N : GUIER and NZWTON (1965) KCS: K1NG-HELE, COOK and SCOTT (1965)
J9
J~a
J18
41
13
19
16 394
13 --666
19 19
121 164
-- 97 -- 189
310 316
19 --95
197
--133
14
19
--
KO: KOZAI (1964) S: SMITH (1963)
854
YOSHIHIDE KOZAI X I O -3
,---'-qo,-9
0.30
/
{ .~ r-.c s 2 . . . .
\
......
'T'---'--'--'--IIIL
~" 1, I
L
T i
l
i
/
O25
! i. Q20
035
/ 0.05
I
;/"'x '
,",/
.
[
; 'V"
0.00
\, --
O,OS
--
0.10
7
i
i
/
q
\/
--045
0~
15 ~
30 ~
45 ~
60 ~
75 ~
90 ~
o.so KO 2 ON KCSI
- ...... ~ .
I t
II I I
I II t
II
I
L
l
I
O,25
O.20
II LI
0.15
0,10
0,05
o.oo
0.05
-
o.ko
-035
~Q2O
Fig. 19.
~
is ~
so ~
45 ~
6o ~
7~
9o0
T h e amplitude, Ae, b a s e d on various values of J2~+1 given in Table VI for n = 14 revolutions per day a n d e = 0.1 as functions of the inclination.
THE EARTH
GRAVITATIONAL
POTENTIAL
855
data themselves. KING-HELE and COOK (1965) and Kozai both used 1961v (i=280.8), 1960t2 (i=4772), 1961 o (i= 66~ and 1961 ecS1 (i=95~.9). However, Kozai's values of f21 are larger than King-Hele and Cook's by 00.0003 for both 1960t2 and 19617c51, although there are not such large systematic differences for the other two satellites. Also, for 196131 (i=380.8) used by both KING-HELE and Cook (1965) and SMITH (1965) the difference for f21 is larger than &.0005. Therefore, it may be necessary to derive orbital data more accurately, and, of course, for more satellites. Table VI presents deviations from the mean values ofd2~ + ~in GUIER and NEWTON's ( 1 9 6 5 ) , K I N G - H E L E , COOK and SCOTT'S (1965), KOZAI'S (1965) and SMITH'S (1963) solutions given in Table III. The values of Ae based on these values of J2j+l in Table VI are plotted in Figures 19 for n = 1 4 revolutions per day and e=0.1 as
J~
15oo ~-
J~
t.__ ~176
t
r
15~
3C
I 4 s~
6 o~
75~
9 0~
Fig. 20a. Fig. 20.
Coefficients o f J2~.1 to c o m p u t e Ae for n ~ 14 revolutions a n d e = 0.1 as functions o f the inclination.
856
YOSHIHIDE KOZA[
J5
J.
Js
J~3
J.
I0001
0o
15 ~
30 ~
45 ~
60 ~
75 ~
Fig. 20b.
functions of inclination. The coefficients of J2j+l to compute these are shown in Figure 20. For other values of the mean motion the coefficients of Jzj+ 1 are derived by multiplying (4j-2)/3-th power of (n/14) by those given in Figure 20, and this factor is given in Table VII. In the upper part of Figure 19 the inclinations of the satellites used for determining Jzj+l are marked. For satellite 1961e61 with inclination 95~ the mark is placed at 84~ Figure 19 shows that there are large discrepancies in A e near i = 10 ~ 30 ~ 45 ~ 60 ~ 70 ~ and 90 ~ These discrepancies are not always due to the lack of satellites near such values of inclinations. For some satellites the data used by different authors also deviate considerably from each other. In Table VIII values of Ae x 103 for the satellites used by at least two authors are shown. In KozAfs column a number of data used, range of Ae, and their mean value are given. The satellite, for which consistent values are given, is 1962flp only. For 1961 o, which all the authors used, Smith's value is larger
THE EARTH GRAVITATIONALPOTENTIAL
857
TABLE VII Factors to be multiplied by coefficients of Jzj- in Figure 20 as functions of the mean motion n
.]3
J5
J7
J9
.711
J13
16 15 14 13 12 11 10 9 8
1.09 1.05 1.00 0.95 0.90 0.85 0.80 0.74 0.69
1.3l 1.15 1.00 0.86 0.73 0.62 0.51 0.41 0.33
1.56 1.26 1.00 0.78 0.60 0.45 0.33 0.23 0.15
1.87 1.38 1.00 0.71 0.49 0.32 0.21 0.13 0.073
2.23 1.51 1.00 0.64 0.40 0.24 0.13 0.056 0.035
2.66 1.66 1.00 0.58 0.32 0.17 0.085 0.039 0.017
t h a n a n y o t h e r ' s v a l u e s . F o r 1961 v K o z a i ' s v a l u e w h i c h is t h e m e a n o f f o u r v a l u e s is l a r g e r t h a n t h e o t h e r s b u t is n o t i n c o n s i s t e n t w i t h S m i t h ' s v a l u e . T h e l a r g e s t d e v i a t i o n is t h a t b e t w e e n K i n g - H e l e , C o o k a n d S c o t t ' s a n d K o z a i ' s v a l u e s f o r 1961c~61. T o i m p r o v e v a l u e s o f 32j+ a it is n e c e s s a r y t o t r a c k n e w s a t e l l i t e s a n d t o e l i m i n a t e discrepancies in the existing data. F o r s a t e l l i t e s w i t h v e r y s m a l l e c c e n t r i c i t y t h e e x p r e s s i o n s (23), e = eo + A e
sin co,
(23)
co = coo + coat + Ao, c o s co,
Amplitudes
Ae •
TABLE VIII 10a for the satellites commonly used by more than one author KO
Satellites
i
GN
KCS
S *
mean
0.452 ~ 0.480
(4)
0.460 +_ 14
61 v
28~
0.434 • 4
59 ~ 1
32.8
0.454 •
0.453 ~ 0.475
(14)
0.466 z52
60 t 2
47.2
0.650 ,2
0.657 ~ 0.663
(4)
0.661 •
62 pit 1
50.1
0.775 •
0.775 ~ 0.782
(3)
0.777 zc4
61 o
66.8
0.251 •
0.250 ~ 0.266
(3)
0.259 ,6
61 a6 1
95.9
0.787 ~ 0.804
(2)
0.791 •
0.258 ,5 0.835 --8
* Number of data used. G N : GUIER and NEWTON (1965) KCS: KING-HEL~, COOK and SCOTT (1965)
0.434 • 16
range
0.275 ,8
S: SMrrH (1963) KO: KOZAI (1964)
858
YOSHIH1DE KOZAI
do not hold, since these expressions are derived by solving equations de -
-
dt
= Ae'co I COS co, (36)
do9 - -
dt
=
0) 1 --
A~'o)
1
sin co,
by assuming that all quantities except A~ are constant and co = o91 9t on the right-hand sides of the equations. However, if e is very small, Ao becomes very large and Ae/e o is not a small quantity. This means that e and A~ on the right-hand side of the equations (36) cannot be regarded as constants any more and co cannot be assumed to be a linear function o f time. However, if we put = e cos co, ~/= e Sill CO, the equations (36) become linear differential equations with constant coefficients with respect to ~ and r/, and they can be solved in the f o r m (24), that is, = e cos co = eo cos (coo + colt), ! r / = e sin co = eo sin (coo + col t) + Ae" 1
(24)
These expressions hold only for small eccentricity, since e 3 is neglected on the righth a n d side.
Ae
Fig. 21. Relations between e0 and e for nearly circular orbit. The coordinates of P, which is a center of circles, are (0, Ae).
THE EARTH GRAVITATIONAL POTENTIAL
859
If values of ~ and r/are plotted in the plane (Figure 21), the curve representing the solution (24) is a circle whose center is at P(0, A~). The radius of this circle is eo and the point Q on the circle is moving with a constant angular velocity coI. If eo is less than A~, the value of co is always between 0 ~ and 180 ~ The satellite, 1964-15A, used by King-Hele, Cook, and Scott has an extremely small value for e. Therefore, I am afraid that the equations (23) should not hold for this satellite, although they used the equations (23).
8. Reference Ellipsoid Once the value of J2 is given, the flattening of the reference ellipsoid, f, can be computed by the following formulas; 3 1 9 15 39 f =~J2+~m+~S 2 +2~8J2rn-5-6 rn2'
m = Rv2/ge,
(37) 27 jz
3
ge=(Gl~I/R 2) 1--,Ua + ~ S z - - m + ~ -
6
47
2-~Jzrnq-~m
)
2 'j
where R is the equatorial radius, v is the rotational angular velocity of the earth,/~, is the relative mass of the atmosphere, and ge is the apparent gravity at the equator. The third-order terms for J2 have been neglected in these equations. Numerical values adopted for the constants in (37) are, J2 = 1.08264 x i0 .3 , R = 6378160 m, GM = 398601 x 10 9 m 3 / s v = 0.000072921ls- ~, lq = 0.000001.
2,
(38)
By using an approximate value of m, g, is computed as g. = 978. 026 gals.
(39)
And then the value of m is obtained as 0.00346778, and the value for the flattening of the reference ellipsoid is derived as f = 0.00335280 = 1/298.258.
(40)
When this reference ellipsoid can be assumed to be geoid, the value of Ja is, of course, 1.08264 x 10 -a, and 4
J4 = - ~ f
2
4
+ ~ m f = - 2.35 x 10 -6 .
(41)
For this geoid the coefficients of zonal harmonics of odd orders and of non-zonal harmonics vanish, and Jzs is of order of J~ and is, therefore, negligibly small except for J2 and J4. Thus, the departure of the actual potential of the earth from that of the
860
YOSHIHIDE KOZAI
reference elliPsoid is expressed as
V' = - G--M-M { ~ J'~\(R-Y r ] Py (sin ]=3
(42)
j pjk (sin q))
_
(Cj,~cos
k2 + S s.,k sin k2
,
j=2k=l
where J4 is the a c t u a l v a l u e o f -/4 m i n u s t h a t g i v e n i n (41) for the reference ellipsoid. A n d the d e p a r t u r e o f geoid f r o m t h e ellipsoid is c o m p u t e d as h = _
R I ~ JJPj(sin q)) j=3
(43)
~ pk(sin~o)(Cj,kCOSk;~+Sj,ksink;O I 9 j=2k=l
T h e e x p r e s s i o n for h b a s e d o n K o Z A F s s o l u t i o n 2 (1964) for z o n a l h a r m o n i c s o f e v e n a n d o d d o r d e r s is w r i t t e n as, h = 0.8 - 18.3 sin ~o - 87.8 sin 2 tp - 119.1 sin 3 ~p + 1042.5 sin 4 ~p + 1191.7 sin 5 ~p - 5074.2 sin 6 ~p - 3636.7 sin 7 r + 12668.0 sin s r + 5230.8 sin 9 ~o - 16676.3 sin t~ q~ - 3556.4 sin i t q0 + 1 0 9 1 3 . 0 sin 12 q) + 926.8 sin 13 q~ - 2791.3 sin ~4 q~ (in m e t e r s ) . ~p h0~ 15 ITI
'
I
,
I
30 o t i '
i
l
i
~
60 ~ i '
90 ~ i
,
I0
I ~
15'" I0
5
5
0
0
-5
\,
-lo
-5 _l-Lo
,,,,,,
-15
~,
-20
-15
',
,
-20
, I ~ I ~1 "i " - -25 60 ~ 90 ~ ~p Fig. 22. Geoid height, h, as a function of geocentric latitude, r Solid line represents geoid height in northern hemisphere, and broken line represents that in southern hemisphere. I
-2500
~
I
L
I
,
30 o
I
861
THE EARTH GRAVITATIONAL POTENTIAL
The values derived by this equation as a function of the geocentric latitude, q~, are plotted in Figure 22. This figure shows that h takes the maximum value, 13.5 meters, at the north pole, and the minimum value, -24.1 meters, at the south pole. Discrepancies in the geoid height computed by values for Jj by different authors are within one meter due to odd-order terms except for KINc-HELE, COOS: and SCOTT'S solution 2 (1965), and are within two meters due to even-order terms except for KIN6-HELE and CooK's solution 2 (1965). Of course, higher-order terms cannot be neglected and if they can be included to compute the geoid height in future, the curve representing the geoid height will become much more complicated than that in Figure 22. 9. Non-Zonal Harmonics
Procedures similar to those for determining coefficients of zonal harmonics cannot be followed for determining those of non-zonal harmonics because of the difference in the periods of perturbations, by which the coefficients are derived from observations. Non-zonal terms in the potential produce mainly periodic perturbations with both short and long periods. However, the periods of the long-periodic perturbations which do not depend on the mean anomaly are not so long for this case as those due to zonal harmonics, and are usually less than one day. Although the amplitudes of the long-periodic perturbations are expected to be about ten times as large as those of the short-periodic terms theoretically, it is sometimes found that the amplitudes for some of the long-periodic terms are smaller than those for the short-periodic terms because of some other factors in the coefficients, and therefore, as it is necessary to use the short-periodic perturbations as well as the long-periodic perturbations to determine coefficients of non-zonal harmonics, it is not possible to detect them from the variations of mean orbital elements derived every two days for this purpose. Since the amplitudes of the perturbations due to non-zonal terms cannot become so large because of their rather short periods, the accuracy in determining coefficients of non-zonal harmonics is considerably lower than that for zonal harmonics, for which the secular terms and the long-periodic perturbations amplified by factor 1/J2 can be used. The non-zonal terms in the potential of the earth are written either as,
GM r
22 )-
R j pjk (sin ~o) (Cj, k cos k2 + Sj k sin k2)
r
j=2k=l or
as,
GM
J p~ (sin ~o) J j k COS k (2 - 2~'k)
r
' j=2k=l
where Cj, k and Sj, k are related to Jj, k and 2j, k as follows:
'
(44)
862
YOSHIHIDE KOZA1
cj,~ = sj,~ cos 1,,~j,~, S j, k = Jj, k sin k2j, k. 1
(45)
The longitude, 2, is measured eastwards, and J2,1 vanishes. The disturbing function due to Jj, k-term is expressed (KAuLA, 1961a) as, J
Uj, k = - - -
Fjk~, (i)
Gjpq (e) Sjkpq (co, l, f2,/~G),
(46)
F p=O
q = -- c~
w h e r e / ~ is the Greenwich sidereal time. By putting
j-k
h-
2
'
if
(j-k)
iseven,
if
(j-k)
isodd,
p' = p,
if
p is equal to or less than
p' = j - p,
if
p is equal to or greater than
j-k-1
h-
2
,
j/2, j/2,
Fi~p(i ), Gjpq(e), and Sjkpq(co, l, ~2, p~) take the following forms: (2j - 20! t! (j -- t)! 22(j-') sinJ-k-2t i
Fjkp(i) = t
c s=O
p - t - c (j s k - _ z t ) ~ '
C
p'--i
Gjv(Zp_j)(e ) = (1 - e2) ~ - j
2d + j - 2p'
d
d=0
sjk,~
2 (2j + 1) (j - k)!
I =t-sj,
C J, k l j - k e v e n kSj-k
+ ) S J ' kkl r
odd
cos {(j - 2p + q) (1 + co) - qco + k (f2 - Pa)}
k odd
where the value of t is summated from 0 to p or h, whichever is less; s is summated from 0 to k; and c is summated over all values making the two binomial coefficients non-zero. Note that Gjp~ is order of e Iql and Fjkp is of sin Ij-k-2pl i. The perturbations in an orbital element due to the non-zonal harmonics are therefore expressed as linear functions of Ci, k and S],k in the following form: J
dE = f2 2 ~ (Pjkvo Cj,k + Qjkpq Sj, k), j=2k=lp,
q
(47)
THE EARTH GRAVITATIONAL POTENTIAL
863
where Pjkpq and Qjkpq have divisor, ( j - 2p + q) (n + coi) - qcol + k (f21 - v). The constant v is the rate of the sidereal time and is 1.00274 revolutions per day. The expressions for the semi-major axis, and, therefore, also for the mean motion do not contain terms with j - 2 p + q = 0, since the disturbing function should be differentiated with respect to the mean anomaly before integration so as to derive the perturbations for the semi-major axis. The expression for the mean anomaly is derived by integrating again the expression for the mean motion, that is, by double integrals, and therefore contains divisor { ( j - 2p + q) (n + 02) - qco~+ k(f2a - v)} 2. For examples IZSAK (1964) gave for the coefficient of C2, 2 in the expression of the mean anomaly for 1959el, Vanguard 2, 4.008 sin {2 (l + co) + 2 (f2 - f t G ) } -
-
-
2.979 sin {2 (f2 - - / ~ ) }
(48)
0.025 sin { - 2 (t - co) + 2 ( a -
and for the coefficient of $3,1 in the expression of the longitude of the ascending node for the same satellite, - 2.255 cos { - co + ( a - ~o)} + 5.481 cos {co + ( e - ~o)} - 1.113 cos { - (co + I) + (O -/~G)} -
(49)
3 . 1 3 5 c o s {(co + 0 + ( e -
KaULA (1963a) gave the coefficient of C3, i in the expression of the eccentricity, for 1960t2, rocket of Echo I, as 1.850 cos (co + ~2 - ~~ + 5.058 cos ( - co + ~2 -/ZG) -- 0.609 COS (-- CO-- 2I + ~2 -- P G) -- 0.001 COS (e) + l + /2 -- #G)
(50)
+ 0.002 COS (-- CO-- I + f2 -- #~). Note that in the expression the amplitude for the first term with short period is larger than that for the second term with a half-day period. By using the same matrix A as in the equation (13) the perturbations due to the non-zonal harmonics in the topocentric coordinates and in the range-rate are written as,
(cod) d6
dp/p dt~/p
= A.dE,
(51)
where dE is a vector consisting of dE in (47), the perturbations in the orbital elements, as components. The right-hand side of this equation is also a linear function of C~, k and Sj, k, which are unknown quantities. By equating the residuals in the observed topocentric coordinates and range-rate after determining mean orbital elements to the left-hand side of Equation (51), a limited number of coefficients, Cj. k and Sj, k, can be solved. This can be done under
864
YOSHIHIDE KOZAI
the assumption that the residuals are entirely due to the limited number of non-zonal harmonics. Even if this assumption were correct, it is necessary to have many observations uniformly distributed with respect to positions in orbit and in time to determine C j, k and Sj, k accurately. Of course, there are many other sources giving systematic errors in the residuals; namely, errors in coordinates of tracking stations, deviations of the air-drag effects from those evaluated theoretically, errors in frequency constant of the radio wave for Doppler observations, and so on. Some of these errors can also be derived as corrections from the observations. However, if perturbations due to the non-zonal harmonics are not taken into account, the true values for mean orbital elements cannot be derived by the differential orbit improvement program, especially, when the observations are not uniformly distributed. For example, suppose that the true mean orbit is a solid line curve drawn in Figure 23 and that the perturbations due to the non-zonal harmonics displace the orbit to the curve described by a broken line. When the observations available are not
./,,, Fig. 23. Bias in orbit determination due to perturbations. Solid line represents the mean orbit, which is distorted to broken line due to the perturbations. The observations are made at *
uniformly distributed, but all the observations are distributed, say, outside the mean orbit as plotted in Figure 23, the differential orbit improvement program provides orbital elements different from those for the mean orbit. Therefore, if the residuals based on these elements are supposed to be due to non-zonal harmonics, the coefficients C~, k and S~, k determined from these residuals are usually less than their actual values. Therefore, it may be necessary to derive a limited number of coefficients, Cj, k and Sj, k, corrections of the orbital elements and of the station coordinates simultaneously from observations, although observations for a long interval of time cannot be used in one file to reduce errors in evaluating air-drag effects to the satellite motion. To reduce errors caused by the air-drag in the determination IZSAK (1962) assigned higher weight to the across-track than to the along-track component of an observation. Even when observations uniformly distributed are available, by tracking one satellite alone C2.2, C4, z, and C6, 2, for example, cannot be separately solved, since their coefficients, Pjkpq, are periodic functions with common periods. In order to solve them separately observations for many satellites with different orbital elements during
THE EARTH GRAVITATIONAL POTENTIAL
865
various periods of time must be used. This means that there are too many unknowns to be solved simultaneously. To avoid this difficulty the method of stepwise least squares is employed by many authors. In the first step, the mean orbital elements for each orbital arc are derived separately by holding values of Cj,k and Sj,~ and of station coordinates fixed. And in the second step the coefficients of the non-zonal harmonics and the corrections of the station coordinates are determined from all the residuals of observations which come from many orbital arcs of many kinds of satellites. These two steps are iterated so that finally satisfactory results are derived. However, even by this procedure systematic
C2, 2
5 i
9 dettreys
3<'10"6 I0 i
Sa.2 X eKoulo
9 Newton
9 Zhongolovitch
9 Kozo i ~5
9 IzsaK
-I0 -6 X io
Fig. 24. Values of Ca, z and $2, ~ determined before 1962 from both satellite motions and gravity measurements. The mark • is placed at the point corresponding to the most probable values now known.
errors might be produced as Figure 23 shows, unless the distribution of observations is uniform. For this reason radio method has its advantages in determining coefficients of non-zonal harmonics. The first determination of C2, 2 and $2,2 from satellite motions was made by IzSAK (1961) with precise Baker-Nunn observations for 1959cd and 1959r/. In Figure 24 values of C2,2 and $2,2 derived by lZSAK (1961), KAULA (1961b), KOZAI (1961) and NEWTON (1962) in 1961 and 1962 from satellite observations are shown with those by JEFFREYS (1959) and ZIqONGOLOVlTCrI(1957) from gravity measurements. As this diagram shows, these values were widely spread at the time. In Figure 24 the most probable values now known are marked. Earlier values are found to be large, probably because there were large errors in reduction of observations.
-2.C
-1,0
0.0
Fig. 25.
-[1
10
@
0
20
X
[3
9
_
-()5
A
F i g . 25a.
O0
10
20
Q
0
(3,33
-05
QO
0.5
@
1'0
(3,1)
@
10 t
A z0 i
Kaula Kozai
o /x
(.1963)
(1963)
(1966)
L -
-IO
9
D A
o
(1962) (4,1)
@ Uotila
Izsak
[]
(1966)
Wagner (1965)
X
@
G u i e r and N e w t o n
Allan
9
"Ld
Andeple (1965)
o
9
9
--
-1.0
00
05
(1965)
Values of Cj, rc and Sj, ~ by various authors. The abscissas represent Cj, ~ and the ordinates S~, ~ in units of 10 -6.
(3,2)
@
(2,2)
LO
N
tn
9
0,5
9- 0 5 f 9 (6.2)
0.0
(5,z)
Q5
f
o~5
l.O
0.0
"~ O5
-(15
9
51 Z O0(7.3)
-O5
-o5
0.0
o9
1'0f [](4'2)9
[]
-0,4
-O
0,0/ ~
.L 00
. (7,4) J
-05,
n&
05
@
Fig. 25b.
-05
[]
9
-10
o.s -o~
(6,5)
1
O0
1.o
O0
(7,6)
D
i
0.5
m
(5,5)
(6,6)
@
o9
05
-05
(4,4)
4 ~176
]05
-05
(5,4)
051 (7,5)
-t -0,5
"] []
o9 (p
9
-0
-qs~9gll
-o.s "I"
9
(6,4) 9
-05
0.0
0.5
[]
Z3 9 LO
(4,3)
(5.33
o.4]1- (6,3)
-110
-0,5
O0
0.5[ 0
0.5
0,5 i
I
[9
(v,t)
--Oil
9
~
0 . 0 ~ 5
o.5
(
J
0
5
0.5
O0
9
(5,t)
(6.1)
5
-05
i
~t
868
YOSHIHIDE KOZAI
These discrepancies have been diminished for the recent values, as the first diagram i n F i g u r e 25 i l l u s t r a t e s . D i a g r a m s e x p r e s s i n g v a l u e s o f C j, k a n d S j, k b y v a r i o u s a u t h o r s a r e s h o w n i n F i g u r e 25, w h e r e t h e a b s c i s s a e r e p r e s e n t Cj, k a n d t h e o r d i n a t e s r e p r e s e n t S j, k i n u n i t s o f 1 0 - 6 . I n t h e s e d i a g r a m s , d a t a f r o m r a d i o o b s e r v a t i o n s o f s a t e l l i t e s a r e denoted by shaded symbols, those from optical observations by open circles, and the gravity data by double circles. ANDERLE (1965) a n d GUIER a n d NEWTON (1965) u s e d D o p p l e r
observations,
and
ALLAN (1966) a n d WAGNER (1964) u s e d d a t a f o r a c c e l e r a t i o n s i n t h e m e a n l o n g i t u d e o f S y n c o m II, t w e n t y - f o u r
h o u r s a t e l l i t e , w h i c h will b e d i s c u s s e d l a t e r . IZSAK (1966),
KAULA ( 1 9 6 3 b ) , a n d KOZAr (1963, u n p u b l i s h e d ) IZSAK (1966) u s e d 2 6 2 4 4 o b s e r v a t i o n s
used optical observations. Although
of eleven satellites, optical data gave smaller
v a l u e s f o r Sj, k a n d Cj, k t h a n r a d i o d a t a , p r o b a b l y concerning
because of the reason explained
F i g u r e 23. I n T a b l e I X r e s u l t s f o r ANDERLE (1965), GUIER a n d NEWTON
(1965) a n d I Z S A I ( ( 1 9 6 6 ) a r e g i v e n , a l t h o u g h ANDERLE a n d GUIER a n d NEWTON g a v e more coefficients than those given here. For non-zonal harmonics
a l s o Cj, k a n d $i, k d o n o t c o n v e r g e r a p i d l y a s / i n c r e a s e s ,
a n d K A U L A (1966) e s t i m a t e d t h a t c o e f f i c i e n t s o f s p h e r i c a l h a r m o n i c s
of order j are
o n a n a v e r a g e o f t h e o r d e r o f 8 x 1 0 - 6 / j 2.
TABLE IX Coefficients of non-zonal harmonics (in units o f 10 -6) ANDERLE (1965)
GtnER and NEWTON (1965)
IZSAK (1966)
j k
G,~
Sj,~
G,~
Sj,~
G,~
Sj,~
2,
2
2.45
-- 1.52
2.38
-- 1.20
2.08
-- 1.25
3, 3, 3,
1 2 3
2.15 0.97 0.57
0.28 -- 0.91 1.65
1.84 1.22 0.66
0.21 -- 0.68 0.98
1.60 0.38 -- 0.17
-- 0.04 -- 0.80 1.40
4, 4, 4, 4,
1 2 3 4
-- 0.50 0.27 1.00 -- 0.47
-- 0.58 0.67 -- 0.17 0.47
-- 0.56 0.42 0.85 -- 0.21
-- 0.44 0.44 0.01 0.19
-- 0.38 0.20 0.69 -- 0.11
-- 0.40 0.58 -- 0.10 0.43
5, 5, 5, 5, 5,
1 2 3 4 5
0.03 0.61 -- 0.30 -- 0.51 0.20
-- 0.12 -- 0.31 -- 0.12 0.13 -- 0.41
0.14 0.27 0.09 -- 0.49 -- 0.03
-- 0.17 -- 0.34 0.10 -- 0.26 -- 0.67
-- 0.14 0.24 -- 0.67 -- 0.13 0.08
-- 0.04 -- 0.27 0.05 0.16 -- 0.41
6, 6, 6, 6, 6, 6,
1 2 3 4 5 6
-- 0.09 0.16 -- 0.02 -- 0.26 -- 0.12 -- 0.43
------
0.19 0.48 0.14 0.26 0.74 0.43
0.00 -- 0.16 0.53 0.31 -- 0.18 0.01
0.10 -- 0.15 0.05 -- 0.51 -- 0.51 -- 0.23
-- 0.02 0.05 0.05 0.07 -- 0.28 -- 0.12
0.12 -- 0.23 0.00 -- 0.39 -- 0.38 -- 0.59
869
THE EARTH GRAVITATIONAL POTENTIAL
10. Resonance Effects As was described in the previous section, the disturbing function contains periodic terms with argument (j-2p + q)(l+ co)-qco +k(O-#c) due to non-zonal spherical harmonics of order j and degree k, that is, Jj.,~-term. Since the coefficients of these terms contain e Lql as factors, the term corresponding to q = 0 takes the largest amplitude, when the eccentricity is very small. Therefore, the most important term has the argument, = (j - 2p) (l + co) + k (f2 - Pc), (52) where p takes a value between 0 and j. When the time derivative of y, that is,
(53)
= (j - 2p) (n + col) + k (f21 - v),
takes a very small value, the perturbations due to this term become very large, especially in the mean longitude, as the expression of the mean longitude is derived by integrating the disturbing force function twice with respect to time, and, therefore, has divisor ~z. For twenty-four hour satellite the values of ~, for which j - 2 p = k, takes an extremely small value, since n is equal to v and the values of cot and O 1 are very small because of the large value of the semi-major axis (6.61 earth radii). The perturbations due to the terms, for which j-2p=k, in the semi-major axis for this satellite are regarded as linear secular terms, and those in the mean longitude as secular accelerations for a limited interval of time. In Figure 26 the drifts of the values of the semi-
-64
I
I
I
I
_e4
I
t
o
- 62
-62 e e
-60
-60
-58
--58
-56
-56
O0 O 0000
O0 4 7 1 7 4 krn
-54
x )C
4 2 1 7 2 km
X
X
7o
X
72 70
X
x X
X XXx
66
x
X
I
Segt QO
68
X
x
x x
x I
I
OctQO 1963
Fig. 26.
x
X
68
Nov.QO
66 I
Dec.O.O
I
I
Jan.O.O
Feb.O.O |964
Drift of the geographical longitude (.) where the satellite crosses the ascending node and of the semi-major axis ( x ) for Syncom II.
870
YOSHIHIDE KOZAI
m a j o r axis and o f the geographical longitude at the point where the satellite crosses the ascending node are shown for S y n c o m II (1963-31A) for six months (WAGNER, 1964). The data were c o m p u t e d at the G o d d a r d Space Flight Center of N A S A . The drift was stopped on N o v e m b e r 12, 1963 by firing o n - b o a r d jets. The inclination o f S y n c o m II is 33 ~ and the eccentricity is less than 2 x 10 -3. The most important term causing the drift is J2, 2 (P = 0), and the values of C2, 2 and $2,2 by WA~NER (1964) plotted in Figure 25 were derived by the drift curves in Figure 26. ALLAN (1966) derived values of Cj, k and Sj, k for (2, 2), (3, 3) and (4, 4) by the data o f S y n c o m II for five drift periods between August 18, 1963 and September 29, 1964. These values were also plotted in the corresponding diagrams of Figure 25, although ALLAN t h o u g h t that the values for (3, 3) and (4, 4) were not so accurate. TABLE X Accelerations in longitude of twenty-four hour satellites (degrees/day 2) • 10-3 (KAULA, 1966) Syncom 1I (i = 33~
Observed Calculated ANDERLE(1965) GUIER and NEWTON (1965) IZSAK(1966) ALLAN(1966)
Syncom II1 (i = 0 ~
1
2
3
1
2
__ 55~ ~ __ 58~
__ 59~ ~ __ 63~ ~
245~ ~ 197~
180~ ~ 178~
179~ 174~
- 1.27
- 1.32
1.13
1.20
1.07
-- 1.26
-- 1.26
1.09
1.11
1.03
-- 1.25 -- 1.05 -- 1.28
-- 1.22 -- 1.06 -- 1.30
1.03 0.98 1.09
0.82 0.92 0.95
0.74 0.85 0.88
KAULA (1966) c o m p a r e d observed accelerations in three periods for Syncom II and in two periods for Syncom I I I ( i = 0 ~ with c o m p u t e d values based on the results by ANDERLE (1965), GUIER and NEWTON (1965), IZSAK (1966) and ALLAN (1966). KAULA'S comparisons are reproduced in Table X, which shows that the most satisfactory agreement b o t h for Syncom II and I I I is found in Anderle's values. Therefore, it m a y be concluded that ANDERLE'S solution (1965) for Cj, k and Sj, k is the best representation of the potential as far as low-order harmonics are concerned. Besides twenty-four h o u r satellite, there is another kind of resonance effect even for close satellites; namely, when n is nearly equal to, say, fourteen revolutions per day, and ~ in (53) takes a small value for j = k = 14 and p = 0 , the amplitude o f the perturbation becomes very large, especially for the mean longitudes. This kind o f perturbations was f o u n d by ANDERLE (1965) for 1961 ol, 1962flpl, and 1963 49B and by YIONOULIS (1965) and GUIER and NEWTON (1965) for 1963 38C and 1963 49B, f r o m observations. The values for C13,13 and S~3,13 derived by them were plotted in Figure 25, although Anderle derived values also for j = 15 and k = 13 and j = 15 and
871
THE EARTH GRAVITATIONAL POTENTIAL
300
&
z~
200
100
a-~
a
oe ar a
~
~a
a
4 ~ . ~~
~,
~z b, o
l~a
-I00
a
a
~ .^
a
a
a a,x,, a~, a
A~
*
~Aa
a 3-
~
a
*
-aA-a~ zx a
~
a
-I00
l&
-300 A~
t 85
86
!
!
!
$7
88
89
t t90
!
91
Fig. 27. Along-track residuals due to .]13, 13 for 1963 49B ( A N D E R L E , 1965).
k = 14. Figure 27 shows the residuals in the along-track component of observations for 1936 49B due to this resonance effect (ANDERLE, 1965). The period of this perturbation is 2.5 days. When more satellite with such resonance effects are tracked, non-zonal terms of high orders will be determined accurately. 11.
Geoid
map
In Figure 28 three geoid maps giving contours of equal geoid height by using ANDERLE'S (1966), GUIER and NEWTON'S 0965) and lZSAK'S (1966) values of Cj, k and Sj, k are given. The unit in these maps is meters. The coordinates where maxima and minima of geoid height are located and the geoid heights there are presented in Table XI. As far as locations of maximum and minimum points are concerned, agreements for the three maps are rather satisfactory, although the heights themselves deviate from each other by twenty meters. In Figure 29 a geoid map based on astro-geodetic data (FISCHER,1960; KAULA,1963c) is given for comparison, although the contours can be drawn only on continents. Since both zonal and non-zonal harmonics of higher orders have been neglected in the determination from satellite observations, fine structures due to higher-order harmonics in the geoid map do not appear in Figure 28. In fact, deflections of vertical lines computed by these maps are at most ten seconds of arc, even though we know
872
YOSHIHIDE KOZAI
that from gravity measurements deflections of vertical lines become as large as a minute of arc somewhere. Therefore, the geoid maps in Figure 28 are regarded as maps for a kind of averaged geoid, and the real geoid map can be derived by connecting geodetic datum points which have not been linked by triangulation. This connection can be made by also satellite observations. As has already been mentioned, corrections of tracking station coordinates are determined simultaneously when coefficients for non-zonal harmonics in the potential are determined from satellite observations. The corrections were derived and published by ANDERLE(1966) and IZSAK (1966). Also VEIS (1965) computed corrections for the coordinates of tracking stations without determining coefficients in the potential. KAULA (1966) computed corrections for the coordinates of North American, European, and Tokyo datum points from those of the tracking stations belonging to the these datum points, and the results are given in Table XII. According to this table agreements in latitude direction are rather satisfactory, whereas those in longitude direction are not, because corrections for the mean longitude of satellites for low inclination satellite and for the longitude of the ascending node for high inclination satellites and those for the longitude of the station cannot be determined separately.
=,,=
W a
,..I
- 180
-150
-120
-90
-60
-30
0
LONGITUDE
30
60
90
120
150
t~,O
{DEG.)
Fig. 28a. Fig. 28. Maps of geoids by ANDERLE(a) (1966), GUIER and NEWTON(b) (1965), and IZSAK(C) (1966). The unit is in meters.
873
THE EARTH GRAVITATIONAL POTENTIAL
~o:
0
o
~
....
o
17t:11: 7:~
' ~ + i ~
Fig.
28b.
+90* +80 ~
~
-
~
~
a.
~
~
-~
.-"
:"
.
§ 60 ~
+40 =
0o
_~o~ -80 -90
::: ......:t.:::Z:;ir:~!2:=:.~
~
"
~-
.
.
_
~
-
..:
:
"I:Z2 :~
_..
~
"s
-160 ~ -140"
-120" - I 0 0 ~ - 8 0 ~
-60"
-400
-20 ~
O*
Fig.
28c.
+20 ~
+40 ~ +60 ~
"80 ~
* I 0 0 ~ +120~
+ 1 4 0 ~ e160~
+180~
YOSHIHIDE KOZAI
874
TABLE Maxima Location
England
South Africa
and minima
ANDERLE ( 1 9 6 6 )
2
--
55
h ~ --
--
67 m
43 m
15 ~
50 ~
45
~o
5
India
h .1.
75 m 75 ~
--
25
--
95m
,~
- - 115 ~
-- 130 ~ --
South Pacific
~o
--
--
North
h )~
-- 75m - - 175 ~
West Atlantic
~o h
F i g . 29.
35m 70 ~ 45m
Geoid heights based on Astro-Geodetic
70
10 63 m 75 ~ 5
--
73m
-- 120 ~ 20 --
46m 115 ~
--
40
--
42m
45 ---
30 --
59m 180 ~
-- 45m - - 175 ~
40 ---
--
25
-- 70m -- 170 ~
~o
77m
50 49 m 150 ~
5 --
15
70
50 ~ --
60 m 75 ~
10
h
h 2
50
39 m 140 ~
h )~
Pacific
15 ~
65 m
45 m 140 ~
~o
--
50
h 2
(o
5~
lZSAK ( 1 9 6 6 )
55
Japan
East Pacific
G U I E R a n d NEWTON ( 1 9 6 5 )
15 ~
~o
~o
XI
in geoid heights (in meters)
22m 75 ~
--
30 --
Data
70m
50 ~ 10j
--
23m
(FISCHER, 1 9 6 0 ; K A U L A , 1963C).
875
THE EARTH GRAVITATIONAL POTENTIAL
TABLE XII Differences of satellite from terrestrial datum positions. Terrestrial positions by KAULA(1961C) Datum North America
AR
R cos ~0.A2 R.A~o AR
Europe
Rcos ~0.A2
IZSAK(1966)
VEIS(1965)
-- 3 m 14 11 -- 41 -- 69
33 m 25 19 -- 49 -- 34
R. A(o
Tokyo
2
3
20 34
15 -- 109
-- 8
9
R Rcos q~.A2 R.A~
ANDERLE (1966) 6m
--15 -- 6
These discrepancies will be diminished by observing satellites simultaneously from more than two stations belonging to different datum points. 12. Values of G M
The value of the geocentric gravitation constant, G M , can be determined by measuring the mean motion and the semi-major axis independently of observations and by Kepler's third law. KAULA (1963b) derived from Baker-Nunn observations G M as (3.986013 + 30)x 1014 m3/sec 2, whereas KAULA (1961C) obtained the value 3.986016 + 28 from terrestrial data. However, the most accurate determination of G M was made by tracking Ranger flights (SJOGREN et al., 1964; WOLLENrtAUPTet al., 1964). The mean value of their results is G M = 3.986009 x 1014 m3/sec 2 ,
(54)
+_6 whereas the International Astronomical Union in 1964 adopted the value of 3.98 603 as G M . By tracking Ranger flights the value of the lunar mass in units of the mass of the earth was found to be 1/81.303. 13. T i d a l Effects
In the preceding sections it was assumed that the shape of the earth or geoid is invariable, that is, the potential of the earth is independent of time. However, we know that the sun and the moon are deforming the earth by their tidal forces not only in the ocean but also in the body. The correction terms in the disturbing function due to the deformations of the earth are written as (KAULA, 1962; KOZAI, 1965),
-
oR5 {m: e2 (cos/~| r 3 k2
m,
}
+ ~ 1'2 (cos p~) , r)
(55)
876
YOSHIHIDE K O Z A I
where R is the earth radius assumed to be constant, m o and m) are the masses, r o and r~ are the geocentric distances, and fie and fl~ are the geocentric angles of radius vector of the satellite to the directions of the sun and the moon, respectively. In (55) P2 is Legendre polynomial of order two and kz is Love's number. Because of the time lag of the earth tides due to the friction, ro, r~, flo, and fl~ should be corrected by taking into account the lag of the tides (KoZAI, 1965). The disturbing factors due to the terms in (55) are k2(no/n) 2 (R/a) s and k2rn ~ (n)/n) 2 (R/a) 5, and are kz(R/a) 5 times as large as those due to the solar and lunar
89.915
0
I
I
I
89.910
I
1
I
a
g, v
0 I,-
89.905
oee e coo
Z -J U Z
00
89.9075; 89.9050
Fig. 30.
. 150
180
o ' . "%g~ 210
240
..
%'.ee~ ~
270 300 (deg,ees
330
360
(a) Inclination for satellite 1963 38C, plotted as functions of 20 - f2. (b) Residuals after correcting known perturbation.
gravitational forces, respectively. The value of k 2 is nearly 0.3. When short-periodic terms are neglected, perturbations do not appear in the expressions of the semi-major axis and the eccentricity. NEWTON (1965) tried to detect the perturbations due to the solar tide in the inclination for two satellites, 1963 38C and 1963 49B, with nearly circular and nearly polar orbits. He plotted the residuals for the inclinations during October 1, 1963 and March 17, 1964 for 1963 38C and during January 22 and December 14, 1964 for 1963 49B against values of 2 o - O (2 o being the mean longitude of the sun) as 2 (2 o -f2) is the argument of the most important perturbation terms. These diagrams are reproduced in Figures 30a and 31a. The longitudes of the ascending nodes for the two satellites move very slowly, and were 11~ and 101~ during periods for 1963 38C and 1963 49B, respectively.
THE EARTH GRAVITATIONAL POTENTIAL
877
The solar gravitational perturbations 0~ cos 2 (2 o - O) should be subtracted from these values. Also perturbations due to the motion of the equator referred to the inertial system of coordinates by the precession and nutation (KozAt, 1960) should be subtracted, since these perturbations have the same period. After these corrections were made, residuals plotted in Figures 30b and 31b still remained systematically and may have been caused by the tidal deformation of the earth. Newton computed from these diagrams values of k 2 as,
k2 = 0.31 _ 0.03, for 1963 38B, k2 =0.33 + 0.04,
for
1963 49C.
These values of k2 are very reasonable, although he derived negative values as the time lag of the tides.
o' 89:960
--
'
I
'
'
9 9
89.955
'
'
I
'
'
Q IQ
--
9o Q
~,
I I
9
9
me
e~
oee
9
o~ oo
89.950 o
j-
_
9 I~176
z
u"i
89.960
Z
9 9 9 ~176176 89.955
I 180
~
I 270
oe%e oOee% ~ I
n
I 0
~. | __ ~e'~
m
n
[ 90
I
n !80
(degrees)
Fig. 31. (a) Inclination for satellite 1963 49B, plotted as functions of 20 -- f2. (b) Residuals after correcting known perturbations.
KozAI (1965) tried to detect also the perturbations in both the inclinations and the longitudes of the ascending nodes for 1959~1, 195911, and 196012, but could not derive consistent results. Besides the perturbations due to the tides, KozAI (1966) derived variations with one year period in the longitude of the ascending node for 1960t2. He thought that the variations were due to change in J2 caused by change in moments of inertia of the earth. And if this idea is correct, the amplitude of the variation of Jz could be computed as 2 x 10 -9, from the data. However, more data are necessary before we can derive any definite conclusion. The computations in this paper were made b y using HITAC 5020 at the Computer Center, University of Tokyo, with assistances by Mr. Y. Hatanaka and Miss C. Takahashi.
878
YOSH1HIDEKOZAI
References ALLAN, R. R. : 1966, I~n Proceedings of the Second International Symposium 'The Use of Artificial Satellites for Geodesy' (in press). ANDERLE, R. J.: 1965, J. Geophys. Res. 70, 2453. ANDERLE, R. J. : 1966, In Trajectories of Artificial Celestial Bodies as Determined from Observations, Springer (in press). BAILIE, A. and BRYANT,R. : 1960, Astron. J. 65, 451. BATRAKOV,Yu. V. : 1963, In Dynamics of Satellites (ed. by M. ROY), Springer, 74. BERROTH, A. and HOEMANN,W. : 1960, Kosmische Geodiisie, G. Braun, Karlsruhe. BROUWER, D.: 1959, Astron. J. 64, 378. COOK, G. E.: 1962, Geophys. or. Roy. Astr. Soc. 6, 271. COOK, A. H.: 1963, Space Sci. Rev. 2, 355. CooK, A. H.: 1965, Geophys. J. Roy. Astr. Soc. 10, 181. COOK, G. E. and PLIMMER,R. N. A. : 1960, Proc. Roy. Soc. A 258, 516. ESCOBAL, P. R. : 1965, Methods of Orbit Determination, John Wiley, New York. FISCHER, I. : 1960, 9". Geophys. Res. 65, 2067. GAPOSCHKIN,E. M. : 1964, Smithsonian Astrophys. Obs. Special Rept. No. 161. GUIER, W. H. and NEWTON, R. R.: 1965, J. Geophys. Res. 70, 4613. IZSAK, I. G.: 1961, Astron. J. 66, 226. IZSAK, f. G. : 1962, In Space Age Astronomy (ed. by A. J. DEUTSCH and W. B. KLEMPERER), Academic Press, New York, 151. IZSAK, I. G. : 1964, J. Geophys. Res. 69, 2621. rZSAK, ]~. G." 1966, In Trajectories of Artificial Celestial Bodies as Determined from Observations, Springer (in press). JEFFREYS,H.." 1959, The Earth, 4th edition, Cambridge Univ. Press. KAULA, W. M. : 1961a, Geophys. J. Roy. Astr. Soc. 5, 104. KAULA, W. M.: 1961b, In Space Research (ed. by H. C. VAN DE HULST), North Holland Publishing Co., Amsterdam, 360. KAULA, W. M.: 1961c, d. Geophys. Res. 66, 1799. KAULA, W. M. : 1962, Advance in Geophysics 9, 191. KAULA, W. M. : 1963a, J. Geophys. Res. 68, 473. KAULA, W. M." 1963b, J. Geophys. Res: 68, 5183. KAULA, W. M. : 1963C, In The Use of Artificial Satellites for Geodesy (ed. by G. VEIS), North-Holland Publishing Co., Amsterdam, 335. KAULA, W. M. : 1966, In Proceedings of the Second International Symposium, 'The Use of Artificial Satellites for Geodesy' (in press). KING-HELE, D. G.: 1958, Proc. Roy. Soc. A 247, 49. KING-HELE, D. G. and COOK, G. E.: 1965, Geophys. J. Roy. Astron. Soc. 10, 17. KING-HELE, D. G., COOK, G. E., and SCOTT,D. W.: 1965, J. Plan. Space Sci. 13, 1213. KOZAI, Y. : 1959a, Smithsonian Astrophys. Obs. Special Rept. No. 22, 10. KOZAI, Y.: 1959b, Astron. J. 64, 367. KOZAI, Y.: 1960, Astron. J. 65, 62l. KOZAI, Y.: 1961, Astron. J. 66, 355. KOZAI, Y.: 1962, Astron. J. 67, 446. KOZAI, Y.: 1964, Publ. Astron. Soc. Japan 16, 263. KOZAI, Y. : 1965, Publ. Astron. Soc. Japan 17, 395. KOZAI, Y.: 1966, In Trajectories of Artificial Celestial Bodies as Determined from Observations, Springer (in press). MERSON, R. H.: 1961, Geophys. J. Roy. Astron. Soc. 5, 17. MERSON, R. H. : 1963, In Dynamics of Satellites (ed. by M. ROY), Springer, 83. MUELLER, I. I. : 1964, Introduction to Satellite Geodesy, Frederick Ungar Publishing Co., New York. MUSEN, P. : 1959, J. Geophys. Res. 64, 2271. NEWTON, R. R.: 1962, d. Geophys. Res. 67, 415. NEWTON, R. R.: 1964, Science, 144, 803. NEWTON, R. R.: 1965, J. Geophys. Res. 70, 5983. O'KEEEE, J. A., ECKELS,A., and SQUIRES,R. K.: 1959, Astron. J. 64, 245.
THE EARTHGRAVITATIONALPOTENTIAL
879
SHAPIRO, t. L : 1963, In Dynamics of Satellites (ed. by M. ROY), Springer, 257. SIOGREN, W. L., CHtJRKENDALL,D. W., HAMILTON,T. W., KIRHOrER, W. E., LtU, A. S., TRASK,D. W., WINNEBERGER,R. A., and WOLLENHAUPT,W. R. : 1964, Jet Propulsion Lab. Technical Rept. No. 32-605. SMITH, D . E. : 1963, J. Plan. Space Sci. 11, 789. SMITH, D. E.: 1965, J. Plan. Space Sci. 13, 1151. SOCHILINA,A. S. : 1963, In Dynamics of Satellites (ed. by M. ROY), Springer, 202. STERNE, T. : 1959, J. Amer. Rocket Soc. 29, 777. UOTILA, U. A. : 1962, Publ. Isostatic Institute of lAG, Helsinki, No. 39. VEIS, G. and MOORE, C. H. : 1960, Seminar Proceedings, Tracking Programs and Orbit Determination, Jet Propulsion Lab. 165. VEIS, G.: 1963, Space Sci. Rev. 2, 250. VEtS, G. : 1965, In Space Research 5 (ed. by P. MULLER),North-Holland Publishing Co., Amsterdam, 849. WAGNER, C. A. : 1965, NASA Technical Note, D-2759. WOLLENHAUPT,W. R., TRASK,D. W., SJOGREN,W. t;, PIAGGI, E. G., CHURKENDALL,D. W., WINNEBERGER,R. A., L1u, A. S., and BERMAN,A. L.: 1964, Jet Propulsion Lab. Technical Rept. No. 32-694. YIONOLrLIS,S. M.: 1965, J. Geophys. Res. 70, 5991. ZHONGOLOVITCH,L D. : 1957, Bull. Institute Theoretical Astron. Leningrad 6, 505. ZHONGOLOVITCH,I. D. : 1960, Bull. Institute Theoretical Astron. Leningrad 7, 743. ZHONGOLOVITCH,L D. and PELLINEN, L. P. : 1962, Bull. Institute Theoretical Astron. Leningrad 8, 381.