Vinti

I. INTRODUCTION

Recent development of altimeters with a capability of measuring the distance from an artificial satellite to a station position within one meter suggests the use of satellite orbit as a base line for surveying the earth or the oceans. There then arises the need of accurate orbital prediction with comparable accuracy and this leads further to an effort of representing the gravitational potential of the earth by a rapidly converging series, so that a satisfactory accuracy can be obtained through an efficient computation.

As is known, the gravitational potential of the earth can be expressed by using spherical harmonic expansion as:

\[U = \frac{\mu}{r} \left[ 1 + \sum_{n=2}^{\infty} \sum_{m=0}^{n} \left( \frac{a}{r} \right)^{n} A_{nm} \cos m(\emptyset - \overline{\emptyset}_{nm}) \overline{P}_{nm} (\cos \theta) \right]\] where \(\overline{A}_{nm}\) is the harmonic coefficients. Written in this form, the spherical harmonic expansion representation of the earth requires a large number of terms; its normalized coefficients decrease rapidly in magnitude from (2,2) to about (8,8), and thereafter they tend to oscillate with very slowly decreasing amplitude (Kohnlein, 1967). A possible geopotential improvement by the use of satellite altimetry was suggested by Lundquist and Giacaglia (1968). These authors obtain a most compact representation for the gravitational potential of the earth by an alternate set of function with the property that only relatively few terms of this set provide the significant contribution to the potential in any region.

Another representation of the gravitational potential of the earth was proposed by Vinti, by considering the potential of a finite body as being divided into two parts. The first part, \(V_{D}\), is represented by a summation of series that can be efficiently and accurately evaluated. One possible representation is the summation of a finite number of spherical harmonics, whose coefficients can be accurately obtained from satellite measurements by standard techniques. Another alternative is representing \(V_{\rm D}\) by the spheroidal potential derived by Vinti (1959, 1961, 1965). The coefficients of the spheroidal harmonic expansion of the potential are rapidly converging, so that truncation after a few terms would result in a high degree of accuracy. The second part, \(\delta V\), which accounts for the remainder of the potential, is represented as a surface distribution \(\sigma\) of mass density of variable sign over a sphere S just clearing the earth. The intent here is to find a representation of this surface distribution which reduces the number of parameters significantly as compared with a spherical harmonic expansion of the same accuracy.

Originally, Vinti proposed to obtain \(\sigma\) from worldwide measurements of surface gravity. The quantity of itself is essentially a measure of the density anomalies near the surface of the earth. It will be quite a while before the necessary gravity data become available since terrestrial gravity data are required at various unsurveyed areas, although by use of satellite altimetry, it is possible to obtain the gravity data over the ocean. It is then proposed to try an interim approach, namely by trying to find o over the sphere just clearing the earth from satellite determinations of the higher harmonic coefficients from various published sources (e.g. Kaula, 1966, A P L 1965, S A O 1967). By comparing contour maps of o obtained from the use of various data, some degree of confidence on such representation can be established. In the present report, the results of the procedure just outlined are presented and discussed. Further discussion on the representation of the earth's gravitational disturbance potential by means of surface density distribution can be found in a recent publication by Vinti (1971).

II. REPRESENTATION OF σ IN SPHERICAL HARMONICS

The gravitational potential V of the earth can be expressed in spherical harmonics as:

\[V = -\frac{\mu}{r} \left\{ 1 - \sum_{n=2}^{\infty} \left( \frac{r_e}{r} \right)^n \left[ J_n P_n(\cos \theta) - \sum_{m=1}^n P_n^m(\cos \theta) \right] \right\}\]

\[(C_{nm}\cos m \phi + S_{nm}\sin m \phi)]\] (1)

where:

\(\text{[rumus tidak dapat ditampilkan dengan baik — lihat PDF asli]}\)

\(C_{nm}\), \(S_{nm}\) - tesseral harmonic coefficients

Figure 1: Coordinate System

It should be remarked that the convergence of (1) at the surface of the earth is not insured (Morritz, 1962). Let the gravitational potential of the earth, following Vinti, be represented in the form:

\[V = V_p + \delta V \tag{2}\] where:

\(V_p\) = some reference potential and

\[\delta V = - \int_{\text{Sphere S}} \frac{G\sigma(\theta', \emptyset') \ d \ S'}{|\underline{r} - \underline{r}'|}\] of radius a (3)

