Add normal_gravitational_potential, normal_gravity_potential, and centrifugal_potential

[MarkWieczorek]

This PR makes several additions in order to compute the normal gravity/gravitational potential on or above the ellipsoid.

• Adds normal_gravity_potential(), normal_gravitational_potential() and centrifugal_potential() methods to both the Ellipsoid and Sphere class.
• Adds two conversion routines, to and from geodetic and ellipsoidal-harmonic coordinates: Ellipsoid.geodetic_to_ellipsoidal_harmonic() and Ellipsoid.ellipsoidal_harmonic_to_geodetic().
• Adds tests to ensure that normal_gravity_potential() = normal_gravitational_potential() + centrifugal_potential().
• Adds test for the coordinate transforms. See below: these are surprisingly not as accurate as I thought they would be.

Notes

1. The centrifugal potential calculation needs the perpendicular distance to the rotation axis. For this, I used the definition of the prime radius of curvature and geodetic latitude which gives: $(N(\phi) + h) \cos\phi$, where $\phi$ is geodetic latitude and $h$ is ellipsoidal height. There are probably other ways that this could be calculated. A separate PR will implement this for triaxial ellipsoids, as we will need to change how we handle longitude_semimajor_axis.

2. The geodetic to ellipsoidal harmonic coordinate transforms work, but I am worried about the precision. The conversion from geodetic to ellipsoidal uses the same code as in the normal_gravity method, which is based on the equations in Lakshmanan (1991). The inverse is just a a couple atans for the latitude, and the ellipsoidal height is computed as the difference of the prime_vertical_radius of the two ellipsoids. Even if I can understand why the height might lose precision this way, the latitude shouldn't.

Here are a couple of examples:

In [2]: longitude, reduced_latitude, u = boule.WGS84.geodetic_to_ellipsoidal_harmonic(0, 45., height=0)

In [3]: boule.WGS84.ellipsoidal_harmonic_to_geodetic(longitude, reduced_latitude, u)
Out[3]: (0, 45.0, 0.0)

In [4]: longitude, reduced_latitude, u = boule.WGS84.geodetic_to_ellipsoidal_harmonic(0, 45., height=1.e-8)

In [5]: boule.WGS84.ellipsoidal_harmonic_to_geodetic(longitude, reduced_latitude, u)
Out[5]: (0, 44.99999999999992, 9.930249518419041e-09)

In [6]: longitude, reduced_latitude, u = boule.WGS84.geodetic_to_ellipsoidal_harmonic(0, 45., height=1)

In [7]: boule.WGS84.ellipsoidal_harmonic_to_geodetic(longitude, reduced_latitude, u)
Out[7]: (0, 44.999999969779665, 0.9999999988228683)

In [8]: longitude, reduced_latitude, u = boule.WGS84.geodetic_to_ellipsoidal_harmonic(0, 45., height=1000)

In [9]: boule.WGS84.ellipsoidal_harmonic_to_geodetic(longitude, reduced_latitude, u)
Out[9]: (0, 44.99996978922941, 1000.0002619091734)

In [10]: longitude, reduced_latitude, u = boule.WGS84.geodetic_to_ellipsoidal_harmonic(0, 45., height=10000)

In [11]: boule.WGS84.ellipsoidal_harmonic_to_geodetic(longitude, reduced_latitude, u)
Out[11]: (0, 44.999698744381085, 10000.026154076979)

Maybe this is good enough. I suspect that the loss of precision comes from the Lakshmanan equations (which have some subtractions of similarly sized numbers). Perhaps there is nothing we can do about it. Nevertheless, I note that if this is a problem, then it will also affect the normal gravity routine that uses the same method.

Relevant issues/PRs: Closes #151