Geopotential model

In geophysics, a geopotential model is the theoretical analysis of measuring and calculating the effects of Earth's gravitational field.

Newton's law

Diagram of two masses attracting one another

Newton's law of universal gravitation states that the gravitational force F acting between two point masses m1 and m2 with centre of mass separation r is given by

\mathbf{F} = - G \frac{m_1 m_2}{r^2}\mathbf{\hat{r}}

where G is the gravitational constant and is the radial unit vector. For an object of continuous mass distribution, each mass element dm can be treated as a point mass, so the volume integral over the extent of the object gives:


\mathbf{\bar{F}} = - G \int\limits_V \frac{\rho }{r^2}\mathbf{\hat{r}}\,dx\,dy\,dz

 

 

 

 

(1)

with corresponding gravitational potential


u\ =\ -\int\limits_V \rho \frac{G}{r}\, dx\,dy\,dz

 

 

 

 

(2)

where ρ = ρ(x, y, z) is the mass density at the volume element and of the direction from the volume element to the point mass.

The case of a homogeneous sphere

In the special case of a sphere with a spherically symmetric mass density then ρ = ρ(s), i.e. density depends only on the radial distance

s = \sqrt{x^2+y^2+z^2} \,.

These integrals can be evaluated analytically. This is the shell theorem saying that in this case:


\bar{F} = -\frac{GM}{R^2}\ \hat{r}

 

 

 

 

(3)

with corresponding potential


u = -\frac{GM}{r}

 

 

 

 

(4)

where M = ∫Vρ(s)dxdydz is the total mass of the sphere.

The deviations of Earth's gravitational field from that of a homogeneous sphere

In reality, Earth is not exactly spherical, mainly because of its rotation around the polar axis that makes its shape slightly oblate. If this shape would have been perfectly known together with the exact mass density ρ = ρ(x, y, z) the integrals (1) and (2) could have been evaluated with numerical methods to find a more accurate model for Earth's gravitational field, but the situation is in fact the opposite. By observing the orbits of spacecraft and the Moon, Earth's gravitational field can be determined quite accurately and the best estimate of Earth's mass is obtained by dividing the product GM as determined from the analysis of spacecraft orbit with a value for G determined to a lower relative accuracy using other physical methods.

From the defining equations (1) and (2) it is clear (taking the partial derivatives of the integrand) that outside the body in empty space the following differential equations are valid for the field caused by the body:


\frac{\partial F_x }{\partial x} + \frac{\partial F_y }{\partial y} + \frac{\partial F_z }{\partial z} = 0

 

 

 

 

(5)


\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} = 0

 

 

 

 

(6)

Functions of the form \phi\ =\ R(r)\ \Theta(\theta)\ \Phi(\varphi) where (r, θ, φ) are the spherical coordinates which satisfy the partial differential equation (6) (the Laplace equation) are called spherical harmonic function.

They take the forms:


\begin{align}
g(x,y,z) & =\frac{1}{r^{n+1}} P^m_n(\sin \theta) \cos m\varphi \,,&\quad 0 \le m \le n \,,&\quad n=0,1,2,\dots \\
h(x,y,z) & =\frac{1}{r^{n+1}}  P^m_n(\sin \theta) \sin m\varphi \,,&\quad 1 \le m \le n \,,&\quad n=1,2,\dots
\end{align}

 

 

 

 

(7)

where spherical coordinates (r, θ, φ) are used, given here in terms of cartesian (x, y, z) for reference:


\begin{align}
&x = r \cos \theta \cos \varphi \\
&y = r \cos \theta \sin \varphi \\
&z = r \sin \theta\,,
\end{align}

 

 

 

 

(8)

also P0n are the Legendre polynomials and Pmn for 1 ≤ mn are the associated Legendre functions.

The first spherical harmonics with n = 0,1,2,3 are presented in the table below.

n Spherical harmonics
  0   \frac{1}{r}
  1   \frac{1}{r^2} P^0_1(\sin\theta) = \frac{1}{r^2} \sin\theta
\frac{1}{r^2} P^1_1(\sin\theta) \cos\varphi= \frac{1}{r^2} \cos\theta \cos\varphi
\frac{1}{r^2} P^1_1(\sin\theta) \sin\varphi= \frac{1}{r^2} \cos\theta \sin\varphi
  2   \frac{1}{r^3} P^0_2(\sin\theta) = \frac{1}{r^3} \frac{1}{2} (3\sin^2\theta -1)
