Small flattening ellipsoids

[MarkWieczorek]

This PR extends the Ellipsoid class to allow the use of flattenings that are close to, and including, zero. The textbook equations used in the previous version ran into numerical instabilities when the flattening was lower than about $10^{-7}$. This was a result of either divide by 0 errors when the linear eccentricity approached zero, or a result of square roots of negative numbers (like -1.e-30) when computing the angle $\beta^{\prime}$.

This PR uses Taylor series expansions about $f=0$ that do not suffer from these numerical problems. The equations and their applicability are presented in #194. In general, for each routine that is concerned, the flattening is checked, and if it is smaller than a given value, the small flattening equations are used instead of the textbook equations.

The following routines have been modified:

• Ellipsoid.reference_normal_gravity_potential
• Ellipsoid.normal_gravitational_potential
• Ellipsoid.normal_gravity_potential
• Ellipsoid.normal_gravity
• Ellipsoid.gravity_equator
• Ellipsoid.gravity_pole
• Ellipsoid.geodetic_to_ellipsoidal_harmonic

The web documentation for the normal gravity routines has been cleaned up a little to better explain the difference between an ellipsoid with zero flattening and a sphere.

Note that the code in Ellipsoid.geodetic_to_ellipsoidal_harmonic is largely repeated in Ellipsoid.normal_gravity. It might be preferable to simply call geodetic_to_ellipsoidal_harmonic within normal_gravity for clarity and brevity.

Note also that this PR is based on this other (not yet merged) PR: #187. It it probably better to address #187 first, given the number of changes made in both of these.

Relevant issues/PRs:
Relevant to #194

[MarkWieczorek]

After fixing a bunch of merge conflicts, this PR is ready to be reviewed. It should be trivial to verify that the case when the flatting is large corresponds to the old code. The only thing to do is to very the small flattening limits that are presented in #194.

[santisoler]

Thanks @MarkWieczorek for all the hard work behind this. Solving numerical instabilities with small numbers and juggling with machine precision requires a lot of effort.

I just pushed one small change to this branch that moves the caution in the Ellipsoid docstring away from the Parameters section and puts it above it. I thought it would be easier to push the change than to suggest it in the review.

I also left two minor suggestions down below.

I think this PR is almost ready. The one thing I would add is more tests for the new cases that were introduced. Right now the 100% coverage is achieved because of the test function that ensures that no warning is being raised, but that tests is not actually ensuring that the result we get out of the new lines is the one we expect.

One way to solve this would be to run all the tests we have for the regular ellipsoids (large flattening) with another ellipsoid with zero flattening. We can quickly do this by adding a new ellipsoid with small flattening to the ELLIPSOIDS list in test_ellipsoid.py.

Moreover, I would also like to add some sort of test that ensures that the physical properties are continuous for an ellipsoid with a flattening slightly below the "small flattening" threshold and for an ellipsoid with a flattening slightly above the threshold. For example, I would expect the gravity_equator for an ellipsoid with flattening=1.24e-5 to be quite close to the one with a flattening=1.26e-5. I can totally help with this if I'm asking too much!

Let me know what do you think!

[MarkWieczorek]

The one thing that is annoying about this PR is that the flattening cutoff is somewhat arbitrary. I don't think that we can do much about that, but what we could do is check that the results for using a flattening of flattening_cutoff +/- small_value give the same result within some predefined relative precision (like 1.e-8).

If we really wanted to do even better, we could add some higher order terms to some of the small flattening limits, which should improve the accuracy if we choose the cutoff value incorrectly. I really don't want to go through all the derivations from scratch, but the Physical geodesy book does give higher order terms for some of the equations (but not all) in the associated issue : $U_0$, $V(u, \beta)$, $\gamma_a$, $\gamma_b$. I don't think that they give an expression for $U(u, \beta)$ and the two terms for the normal gravity above the ellipsoid and $\beta^{\prime}$, but I'll look again.

In any case, my main motivation for the PR was simply to be able to use a flattening of 0.0.

[santisoler]

I understand the motivation.

I wouldn't stress too much on deriving the higher order terms until we are sure we need those. That's why I suggested to focus on adding more tests that ensure that the approximations we have already are accurate enough.

For example, if we run this little script:

import numpy as np
import boule as bl
import matplotlib.pyplot as plt

kwargs = dict(
    name="foo",
    semimajor_axis=bl.WGS84.semimajor_axis,
    geocentric_grav_const=bl.WGS84.geocentric_grav_const,
    angular_velocity=bl.WGS84.angular_velocity,
)

threshold = 1.25e-5
half_range = 0.25e-5
flattenings = np.linspace(threshold - half_range, threshold + half_range, 501)
ellipsoids = (bl.Ellipsoid(flattening=f, **kwargs) for f in flattenings)
gravity_equator = np.array([e.gravity_equator for e in ellipsoids])

plt.plot(flattenings, gravity_equator)
plt.axvline(x=threshold, linestyle="--", color="grey")
plt.xlabel("Flattening")
plt.ylabel("Gravity on the equator [m/s2]")
plt.show()

We get this:

Where the grey dotted line represents the threshold. I can still see some instabilities on flattenings larger than the threshold, so maybe we would need to increase it and see if the approximation is still valid. If not, I agree that maybe we should explore higher order terms.