We previously had a discussion about how to define normal gravity on a sphere. At the time, I think that I was probably a little confused, but now that I see how Boule works and defines things, and after have looked through Chapter 2 of Physical Geodesy again, I think that the way we treat the Sphere class is not entirely consistent with the Ellipsoid class.
Start with the properties of a reference ellipsoid, as found in the boule documentation:
- The gravity potential (gravitation + centrifugal) is constant on the surface.
- The internal density structure is unspecified but must lead to a constant potential at the surface.
However, for a sphere, the assumptions are different:
- The internal density structure is unspecified but must be either homogeneous or vary radially (e.g., in concentric layers).
- A constant gravity potential on the surface of a rotating sphere is not possible (when the internal density structure varies only as a function of radius).
Because of this, if you take an Ellipsoid, and gradually decrease the flattening to zero, you will not asymptotically approach the results of a sphere. This is because the gravity potential is not constant on the surface of the sphere, but it is on an ellipsoid.
Now, I agree that the assumptions for the Sphere class are correct if we assume that the body is a fluid in hydrostatic equilibrium. However, we can add arbitrary density anomalies in the mantle in order to generate a gravity potential (gravitation+centrifugal) at the surface that is constant, just like in the case for Ellipsoid. The mathematical problem is that if you put flattening=0 in the Ellipsoid equations, you will get divide by zero errors, but as I will show below, these can be avoided by taking the limit where the semiminor axis approaches the semimajor axis.
What I propose is the following: use the same assumptions for Sphere as for Ellipsoid, and then with the equations in Ellipsoid, find the asymptotic limits as the flattening goes to zero. I note that the following equations could probably have been more easily derived starting with spherical harmonics instead of ellipsoidal harmonics, but I will follow the geodesy tradition of doing things the hard way!
reference normal gravity potential
Start with eq. 2-123,
$$
U_0 = \frac{GM}{E} \arctan{\frac{E}{b}}+ \frac{1}{3} \omega^2 a^2
$$
and use the small-angle approximation (only the first two terms are necessary here)
$$
\arctan{x}\simeq x -\frac{x^3}{3} + \frac{x^5}{5} + O(7)
$$
When $b$ approaches $a$ and $E$ goes towards zero, we have
$$
U_0 \simeq GM \left(\frac{1}{b} - \frac{E^2}{3 b^3}\right) + \frac{1}{3} \omega^2 a^2
$$
For the case of a sphere with $b=a$, the reference potential is
$$
U_0 = \frac{GM}{a} + \frac{1}{3} \omega^2 a^2
$$
normal gravitation potential
Eq. 2-124 for the normal gravitation potential (where $\beta$ is the reduced latitude)
$$
V(u, \beta) = \frac{GM}{E} \arctan{\frac{E}{u}} + \frac{1}{3} \omega^2 a^2 \frac{q}{q_0} \left[\frac{3}{2} \sin^2 \beta - \frac{1}{2} \right]
$$
has the offending term
$$
\frac{q}{q_0} = \frac{ \frac{1}{2}\left[\left(1+3\frac{u^2}{E^2} \right) \arctan{\frac{E}{u}} - 3 \frac{u}{E}\right]}{\frac{1}{2}\left[\left(1+3\frac{b^2}{E^2} \right) \arctan{\frac{E}{b}} - 3 \frac{b}{E}\right]}
$$
Using the higher-order small angle approximation of the arctan function above (all three terms are necessary), we find that
$$
q \simeq \frac{2}{15} \frac{E^3}{u^3}
$$
and
$$
\frac{q}{q_0} \simeq \frac{b^3}{u^3}
$$
As $b$ approaches $a$ the normal gravity potential approaches
$$
V(u, \beta) \simeq GM \left(\frac{1}{u} - \frac{E^2}{3 u^3}\right)+ \frac{1}{3} \omega^2 a^2 \frac{b^3}{u^3} \left[\frac{3}{2} \sin^2 \beta - \frac{1}{2} \right]
$$
For the case of a sphere, with $a=b$, $u=r$ and $\beta =\phi$ we have
$$
\frac{q}{q_0} \simeq \frac{a^3}{r^3}
$$
the normal gravity potential 2-124 is
$$
V(r, \phi) = \frac{GM}{r} + \frac{1}{3} \omega^2 a^2 \frac{a^3}{r^3} \left[\frac{3}{2} \sin^2 \phi - \frac{1}{2} \right]
$$
where $\phi$ is the spherical geocentric latitude and $r$ is geocentric spherical radius.
normal gravity potential
The normal gravity potential
$$
U(u, \beta) = \frac{GM}{E} \arctan{\frac{E}{u}} + \frac{1}{2} \omega^2 a^2 \frac{q}{q_0} \left[\sin^2 \beta - \frac{1}{3} \right] + \frac{1}{2} \omega^2 \left(u^2 + E^2\right) \cos^2 \beta
$$
also has the offending $q/q_0$ term. Substituting the above approximation into 2-126, as $b$ approaches $a$, the normal gravity potential goes to
$$
U(u, \beta) \simeq GM \left(\frac{1}{u} - \frac{E^2}{3 u^3}\right) + \frac{1}{2} \omega^2 a^2 \frac{b^3}{u^3} \left[\sin^2 \beta - \frac{1}{3} \right] + \frac{1}{2} \omega^2 \left(u^2 + E^2\right) \cos^2 \beta
$$
For the case of a sphere with $a=b$, $u=r$, and $\beta=\phi$, the normal gravity potential is
$$
U(r, \phi) = \frac{GM}{r} + \frac{1}{2} \omega^2 a^2 \frac{a^3}{r^3} \left[\sin^2 \phi - \frac{1}{3} \right] + \frac{1}{2} \omega^2 r^2 \cos^2 \phi
$$
normal gravity
The normal gravity on the ellipsoid is given by 2-146, which is Somilgiana's formula. For this we need the normal gravity at the pole and equator. Each of these terms has the offending ratio
$$
\frac{E q^{\prime}_0}{b q_0}
$$
where
$$q^{\prime}_0 = 3 \left(1 + \frac{b^2}{E^2}\right) \left(1 - \frac{b}{E} \arctan \frac{E}{b} \right) -1 $$
Using the small-angle approximation (all three terms are necessary), we find that
$$
q^{\prime}_0 \simeq \frac{2}{5} \frac{E^2}{b^2}
$$
and that
$$
\frac{E q^{\prime}_0}{b q_0} \simeq 3
$$
The gravity at the equator, eq 2-141, is
$$
\gamma_a \simeq \frac{GM}{a b} \left( 1 - \frac{3}{2} m \right)
$$
and the gravity at the pole, eq 2-142, is
$$
\gamma_b \simeq \frac{GM}{a^2} \left( 1 + m \right)
$$
For the case of a sphere with $a=b$, the above two equations reduce to
$$
\gamma_a = \frac{GM}{a^2} \left( 1 - \frac{3}{2} m \right)
$$
and
$$
\gamma_b = \frac{GM}{a^2} \left( 1 + m \right)
$$
which yields for the normal gravity
$$\gamma = \gamma_a \cos^2 \phi + \gamma_b \sin^2 \phi$$
When computing the normal gravity above the ellipsoid, we use the equations in the appendix of Li and Götze (2001). These equations have two offending terms when the flattening goes to zero. One is nearly the same as computed above:
$$ \frac{E q^{\prime}}{q_0} = 3 \frac{b^3}{\left(b^{\prime}\right)^2} $$
The second involves the term concerning $\cos \beta'$, which can be rewritten as
$$
\cos \beta' = \left[\frac{1}{2} + \frac{R}{2} - \frac{R}{2}\left(1 + \frac{1}{R^2} - 2 \frac{D}{R^2} \right) \right]^{1/2}
$$
where $R=r^{\prime \prime 2} / E^2$ and $D=d^{\prime \prime 2} / E^2$. Using the approximation
$$(1+x)^{1/2} = 1 + \frac{x}{2} - \frac{x^2}{8} $$
for the inner square root in the above expression, in the limit where $E$ approaches zero, we find
$$
\cos \beta' = \left[ \frac{1}{2} + \frac{1}{2} \frac{d^{\prime\prime 2}}{ r^{\prime\prime 2}} + \frac{E^2}{4} \left( \frac{d^{\prime\prime 4}}{r^{\prime\prime 6}} - \frac{1}{r^{\prime\prime 2}} \right) \right]^{1/2}
$$
![dependabot[bot] avatar](https://avatars.githubusercontent.com/in/29110?v=4 "dependabot[bot]")
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent