Add total_field_anomaly function (#510)

Add new function that computes the total field anomaly from the magnetic
field components and the inclination and declination of the IGRF.
Include the function in the API index and add tests for it.

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)