where:

• \(\underline{\underline{r}}\) field point radial distance from the center of the earth.
• \(\underline{\underline{r}}\) source point radial distance from the center of the earth.

ï

The choice of a reference potential \(V_p\) is particularly dictated by the accuracy and effectiveness of the computation, and is at our disposal. One could, for example, represent \(V_p\) as the sum of a finite number of spherical harmonics, whose coefficients can be accurately obtained from satellite measurements. The intent here is that \(V_p\) should be represented by a sum of finite terms in a series expansion, which is rapidly converging.

In a recent development, Vinti (1959, 1961, 1965, 1966) has found the coefficients of a spheroidal potential of the earth's gravitational field. The gravitational potential \(\rm V_S\) outside of the earth, expanded in spherical harmonics, can then be represented in oblate spheroidal coordinate as:

\[V_{S} = -\frac{\mu \rho}{\rho^2 - cn^2} \tag{4a}\] or in spherical harmonics as:

\[V_{s} = -\frac{\mu}{r} \left[ 1 - \sum_{n=2}^{\infty} \left( \frac{r_{e}}{r} \right)^{n} J_{n}^{x} P_{n}(\cos \theta) \right]\] (4)

If here the earth is regarded as being symmetrical with respect to the equatorial plane, all the odd harmonics drop out, that is, \(J_1 = J_3 = J_5 \ldots = 0\). Deriving \(V_s\) as a solution of Laplace equation in oblate spheroidal coordinates, with axial symmetry and leading to exact separability of the Hamilton Jacobi's equation for the motion of a satellite, the following relations for \(J_n^x\) were obtained:

\[J_n^x = (-1)^{\frac{n}{2}} + 1 \quad J_2^{\frac{n}{2}}\] for \(n = 2, 4 \dots\) (even)
\(J_n^x = 0\) for \(n = 1, 3, 5 \dots\) (odd)

In the derivation, the center of mass was used as the origin.

Clearly \(J_n^x\) diminishes rapidly with increasing n, since \(J_2\stackrel{\sim}{=} 1.08 \times 10^{-3}\); consequently \(V_s\) is rapidly converging and suitable for being used as the reference potential \(V_p\). A remark is in order at this point.

The disturbance potential or the "potential tail", \(\delta V\), is represented as a surface integral of a fictitous surface density. This idea is suggested from the dependency of the higher harmonic coefficients on the situation near the crust of the earth, as can be concluded from the following relations:

\[J_{n} = -\frac{1}{M} \int \int \int \left(\frac{\mathbf{r}}{\mathbf{r}_{o}}\right)^{n} P_{n}(\cos \theta) \rho d \tau\] (6)

\[C_{n,m} + iS_{n,m} = \frac{2}{M} \frac{(n-m)!}{(n+m)!} \iiint_{earth} (\frac{r}{r_e})^n P_n^m(\cos \theta) e^{im\phi} \rho d \tau \qquad (7)\] where:

1 - mass of the earth

ρ - density of the earth

d T - element of volume

We can see that for large values of n the main contributions to \(J_n\), \(C_{n,m}\) and \(S_{n,m}\) come from the crust since the factor \((\frac{r}{r_e})^n\) in the integrand changes slowly there, hence the idea of surface density distribution to represent this "potential tail".

It is customary to write:

\[V = -\frac{\mu}{r} \left[1 + \sum_{n=2}^{\infty} \sum_{m=0}^{n} \left(\frac{r_e}{r}\right)^n P_n^m(\cos \theta) \left(C_{nm} \cos m \phi\right) + C_{nm} \right]\]

\[S_{nm} \sin m \phi)\] (8)

so that:

\[J_{n} = -C_{no}\] (9)

Thus:

\[J_n^x = C_{20}^{\frac{n}{2}} \text{ for n even}\] (5a)

Letting \(V_p = V_s\), from (4) and (1) we obtain:

\[\delta V = -\frac{\mu}{r} \left[ \sum_{n=2}^{\infty} \left( \frac{r_e}{r} \right)^n J_n^{XX} P_n(\cos \theta) + \sum_{n=2}^{\infty} \left( \frac{r_e}{r} \right)^n \sum_{m=0}^{n-1} P_n^m(\cos \theta) (C_{nm} \cos m \phi + S_{nm} \sin m \phi) \right]\](10)

