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

[VascoSch92]

This is a bit strange because in the docstring they have other names :)

[MarkWieczorek]

I could change this to just u and reduced_latitude. The only reason I did it this way was that these are the names of the parameters in the paper we cite, and they are also the names of the variables used in the normal gravity routine. When writing this, I originally used u and then changed it at the end, so I don't have any preference either way...

[VascoSch92]

it makes sense. Perhaps change just the docstring so you don't have surprise ;-) but it is a nit

[MarkWieczorek]

This PR was incorporated into #187