Add computation of gravitational tensor components for point sources (f…

…atiando#288)

Add new kernel functions to compute gravity tensor components generated by
point sources. Add test functions for the new feature: check that the diagonal
elements satisfy the Laplace equation, compare all components against finite
difference computations from the gravity acceleration. Add test class for
checking the symmetry of tensor components. Refactor old test functions for
point gravity: merge some functions into single ones through pytest
parametrizations. Avoid using "gradient" for specifying the gravity
acceleration vector: the "gravity gradient" is usually used to refer to the
tensor.

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,15 +371,105 @@ 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)
     )
     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