Refactor code to have a single `Ellipsoid` class #616

santisoler · 2025-12-01T19:55:37Z

Refactor the ellipsoid code to have a single Ellipsoid class that can represent any type of ellipsoid. Ditch the specific classes for each ellipsoid type and the create_ellipsoid function. Update the forward modelling code for gravity and magnetics. Make use of an extra rotation matrix that aligns the x, y, z local coordinate system to the ellipsoid's semiaxes in decreasing order, so a >= b >= c. Modify the analytic solutions for oblate ellipsoids in such way that they are defined by a == b > c (instead of a < b == c as Takahashi et al. and Clark et al. do). Make rotation angles and center as optional arguments for the Ellipsoid class. Fix numerical instabilities of triaxial ellipsoids when they approximate prolate or oblate ellipsoids (when two semiaxes lengths are close enough to each other). Update the tests and add a few more for the new bits of code.

Relevant issues/PRs:

Part of #594

Start implementing gravity field with a permutation matrix to sort out the semiaxes.

Update the elliptical integrals for oblates so they assume that `a == b > c`. Simplify the permutation matrix since now we don't need to take care of oblates as a special case.

Remove quite a lot of tests since we now have a single ellipsoid class.

It makes it more explicit about what's happening with the semiaxes. The permutation matrix should match that sorting, but we can cover that in tests.

Update the magnetic code to use the single `Ellipsoid` class, make use of the permutation matrix, and update the internal demagnetization tensor equations for the oblate ellipsoid defined as `a == b > c`.

Do this to differentiate it from just the rotation matrix, since it also includes the permutation one.

santisoler · 2025-12-01T20:02:48Z

TODO:

add tests for the ~~get_permutation_matrix~~ get_semiaxes_rotation_matrix function (fixed a bug on it)
test gravity and magnetic fields of triaxial ellipsoids vs oblates and prolates (when two of the semiaxes are similar to each other).

The permutation matrix used wasn't right: we need to apply a proper rotation matrix (with 90 degrees) to ensure that the local coordinate system is a right-handed one.

Test ellipsoids defined with semiaxes passed in arbitrary orders against equivalent ellipsoids defined with ``a >= b >= c``, but with necessary rotations.

After changing the definition of the oblate to `a == b > c`, we need to change the condition for oblates in such function.

Fix the numerical instabilities of gravity fields produced by triaxial ellipsoids that approximate a prolate and an oblate ellipsoid. Fallback to the prolate and oblate solutions when two of the three semiaxes are close enough to each other.

Fallback to the prolate and oblate solutions for the internal demagnetization tensor when triaxial has two semiaxes close enough to each other.

harmonica/_forward/ellipsoids/utils.py

harmonica/tests/ellipsoids/test_utils.py

harmonica/_forward/ellipsoids/gravity.py

Since we merged all classes to a single one, we don't actually need an `Ellipsoid` type class. It's rare that users will pass an object that is not an instance of `hm.Ellipsoid` or a child of it.

Add sanity checks to some of the private functions in which the semiaxes are assumed to be sorted such as `a >= b >= c`. Most of these errors are never going to be accessed by code ran by users, but they might catch bugs when refactoring. Adding those checks there since it's hard to test the behaviour from outside.

Split the computation of lambda for oblates and prolates. Since we changed the definition of oblate ellipsoids to `a = b > c`, the solution of lambda through the second degree equation is different than the one for prolates.

santisoler added 13 commits November 27, 2025 13:56

Start drafting a single Ellipsoid class

da1c361

Start implementing gravity field with a permutation matrix to sort out the semiaxes.

Rename variable holding rotation and permutation matrix

963292a

Fix get_permutation_matrix for spheres

b1d369e

Set angles and center as optional arguments

22c8547

Update elliptical integrals to oblates as a == b > c

abb9379

Update the elliptical integrals for oblates so they assume that `a == b > c`. Simplify the permutation matrix since now we don't need to take care of oblates as a special case.

Update gravity tests

322edb4

Update test_ellipsoids.py

8ffc88c

Remove quite a lot of tests since we now have a single ellipsoid class.

Update docstrings of get_permutation_matrix

2be3c55

Use sorted to sort semiaxes in gravity function

d430c7f

It makes it more explicit about what's happening with the semiaxes. The permutation matrix should match that sorting, but we can cover that in tests.

Update magnetic code to Ellipsoid class and oblates as a == b > c

a1cdfc8

Update the magnetic code to use the single `Ellipsoid` class, make use of the permutation matrix, and update the internal demagnetization tensor equations for the oblate ellipsoid defined as `a == b > c`.

Remove unused import

5a8211d

Get rid of mentions of the old ellipsoid classes

9ceb552

Rename matrix to transformation

6084607

Do this to differentiate it from just the rotation matrix, since it also includes the permutation one.

santisoler mentioned this pull request Dec 1, 2025

Rotate ellipsoids in create_ellipsoid based on semiaxes order #613

Closed

santisoler added 15 commits December 1, 2025 16:55

Add test comparing gravity of triaxial vs oblate and prolate

acca7b7

Add test comparing mag field of triaxials vs oblate and prolate

5383d85

Fix comment in test

ff85ce9

Fix comments in tests

83998cc

Merge branch 'main' into single-ellipsoid-class

1848ebf

Fix bug in permutation matrix

bd92efd

The permutation matrix used wasn't right: we need to apply a proper rotation matrix (with 90 degrees) to ensure that the local coordinate system is a right-handed one.

Test semiaxes in any order with equivalent ellipsoids

3efba9f

Test ellipsoids defined with semiaxes passed in arbitrary orders against equivalent ellipsoids defined with ``a >= b >= c``, but with necessary rotations.

Update tests for __str__ method

2992345

Add test to get_semiaxes_rotation_matrix function

1a25659

Fix is_almost_a_sphere function

65d907d

After changing the definition of the oblate to `a == b > c`, we need to change the condition for oblates in such function.

Add test for is_almost_a_sphere

c723ca5

Add test for numerical instabilities on gravity fields

6549a17

Fix numerical instabilities in internal demagnetization tensor

528ffeb

Fallback to the prolate and oblate solutions for the internal demagnetization tensor when triaxial has two semiaxes close enough to each other.

Add tests for the magnetic numerical instabilities

384df38