Repository: fatiando/harmonica
Pull request #119 (2019-11-04 16:49:50)

Add horizontal Cartesian components for point masses
Implement horizontal Cartesian components `g_northing`and `g_easting` for point masses. This PR addresses some of #85 . Components in spherical coordinate system will be implemented in the future. 


**Reminders**

- [x] Run `make format` and `make check` to make sure the code follows the style guide.
- [x] Add tests for new features or tests that would have caught the bug that you're fixing.
- [x] Add new public functions/methods/classes to `doc/api/index.rst`.
- [x] Write detailed docstrings for all functions/methods.
- [x] If adding new functionality, add an example to the docstring, gallery, and/or tutorials.


leouieda: @birocoles notice that your branch is out of date with respect to the master branch. You might want to merge in the changes before proceeding (you can do this locally or through Github)
birocoles: Hi @santisoler! I may have done something wrong when I tried to pull things from master. Actually, I think the sanity checks are important. Whatever. I will bring this sanity check and the related test back to the code.
birocoles: @santisoler , I think everything is working now.
birocoles: I have just fixed some problems with too long lines. What I have to do now?
santisoler: Thanks for tackling this @birocoles! :1st_place_medal: 

I think this is good to be merged. If you agree, please remove the WIP from the PR title and ask @leouieda for a review so we can merge it.
Thanks again for the hard work!
birocoles: Thanks @santisoler 

Patch:
diff --git a/harmonica/forward/point_mass.py b/harmonica/forward/point_mass.py
index 3aaa428..caa8ba2 100644
--- a/harmonica/forward/point_mass.py
+++ b/harmonica/forward/point_mass.py
@@ -44,17 +44,27 @@ def point_mass_gravity(
 
     .. math::
 
-        \vec{g} = \nabla V.
+        \vec{g} = \nabla V
 
-    Therefore, the :math:`z` component of :math:`\vec{g}` at the point
-    :math:`P` can be computed as (remember that :math:`z` points upwards):
+    and has components :math:`g_{northing}(P)`, :math:`g_{easting}(P)` and
+    :math:`g_{upward}(P)` given by
 
     .. math::
 
-        g_{up}(P) = - \frac{G m}{l^3} (z - z_p).
+        g_{northing}(P) = - \frac{G m}{l^3} (x - x_p),
+
+    .. math::
+
+        g_{easting}(P) = - \frac{G m}{l^3} (y - y_p)
+
+    and
+
+    .. math::
+
+        g_{upward}(P) = - \frac{G m}{l^3} (z - z_p).
 
     We define the downward component of the gravitational acceleration as the
-    opposite of :math:`g_{up}`:
+    opposite of :math:`g_{upward}` (remember that :math:`z` points upwards):
 
     .. math::
 
@@ -138,6 +148,8 @@ def point_mass_gravity(
 
         - Gravitational potential: ``potential``
         - Downward acceleration: ``g_z``
+        - Northing acceleration: ``g_northing``
+        - Easting acceleration: ``g_easting``
 
     coordinate_system : str (optional)
         Coordinate system of the coordinates of the computation points and the
@@ -166,6 +178,8 @@ def point_mass_gravity(
         "cartesian": {
             "potential": kernel_potential_cartesian,
             "g_z": kernel_g_z_cartesian,
+            "g_northing": kernel_g_northing_cartesian,
+            "g_easting": kernel_g_easting_cartesian,
         },
         "spherical": {
             "potential": kernel_potential_spherical,
@@ -195,7 +209,7 @@ def point_mass_gravity(
     )
     result *= GRAVITATIONAL_CONST
     # Convert to more convenient units
-    if field == "g_z":
+    if field in ("g_easting", "g_northing", "g_z"):
         result *= 1e5  # SI to mGal
     return result.reshape(cast.shape)
 
@@ -242,9 +256,10 @@ def kernel_potential_cartesian(
     """
     Kernel function for gravitational potential field in Cartesian coordinates
     """
-    return 1 / distance_cartesian(
+    distance = distance_cartesian(
         (easting, northing, upward), (easting_p, northing_p, upward_p)
     )
+    return 1 / distance
 
 
 @jit(nopython=True)
@@ -253,11 +268,43 @@ def kernel_g_z_cartesian(easting, northing, upward, easting_p, northing_p, upwar
     Kernel for downward component of gravitational gradient in Cartesian coords
     """
     distance = distance_cartesian(
-        [easting, northing, upward], [easting_p, northing_p, upward_p]
+        (easting, northing, upward), (easting_p, northing_p, upward_p)
     )
+    # Remember that the ``g_z`` field returns the downward component of the
+    # gravitational acceleration. As a consequence, it is multiplied by -1.
+    # Notice that the ``g_z`` does not have the minus signal observed at the
+    # compoents ``g_northing`` and ``g_easting``.
     return (upward - upward_p) / distance ** 3
 
 
+@jit(nopython=True)
+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
+    """
+    distance = distance_cartesian(
+        (easting, northing, upward), (easting_p, northing_p, upward_p)
+    )
+    return -(northing - northing_p) / distance ** 3
+
+
+@jit(nopython=True)
+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
+    """
+    distance = distance_cartesian(
+        (easting, northing, upward), (easting_p, northing_p, upward_p)
+    )
+    return -(easting - easting_p) / distance ** 3
+
+
 @jit(nopython=True)
 def jit_point_mass_spherical(
     longitude, latitude, radius, longitude_p, latitude_p, radius_p, masses, out, kernel


Test patch:
diff --git a/harmonica/tests/test_point_mass.py b/harmonica/tests/test_point_mass.py
index 3169335..9c875ff 100644
--- a/harmonica/tests/test_point_mass.py
+++ b/harmonica/tests/test_point_mass.py
@@ -7,6 +7,7 @@ import pytest
 
 from ..constants import GRAVITATIONAL_CONST
 from ..forward.point_mass import point_mass_gravity
+from ..forward.utils import distance_cartesian
 
 
 def test_invalid_coordinate_system():
@@ -342,6 +343,231 @@ def test_g_z_symmetry():
     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_mass_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_mass_gravity(
+        coordinates, point_mass, mass, "potential", "cartesian"
+    )
+    # 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)
+
+    # 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)
+
+
+@pytest.mark.use_numba
+def test_g_z_sign():
+    """
+    Test if g_z field of a positive point mass has 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
+    # bellow 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_mass_gravity(coordinates, point_mass, mass, "g_z", "cartesian")
+    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_mass_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():
+    """
+    Test the relative error in computing the g_northing component
+    """
+    # 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_mass_gravity(
+        coordinates_p, point_mass, mass, "g_northing", "cartesian"
+    )
+    # 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_mass_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)
+
+
+@pytest.mark.use_numba
+def test_g_northing_sign():
+    """
+    Test if g_northing field of a positive point mass has the correct sign
+    """
+    # 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_mass_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])
+
+
+@pytest.mark.use_numba
+def test_g_easting_symmetry():
+    """
+    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_mass_gravity(
+        coordinates, point_mass, masses, "g_easting", "cartesian"
+    )
+    npt.assert_allclose(results[0], -results[1])
+
+
+@pytest.mark.use_numba
+def test_g_easting_relative_error():
+    """
+    Test the relative error in computing the g_easting component
+    """
+    # 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_mass_gravity(
+        coordinates_p, point_mass, mass, "g_easting", "cartesian"
+    )
+    # 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_mass_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)
+
+
+@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]
+    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_mass_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])
+
+
 # ---------------------------
 # Spherical coordinates tests
 # ---------------------------