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

[MarkWieczorek]

I added a centrifugal_potential method for the triaxial ellipsoid now that PR with the semimajor_axis_longitude attribute has been accepted.

[MarkWieczorek]

Note: I tested three different analytic solutions for going between geodetic and ellipsoidal harmonic coordinates, and they all gave the same precision. The alternative forms can be found in the manuscript "Closed-form transformation between geodetic and ellipsoidal coordinates" by Featherstone and Claessens (2008, http://link.springer.com/10.1007/s11200-008-0002-6). I think that the precision is good enough, and I suspect that the problem comes down to the fact that you need to solve a quartic equation.

[MarkWieczorek]

The last comit precomputes some variables that are used repeatedly in the normal gravity computations (such as self.linear_eccentricy which is not a constant, but requires a computation). It also removes the redundant code for computing the ellipsoidal harmonic coordinates, as this is done in a separate dedicated routine. Doing so makes the code for normal_gravity, normal_gravitational_potential and normal_gravity_potential more similar. If we need to swap out geodetic_to_ellipsoidal_harmonic for a more accurate routine later, this will make our lives easier.

I don't plan on making any further changes to this PR, unless requested.

[leouieda]

@MarkWieczorek this is truly fantastic! Thank you for all your work on this and for the pointer to the paper. I'm reviewing this by pieces whenever I can find some spare time. I'm eager to add this to Boule since it would allow us to work with the disturbing potential instead of geoid heights, which I find more intuitive.