\frac{1}{r^3} P^1_2(\sin\theta) \cos\varphi = \frac{1}{r^3} 3 \sin\theta \cos\theta\ \cos\varphi
\frac{1}{r^3} P^1_2(\sin\theta) \sin\varphi = \frac{1}{r^3} 3 \sin\theta \cos\theta \sin\varphi
\frac{1}{r^3} P^2_2(\sin\theta) \cos2\varphi = \frac{1}{r^3} 3 \cos^2 \theta\ \cos2\varphi
\frac{1}{r^3} P^2_2(\sin\theta) \sin2\varphi = \frac{1}{r^3} 3 \cos^2 \theta \sin 2\varphi
  3   \frac{1}{r^4} P^0_3(\sin\theta) = \frac{1}{r^4} \frac{1}{2} \sin\theta\ (5\sin^2\theta -3)
\frac{1}{r^4} P^1_3(\sin\theta) \cos\varphi = \frac{1}{r^4} \frac{3}{2}\ (5\ \sin^2\theta - 1) \cos\theta \cos\varphi
\frac{1}{r^4} P^1_3(\sin\theta) \sin\varphi = \frac{1}{r^4} \frac{3}{2}\ (5\ \sin^2\theta - 1) \cos\theta \sin\varphi
\frac{1}{r^4} P^2_3(\sin\theta) \cos 2\varphi = \frac{1}{r^4} 15 \sin\theta \cos^2 \theta \cos 2\varphi
\frac{1}{r^4} P^2_3(\sin\theta) \sin 2\varphi = \frac{1}{r^4} 15 \sin\theta \cos^2 \theta \sin 2\varphi
\frac{1}{r^4} P^3_3(\sin\theta) \cos 3\varphi = \frac{1}{r^4} 15 \cos^3 \theta \cos 3\varphi
\frac{1}{r^4} P^3_3(\sin\theta) \sin 3\varphi = \frac{1}{r^4} 15 \cos^3 \theta \sin 3\varphi

The model for Earth's gravitational potential is a sum


u = -\frac{\mu }{r} + \sum_{n=2}^{N_z}  \frac{J_n P^0_n(\sin\theta) }{r^{n+1}} + \sum_{n=2}^{N_t} \sum_{m=1}^n \frac{ P^m_n(\sin\theta) (C_n^m \cos m\varphi + S_n^m \sin m\varphi)}{r^{n+1}}

 

 

 

 

(9)

where \mu=GM and the coordinates (8) are relative the standard geodetic reference system extended into space with origin in the center of the reference ellipsoid and with z-axis in the direction of the polar axis.

The zonal terms refer to terms of the form:

\frac{P^0_n(\sin\theta) }{r^{n+1}} \quad n=0,1,2,\dots

and the tesseral terms terms refer to terms of the form:

\frac{ P^m_n(\sin\theta) \cos m\varphi}{r^{n+1}}\,, \quad 1 \le m \le n \quad n=1,2,\dots
\frac{ P^m_n(\sin\theta) \sin m\varphi}{r^{n+1}}

The zonal and tesseral terms for n = 1 are left out in (9).

The different coefficients Jn, Cnm, Snm, are then given the values for which the best possible agreement between the computed and the observed spacecraft orbits is obtained.

As P0n(x) = −P0n(−x) non-zero coefficients Jn for odd n correspond to a lack of symmetry "north–south" relative the equatorial plane for the mass distribution of Earth. Non-zero coefficients Cnm, Snm correspond to a lack of rotational symmetry around the polar axis for the mass distribution of Earth, i.e. to a "tri-axiality" of Earth.

For large values of n the coefficients above (that are divided by r(n + 1) in (9)) take very large values when for example kilometers and seconds are used as units. In the literature it is common to introduce some arbitrary "reference radius" R close to Earth's radius and to work with the dimensionless coefficients

\tilde{J_n} = -\frac{J_n}{\mu\ R^n}
\tilde{C_{n}^m} = -\frac{C_{n}^m}{\mu\ R^n}
\tilde{S_{n}^m} = -\frac{S_{n}^m}{\mu\ R^n}

and to write the potential as