where:

\[J_n^{XX} = C_{no}\], for n odd \(J_n^{XX} = C_{no} - (C_{20})^{\frac{n}{2}}\), for n even (11)

Appendix A describes the inversion of (3), if \(\delta V\) is given by (10). The result is given by:

\[G\sigma = \frac{\mu}{4\pi a^2} \left\{ \sum_{n=2}^{\infty} (2n+1) \left(\frac{r_e}{a}\right)^n \left[ J_n^{XX} P_n(\cos \emptyset) + \sum_{m=0}^{\infty} P_n^m(\cos \theta) \left( C_{nm} \cos m \emptyset + S_{nm} \sin m \emptyset \right) \right] \right\}\](12)

This is the relation that we seek in determining the surface density distribution from satellite determination of the spherical harmonic coefficients.

It can be added, 'that if \(\sigma\) (\(\theta\),\(\phi\)) can be determined from other methods, the higher harmonic coefficients of the potential can be obtained by inverting equation (8). Using the orthogonality properties of the spherical harmonics, the following results can be readily derived:

\[J_{n}^{xx} = \frac{a^{2}}{\mu} \left(\frac{a}{r_{e}}\right)^{n} \int_{0}^{\pi} P_{n} \cos(\theta) \sin \theta d\theta \int_{0}^{2\pi} G\sigma(\theta, \phi) d\phi\] (13)

\[C_{nm} + i S_{nm} = \frac{2a^2}{\mu} \left(\frac{a}{r_e}\right)^n \frac{(n-m)!}{(n+m)!} \int_0^{\pi} P_n^m(\cos \theta) \sin \theta d\theta\] \[\int_0^{2\pi} G\sigma(\theta, \phi) e^{im\phi} d\phi \qquad (14)\] \[m = 1, 2, \dots, n\]

III. EARTH PARAMETERS AND POTENTIAL COEFFICIENTS

The data to be used in this study are the harmonic coefficients recently determined by various authors. The semimajor axis, gravitational constant and flattening consistent with each set of coefficients were also used. The main purpose of using these different sets of earth parameters and potential coefficients is two fold: to assess the qualitative

consistency of the resulting contour plot of the surface denslty dlstributlon and to investigate the feaslbility of the proposed representatlon of surface density.

The data sets used were:

• (a) seodata (L967)
• (b) Rapp data (1) (L967)
• Cc) napp data (2) (1968)
• (d)ALPdara (1965)

Data set Ca) ls based on the satel-l-lte deterrnlnatlon by Kozal of the zonaL harmonics up to J20, and the determination by KohnleLn of the tesseral ha::nonLcs complete to degree 15, whlch are further based on Gaposhkinrs satel-lite data and Kaulars gravlty data.

