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 🙏