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.

[santisoler]

This looks amazing @MarkWieczorek! I think it's in really good shape to be merged. I just left one minor suggestion. Let me know what do you think.

One other thing I noticed is that the normal_gravity, normal_gravity_potential and normal_gravitational_potential defines the q and q_0 quantities. Does it worth having private methods for those? Just raising the question to see if that would reduce repeated code, but I don't have strong feelings about it.

Anyways, thanks again for the hard work!

[santisoler]

Seeing this reminded me that we can use pytest.warns to check if a code is raising warnings. For example, these lines could be replaced by:

    msg = (
        "Formulas used are valid for points outside the ellipsoid."
        "Height must be greater than or equal to zero."
    )
    with pytest.warns(UserWarning, match=msg):
        ellipsoid.normal_gravity_potential(latitude, height=-10)
    with pytest.warns(UserWarning, match=msg):
        ellipsoid.normal_gravitational_potential(latitude, height=-10)

With pytest.warns we can also check that the text of the raised warning to ensure that the code raised the right warning.

We don't need to change these lines here, I'll open a different PR for this.

[MarkWieczorek]

One other thing I noticed is that the normal_gravity, normal_gravity_potential and normal_gravitational_potential defines the q and q_0 quantities. Does it worth having private methods for those? Just raising the question to see if that would reduce repeated code, but I don't have strong feelings about it.

That's a good point. I think that we could create two auxilliary functions for q and q_0. I don't know if we need these for anything besides these three methods, but it would probably clean up the code a bit. Nevertheless, could we do this after #197 ? Changing this now would complicate that PR.

[MarkWieczorek]

The docs are failing with this error. I don't know how to resolve this...

Run conda install --quiet --file requirements-full.txt python==$PYTHON
  conda install --quiet --file requirements-full.txt python==$PYTHON
  shell: /usr/bin/bash -l {0}
  env:
    REQUIREMENTS: env/requirements-docs.txt env/requirements-build.txt
    PYTHON: 3.12
    INPUT_RUN_POST: true
    CONDA_PKGS_DIR: /home/runner/conda_pkgs_dir
Channels:
 - conda-forge
 - defaults
Platform: linux-64
Collecting package metadata (repodata.json): ...working... done
Solving environment: ...working... failed

LibMambaUnsatisfiableError: Encountered problems while solving:
  - package build-0.10.0-py310h06a4308_0 requires python >=3.10,<3.11.0a0, but none of the providers can be installed

Could not solve for environment specs
The following packages are incompatible
├─[ build ](https://github.com/fatiando/boule/actions/runs/10585105159/job/29330939334?pr=187#logs)is installable with the potential options
│  ├─ build 0.10.0 would require
│  │  └─ python >=3.10,<3.11.0a0 , which can be installed;
│  ├─ build 0.10.0 would require
│  │  └─ python >=3.11,<3.12.0a0 , which can be installed;
│  ├─ build 0.10.0 would require
│  │  └─ python >=3.8,<3.9.0a0 , which can be installed;
│  └─ build 0.10.0 would require
│     └─ python >=3.9,<3.10.0a0 , which can be installed;
└─ python 3.12  is not installable because it conflicts with any installable versions previously reported.

Error: Process completed with exit code 1.

[santisoler]

One other thing I noticed is that the normal_gravity, normal_gravity_potential and normal_gravitational_potential defines the q and q_0 quantities. Does it worth having private methods for those? Just raising the question to see if that would reduce repeated code, but I don't have strong feelings about it.

That's a good point. I think that we could create two auxilliary functions for q and q_0. I don't know if we need these for anything besides these three methods, but it would probably clean up the code a bit. Nevertheless, could we do this after #197 ? Changing this now would complicate that PR.

Yes, totally. I just mentioned after I noticed it, but I don't want that to block this or other PRs. We can totally leave it for later.

The docs are failing with this error. I don't know how to resolve this...

This is unrelated to your PR. The reason behind it is because we were installing build from conda-forge, while the correct package for it should be python-build. I've already fixed that in #207. Updating the branch should solve it.

[leouieda]

Merged! Sorry for the long delay and many many thanks for implementing all of this @MarkWieczorek!