u = -\frac{\mu }{r} \left(1 + \sum_{n=2}^{N_z}  \frac{\tilde{J_n} P^0_n(\sin\theta) }{{(\frac{r}{R})}^n} + \sum_{n=2}^{N_t} \sum_{m=1}^n \frac{ P^m_n(\sin\theta) (\tilde{C_{n}^m} \cos m\varphi + \tilde{S_{n}^m} \sin m\varphi)}{{(\frac{r}{R})}^n}\right)

 

 

 

 

(10)

The dominating term (after the term −μ/r) in (9) is the "J2 term":

u = \frac{J_2\ P^0_2(\sin\theta)}{r^3} = J_2 \frac{1}{r^3} \frac{1}{2} (3\sin^2\theta -1) = J_2 \frac{1}{r^5} \frac{1}{2} (3 z^2 -r^2)

Relative the coordinate system


\begin{align}
&\hat{\varphi}=-\sin \varphi \hat{x} + \cos \varphi \hat{y} \\
&\hat{\theta}=-\sin \theta\ (\cos \varphi\ \hat{x} + \sin \varphi \hat{y})+ \cos\theta \hat{z} \\
&\hat{r}= \cos \theta\ (\cos \varphi\ \hat{x}\ +\ \sin \varphi\ \hat{y})+\ \sin\theta\ \hat{z}
\end{align}

 

 

 

 

(11)

Figure 1: The unit vectors \hat{\varphi}\ ,\ \hat{\theta}\ ,\ \hat{r}

illustrated in figure 1 the components of the force caused by the "J2 term" are


\begin{align}
&F_\theta = -\frac{1}{r}\ \frac{\partial u }{\partial \theta} = -J_2\ \frac{1}{r^4} 3\ \cos\theta\ \sin\theta \\
&F_r = -\frac{\partial u }{\partial r} = J_2\ \frac{1}{r^4} \frac{3}{2}\ \left(3\sin^2\theta\ -\ 1\right)
\end{align}

 

 

 

 

(12)

In the rectangular coordinate system (x, y, z) with unit vectors (x̂ ŷ ẑ) the force components are:


\begin{align}
&F_x = -\frac{\partial u }{\partial x} = J_2\ \frac{x}{r^7} \left(6z^2 - \frac{3}{2} (x^2 + y^2)\right) \\
&F_y = -\frac{\partial u }{\partial y} = J_2\ \frac{y}{r^7} \left(6z^2 - \frac{3}{2} (x^2 + y^2)\right) \\
&F_z = -\frac{\partial u }{\partial z} = J_2\ \frac{z}{r^7} \left(3z^2 - \frac{9}{2} (x^2 + y^2)\right)
\end{align}

 

 

 

 

(13)

The components of the force corresponding to the "J3 term"

u = \frac{J_3 P^0_3(\sin\theta) }{r^4} = J_3 \frac{1}{r^4} \frac{1}{2} \sin\theta (5\sin^2\theta -3) = J_3 \frac{1}{r^7} \frac{1}{2} z (5 z^2 - 3 r^2)

are


\begin{align}
&F_\theta = -\frac{1}{r} \frac{\partial u }{\partial \theta} = -J_3 \frac{1}{r^5} \frac{3}{2} \cos\theta \left(5 \sin^2\theta -1\right) \\
&F_r = -\frac{\partial u }{\partial r} = J_3 \frac{1}{r^5} 2 \sin\theta \left(5\sin^2\theta - 3\right)
\end{align}

 

 

 

 

(14)

and


\begin{align}
&F_x = -\frac{\partial u }{\partial x} = J_3 \frac{x z}{r^9} \left(10 z^2 - \frac{15}{2} (x^2 + y^2)\right) \\
&F_y = -\frac{\partial u }{\partial y} = J_3 \frac{y z}{r^9} \left(10 z^2 - \frac{15}{2} (x^2 + y^2)\right) \\
&F_z = -\frac{\partial u }{\partial z} = J_3 \frac{1}{r^9} \left(4 z^2\ \left( z^2 - 3 (x^2 + y^2)\right) + \frac{3}{2} (x^2 + y^2)^2\right)
\end{align}

 

 

 

 

(15)

The exact numerical values for the coefficients deviate (somewhat) between different Earth models but for the lowest coefficients they all agree almost exactly.

