Forward modelling of ellipsoidal bodies for gravity and magnetic fields -

[KellySavanna]

This repo produces the magnetic and gravitational forward models of ellipsoidal bodies, much like the current Harmonica forward modelling packages. This was achieved with use of magnetic and gravitational ellipsoid derivations outlined in:

• Clark, D. A., Saul, S. J., & Emerson, D. W. (1986). Magnetic and gravity anomalies of a triaxial ellipsoid. Exploration Geophysics, 17(4), 189–200. https://doi.org/10.1071/EG986189
• Takahashi, D., & Oliveira Jr., V. C. (2017). Ellipsoids (v1.0): 3-D magnetic modelling of ellipsoidal bodies. Geoscientific Model Development, 10(9), 3591–3608. https://doi.org/10.5194/gmd-10-3591-2017.

There are three main .py files included: ellipsoid_magnetics.py for producing the total magnetic field due to ellipsoidal bodies, ellipsoid_gravity.py for the equivalent gravity fields, and create_ellipsoids.py to create the input ellipsoid for these functions. Tests for the code are included in /tests.

The user can specify the type of ellipsoid (triaxial: a > b > c, prolate: a > b = c, oblate a < b = c where a is the semi-major axis), the rotation components of the ellipsoid in intrinsic Euler angles (yaw, pitch and roll), and the centre of the ellipsoid in their given coordinate system. It also allows for multiple ellipsoidal bodies to be modelled, and produce and produce an array of the total anomaly as three components (be, bn, bu) or (ge, gn, gu) as (easting, northing, upward) magnetic and gravity components, respectfully. Users can also specify which component of the field to return (easting or northing or upward), and relevant input parameters such as susceptibility (isotropic or anisotropic), density and regional field $H_0$.

This code would be a useful addition to Harmonica, particularly in forward modelling for applications such as ore bodies.

! Work in progress !

Current issues in progress:

• test_magnetics_ellipsoids.py has function test_mag_ext_int_boundary. The test is to check continuity across a the internal-external boundary, which currently fails unless the tolerance is greater than 1e2. The max rel difference is given as about 4%, which remains true for several tested magnitudes of tested radius. This is currently being looked into as to whether it is a instability/rounding/precision error or a bug in the code.
• There are some areas which have potential for optimisation, e.g. currently the ellipsoid_magnetics function uses for-loops to avoid having to mask with an N x 3 x 3 matrix.
• Currently the functions produce total field, but perhaps a user specified option to choose total field or just the anomaly would be beneficial to the interface?

Open to suggestions! Thank you 🙏

[leouieda]

@KellySavanna thank you for the contribution and welcome to Fatiando! Great to have you here ❤️

I have enabled GitHub Actions for your pull request. Please take a look at the logs by clicking on each failed job to see what went wrong. Don't worry, this is very normal and we all get failing jobs almost every time.

I'll have a closer look at the code later but for now I can say that we don't need an option to return the total-field anomaly. The anomaly can already be calculated from the B field using https://www.fatiando.org/harmonica/latest/api/generated/harmonica.total_field_anomaly.html#harmonica.total_field_anomaly

Let me know if you need help with anything in particular!

[KellySavanna]

Thank you! Hopefully most of these issues have now been solved.

[nwilliams-kobold]

Great contribution @KellySavanna. I've been trying to reconcile the Takahashi and Clarke formulations, but now I don't have to! It would be super helpful to have a get_ellipsoid helper function that accepts 3 arbitrary axis dimensions, and then returns the appropriate class. It will commonly be the case that a user has an arbitrary ellipsoid, and just wants to simulate it, without figuring out whether it is prolate, oblate, triaxial. It would be even better if the helper function also handled the case of a == b == c by returning a spherical class (which may just be a wrapper of the existing point source and dipole classes). Is that doable?

[santisoler]

Hi @nwilliams-kobold. Thanks for your inputs on this PR.

We do plan to include a constructor function for ellipsoids that would work like you describe it.

@KellySavanna also started adding support for remnant magnetization. We were discussing whether the magnetization vector that is passed to the function should be defined in the global coordinate system (easting, northing, upward), or aligned with the ellipsoid's axes (the local coordinate system). The former would define a constant magnetization vector, regardless of the orientation of the ellipsoid. The latter would make the magnetization vector to rotate with the ellipsoid.

We are more inclined to define the magnetization vector in terms of its easting, northing and updated components, but we would like to hear your thoughts considering your potential applications for this code.

[nwilliams-kobold]

Super cool!

I definitely suggest defining the direction vector and ellipsoid rotation independently. We know they are degenerately-coupled (we can't resolve both at the same time with constraint), but we always hope we can find independent control of one or the other: e.g., paleomag measurements to constrain the vector, or positional or orientation constraints to constrain the ellipsoid shape and orientation. Thus allowing independent and natural inclusion of either directly (vector as east-north-up, and rotation as relative to east-north-up), without having to constantly rotate reference frames is very helpful, and less prone to error.

[santisoler]

@KellySavanna, I started pushing a few changes to your PR. Remember to pull them if you want to continue working on them. I'll probably keep working more on it this week. Thanks! I appreciate how neatly the code is, super easy to follow!

[santisoler]

@KellySavanna, congrats for this! Amazing to see it merged! 🚀 🏅

I'll keep working on the ellipsoids code in separate PRs. I have some ideas for it that I started listing in #594. If at any point you want to (and can) resume working on it, feel free to leave a comment and pick something to do.

One more thing, feel free to open a PR to add your name to the AUTHORS.md file. You can find more information about this file and authorship in Fatiando in our Autorship Guidelines.

[KellySavanna]

@santisoler thank you so much for all your help and support on this project!

[lheagy]

🚀 So excited to see this merged! A big congrats to @KellySavanna!