fatiando
Rotate ellipsoids in `create_ellipsoid` based on semiaxes order
Make the create_ellipsoid function to rotate the ellipsoid 90 degrees on the right angles when semiaxes lengths are provided in arbitrary order. An ellipsoid defined by a, b, c with no initial rotation angles passed to the function will be oriented so its semiaxes match with easting, norting, upward directions. Any rotation angles passed to the function will be added to those 90 degrees ones.
Relevant issues/PRs:
Part of #594
Before I continue working on this PR, I'd appreciate some feedback. What would be the intuitive output of calling the create_ellipsoid function as follows?
ellipsoid = hm.create_ellipsoid(1.0, 2.0, 3.0, center=(0, 0, 0))
Option 1: Sort the semiaxes, don't rotate
A triaxial ellipsoid that has a=3.0, b=2.0, c=1.0 with all angles set to zero. So its dimensions are 3.0 in easting, 2.0 in northing, 1.0 in upward.
Option 2: Rotate the ellipsoid so semiaxes order align with enu system
A triaxial ellipsoid that has a=3.0, b=2.0, c=1.0 with yaw=0, pitch=90 and roll=0. This means that its dimensions are 1.0 in easting, 2.0 in northing, 3.0 in upward.
If Option 2 is preferred
If we go with option 2, I'll move forward with this PR. Just one thing to consider, if we also pass some rotation angles to it, like:
ellipsoid = hm.create_ellipsoid(1.0, 2.0, 3.0, pitch=30, center=(0, 0, 0))
Is it acceptable that now ellipsoid.pitch == 120.0? The 30 degrees will get added to the 90 degrees needed to align the passed semiaxes lengths to the coordinate system.
cc @KellySavanna, @leouieda, @nwilliams-kobold
@santisoler Is there an explicit need to do either? Isn't there also:
Option 0: Construct directly without sorting or rotating
where
ellipsoid = hm.create_ellipsoid(1.0, 2.0, 3.0, center=(0, 0, 0))
implies no rotations or transforms and would simply give an ellipsoid centered at (0,0,0) with a=1->east, b=2->north, c=3->up (math convention), regardless of axis lengths. I'd equally be happy with defining the convention as a=1->north, b=2->east, c=3->up (geology convention). But either way if you don't specify a rotation you get exactly what you asked for: those axis length with no transforms applied. No double guessing what the transform will do (like in your Option 2). If the user wants to sort the axes, then they can do so before the call.
When we create an ellipsoid we specify the major axis length (a), and then the other axes are defined as ratios (b/a and c/b) such that we always know they are sorted and decreasing in size.
That's a good point. From the point of view of the definition of the ellipsoid itself, it wouldn't be a problem to define them like that. The tricky part comes with the analytic solutions for the gravity and magnetic fields. The equations are usually written in terms of specific definitions for the ellipsoids:
- Triaxial:
a > b > c - Prolate:
a > b = c - Oblate:
a < b = c
We decided to enforce these conditions in the ellipsoid classes, so the forward modelling functions don't have to sort out ordering of the semiaxes plus apply the (extra) required rotations.
One minor detail I haven't explicitly mentioned. The a, b, c semiaxes are aligned with the x, y, z local coordinate system, not necessarily with easting, northing, upward. In case the ellipsoid doesn't have any rotations, then a -> easting, b -> northing, c -> upward. But, for example for a pitch=90, then a->upward, b->northing, c->easting.
These definitions have a few other benefits. Like reducing the multiple definitions of ellipsoids. For example, any prolate ellipsoid with a->east, b->north, c->upward with constant a, b, and c will be defined as follows:
ProlateEllipsoid(a, b, c, yaw=180*n, pitch=90*m, ...)
with any n and m integers. And these are unique ways of defining that prolate ellipsoid.
If we allow passing a, b, c in any order (like a = b < c), then these two ellipsoids are going to represent the same prolates:
ProlateEllipsoid(c, b, a, yaw=0, pitch=0, ...)
ProlateEllipsoid(a, b, c, yaw=0, pitch=90, ...)
Moreover, by sticking with those definitions we know beforehand which angles will generate invariant rotations for the oblate and prolate (roll in both cases). And that helps to simplify the interface for them: the prolate and oblate have always roll=0.0 and cannot take roll angles.
where
ellipsoid = hm.create_ellipsoid(1.0, 2.0, 3.0, center=(0, 0, 0))implies no rotations or transforms and would simply give an ellipsoid centered at (0,0,0) with a=1->east, b=2->north, c=3->up (math convention)
With option 2, the obtained ellipsoid will have a->easting, b-> northing, c->upward, regardless the ordering of a, b and c. The only difference with option 0 is that it's achieved by adding an extra 90 degrees rotation to the pitch, to satisfy the former conditions. So I guess option 2 is closer to option 0 than option 1.
Option 3: Redesign?
Another alternative is to ditch the multiple ellipsoid classes and keep only an Ellipsoid class that takes a, b, and c in any order, the center, and yaw, pitch, roll angles. And let the forward modelling functions to sort out the order of the semiaxes lengths and the required extra rotations. This would require quite a rewrite, but it might be some path to consider.
With such object, we won't even need the create_ellipsoid function. One can just create the ellipsoid by instantiating the class itself.
I started drafting the Option 3 in #616. I have to say that I like it better than having three ellipsoid classes + a factory function so far. A single class introduces some degeneracy, but I think the interface is much better with it.