u {

lt, t

Data set (b) is the surface data obtalned by Rapp (L967), which ls complete onJ-y to the 8th order and degree. Data set (c) fs obtalned by Rapp (L967, L96L, L969) fron a combLnation of satelllte data of the S A O standard earth (l-966), terrestrlal data and model anomalles. Data set (d) ls obtalned by Guier an Nenton. Using orbital determinations and data reslduals from the Doppler tracklng of flve satellltes, they then deduced for the zonaL iharmonlcs of the earthrs gravlty field of odd degree through the ninth, the nonzonal harmonlcs of all degrees fron the second through the eighth, and the sectorlal harmonics of the thlrteenth degree and order.

It ls iroted that fron comparLson wlth terrestrlal gravlty data, Kaula (L966) has chosen Gaposhklnrs solutlon, whlch ls the basls of Kohnl-elnfs harmonlc coefficlents, as one that best agrees with terrestrlal gravlty data. Kohnl-elnfs tesseral harnonlc coefficlents codrprlse the most complete set at the present time. Rapp (14114) data has smaller and more conslstent Ftandard devlation. The earth par€rmeters and potentlal coefflcients for each set of these data are tabulated in table 1 below.

Table I

IV. DISCUSSION OF RESULTS

Figures 1 to 4 show the contour map of \(\emptyset\) = Go, here being multiplied by \(10^6\), in m/sec² using four sets of data of Rapp (8,8), Rapp (14,14), A P L (8,8) and Kohnlein (15,15). Numerical values are tabulated in tables 2 and 3. By incorporating nine higher harmonic coefficients, the contour map displays an increase of the number of extrema. For example, corresponding to data set of Rapp (8,8), there are 21 extrema (11 minima and 10 maxima) and to data sets of Rapp (14,14) there are 45 extrema. In addition, there exist some similarities between the contour map of G, of the gravity anomalies and of the geoidal undulation. The relationship between G and the gravity anomalies can be observed from the following relations:

\[4\pi G\sigma = 2 \delta g_{ra} + \left(\frac{\delta V}{a}\right)_{r=a}\] (15)

where:

\[\delta g_r = \left[\frac{\partial}{\partial r} \left( V - V_p \right) \right]_{r=a} + \omega^2 a \sin^2 \theta\] (16)

which was derived by Vinti in an unpublished note. Now

\[\Delta g = \left[ -\frac{\partial T}{\partial r} - \frac{2}{r} T \right]_{r=a}\] (17)

where:

\[T = V - V_{ell}\] (18)

Essentially \(V_{ell} \stackrel{\sim}{=} V_s\), to a very high degree of approximation (in the order of \(J_3\)). These equations clearly describe a direct relation between GO and \(\Delta g\). Note also that the term \((\frac{\delta V}{a})\) contributes to about 12% of the first term on the right hand side of equation (15).

Obviously, the surface distribution \(\sigma\), which characterises the potential tail, reflects the irregularities of the earth near its crust. S A O data yields rms value of G\(\sigma\) of 16,8 x \(10^{-6}\) m/sec² while Rapp data gives 16,7 x \(10^{-6}\) m/sec². The above values are based on:

\[V_p = V_s + J_3(\frac{r_e}{a})^3 P_3(\cos \theta)\] (19)

and letting \(a = r_e\). \(J_2^x\) was taken to be the average value of the \(J_2\) given by Rapp and S A O. By choosing \(V_p = V_s\), numerical computation yields contour maps that do not differ markedly from the ones described above, except near the poles.

The contour maps shown in figures 1 to 2 indicate that all these data sets yield the presence of main patterns that exist in the geoid map, such as the New Guinea high, Indian low, Puerto Rican low, United States West Coast low, Gibraltar high, etc.

To illustrate the merit of the results, consider the following figures. Since \(G = 6.66 \times 10^{-8} \text{ cm}^3/\text{g} - \text{sec}^2\), and the rms value of GO is about \(17 \times 10^{-6} \text{ m/sec}^2\), then one obtains

\[\sigma_{\rm rms} = 25,500 \text{ g/cm}^2\]

If one assumes that the density anomalies which contribute to the higher coefficients of the gravitational potential occur in the crust of thickness b, then the rms value of density anomalies is given by:

\[\Delta \rho = \sigma/b\]

Since the mean value of b = 30 km (Heiskanen and Moritz, 1967), then it follows that:

\[\Delta \rho_{\rm rms} = \frac{\sigma_{\rm rms}}{b} = \frac{25,500}{3 \times 10^6} = 0.0085 \text{ g/cm}^3\]

It can be observed, that \(|\sigma G \times 10^{-6}| < 100\), hence \(\Delta \rho < 0.05\) g/cm³ for b = 30 km, and \(\Delta \rho\) < 0.15 g/cm³ for b = 10 km.

It is of interest to see the contour map near the poles. In this purpose, polar maps were drawn, as shown by figures 3 to 6, which were based on Rapp (14,14) data and S A O (15,15) data. Fairly good agreement between these two data are shown around the North pole, although it is not true around the In general, however, we can conclude that the two South pole. sets of data of Rapp and S A O yield Go distribution with an agreement better than qualitative, and the idea of employing Go as a basis for predicting accurate potential tail seems to be feasible. As a further step, one may be inclined to employ a finite set of lumped mass rather than continuous surface density distribution, leading to the idea of mascons. Some relationship between the mass density distribution of the moon of some spherical surface beneath the moon surface to the location of the mascons was suggested recently by Campbell, O'Leary and Sagan. Such an idea of representing the surface

density irregularities by mascons and using these to predict the potential tail is currently studied by Vinti and Price.

Thus if one could predict the surface density distribution by using satellite altimetry and terrestrial gravity data, more accurate and efficient prediction of the potential tail is possible. Future work will be devoted to evaluation of \(\sigma\) from satellite altimetry and terrestrial data and analytical evaluation of \(\sigma\) to obtain the corrective integral.

V. ACKNOWLEDGEMENT

The author would like to acknowledge Dr. J.P. Vinti, who initiated and supervised the present investigation, Dr. S.J. Madden, for fruitful discussion and Mr. Carl Goodwin for various assistance. This work was performed at the Experimental Astronomy Laboratory (now Measurement System Laboratory) at the Massachusetts Institute of Technology and sponsored by NASA Grant N G R - 22 - 009 - 262 in 1968 - 1969.

APPENDIX A: RELATION BETWEEN G AND V

The potential tail \(\delta V\) is given as a corrective integral of \(\sigma\) as follows:

\[\delta V = -G \int_{S} \frac{\sigma |\theta', \phi'|}{|\underline{r} - \underline{r}'|} d s'\]

Now let:

\[\sigma(\theta', \emptyset') = \sum_{n=0}^{\infty} \sum_{m=0}^{\infty} (A_{nm} \cos m \emptyset' + B_{nm} \sin m \emptyset') P_n^m (\cos \theta') (\frac{r_e}{a})^n\]

Furthermore, \(\frac{1}{|\underline{r} - \underline{r}'|}\) is the generating function of the Legendre polynomial:

\[\frac{1}{\left|\underline{\mathbf{r}}-\underline{\mathbf{r}}'\right|} = \frac{1}{\left|\mathbf{r}\right| \sqrt{1 - \frac{2\mathbf{r}'}{\mathbf{r}}} \cos \psi + (\frac{\mathbf{r}'}{\mathbf{r}})^{2}}\]\[= \frac{1}{\mathbf{r}} \sum_{n=0}^{00} (\frac{\mathbf{r}'}{\mathbf{r}})^{n} P_{n}(\cos \psi)\] for:

\[\left|\frac{r'}{r}\right|\] less than the smaller of \(\left|\frac{r'}{r}\right| \pm \sqrt{\left(\frac{r'}{r}\right)^2 - 1}\)

with r > r', r' being the radius of the sphere (=a) just clearing the earth and r is outside of this sphere. Here \(\psi\) is the angle between the vectors \(\underline{r}'(r', \theta', \emptyset')\) and \(\underline{r}(r, \theta, \emptyset)\). From the addition theorem (see figure A)

\[\begin{split} P_{n}(\cos \psi) &= P_{n}(\cos \theta) \ P_{n}(\cos \theta') + \\ &2 \sum_{m=1}^{n} \frac{(n-m)!}{(n+m)!} P_{n}^{m} (\cos \theta) \ P_{n}^{m}(\cos \theta') \\ &\{\cos m \phi \cos m \phi' + \sin m \phi \sin m \phi'\} \end{split}\]

Furthermore:

\[dS = a^2 \sin \theta' d\theta' d\phi'\]

Figure A

Hence:

\[\delta V = -\frac{a^2 G}{r} \int_0^T \int_0^T \sum_{n=0}^\infty \sum_{m=0}^\infty P_n^m(\cos \theta) (A_{nm} \cos m \phi' + B_{nm} \sin m \phi') \cdot \sum_{n=0}^\infty \left(\frac{a}{r}\right)^n (P_n(\cos \theta) P_n(\cos \theta') + 2 \sum_{m=1}^n \frac{(n-m)!}{(n+m)!} P_n^m(\cos \theta) P_m^m(\cos \theta')\] \[\{\cos m \phi \cos m \phi' + \sin m \phi \sin m \phi'\} \sin \theta' d \theta' d \phi'\]

Using the orthogonality relations:

we then obtain:

\[\delta V = -\frac{a^2 G}{r} \left[ \sum_{n=0}^{\infty} \left( \frac{a}{r} \right)^n \left( \frac{r_e}{a} \right)^n \frac{4\pi}{2n+1} \left\{ P_n(\cos \theta) A_{no} + \sum_{m=0}^{n} P_n^m(\cos \theta) (A_{nm} \cos m \phi + B_{nm} \sin m \phi) \right\} \right]\]

Now:

\[\delta V = -\frac{\mu}{r} \left[ \sum_{n=2}^{\infty} \left( \frac{a}{r} \right)^n J_n^{xx} P_n(\cos \theta) + \frac{c}{r} \left( \frac{a}{r} \right)^n \sum_{m=0}^{n} P_n^m(\cos \theta) (C_{nm} \cos m \emptyset + S_{nm} \sin m \emptyset) \right]\]

Comparing coefficients of \(P_n^m(\cos\,\theta)\text{, we obtain:}\)

\[\left(\frac{r_e}{a}\right)^n A_{no} = + \frac{\mu}{4\pi a^2 G} (2n + 1) J_n^{xx}\]

\[\left(\frac{r_e}{a}\right)^n (A_{nm} + iB_{nm}) = (C_{nm} + iS_{nm}) \frac{\mu}{4\pi a_G^2} (2n + 1)\]

Hence:

\[\begin{split} \sigma &= + \frac{M}{4\pi a^2} \sum_{n=2}^{00} (2n+1) \left\{ J_n^{XX} P_n(\cos\theta) + P_n^m(\cos\theta) (C_{nm}^{\cos\theta} \cos \theta) + S_{nm}^{\sin\theta} \sin \theta \right\} \left( \frac{r_e}{a} \right)^n \end{split}\]

APPENDIX B: GENERATION OF LEGENDRE FUNCTIONS AND ASSOCIATED LEGENDRE FUNCTIONS

The Legendre function \(P_n(x)\) and the associated Legendre function \(P_n^m(x)\) are defined respectively by:

\[P_n(x) = \frac{1}{2^n n!} \frac{d^n}{dx^n} (1 - x^2)^n\] B.1 and

\[P_{n}^{m}(x) = (1 - x^{2})^{\frac{m}{2}} \frac{d^{n}}{dx^{m}} P_{n}(x)\] \[= \frac{(1 - x^{2})^{\frac{m}{2}}}{2^{n} n!} \frac{d^{n+m}}{dx^{n+m}} (1 - x^{2})^{n}\] B.2 for \[|x| \le 1\]

These functions are of the first kind. The Legendre function \(P_{n}\) may be generated from the following recursion formula:

\[P_n(x) = \frac{2n-1}{n} \times P_{n-1}(x) - \frac{n-1}{n} P_{n-2}(x)\] B.3 with:

\[P_1(x) = x\] \[P_0(x) = 1\]

To generate the associated Legendre functions of the first kind, the following recursion formula may be used:

\[P_{n}^{m}(x) = x P_{n-1}^{m} + (n+m-1)(1-x^{2})^{\frac{1}{2}} P_{n-1}^{m-1}\] \[for 0 \le m < n\] \[|x| \le 1\] while \(P_n^n\) can be calculated from:

\[P_{n}^{n}(x) = \frac{(2n)!}{2^{n}n!} (1 - x^{2})^{\frac{n}{2}}\] \[= 1.3. \dots (2n - 1)(1 - x^{2})^{\frac{n}{2}}\] B.5a or

\[P_n^n(x) = (2n - 1)(1 - x^2)^{\frac{1}{2}} P_{n-1}^{n-1}(x)\] B.5b for n > o

The formula A.5a is a special case of the hypergeometric series representation of \(P_n^m\):

\[P_n^m(x) = \frac{(n+m+1)}{2^m m! \Gamma(n-m+1)} (1-x^2)^{\frac{m}{2}} F(m-n, m+n+1, m+1, \frac{1-x}{2}) \qquad B.6\]

Several of these recursion formula can be found in Erdelyi et al.; the recursion formula B.4 has a property of being bounded everywhere between \(-1 \le x \le 1\).

For each n and m, the recursion formula B.4 may be started if \(P_m^m\) and \(P_{n-m}^o\) are known from B.5 and B.3, respectively, as can be seen from the diagram shown in figure B.

This fact is useful in performing an analysis to estimate the error in computing the Legendre functions.

Figure B

For the tesseral we nay use:

\[(n-m+1)P_{n+1}^{m}(x) = (2n+1)x P_{n}^{m}(x) - (n+m)P_{n-1}^{m}(x)\] and \[P_{n}^{m}(x) = 2(m-1)x(1-x^{2})^{\frac{1}{2}} P_{n}^{m-1}(x) - (n-m+2)(n+m-1)P_{n}^{m-2}(x)\] B.8