Add TriaxialEllipsoid class with geometric parameters

[dabiged]

Skeleton class for a Triaxial Ellipsoid. I have implemented the basic framework and a number of properties from Issue #47 .
More advanced methods raise a NotImplementedError.
I have taken formula from wikipedia and this paper: https://www.sciencedirect.com/science/article/abs/pii/S0019103584711183?via%3Dihub

ToDo:

• Internal Doctests using a generic watermelon triaxial ellipsoid.
• Test within the class that semimajor > semimedium > semiminor.
• Tests to the test directory.
• realization of a triaxial ellipsoid (Vesta / Eros?)
• Examples to the gallery

Submitting this Pull request early to get feedback.

• Should the class be called TriAxialEllipsoid or TriaxialEllipsoid
• Is the current watermelon example suitable?

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.

[dabiged]

I would love some feedback on this before I go much further.
Implementing things like a radius method, or anything more advanced seems to require
elliptic integrals which I do not feel comfortable doing.

[santisoler]

Hi @dabiged! I'm really happy to see an open PR for Triaxial Ellipsoids.

I don't have too much knowledge on this type of ellipsoids, maybe @MarkWieczorek can give a higher quality review. Nevertheless I'll leave some inline comments (probably more related to code than to the underlying math).

[santisoler]

I'm wondering if it's better to write a special method for testing the relations between the semimajor_axis, semimedium_axis and the semiminor_axis. If attrs make it too complex we can leave these tests separeted as they are now.

[dabiged]

This is my first time using the attr toolset. I have no idea how to write this and the attr documentation seems quite sparse. Can you give me any idea where to start?

[santisoler]

We could add a new semimedium_axis validator that checks the semimedium_axis against the other two. Something like:

    @semimedium_axis.validator
    def _check_valid_axes(self, semimedium_axis, value):
        if self.semiminor_axis > value or value > self.semimajor_axis:
            raise ValueError("Invalid set of axes: ....")

This way we can have a single ValueError for any possible combination of invalid set of axes.
Although it's weird to have a semimedium_axis.validator decorator for a method that checks the entire set of axes, I couldn't came out with a better solution.

cc @leouieda

[santisoler]

We could also add another validator that raises a warning if semimedium_axis is equal to semiminor_axis or equal to semimajor_axis and instruct people to use the Ellipsoid class instead.

[santisoler]

Nice! If we add this method on TriaxialEllipsoid we should do the same on Ellipsoid and Sphere, but that's for a separate PR.

[dabiged]

I have added this as issue #75 so I don't forget!

[santisoler]

We don't need to define a radius for the TriaxialEllipsoid. The radius is only suited for the Sphere.

Suggested change

	def radius(self, latitude, longitude):
	r"""
	Return the radius of the Triaxial Ellipsoid at the given latitude and longitude [metres]
	"""
	raise NotImplementedError

[dabiged]

This method was mentioned in the feature request.

I assumed this was a vector running from the centre of the ellipsoid (aka the origin [is it called the geocentric location??]) to a point on the surface at a given latitude and longitude.
For a sphere this is a constant.
For a biaxial ellipsoid this maps to an ellipse along lines of longitude and a constant on lines of latitude.
For a triaxial ellipsoid I think it is this value (source: wikipedia page on ellipsoids).

Should this method be called something else?

[leouieda]

This is akin to the geocentric_radius’s of Ellipsoid` https://www.fatiando.org/boule/dev/api/generated/boule.Ellipsoid.html#boule.Ellipsoid.geocentric_radius

[santisoler]

Sorry for the misunderstanding: I though you were trying to implement a not implemented method in replacement of the radius attribute of the Sphere.
Now that I know what you meant by radius, I completely agree that the class should have this method.
Nevertheless, radius is a very misleading name. The Ellipsoid class has a geocentric_radius method which returns the quantity you describe. That name explicitly says that we are about to return the norm of the radius vector to the surface of the ellipsoid measured from a geocentric coordinate system, so the origin is at the center of the ellipsoid.

In the TriaxialEllipsoid, the geocentric_radius should get the longitude and the latitude, as you already did. But, remember that on every Fatiando library we keep the following order for coordinates: longitude, latitude, height, ... or easting, northing, height, ..., except of very special cases where the reason is well documented.

So feel free to ignore this suggestion and please implement the geocentric_radius as you described.

[dabiged]

Thanks for clearing up the misunderstanding there.
I was completely wrong and didn't see geocentric_radius method.

I have removed the radius method from the latest commit.

[santisoler]

I'm wondering if we can benefit from creating some inheritance, I mean make Ellipsoid a subclass of TriaxialEllipsoid. Even though from a mathematic perspective is true (any biaxial ellipsoid is a triaxial ellipsoid with one pair of axis of equal lenght), maybe it could not bring any advantage from a code design perspective.

cc @leouieda

[leouieda]

@santisoler @dabiged this is bringing up some interesting flaws in our current design. I’m usually hesitant to add a lot of inheritance. It’s rarely a good solution in my experience. Instead, see #76 for a proposal to have a single class. I think this would fix a lot of issues with methods taking unused arguments.

[MarkWieczorek]

I am ok with geocentric_radius, but at the same time, in spherical coordinates (r, theta, phi), r is "radius". Maybe geodesists aren't used to working with spherical coordinates, so having this be explicit is ok.

Code discussed is no longer part of this PR

[leouieda]

Alright, I finished the few things left over from this PR. Removed all of the not implemented methods instead to simplify things. We can add them later when we have the need for such calculations.