Add method for computing prime vertical radius

[santisoler]

Add private method to boule.Ellipsoid for computing the prime vertical radius.
The definition of the prime vertical radius can be seen on Vermeille (2002) and Vajda (2004).
It's usually symbolized with the letter N and depends on the geometry of the ellipsoid and it must be evaluated on a certain latitude angle.
We were already computing it for applying the geodetic to geodetic to spherical coordinate conversion.
It is also needed for computing the Euclidean distance between two points given in geodetic coordinates (see Vajda, 2004), feature that will be added after solving fatiando/harmonica#154 .
In order to prevent rewriting the equation, it's better to have a private function inside boule.Ellipsoid.
The new method takes the precomputed sine of the latitude angle to make it more efficient.

Reminders:

• Run make format and make check to make sure the code follows the style guide.
• Add tests for new features or tests that would have caught the bug that you're fixing.
• Add new public functions/methods/classes to doc/api/index.rst and the base __init__.py file for the package.
• Write detailed docstrings for all functions/classes/methods. It often helps to design better code if you write the docstrings first.
• If adding new functionality, add an example to the docstring, gallery, and/or tutorials.
• Add your full name, affiliation, and ORCID (optional) to the AUTHORS.md file (if you haven't already) in case you'd like to be listed as an author on the Zenodo archive of the next release.

[santisoler]

@leouieda Please take a look at this when you have time.

I decided it's effectively better to have this method. Because the method it's private, I took the liberty of setting the sine of the latitude angle as argument, so it can be efficiently used for coordinates conversions (in Boule) and for calculating geodetic distances (in Harmonica).

[leouieda]

Thanks, @santisoler! Right, this is better than copying the calculation everywhere. I like the thought of just asking the user to do the sine calculation. We might as well make this a public method anyway. Left a few suggestions for the docstring.

Still needs some tests. One way to test this is to calculate on the pole and equator and check against semi-major and semi-minor axis. Another is to try to find pre-computed values at mid latitudes for WGS84 or the like and test against that.

[santisoler]

I like the thought of just asking the user to do the sine calculation.

👍

We might as well make this a public method anyway. Left a few suggestions for the docstring.

Yes, I agree. I even though about making it public after I opened the PR yesterday.

Still needs some tests. One way to test this is to calculate on the pole and equator and check against semi-major and semi-minor axis. Another is to try to find pre-computed values at mid latitudes for WGS84 or the like and test against that.

Yes, you're right about that. Even though the coordinate conversions tests are actually using it. we might want to add specific tests for this new method.
One little thing, the prime vertical radius at the poles is not equal to the semi-minor axis.
On Wikipedia you can find an alternative way of expressing the prime vertical radius:

If latitude is 90 or -90, then the cosine is 0, therefore the prime vertical radius will be equal to:

That's why, when converting the coordinates to Cartesian and specially computing the z coordinate, the prime vertical radius is also multiplied by (1 - e^2):

I'll add a test function to compare the prime vertical radius on the equator and the poles.
I'll see if I can find precomputed values for other latitudes.

[leouieda]

@santisoler looks good to go! I fixed the issues with the Mac builds. Merging and releasing 0.2.0