diff --git a/harmonica/forward/point.py b/harmonica/forward/point.py
index 4bc2991be..c328b89dc 100644
--- a/harmonica/forward/point.py
+++ b/harmonica/forward/point.py
@@ -164,6 +164,13 @@ def point_gravity(
         - Downward acceleration: ``g_z``
         - Northing acceleration: ``g_northing``
         - Easting acceleration: ``g_easting``
+        - Tensor components:
+            - ``g_ee``
+            - ``g_nn``
+            - ``g_zz``
+            - ``g_en``
+            - ``g_ez``
+            - ``g_nz``
 
     coordinate_system : str (optional)
         Coordinate system of the coordinates of the computation points and the
@@ -214,6 +221,9 @@ def point_gravity(
     # Convert to more convenient units
     if field in ("g_easting", "g_northing", "g_z"):
         result *= 1e5  # SI to mGal
+    tensors = ("g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz", "g_ne", "g_ze", "g_zn")
+    if field in tensors:
+        result *= 1e9  # SI to Eotvos
     return result.reshape(cast.shape)
 
 
@@ -275,6 +285,17 @@ def get_kernel(coordinate_system, field):
             "g_z": kernel_g_z_cartesian,
             "g_northing": kernel_g_northing_cartesian,
             "g_easting": kernel_g_easting_cartesian,
+            # diagonal tensor components
+            "g_ee": kernel_g_ee_cartesian,
+            "g_nn": kernel_g_nn_cartesian,
+            "g_zz": kernel_g_zz_cartesian,
+            # non-diagonal tensor components
+            "g_en": kernel_g_en_cartesian,
+            "g_ez": kernel_g_ez_cartesian,
+            "g_nz": kernel_g_nz_cartesian,
+            "g_ne": kernel_g_en_cartesian,
+            "g_ze": kernel_g_ez_cartesian,
+            "g_zn": kernel_g_nz_cartesian,
         },
         "spherical": {
             "potential": kernel_potential_spherical,
@@ -291,6 +312,11 @@ def get_kernel(coordinate_system, field):
     return kernel
 
 
+# ------------------------------------------
+# Kernel functions for Cartesian coordinates
+# ------------------------------------------
+
+
 @jit(nopython=True)
 def kernel_potential_cartesian(
     easting, northing, upward, easting_p, northing_p, upward_p
@@ -304,10 +330,16 @@ def kernel_potential_cartesian(
     return 1 / distance
 
 
+#  Acceleration components
+#  -----------------------
+
+
 @jit(nopython=True)
 def kernel_g_z_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
     """
-    Kernel for downward component of gravitational gradient in Cartesian coords
+    Kernel for downward component of gravitational acceleration
+
+    Use Cartesian coords.
     """
     distance = distance_cartesian(
         (easting, northing, upward), (easting_p, northing_p, upward_p)
@@ -324,8 +356,9 @@ def kernel_g_northing_cartesian(
     easting, northing, upward, easting_p, northing_p, upward_p
 ):
     """
-    Kernel function for northing component of gravitational gradient in
-    Cartesian coordinates
+    Kernel function for northing component of gravitational acceleration
+
+    Use Cartesian coordinates
     """
     distance = distance_cartesian(
         (easting, northing, upward), (easting_p, northing_p, upward_p)
@@ -338,8 +371,9 @@ def kernel_g_easting_cartesian(
     easting, northing, upward, easting_p, northing_p, upward_p
 ):
     """
-    Kernel function for easting component of gravitational gradient in
-    Cartesian coordinates
+    Kernel function for easting component of gravitational acceleration
+
+    Use Cartesian coordinates
     """
     distance = distance_cartesian(
         (easting, northing, upward), (easting_p, northing_p, upward_p)
@@ -347,6 +381,95 @@ def kernel_g_easting_cartesian(
     return -(easting - easting_p) / distance**3
 
 
+#  Tensor components
+#  -----------------
+
+
+@jit(nopython=True)
+def kernel_g_ee_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
+    """
+    Kernel function for g_ee component of gravitational tensor
+
+    Use Cartesian coordinates
+    """
+    distance = distance_cartesian(
+        (easting, northing, upward), (easting_p, northing_p, upward_p)
+    )
+    return 3 * (easting - easting_p) ** 2 / distance**5 - 1 / distance**3
+
+
+@jit(nopython=True)
+def kernel_g_nn_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
+    """
+    Kernel function for g_nn component of gravitational tensor
+
+    Use Cartesian coordinates
+    """
+    distance = distance_cartesian(
+        (easting, northing, upward), (easting_p, northing_p, upward_p)
+    )
+    return 3 * (northing - northing_p) ** 2 / distance**5 - 1 / distance**3
+
+
+@jit(nopython=True)
+def kernel_g_zz_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
+    """
+    Kernel function for g_zz component of gravitational tensor
+
+    Use Cartesian coordinates
+    """
+    distance = distance_cartesian(
+        (easting, northing, upward), (easting_p, northing_p, upward_p)
+    )
+    return 3 * (upward - upward_p) ** 2 / distance**5 - 1 / distance**3
+
+
+@jit(nopython=True)
+def kernel_g_en_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
+    """
+    Kernel function for g_en component of gravitational tensor
+
+    Use Cartesian coordinates
+    """
+    distance = distance_cartesian(
+        (easting, northing, upward), (easting_p, northing_p, upward_p)
+    )
+    return 3 * (easting - easting_p) * (northing - northing_p) / distance**5
+
+
+@jit(nopython=True)
+def kernel_g_ez_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
+    """
+    Kernel function for g_ez component of gravitational tensor
+
+    Use Cartesian coordinates
+    """
+    distance = distance_cartesian(
+        (easting, northing, upward), (easting_p, northing_p, upward_p)
+    )
+    # Add a minus sign to account that the z axis points downwards.
+    return -3 * (easting - easting_p) * (upward - upward_p) / distance**5
+
+
+@jit(nopython=True)
+def kernel_g_nz_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
+    """
+    Kernel function for g_nz component of gravitational tensor
+
+    Use Cartesian coordinates
+    """
+    distance = distance_cartesian(
+        (easting, northing, upward), (easting_p, northing_p, upward_p)
+    )
+    # Add a minus sign to account that the z axis points downwards.
+    return -3 * (northing - northing_p) * (upward - upward_p) / distance**5
+
+
+# ------------------------------------------
+# Kernel functions for Cartesian coordinates
+# ------------------------------------------
+
+
 @jit(nopython=True)
 def kernel_potential_spherical(
     longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
@@ -360,12 +483,18 @@ def kernel_potential_spherical(
     return 1 / distance
 
 
+#  Acceleration components
+#  -------------------
+
+
 @jit(nopython=True)
 def kernel_g_z_spherical(
     longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
 ):
     """
-    Kernel for downward component of gravitational gradient in spherical coords
+    Kernel for downward component of gravitational acceleration
+
+    Use spherical coordinates
     """
     distance, cospsi, _ = distance_spherical_core(
         longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
diff --git a/harmonica/tests/test_point_gravity.py b/harmonica/tests/test_point_gravity.py
index a71883a8e..128d66841 100644
--- a/harmonica/tests/test_point_gravity.py
+++ b/harmonica/tests/test_point_gravity.py
@@ -150,7 +150,7 @@ def test_point_mass_gravity_deprecated(point_mass, sample_coordinates_potential)
 
 
 @pytest.mark.use_numba
-def test_potential_cartesian_symmetry():
+def test_potential_symmetry_cartesian():
     """
     Test if potential field of a point mass has symmetry in Cartesian coords
     """
@@ -176,252 +176,461 @@ def test_potential_cartesian_symmetry():
 
 
 @pytest.mark.use_numba
-def test_g_z_symmetry():
+@pytest.mark.parametrize("field", ("g_northing", "g_easting", "g_z"))
+def test_acceleration_symmetry_cartesian(field):
     """
-    Test if g_z field of a point mass has symmetry in Cartesian coordinates
+    Test if the acceleration components verify the expected symmetry
+
+    Use Cartesian coordinates
     """
     # Define a single point mass
     point_mass = [1.1, 1.2, 1.3]
     masses = [2670]
-    # Define a pair of computation points above and below the point mass
+    # Define a pair of computation points at each side of the point mass along
+    # the direction given by the field parameter
     distance = 3.3
     easting = point_mass[0] * np.ones(2)
     northing = point_mass[1] * np.ones(2)
     upward = point_mass[2] * np.ones(2)
-    upward[0] += distance
-    upward[1] -= distance
+    if field == "g_northing":
+        northing[0] += distance
+        northing[1] -= distance
+    elif field == "g_easting":
+        easting[0] += distance
+        easting[1] -= distance
+    elif field == "g_z":
+        upward[0] += distance
+        upward[1] -= distance
     coordinates = [easting, northing, upward]
-    # Compute g_z gravity field on each computation point
-    results = point_gravity(coordinates, point_mass, masses, "g_z", "cartesian")
+    # Compute gravity acceleration component on each computation point
+    results = point_gravity(coordinates, point_mass, masses, field, "cartesian")
     npt.assert_allclose(results[0], -results[1])
 
 
-@pytest.mark.use_numba
-def test_g_z_relative_error():
-    """
-    Test the relative error in computing the g_z component
-    """
-    # Define a single point mass
-    point_mass = (1, -67, -300.7)
-    mass = 250
-    coordinates_p = (0, -39, -13)
-    # Compute the z component
-    exact_deriv = point_gravity(coordinates_p, point_mass, mass, "g_z", "cartesian")
-    # Compute the numerical derivative of potential
-    delta = 0.1
-    easting = np.zeros(2) + coordinates_p[0]
-    northing = np.zeros(2) + coordinates_p[1]
-    upward = np.array([coordinates_p[2] - delta, coordinates_p[2] + delta])
-    coordinates = (easting, northing, upward)
-    potential = point_gravity(coordinates, point_mass, mass, "potential", "cartesian")
+def acceleration_finite_differences(coordinates, point, mass, field, delta=0.05):
+    """
+    Compute acceleration components through finite differences
+
+    Parameters
+    ----------
+    coordinates : tuple
+        The coordinates of the computation point where the approximated
+        acceleration components will be computed.
+    point : tuple
+        The coordinates of the point source.
+    mass : float
+        Mass of the point source.
+    field : str
+        Acceleration component that needs to be approximated ("g_easting",
+        "g_northing", "g_z").
+    delta : float
+        Distance use to compute the finite difference in meters.
+
+    Returns
+    -------
+    finite_diff : float
+        Approximation of the acceleration component.
+    error : float
+        Relative error of the approximation (unitless).
+    """
+    # Build a two computation points slightly shifted from the original
+    # computation point by a small delta
+    coordinates_pair = tuple([coord, coord] for coord in coordinates)
+    if field == "g_easting":
+        index = 0
+    elif field == "g_northing":
+        index = 1
+    elif field == "g_z":
+        index = 2
+    coordinates_pair[index][0] -= delta
+    coordinates_pair[index][1] += delta
+    # Compute the potential on both points
+    potential = point_gravity(
+        coordinates_pair, point, mass, field="potential", coordinate_system="cartesian"
+    )
+    # Compute the difference between the two values
+    finite_diff = (potential[1] - potential[0]) / (2 * delta)
+    # Convert to mGal
+    finite_diff *= 1e5
     # Remember that the ``g_z`` field returns the downward component of the
     # gravitational acceleration. As a consequence, the numerical
-    # derivativative is multiplied by -1.
-    approximated_deriv = -1e5 * (potential[1] - potential[0]) / (2.0 * delta)
-
-    # Compute the relative error
-    relative_error = np.abs((approximated_deriv - exact_deriv) / exact_deriv)
+    # derivative is multiplied by -1.
+    if field == "g_z":
+        finite_diff *= -1
+    # Compute the bounding error of the approximation
+    distance = distance_cartesian(coordinates, point)
+    relative_error = 3 / 2 * (delta / distance) ** 2
+    return finite_diff, relative_error
 
-    # Bound value
-    distance = distance_cartesian(coordinates_p, point_mass)
-    bound_value = 1.5 * (delta / distance) ** 2
 
+@pytest.mark.use_numba
+@pytest.mark.parametrize("field", ("g_northing", "g_easting", "g_z"))
+@pytest.mark.parametrize(
+    "coordinates, point, mass",
+    (
+        [(0, -39, -13), (1, -67, -300.7), 250],
+        [(-3, 24, -10), (20, 54, -500.7), 200],
+    ),
+    ids=["set1", "set2"],
+)
+def test_acceleration_finite_diff_cartesian(coordinates, point, mass, field):
+    """
+    Test acceleration components against a finite difference of the potential
+    """
+    # Compute the z component
+    result = point_gravity(coordinates, point, mass, field, "cartesian")
+    # Compute the derivative of potential through finite differences
+    finite_diff, relative_error = acceleration_finite_differences(
+        coordinates, point, mass, field
+    )
     # Compare the results
-    npt.assert_array_less(relative_error, bound_value)
+    npt.assert_allclose(result, finite_diff, rtol=relative_error)
 
 
 @pytest.mark.use_numba
-def test_g_z_sign():
+@pytest.mark.parametrize("field", ("g_northing", "g_easting", "g_z"))
+def test_acceleration_sign(field):
     """
-    Test if g_z field of a positive point mass has the correct sign
+    Test if acceleration components have the correct sign
     """
     # Define a single point mass
     point_mass = [-10, 100.2, -300.7]
     mass = [2670]
-    # Define three computation points located above, at the same depth and
-    # below the point mass
-    easting = np.zeros(3)
-    northing = np.zeros(3) + 52.3
-    upward = np.array([100.11, -300.7, -400])
-    coordinates = [easting, northing, upward]
-    # Compute g_z gravity field on each computation point
-    results = point_gravity(coordinates, point_mass, mass, "g_z", "cartesian")
+    # Define computation points
+    coordinates = [np.zeros(3) for i in range(3)]
+    if field == "g_easting":
+        coordinates[0] = np.array([-150.7, -10, 79])
+    elif field == "g_northing":
+        coordinates[1] = np.array([0, 100.2, 210.7])
+    elif field == "g_z":
+        coordinates[2] = np.array([100.11, -300.7, -400])
+    # Compute acceleration component
+    results = point_gravity(coordinates, point_mass, mass, field, "cartesian")
+    # Check if the sign of the results is right
     assert np.sign(mass) == np.sign(results[0])
     npt.assert_allclose(results[1], 0)
     assert np.sign(mass) == -np.sign(results[2])
 
 
-@pytest.mark.use_numba
-def test_g_northing_symmetry():
-    """
-    Test if g_northing field of a point mass has symmetry in Cartesian
-    coordinates
-    """
-    # Define a single point mass
-    point_mass = [-7.9, 25, -130]
-    masses = [2670]
-    # Define a pair of computation points northward and southward the point
-    # mass
-    distance = 6.1
-    easting = point_mass[0] + np.zeros(2)
-    northing = point_mass[1] + np.zeros(2)
-    upward = point_mass[2] + np.zeros(2)
-    northing[0] += distance
-    northing[1] -= distance
-    coordinates = [easting, northing, upward]
-    # Compute g_northing gravity field on each computation point
-    results = point_gravity(coordinates, point_mass, masses, "g_northing", "cartesian")
-    npt.assert_allclose(results[0], -results[1])
-
-
-@pytest.mark.use_numba
-def test_g_northing_relative_error():
+@run_only_with_numba
+@pytest.mark.parametrize(
+    "field",
+    # fmt: off
+    (
+        "potential", "g_z", "g_northing", "g_easting",
+        "g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz"
+    ),
+    # fmt: on
+)
+def test_point_mass_cartesian_parallel(field):
     """
-    Test the relative error in computing the g_northing component
+    Check if parallel and serial runs return the same result
     """
-    # Define a single point mass
-    point_mass = (1, -67, -300.7)
-    mass = 250
-    coordinates_p = (0, -39, -13)
-    # Compute the northing component
-    exact_deriv = point_gravity(
-        coordinates_p, point_mass, mass, "g_northing", "cartesian"
+    region = (2e3, 10e3, -3e3, 5e3)
+    points = vd.scatter_points(region, size=30, extra_coords=-1e3, random_state=0)
+    masses = np.arange(points[0].size)
+    coordinates = vd.grid_coordinates(region=region, spacing=1e3, extra_coords=0)
+    result_serial = point_gravity(
+        coordinates, points, masses, field=field, parallel=False
     )
-    # Compute the numerical derivative of potential
-    delta = 0.1
-    easting = np.zeros(2) + coordinates_p[0]
-    northing = np.array([coordinates_p[1] - delta, coordinates_p[1] + delta])
-    upward = np.zeros(2) + coordinates_p[2]
-    coordinates = (easting, northing, upward)
-    potential = point_gravity(coordinates, point_mass, mass, "potential", "cartesian")
-    approximated_deriv = 1e5 * (potential[1] - potential[0]) / (2.0 * delta)
-
-    # Compute the relative error
-    relative_error = np.abs((approximated_deriv - exact_deriv) / exact_deriv)
-
-    # Bound value
-    distance = distance_cartesian(coordinates_p, point_mass)
-    bound_value = 1.5 * (delta / distance) ** 2
-
-    # Compare the results
-    npt.assert_array_less(relative_error, bound_value)
+    result_parallel = point_gravity(
+        coordinates, points, masses, field=field, parallel=True
+    )
+    npt.assert_allclose(result_serial, result_parallel)
 
 
 @pytest.mark.use_numba
-def test_g_northing_sign():
-    """
-    Test if g_northing field of a positive point mass has the correct sign
+def test_laplace_equation_cartesian():
     """
-    # Define a single point mass
-    point_mass = [-10, 100.2, -300.7]
-    mass = [2670]
-    # Define three computation points located above the point mass, along the
-    # north axis
-    easting = np.zeros(3)
-    northing = np.array([0, 100.2, 210.7])
-    upward = np.zeros(3)
-    coordinates = [easting, northing, upward]
-    # Compute g_northing gravity field on each computation point
-    results = point_gravity(coordinates, point_mass, mass, "g_northing", "cartesian")
-    assert np.sign(mass) == np.sign(results[0])
-    npt.assert_allclose(results[1], 0)
-    assert np.sign(mass) == -np.sign(results[2])
-
+    Check if the diagonal components of the tensor satisfy Laplace equation
 
-@pytest.mark.use_numba
-def test_g_easting_symmetry():
+    Use Cartesian coordinates.
     """
-    Test if g_easting field of a point mass has symmetry in Cartesian
-    coordinates
-    """
-    # Define a single point mass
-    point_mass = [191, -5, 0]
-    masses = [2670]
-    # Define a pair of computation points northward and southward the point
-    # mass
-    distance = 4.6
-    easting = point_mass[0] + np.zeros(2)
-    northing = point_mass[1] + np.zeros(2)
-    upward = point_mass[2] + np.zeros(2)
-    easting[0] += distance
-    easting[1] -= distance
-    coordinates = [easting, northing, upward]
-    # Compute g_easting gravity field on each computation point
-    results = point_gravity(coordinates, point_mass, masses, "g_easting", "cartesian")
-    npt.assert_allclose(results[0], -results[1])
+    region = (2e3, 10e3, -3e3, 5e3)
+    points = vd.scatter_points(region, size=30, extra_coords=-1e3, random_state=0)
+    masses = np.arange(points[0].size)
+    coordinates = vd.grid_coordinates(region=region, spacing=1e3, extra_coords=0)
+    g_ee = point_gravity(coordinates, points, masses, field="g_ee")
+    g_nn = point_gravity(coordinates, points, masses, field="g_nn")
+    g_zz = point_gravity(coordinates, points, masses, field="g_zz")
+    # Check if the Laplacian of the gravitational field is close to zero
+    npt.assert_allclose(g_ee + g_nn, -g_zz)
 
 
 @pytest.mark.use_numba
-def test_g_easting_relative_error():
+@pytest.mark.parametrize(
+    "field, flipped_field", [("g_en", "g_ne"), ("g_ez", "g_ze"), ("g_nz", "g_zn")]
+)
+def test_tensor_non_diagonal_components(field, flipped_field):
     """
-    Test the relative error in computing the g_easting component
+    Check if function computes g_xy as the same as g_yx
     """
-    # Define a single point mass
-    point_mass = (20, 54, -500.7)
-    mass = 200
-    coordinates_p = (-3, 24, -10)
-    # Compute the easting component
-    exact_deriv = point_gravity(
-        coordinates_p, point_mass, mass, "g_easting", "cartesian"
+    region = (2e3, 10e3, -3e3, 5e3)
+    points = vd.scatter_points(region, size=30, extra_coords=-1e3, random_state=0)
+    masses = np.arange(points[0].size)
+    coordinates = vd.grid_coordinates(region=region, spacing=1e3, extra_coords=0)
+    npt.assert_allclose(
+        point_gravity(coordinates, points, masses, field=field),
+        point_gravity(coordinates, points, masses, field=flipped_field),
     )
-    # Compute the numerical derivative of potential
-    delta = 0.1
-    easting = np.array([coordinates_p[0] - delta, coordinates_p[0] + delta])
-    northing = np.zeros(2) + coordinates_p[1]
-    upward = np.zeros(2) + coordinates_p[2]
-    coordinates = (easting, northing, upward)
-    potential = point_gravity(coordinates, point_mass, mass, "potential", "cartesian")
-    approximated_deriv = 1e5 * (potential[1] - potential[0]) / (2.0 * delta)
-
-    # Compute the relative error
-    relative_error = np.abs((approximated_deriv - exact_deriv) / exact_deriv)
-
-    # Bound value
-    distance = distance_cartesian(coordinates_p, point_mass)
-    bound_value = 1.5 * (delta / distance) ** 2
 
-    # Compare the results
-    npt.assert_array_less(relative_error, bound_value)
 
+class TestTensorSymmetryCartesian:
 
-@pytest.mark.use_numba
-def test_g_easting_sign():
-    """
-    Test if g_easting field of a positive point mass has the correct sign
-    """
-    # Define a single point mass
-    point_mass = [-10, 100.2, -300.7]
+    # Define sample point source and its mass
+    point = [1.1, 1.2, 1.3]
     mass = [2670]
-    # Define three computation points located above the point mass, along the
-    # east axis
-    easting = np.array([-150.7, -10, 79])
-    northing = np.zeros(3)
-    upward = np.zeros(3)
-    coordinates = [easting, northing, upward]
-    # Compute g_easting gravity field on each computation point
-    results = point_gravity(coordinates, point_mass, mass, "g_easting", "cartesian")
-    assert np.sign(mass) == np.sign(results[0])
-    npt.assert_allclose(results[1], 0)
-    assert np.sign(mass) == -np.sign(results[2])
 
+    # Define tensor components tuple
+    diagonal_fields = ["g_ee", "g_nn", "g_zz"]
+    nondiagonal_fields = ["g_en", "g_ez", "g_nz"]
+    tensor_fields = diagonal_fields + nondiagonal_fields
+
+    def mirrored_computation_points(self, directions):
+        """
+        Create mirrored computation points to the point source
+
+        The mirrored computation points will be mirror images along one of the
+        planes given by the tensor component.
+
+        .. code::
+
+                    |
+               *    |    *     m: point source of mass m
+                    |          *: each one of the mirrored computation points
+                    |
+            --------m--------
+                    |
+                    |
+                    |
+                    |
+
+        Parameters
+        ----------
+        directions : str
+            Directions along which the computation points will be mirrored.
+            The computation points will be shifted an equal distance along the
+            first direction and then they'll be mirrored along the second one.
+            For example: ``("n", "z")``
+
+        Returns
+        -------
+        coordinates : tuple
+            Tuple containing the coordinates of the mirrored computation points
+            in the following order: easting, northing, upward.
+        """
+        distance = 3.3
+        easting = self.point[0] * np.ones(2)
+        northing = self.point[1] * np.ones(2)
+        upward = self.point[2] * np.ones(2)
+        direction_shift, direction_mirror = directions[:]
+        # Shift computation points
+        if direction_shift == "n":
+            northing[0] += distance
+            northing[1] += distance
+        elif direction_shift == "e":
+            easting[0] += distance
+            easting[1] += distance
+        elif direction_shift == "z":
+            upward[0] += distance
+            upward[1] += distance
+        # Mirror computation points
+        if direction_mirror == "n":
+            northing[0] += distance
+            northing[1] -= distance
+        elif direction_mirror == "e":
+            easting[0] += distance
+            easting[1] -= distance
+        elif direction_mirror == "z":
+            upward[0] += distance
+            upward[1] -= distance
+        return (easting, northing, upward)
+
+    def opposite_computation_points(self, directions):
+        """
+        Create opposite computation points to the point source
+
+        The opposite computation points will live in the diagonal given by the
+        direction of the tensor component and be equidistant to the point
+        source.
+        For example:
+
+        * if ``directions=("e", "n")``, the two computation points will be
+          along the easting-northing diagonal, equidistant to the point mass.
+        * if ``directions=("n", "z")``, the two computation points will be
+          along the northing-upward diagonal, equidistant to the point mass.
+
+        .. code::
+
+                    |
+                    |    *     m: point source of mass m
+                    |          *: each one of the opposite computation points
+                    |
+            --------m--------
+                    |
+                    |
+               *    |
+                    |
+
+        Parameters
+        ----------
+        directions : tuple
+            Tuple with the directions along which the opposition will be
+            carried out. For example: ``("n", "z")``
+
+        Returns
+        -------
+        coordinates : tuple
+            Tuple containing the coordinates of the opposite computation points
+            in the following order: easting, northing, upward.
+        """
+        distance = 3.3
+        easting = self.point[0] * np.ones(2)
+        northing = self.point[1] * np.ones(2)
+        upward = self.point[2] * np.ones(2)
+        if "n" in directions:
+            northing[0] += distance
+            northing[1] -= distance
+        if "e" in directions:
+            easting[0] += distance
+            easting[1] -= distance
+        if "z" in directions:
+            upward[0] += distance
+            upward[1] -= distance
+        return (easting, northing, upward)
+
+    @pytest.mark.use_numba
+    @pytest.mark.parametrize("field", tensor_fields)
+    def test_opposite(self, field):
+        """
+        Test tensor components symmetry on opposite computation points
+
+        The values of each tensor component should be the same on opposite
+        points.
+        """
+        # Define opposite computation points
+        directions = (field[-2], field[-1])
+        coordinates = self.opposite_computation_points(directions)
+        # Compute gravity tensor component on each computation point
+        results = point_gravity(coordinates, self.point, self.mass, field, "cartesian")
+        # Check for expected symmetry
+        npt.assert_allclose(results[0], results[1])
+
+    @pytest.mark.use_numba
+    @pytest.mark.parametrize("field", nondiagonal_fields)
+    @pytest.mark.parametrize("swap_directions", (False, True), ids=["noswap", "swap"])
+    def test_mirrored_nondiagonals(self, field, swap_directions):
+        """
+        Test nondiagonal tensor components symmetry on mirrored points
+
+        For non-diagonal tensor components, their values should be opposite on
+        mirrored computation points.
+        """
+        # Choose the direction corresponding to the diagonal tensor component
+        directions = (field[-2], field[-1])
+        if swap_directions:
+            directions = directions[::-1]
+        # Define mirrored computation points
+        coordinates = self.mirrored_computation_points(directions)
+        # Compute gravity tensor component on each computation point
+        results = point_gravity(coordinates, self.point, self.mass, field, "cartesian")
+        # Check for expected symmetry
+        npt.assert_allclose(results[0], -results[1])
+
+
+def tensor_finite_differences(coordinates, point, mass, field, delta=0.05):
+    """
+    Compute tensor components through finite differences
+
+    Parameters
+    ----------
+    coordinates : tuple
+        The coordinates of the computation point where the approximated
+        tensor components will be computed.
+    point : tuple
+        The coordinates of the point source.
+    mass : float
+        Mass of the point source.
+    field : str
+        Tensor component that needs to be approximated ("g_ee",
+        "g_nn", "g_zz", "g_en", "g_ez", "g_nz").
+    delta : float
+        Distance use to compute the finite difference in meters.
+
+    Returns
+    -------
+    finite_diff : float
+        Approximation of the tensor component.
+    error : float
+        Relative error of the approximation (unitless).
+    """
+    # Determine the direction along which the finite difference will be
+    # computed
+    direction_i, direction_j = field[-2], field[-1]
+    if direction_i == "e":
+        direction_i = "easting"
+    if direction_i == "n":
+        direction_i = "northing"
+    if direction_j == "e":
+        direction_j = "easting"
+    if direction_j == "n":
+        direction_j = "northing"
+    # Build a two computation points slightly shifted from the original
+    # computation point by a small delta
+    coordinates_pair = tuple([coord, coord] for coord in coordinates)
+    if direction_j == "easting":
+        index = 0
+    elif direction_j == "northing":
+        index = 1
+    elif direction_j == "z":
+        index = 2
+    coordinates_pair[index][0] -= delta
+    coordinates_pair[index][1] += delta
+    # Compute the acceleration on both points
+    acceleration = point_gravity(
+        coordinates_pair,
+        point,
+        mass,
+        field=f"g_{direction_i}",
+        coordinate_system="cartesian",
+    )
+    # Compute the difference between the two values
+    finite_diff = (acceleration[1] - acceleration[0]) / (2 * delta)
+    # Convert to Eotvos
+    finite_diff *= 1e-5  # convert back from mGal to SI
+    finite_diff *= 1e9  # convert to Eotvos
+    # The z axis points downwards, so the finite_diff should be multiplied by
+    # -1
+    if direction_i == "z" or direction_j == "z":
+        finite_diff *= -1
+    # Compute the bounding error of the approximation
+    distance = distance_cartesian(coordinates, point)
+    relative_error = 6 * (delta / distance) ** 2
+    return finite_diff, relative_error
 
-@run_only_with_numba
-def test_point_mass_cartesian_parallel():
-    """
-    Check if parallel and serial runs return the same result
+
+@pytest.mark.use_numba
+@pytest.mark.parametrize("field", ("g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz"))
+@pytest.mark.parametrize(
+    "coordinates, point, mass",
+    (
+        [(0, -39, -13), (1, -67, -300.7), 250],
+        [(-3, 24, -10), (20, 54, -500.7), 200],
+    ),
+    ids=["set1", "set2"],
+)
+def test_tensor_finite_diff_cartesian(coordinates, point, mass, field):
+    """
+    Test tensor components against a finite difference of the acceleration
     """
-    region = (2e3, 10e3, -3e3, 5e3)
-    points = vd.scatter_points(region, size=30, extra_coords=-1e3, random_state=0)
-    masses = np.arange(points[0].size)
-    coordinates = vd.grid_coordinates(region=region, spacing=1e3, extra_coords=0)
-    for field in ("potential", "g_z", "g_northing", "g_easting"):
-        result_serial = point_gravity(
-            coordinates, points, masses, field=field, parallel=False
-        )
-        result_parallel = point_gravity(
-            coordinates, points, masses, field=field, parallel=True
-        )
-        npt.assert_allclose(result_serial, result_parallel)
+    # Compute the z component
+    result = point_gravity(coordinates, point, mass, field, "cartesian")
+    # Compute the derivative of potential through finite differences
+    finite_diff, relative_error = tensor_finite_differences(
+        coordinates, point, mass, field
+    )
+    # Compare the results
+    npt.assert_allclose(result, finite_diff, rtol=relative_error)
 
 
 # ---------------------------
@@ -453,7 +662,8 @@ def test_point_mass_on_origin():
 
 
 @pytest.mark.use_numba
-def test_point_mass_same_radial_direction():
+@pytest.mark.parametrize("field", ("potential", "g_z"))
+def test_point_mass_same_radial_direction(field):
     """
     Check potential and g_z of point mass and computation point on same radius
     """
@@ -475,13 +685,10 @@ def test_point_mass_same_radial_direction():
                     "g_z": GRAVITATIONAL_CONST * mass / height**2 * 1e5,
                 }
                 # Compare results with analytical solutions
-                for field, solution in analytical.items():
-                    npt.assert_allclose(
-                        point_gravity(
-                            coordinates, point_mass, mass, field, "spherical"
-                        ),
-                        solution,
-                    )
+                npt.assert_allclose(
+                    point_gravity(coordinates, point_mass, mass, field, "spherical"),
+                    analytical[field],
+                )
 
 
 @pytest.mark.use_numba