MarkWieczorek
commented
7 months ago And here are some images showing how the small angle-limit formulas compare with the standard equations for the parameters
- $\frac{GM}{E} \arctan(E/b)$
- $q_0$
- $q^{\prime}_0$
- $E q^{\prime}_0 / (q_0 b)$
For the case of Earth, you can see the approximate formulas are good for flattenings less than about $10^{-3}$, and that the standard formulas run into numerical issues below about $10^{-7}$. The results are nearly the same Ceres, which is considerably smaller.
Inspecting how the standard equations degrade in precision with decreasing flattening (for bodies the size of Earth and Ceres), the small-angle equations should be used in the following situations:
- $\frac{1}{E} \arctan(E/b)$ : $f< 4 \times 10^{-16}$
- $q_0$ : $f< 1.5 \times 10^{-5}$
- $q^{\prime}_0$ : $f< 1 \times 10^{-5}$
- $E q^{\prime}_0 / (q_0 b)$ : $f< 1.25 \times 10^{-5}$
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:
However, for a sphere, the assumptions are different:
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=0in 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!
In Boule.Ellipsoid, numerical issues are encountered when the flattening goes to zero. Here, I provide equations that are more precise in the limit where the flattening goes to zero. In fact, using these new equations, one can even use a flattening of zero in the ellipsoid class. I note that a flattening of zero in the ellipsoid class is not equivalent to a Boule.Sphere class, and I have updated the documentation to make this a little more clear.
Though I have derived these equations myself, I provide equation numbers in the Physical Geodesy book where the equations can be verified.
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 $$
This compares with eq. 2-185 (but which gives an extra fourth order term not included here).
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 gravitational 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] $$
This compares with eq. 2-158 (though they give higher order terms as well).
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 gravitational 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 (eq. 2-126)
$$ 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 on ellipsoid
The normal gravity on the ellipsoid is given by 2-145 and 2-146, which is Somilgiana's formula:
$$ \tan \beta = \frac{b}{a} \tan \phi $$
$$ \gamma = \frac{a \gamma_a \cos^2 \phi + b \gamma_b \sin^2 \phi}{\sqrt{a^2 \cos^2 \phi + b^2 \sin^2 \phi}} $$
For this we need the normal gravity at the pole and equator (2-141, 2-142):
$$ \gamma_a = \frac{GM}{ab} \left(1-m- \frac{m}{6}\frac{e^{\prime} q^{\prime}_0}{q_0} \right) $$
$$ \gamma_b = \frac{GM}{a^2} \left(1 + \frac{m}{3}\frac{e^{\prime} q^{\prime}_0}{q_0} \right) $$
where
$$ e^{\prime} = \frac{E}{b} $$
and
$$ m= \frac{\omega^2 a^2 b}{GM} $$
Each of these terms has the offending ratio
$$ \frac{E q^{\prime}_0}{b q_0} $$
where (eq. 2-133)
$$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) $$
The above two equations compare with 2-186 (which includes higher order terms).
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$$
normal gravity above the ellipsoid
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} \simeq 3 \frac{b^3}{\left(b^{\prime}\right)^2} $$
The second involves the term concerning $\cos \beta'$ (eq. A2 of Li and Götze), which can be rewritten as
$$ \cos \beta' = \left[\frac{1}{2} + \frac{R}{2} - \frac{R}{2}\sqrt{1 + \frac{1}{R^2} - 2 \frac{D}{R^2}} \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} \simeq 1 + \frac{x}{2} - \frac{x^2}{8} $$
for the inner square root in the above expression, in the limit where $E$ approaches zero (ignoring terms in $E^4$ and $E^6$ within the square root), 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} $$