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: https://github.com/fatiando/boule/pull/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