For JGM-3 the values are:

μ = 398600.440 km3⋅s−2
J2 = 1.7555 × 1010 km5⋅s−2
J3 = −2.619 × 1011 km6⋅s−2

With a "reference radius" R of 6378.1363 km corresponding dimensionless parameters are

\tilde{J_2} = 1.0826 \times 10^{-3}
\tilde{J_3} = -2.532 \times 10^{-6}

For example, at a radius of 6600 km (about 200 km above Earth's surface) J3/(J2r) is about 0.002, i.e the correction to the "J2 force" from the "J3 term" is in the order of 2 permille. The negative value of J3 implies that for a point mass in Earth's equatorial plane the gravitational force is tilted slightly towards the south due to the lack of symmetry for the mass distribution of Earth's "north–south".

Recursive algorithms used for the numerical propagation of spacecraft orbits

Spacecraft orbits are computed by the numerical integration of the equation of motion. For this the gravitational force, i.e. the gradient of the potential, must be computed. Efficient recursive algorithms have been designed to compute the gravitational force for any N_z and N_t and such algorithms are used in standard orbit propagation software

Available models

The earliest Earth models in general use by NASA and ESRO/ESA were the "Goddard Earth Models" developed by Goddard Space Flight Center denoted "GEM-1", "GEM-2", "GEM-3", and so on. Later the "Joint Earth Gravity Models" denoted "JGM-1", "JGM-2", "JGM-3" developed by Goddard Space Flight Center in cooperation with universities and private companies became available. The newer models generally provided higher order terms than their precursors. The EGM96 uses Nz = Nt = 360 resulting in 130317 coefficients.

For a normal Earth satellite for which an orbit determination/prediction accuracy of a few meters is sufficient the "JGM-3" truncated to Nz = Nt = 36 (1365 coefficients) is usually sufficient. Inaccuracies from the modeling of the air-drag and to a lesser extent the solar radiation pressure will exceed the inaccuracies caused by the gravitation modeling errors.

Spherical harmonics

Main article: Spherical harmonics

The following is a compact account of the spherical harmonics used to model Earth's gravitational field. The spherical harmonics are derived from the approach of looking for harmonic functions of the form


\phi\ =\ R(r)\ \Theta(\theta)\ \Phi(\varphi)

 

 

 

 

(16)

where (r, θ, φ) are the spherical coordinates defined by the equations (8). By straightforward calculations one gets that for any function f


\frac{\partial^2 f}{\partial x^2}\ +\ \frac{\partial^2 f}{\partial y^2}\ +\ \frac{\partial^2 f}{\partial z^2}\ =\ {1 \over r^2}{\partial \over \partial r}\left(r^2 {\partial f \over \partial r}\right) 
  + {1 \over r^2\cos\theta}{\partial \over \partial \theta}\left(\cos\theta {\partial f \over \partial \theta}\right) 
  + {1 \over r^2\cos^2\theta}{\partial^2 f \over \partial \varphi^2}

 

 

 

 

(17)

Introducing the expression (16) in (17) one gets that


\frac{r^2}{\phi}\left(\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2}+\frac{\partial^2 \phi}{\partial z^2}\right)\ =
\frac{1}{R}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right)+\frac{1}{\Theta\cos\theta}\frac{d}{d\theta}\left(\cos\theta \frac{d\Theta}{d\theta}\right)+\frac{1}{\Phi\cos^2\theta}\frac{d^2\Phi}{d\varphi^2}

 

 

 

 

(18)

As the term

\frac{1}{R}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right)

only depends on the variable r and the sum

\frac{1}{\Theta\cos\theta}\frac{d}{d\theta}\left(\cos\theta \frac{d\Theta}{d\theta}\right) + \frac{1}{\Phi\cos^2\theta}\frac{d^2\Phi}{d\varphi^2}

only depends on the variables θ and φ. One gets that φ is harmonic if and only if


\frac{1}{R}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right)\ =\ \lambda

 

 

 

 

(19)

and


\frac{1}{\Theta\cos\theta}\frac{d}{d\theta}\left(\cos\theta \frac{d\Theta}{d\theta}\right) + \frac{1}{\Phi\cos^2\theta}\frac{d^2\Phi}{d\varphi^2}\ =\ -\lambda

 

 

 

 

(20)

for some constant \lambda

