Add computation of gravitational tensor components for point sources

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