Merge branch 'main' into depth_bug_eqs

doc/api/index.rst

@@ -135,3 +135,4 @@ Utilities
 
     magnetic_vec_to_angles
     magnetic_angles_to_vec
+    total_field_anomaly

harmonica/__init__.py

@@ -31,7 +31,7 @@
     total_gradient_amplitude,
     upward_continuation,
 )
-from ._utils import magnetic_angles_to_vec, magnetic_vec_to_angles
+from ._utils import magnetic_angles_to_vec, magnetic_vec_to_angles, total_field_anomaly
 from ._version import __version__
 
 # Append a leading "v" to the generated version by setuptools_scm

harmonica/_utils.py

@@ -148,3 +148,56 @@ def magnetic_vec_to_angles(magnetic_e, magnetic_n, magnetic_u, degrees=True):
         (inclination,) = inclination
         (declination,) = declination
     return intensity, inclination, declination
+
+
+def total_field_anomaly(magnetic_field, inclination, declination):
+    r"""
+    The total field anomaly from the anomalous magnetic field.
+
+    Compute the total field anomaly from the anomalous magnetic field given the
+    regional field direction.
+
+    .. note::
+
+        Inclination is measured positive downward from the horizontal plane and
+        declination is measured with respect to North and it is positive east.
+
+    Parameters
+    ----------
+    magnetic_field : tuple of floats or tuple of arrays
+        Three component vector of the anomalous magnetic field.
+    inclination : float or array
+        Inclination angle of the regional field.
+        It must be in degrees.
+    declination : float or array
+        Declination angle of the regional field.
+        It must be in degrees.
+
+    Returns
+    -------
+    total_field_anomaly : float or array
+        The magnetic total field anomaly in the same units as the
+        ``magnetic_field``.
+
+    Notes
+    -----
+    Given the magnetic field, :math:`\mathbf{B}`, the regional field,
+    :math:`\mathbf{F}`, the total field anomaly can be computed as:
+
+    .. math::
+
+        \Delta T = \mathbf{B} \cdot \mathbf{\hat{F}}
+
+    where :math:`\mathbf{\hat{F}}` is the unit vector in the same direction
+    as the regional field.
+
+    Examples
+    --------
+    >>> tfa = total_field_anomaly([0, 0, -50e3], 90.0, 0.0)
+    >>> print(tfa)
+    50000.0
+    """
+    b_e, b_n, b_u = tuple(np.array(i) for i in magnetic_field)
+    f_e, f_n, f_u = magnetic_angles_to_vec(1, inclination, declination)
+    tfa = b_e * f_e + b_n * f_n + b_u * f_u
+    return tfa

harmonica/tests/test_utils.py

@@ -8,7 +8,7 @@
 import numpy.testing as npt
 import pytest
 
-from .. import magnetic_angles_to_vec, magnetic_vec_to_angles
+from .. import magnetic_angles_to_vec, magnetic_vec_to_angles, total_field_anomaly
 
 VECTORS = [
     [0.5, 0.5, -0.70710678],
@@ -114,3 +114,52 @@ def test_identity(arrays, start_with):
         npt.assert_almost_equal(
             magnetic_angles_to_vec(*angles), (magnetic_e, magnetic_n, magnetic_u)
         )
+
+
+@pytest.mark.parametrize("direction", ("easting", "northing", "upward"))
+def test_tfa(direction):
+    b = [30.0, -40.0, 50.0]
+    if direction == "easting":
+        inc, dec = 0.0, 90.0
+        expected_tfa = b[0]
+    elif direction == "northing":
+        inc, dec = 0.0, 0.0
+        expected_tfa = b[1]
+    else:
+        inc, dec = -90.0, 0.0
+        expected_tfa = b[2]
+    tfa = total_field_anomaly(b, inc, dec)
+    npt.assert_allclose(tfa, expected_tfa)
+
+
+@pytest.mark.parametrize("direction", ("easting", "northing", "upward"))
+def test_tfa_b_as_array(direction):
+    b = [[20, -30, -40, 50], [-60, 70, -80, 10], [21, -31, 41, -51]]
+    if direction == "easting":
+        inc, dec = 0.0, 90.0
+        expected_tfa = b[0]
+    elif direction == "northing":
+        inc, dec = 0.0, 0.0
+        expected_tfa = b[1]
+    else:
+        inc, dec = -90.0, 0.0
+        expected_tfa = b[2]
+    tfa = total_field_anomaly(b, inc, dec)
+    npt.assert_allclose(tfa, expected_tfa)
+
+
+def test_tfa_inc_dec_as_array():
+    be = [20.0, -30.0, -40.0, 50.0]
+    bn = [-60.0, 70.0, -80.0, 10.0]
+    bu = [21.0, -31.0, 41.0, -51.0]
+    b = [be, bn, bu]
+    inc = [0.0, 0.0, 45.0, -90.0]
+    dec = [90.0, 0.0, 45.0, 0.0]
+    tfa = total_field_anomaly(b, inc, dec)
+    expected_tfa = [
+        be[0],
+        bn[1],
+        0.5 * be[2] + 0.5 * bn[2] - 1 / np.sqrt(2) * bu[2],
+        bu[3],
+    ]
+    npt.assert_allclose(tfa, expected_tfa)