From (20) then follows that

\frac{1}{\Theta}\ \cos\theta\ \frac{d}{d\theta}\left(\cos\theta \frac{d\Theta}{d\theta}\right)\ + \lambda\ \cos^2\theta\ +\ \frac{1}{\Phi}\frac{d^2\Phi}{d\varphi^2}\ =\ 0

The first two terms only depend on the variable \theta and the third only on the variable \varphi.

From the definition of φ as a spherical coordinate it is clear that Φ(φ) must be periodic with the period 2π and one must therefore have that


\frac{1}{\Phi}\frac{d^2\Phi}{d\varphi^2}\ =\ -m^2

 

 

 

 

(21)

and


\frac{1}{\Theta}\ \cos\theta\ \frac{d}{d\theta}\left(\cos\theta \frac{d\Theta}{d\theta}\right)\ + \lambda\ \cos^2\theta=\ m^2

 

 

 

 

(22)

for some integer m as the family of solutions to (21) then are


\Phi(\varphi)\ =a\ \cos m\varphi\ +\ b\ \sin m\varphi

 

 

 

 

(23)

With the variable substitution

x=\sin \theta

equation (22) takes the form


\frac{d}{dx}\left((1-x^2) \frac{d\Theta}{dx}\right)+\left(\lambda -\frac{m^2}{1-x^2} \right)\Theta=0

 

 

 

 

(24)

From (19) follows that in order to have a solution \phi with

R(r) = \frac{1}{r^{n+1}}

one must have that

\lambda = n (n+1)

If Pn(x) is a solution to the differential equation


\frac{d}{dx}\left((1-x^2)\ \frac{dP_n}{dx}\right)\ +\ n(n+1)\ P_n\ =\ 0

 

 

 

 

(25)

one therefore has that the potential corresponding to m = 0

\phi = \frac{1}{r^{n+1}}\ P_n(\sin\theta)

which is rotational symmetric around the z-axis is an harmonic function

If P_{n}^{m}(x) is a solution to the differential equation


\frac{d}{dx}\left((1-x^2)\ \frac{dP_{n}^{m}}{dx}\right)\ +\ \left(n(n+1) -\frac{m^2}{1-x^2} \right)\ P_{n}^{m}\ =\ 0

 

 

 

 

(26)

with m ≥ 1 one has the potential


\phi = \frac{1}{r^{n+1}}\ P_{n}^{m}(\sin\theta)\ (a\ \cos m\varphi\ +\ b\ \sin m\varphi)

 

 

 

 

(27)

where a and b are arbitrary constants is a harmonic function that depends on φ and therefore is not rotational symmetric around the z-axis

The differential equation (25) is the Legendre differential equation for which the Legendre polynomials defined


\begin{align}
& P_{0}(x) = 1 \\
& P_{n}(x)=\frac{1}{2^n n!}\ \frac{d^n(x^2-1)^n}{dx^n} \quad n \ge 1 \\
\end{align}

 

 

 

 

(28)

are the solutions.

The arbitrary factor 1/(2nn!) is selected to make Pn(−1)=−1 and Pn(1) = 1 for odd n and Pn(−1) = Pn(1) = 1 for even n.

The first six Legendre polynomials are:


\begin{align}
& P_{0}(x) = 1 \\
& P_{1}(x) = x \\
& P_{2}(x)=\frac{1}{2} \left(3x^2-1\right) \\
& P_{3}(x)=\frac{1}{2} \left(5x^3-3x\right) \\
& P_{4}(x)=\frac{1}{8} \left(35x^4-30x^2+3\right) \\
& P_{5}(x)=\frac{1}{8} \left(63x^5-70x^3+15x\right) \\
\end{align}

 

 

 

 

(29)

The solutions to differential equation (26) are the associated Legendre functions


  P_{n}^{m}(x)\ = (1-x^2)^{\frac{m}{2}}\ \frac{d^m P_n}{dx^m} \quad 1 \le m \le n

 

 

 

 

(30)

One therefore has that


  P_{n}^{m}(\sin\theta)=\cos^m \theta\ \frac{d^n P_n}{dx^n} (\sin\theta)

References

External links

This article is issued from Wikipedia - version of the Wednesday, May 04, 2016. The text is available under the Creative Commons Attribution/Share Alike but additional terms may apply for the media files.