diff --git a/doc/api/index.rst b/doc/api/index.rst
index cf87e025b..cca7df4c4 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -125,4 +125,5 @@ Utilities
 .. autosummary::
    :toctree: generated/
 
-    test
+    magnetic_vec_to_angles
+    magnetic_angles_to_vec
diff --git a/harmonica/__init__.py b/harmonica/__init__.py
index 62f6901be..97de1f8b5 100644
--- a/harmonica/__init__.py
+++ b/harmonica/__init__.py
@@ -29,6 +29,7 @@
     reduction_to_pole,
     upward_continuation,
 )
+from ._utils import magnetic_angles_to_vec, magnetic_vec_to_angles
 from ._version import __version__
 
 
diff --git a/harmonica/_utils.py b/harmonica/_utils.py
new file mode 100644
index 000000000..f5107bea3
--- /dev/null
+++ b/harmonica/_utils.py
@@ -0,0 +1,150 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+import numpy as np
+
+
+def magnetic_angles_to_vec(intensity, inclination, declination):
+    """
+    Convert magnetic field angles to magnetic field vector
+
+    Convert intensity, inclination and declination angles of the magnetic field
+    to a 3-component magnetic vector.
+
+    .. note::
+
+        Inclination is measured positive downward from the horizontal plane,
+        and declination is measured with respect to north, where positive
+        angles indicate east.
+
+    Parameters
+    ----------
+    intensity: float or array
+        Intensity (norm) of the magnetic vector in A/m.
+    inclination : float or array
+        Inclination angle of the magnetic vector in degree.
+        It must be in ``degrees``.
+    declination : float or array
+        Declination angle of the magnetic vector.
+        It must be in ``degrees``.
+
+    Returns
+    -------
+    magnetic_e : float or array
+        Easting component of the magnetic vector.
+    magnetic_n : float or array
+        Northing component of the magnetic vector.
+    magnetic_u : float or array
+        Upward component of the magnetic vector.
+
+    Examples
+    --------
+    >>> mag_e, mag_n, mag_u = magnetic_angles_to_vec(3.0, 45.0, 45.0)
+    >>> print(mag_e, mag_n, mag_u)
+    1.5 1.5 -2.121
+    """
+    # Transform to radians
+    inc_rad = np.radians(inclination)
+    dec_rad = np.radians(declination)
+    # Calculate the 3 components
+    magnetic_e = intensity * np.cos(inc_rad) * np.sin(dec_rad)
+    magnetic_n = intensity * np.cos(inc_rad) * np.cos(dec_rad)
+    magnetic_u = -intensity * np.sin(inc_rad)
+    return magnetic_e, magnetic_n, magnetic_u
+
+
+def magnetic_vec_to_angles(magnetic_e, magnetic_n, magnetic_u, degrees=True):
+    r"""
+    Convert magnetic field vector to magnetic field angles
+
+    Convert the 3-component magnetic vector to intensity, and inclination and
+    declination angles.
+
+    .. 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_e : float or array
+        Easting component of the magnetic vector.
+    magnetic_n : float or array
+        Northing component of the magnetic vector.
+    magnetic_u : float or array
+        Upward component of the magnetic vector.
+    degrees : bool (optional)
+        If True, the angles are returned in degrees.
+        If False, the angles are returned in radians.
+        Default True.
+
+    Returns
+    -------
+    intensity: float or array
+        Intensity of the magnetic vector.
+    inclination : float or array
+        Inclination angle of the magnetic vector.
+        If ``degrees`` is True, then the angle is returned in degree, else it's
+        returned in radians.
+    declination : float or array
+        Declination angle of the magnetic vector.
+        If ``degrees`` is True, then the angle is returned in degrees, else
+        it's returned in radians.
+
+    Notes
+    -----
+    The intensity of the magnetic vector is calculated as:
+
+    .. math::
+
+        T = \sqrt{B_e^2 + B_n^2 + B_u^2}
+
+    where :math:`B_e`, :math:`B_n`, :math:`B_u` are the easting, northing and
+    upward components of the magnetic vector, respectively.
+
+    The inclination angle is defined as the angle between the magnetic field
+    vector and the horizontal plane:
+
+    .. math::
+
+        I = \arctan \frac{- B_u}{\sqrt{B_e^2 + B_n^2}}
+
+    And the declination angle is defined as the azimuth of the projection of
+    the magnetic field vector onto the horizontal plane (starting from the
+    northing direction, positive to the east and negative to the west):
+
+    .. math::
+
+        D = \arcsin \frac{B_e}{\sqrt{B_e^2 + B_n^2}}
+
+    Examples
+    --------
+    >>> intensity, inc, dec = magnetic_vec_to_angles(1.5, 1.5, -2.12132)
+    >>> print(intensity, inc, dec)
+    3.0 45.0 45.0
+    """
+    # Compute the intensity as a norm
+    intensity = np.sqrt(magnetic_e**2 + magnetic_n**2 + magnetic_u**2)
+    # Compute the horizontal component of the magnetic vector
+    horizontal_component = np.atleast_1d(np.sqrt(magnetic_e**2 + magnetic_n**2))
+    # Mask the values equal to zero in the horizontal component
+    horizontal_component = np.ma.masked_values(horizontal_component, 0.0)
+    # Calculate the inclination and declination using the mask
+    inclination = np.arctan(-magnetic_u / horizontal_component)
+    declination = np.arcsin(magnetic_e / horizontal_component)
+    # Fill the masked values
+    inclination = inclination.filled(-np.sign(magnetic_u) * np.pi / 2)
+    declination = declination.filled(0)
+    # Convert to degree if needed
+    if degrees:
+        inclination = np.degrees(inclination)
+        declination = np.degrees(declination)
+    # Cast to floats if all components are floats
+    if intensity.ndim == 0:
+        (inclination,) = inclination
+        (declination,) = declination
+    return intensity, inclination, declination
diff --git a/harmonica/tests/test_utils.py b/harmonica/tests/test_utils.py
new file mode 100644
index 000000000..8c0808204
--- /dev/null
+++ b/harmonica/tests/test_utils.py
@@ -0,0 +1,116 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+import numpy as np
+import numpy.testing as npt
+import pytest
+
+from .. import magnetic_angles_to_vec, magnetic_vec_to_angles
+
+VECTORS = [
+    [0.5, 0.5, -0.70710678],
+    [0.5, 0.5, 0.70710678],
+    [-0.5, 0.5, -0.70710678],
+    [0, 0, -1],  # Over -z axis
+    [1, 0, 0],  # Over east (y) axis
+    [0, 1, 0],  # Over north (x) axis
+]
+
+ANGLES = [[1, 45, 45], [1, -45, 45.0], [1, 45, -45], [1, 90, 0], [1, 0, 90], [1, 0, 0]]
+
+
+@pytest.mark.parametrize("angles, vector", [(a, v) for a, v in zip(ANGLES, VECTORS)])
+def test_magnetic_ang_to_vec_float(angles, vector):
+    """
+    Check if the function returns the expected values for a given intensity
+    inclination and declination as float
+    """
+    intensity, inclination, declination = angles
+    magnetic_e, magnetic_n, magnetic_u = vector
+    npt.assert_almost_equal(
+        magnetic_angles_to_vec(intensity, inclination, declination),
+        (magnetic_e, magnetic_n, magnetic_u),
+    )
+
+
+@pytest.mark.parametrize("degrees", (False, True), ids=("radians", "degrees"))
+@pytest.mark.parametrize("angles, vector", [(a, v) for a, v in zip(ANGLES, VECTORS)])
+def test_magnetic_vec_to_angles_float(angles, vector, degrees):
+    """
+    Check if the function returns the expected values for a given magnetic
+    vector as float
+    """
+    intensity, inclination, declination = angles
+    magnetic_e, magnetic_n, magnetic_u = vector
+    if not degrees:
+        inclination, declination = np.radians(inclination), np.radians(declination)
+    npt.assert_allclose(
+        magnetic_vec_to_angles(magnetic_e, magnetic_n, magnetic_u, degrees=degrees),
+        (intensity, inclination, declination),
+    )
+
+
+@pytest.fixture(name="arrays", params=["single-element", "multi-element"])
+def angles_vectors_as_arrays(request):
+    """
+    Generate magnetic angles and vectors as arrays
+    """
+    if request.param == "single-element":
+        # Generate arrays with a single element
+        intensity, inclination, declination = tuple(np.atleast_1d(i) for i in ANGLES[0])
+        magnetic_e, magnetic_n, magnetic_u = tuple(np.atleast_1d(i) for i in VECTORS[0])
+    else:
+        intensity, inclination, declination = np.vstack(ANGLES).T
+        magnetic_e, magnetic_n, magnetic_u = np.vstack(VECTORS).T
+    return (intensity, inclination, declination), (magnetic_e, magnetic_n, magnetic_u)
+
+
+def test_magnetic_ang_to_vec_array(arrays):
+    """
+    Check if the function returns the expected values for a given intensity,
+    inclination and declination a array
+    """
+    intensity, inclination, declination = arrays[0]
+    magnetic_e, magnetic_n, magnetic_u = arrays[1]
+    npt.assert_almost_equal(
+        magnetic_angles_to_vec(intensity, inclination, declination),
+        (magnetic_e, magnetic_n, magnetic_u),
+    )
+
+
+@pytest.mark.parametrize("degrees", (False, True), ids=("radians", "degrees"))
+def test_magnetic_vec_to_angles_array(arrays, degrees):
+    """
+    Check if the function returns the expected values for the given magnetic
+    vector as arrays
+    """
+    intensity, inclination, declination = arrays[0]
+    magnetic_e, magnetic_n, magnetic_u = arrays[1]
+    if not degrees:
+        inclination, declination = np.radians(inclination), np.radians(declination)
+    npt.assert_allclose(
+        magnetic_vec_to_angles(magnetic_e, magnetic_n, magnetic_u, degrees=degrees),
+        (intensity, inclination, declination),
+    )
+
+
+@pytest.mark.parametrize("start_with", ("angles", "vectors"))
+def test_identity(arrays, start_with):
+    """
+    Check if applying both conversions return the original set of vectors
+    """
+    if start_with == "angles":
+        intensity, inclination, declination = arrays[0]
+        vector = magnetic_angles_to_vec(intensity, inclination, declination)
+        npt.assert_almost_equal(
+            magnetic_vec_to_angles(*vector), (intensity, inclination, declination)
+        )
+    else:
+        magnetic_e, magnetic_n, magnetic_u = arrays[1]
+        angles = magnetic_vec_to_angles(magnetic_e, magnetic_n, magnetic_u)
+        npt.assert_almost_equal(
+            magnetic_angles_to_vec(*angles), (magnetic_e, magnetic_n, magnetic_u)
+        )
diff --git a/pyproject.toml b/pyproject.toml
index 6a5b0e2aa..edf1519ec 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -10,6 +10,7 @@ local_scheme =  "node-and-date"
 write_to =  "harmonica/_version_generated.py"
 
 [tool.pytest.ini_options]
+doctest_optionflags = "NUMBER"
 markers = [
     "use_numba: mark test functions that call Numba jitted functions"
 ]