Add conversion routines between ellipsoidal-harmonic and geodetic coordinates

[MarkWieczorek]

This PR adds two conversion routines, to and from geodetic and ellipsoidal-harmonic coordinates:

• Ellipsoid.geodetic_to_ellipsoidal_harmonic()
• Ellipsoid.ellipsoidal_harmonic_to_geodetic()

These will be required when computing the normal gravitational potential above the reference ellipsoid in issue #151.

The routines 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, but if this is the case, then it affects the Normal gravity routine.

Relevant issues/PRs: Related to #151

[MarkWieczorek]

This PR was incorporated into #187