Repository: fatiando/harmonica (python)

Issue:
Function to compute the total field anomaly from magnetic field components
**Description of the desired feature:**

As discussed in today's meeting, it would be nice to include a function to compute the total field anomaly from the magnetic field components and with inclination and declination of the regional field. A function like such would be necessary to implement the ideas described in #426.

The function could be called `total_field_anomaly` and look something like:
```python
def total_field_anomaly(magnetic_field, inc, dec):
    b_e, b_n, b_u = magnetic_field
    ...
    tfa = ...
    return tfa
```

In order to get this done we should:
1. Implement the function in `harmonica/_utils.py`.
2. Make sure we are importing it in `harmonica/__init__.py` so it's accesible.
3. Add the new function to `doc/api/index.rst`.
4. Write tests for the function.


**Are you willing to help implement and maintain this feature?** 

From what we talked in today's meeting, @indiauppal is willing to take the lead for this 🚀. I'll let her assign herself to this issue if she wants to.


Pull request description:
Add `total_field_anomaly` function
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.

**Relevant issues/PRs:**
Fixes #506 

Patch:
diff --git a/doc/api/index.rst b/doc/api/index.rst
index 45bacd1..8a6770f 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -135,3 +135,4 @@ Utilities
 
     magnetic_vec_to_angles
     magnetic_angles_to_vec
+    total_field_anomaly
diff --git a/harmonica/__init__.py b/harmonica/__init__.py
index 6b73422..d578286 100644
--- a/harmonica/__init__.py
+++ b/harmonica/__init__.py
@@ -31,7 +31,7 @@ from ._transformations import (
     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
diff --git a/harmonica/_utils.py b/harmonica/_utils.py
index f5107be..a275bf5 100644
--- a/harmonica/_utils.py
+++ b/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


Test patch:
diff --git a/harmonica/tests/test_utils.py b/harmonica/tests/test_utils.py
index 8c08082..f97e0d3 100644
--- a/harmonica/tests/test_utils.py
+++ b/harmonica/tests/test_utils.py
@@ -8,7 +8,7 @@ import numpy as np
 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)