Improvements to ellipsoids

[santisoler]

Description of the desired feature:

Listing here some ideas to improve ellipsoids code:

• [x] Move the physical properties of the ellipsoids to the objects. Remove them from the functions calls (#597).
• [x] Expose the rotation matrix through a property in ellipsoid classes (#597).
• [x] Optimize the magnetic forward by replacing for loop over observation points: build all demagnetization tensors in vectorized fashion (through a 3D Numpy array) and use them to compute the dot products (#603).
• [x] Optimize computation of lambda for prolate and oblate: we can solve the second-order equation rather than the third-order one (#604).
• [x] Add a Sphere class that represents an homogeneous sphere as a source (#602).
  • Use that class in tests, when comparing ellipsoids to spheres.
• [x] Add a constructor function (e.g. ~~get_ellipsoid~~ create_ellipsoid) that takes semiaxes lengths and returns and object of the right class (with semiaxes in the right order). (#606)
• [x] Fallback to sphere forward modelling if semiaxes lenghts are almost equal (e.g. a/b == 1.0000000001), to avoid numerical instabilities. Optionally: raise a warning to let users know about it (#611).
• [x] Add method to export ellipsoids to PyVista (#614).
• [ ] Add a __str__ method to ellipsoids, so we can print them and their properties in a nice looking way (#615).
• [ ] Solve numerical instabilities when triaxial ellipsoids approximate an oblate or a prolate. Basically when two of the three semiaxes are close to each other.
  • From my experiments, the triaxial solutions are quite robust. The difference between two semiaxes must be really small (~1e-14) for the instabilities to appear.
  • Falliing back to prolate is easy, just check if b ~ c and then use the prolate analytic solutions.
  • Falling back to oblate is tricky because oblates are defined as a < b = c, while a triaxial that approximates an oblate is is a > b ~ c. Using the analytic solutions for the oblate as they are now is tricky (because we would need to apply a rotation as well). Maybe something to open an issue for.
• [ ] Add ellipsoids to the User Guide (same as the other source types).
• [ ] Possible optimization: there's no need to compute the internal demagnetization tensor if susceptibility is zero (or None) and if no observation point lies inside the ellipsoid. Maybe we can avoid computing those expensive functions all together.

Are you willing to help implement and maintain this feature?

Sure! I'll be continuing the work that @KellySavanna started in #568. Everyone is more than welcomed to help!

[leouieda]

Looks great! How about make_ellipsoid or create_ellipsoid instead of get_?

[santisoler]

I agree! Much better than get.

[johnweis0480]

Is there any plan to implement gradiometry data? It would be really useful for testing!

[santisoler]

Hi @johnweis0480! That's a nice thing we could add. We already have the math implemented for the magnetics, so the gravity gradiometry wouldn't be too hard to add.

We might want to think about the interface for it. But we can add it to the whishlist, for sure.