Improvements to ellipsoids

[santisoler]

Description of the desired feature:

Listing here some ideas to improve ellipsoids code:

• Move the physical properties of the ellipsoids to the objects. Remove them from the functions calls (Move physical properties to ellipsoid classes #597).
• Expose the rotation matrix through a property in ellipsoid classes (Move physical properties to ellipsoid classes #597).
• 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 (Optimize ellipsoid's magnetic forward modelling #603).
• Optimize computation of lambda for prolate and oblate: we can solve the second-order equation rather than the third-order one (Optimize lambda for prolate and oblate ellipsoids #604).
• Add a Sphere class that represents an homogeneous sphere as a source (Add a new Sphere class as a special case of ellipsoids #602).
  • Use that class in tests, when comparing ellipsoids to spheres.
• 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). (Add a new create_ellipsoid function #606)
• 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 (Solve numerical instabilities of almost spheroids #611).
• Add method to export ellipsoids to PyVista (Add method to export ellipsoids as PyVista objects #614).
• Add a __str__ method to ellipsoids, so we can print them and their properties in a nice looking way (Add __str__ method to ellipsoid classes #615).
• Refactor code and merge all ellipsoid classes into a single Ellipsoid class (Refactor code to have a single Ellipsoid class #616).
  • Ditch the create_ellipsoid function since it's not longer needed.
  • Don't impose order in semiaxes lenghts nor extra rotation angles.
  • Apply extra rotation before forward modelling, so local coordinate system x, y, z are aligned with a, b, c where a >= b >= c, for any ellipsoid type.
• Solve numerical instabilities when triaxial ellipsoids approximate an oblate or a prolate. Basically when two of the three semiaxes are close to each other (tackled by Refactor code to have a single Ellipsoid class #616).
  • 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. This was tackled by Refactor code to have a single Ellipsoid class #616. After standardizing the definition of the ellipsoids, it was easy to fallback to oblate solutions when triaxial is close enough to an oblate.
• Standardize the external field input as the magnetic field components in nT, instead of magnitude, inclination and declination. Follow a similar scheme as the one used in other functions (Modify inducing field argument in ellipsoids mag forward #622).
• 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.
• Add a __repr__ method to Ellipsoid, so we can print them in a single line when showed in error messages (Add Ellipsoid.__repr__() method #627)

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.