Repository: fatiando/harmonica

Issue:
Point masses gravity computation on Cartesian coordinates
**Description of the desired feature**

The most simple geometry that generates gravitational fields are point masses.
We should make a similar implementation that the one proposed on #47 but for Cartesian coordinates.

**Are you willing to help implement and maintain this feature?** Yes/No

Discussion:
santisoler: @leouieda Could you take a quick look at this. Specially to the design of the extended `point_mass_gravity()` function. I'm not quite sure about the dictionaries I defined, so I wanted to have a second opinion.

I added some functions to compute distances in Cartesian coordinates inside `forward/point_mass.py`. I intend to move them as proposed in #69.

If you can't find time to review it, let me know. I think the main issues that can live here are related with design, so maybe I could merge it and we can discuss it afterwards when you have time.
leouieda: Looks good! It would be good to have a gallery example as well to visually verify that the results are correct. Maybe one with 2 masses at different depths with opposite densities?
santisoler: Thanks for the quick review @leouieda. I've added a simple example as you proposed. If you like it, feel free to merge this.
santisoler: @leouieda Would you mind taking a final look at this? If you agree with the new variable name and don't have any other suggestion, I think this is ready to be merged.
leouieda: @santisoler this is great! I love that docstring! Made a few corrections to the text but once these are done feel free to merge!
santisoler: Thanks for the review @leouieda! I like how it looks, so I'll merge it right away!

Patch:
diff --git a/doc/references.rst b/doc/references.rst
index 15b82b5..ba43242 100644
--- a/doc/references.rst
+++ b/doc/references.rst
@@ -3,7 +3,9 @@ References
 
 .. [AmanteEakins2009] Amante, C. and B.W. Eakins, 2009. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS NGDC-24. National Geophysical Data Center, NOAA. doi:`10.7289/V5C8276M <https://doi.org/10.7289/V5C8276M>`__
 .. [BarthelmesKohler2016] Barthelmes, F. and Kohler, W. (2016), International Centre for Global Earth Models (ICGEM), in: Drewes, H., Kuglitsch, F., Adam, J. et al., The Geodesists Handbook 2016, Journal of Geodesy (2016), 90(10), pp 907-1205, doi:`10.1007/s00190-016-0948-z <https://doi.org/10.1007/s00190-016-0948-z>`__
+.. [Blakely1995] Blakely, R. (1995). Potential Theory in Gravity and Magnetic Applications. Cambridge: Cambridge University Press. doi:`10.1017/CBO9780511549816 <https://doi.org/10.1017/CBO9780511549816>`__
 .. [Forste_etal2014] Förste, Christoph; Bruinsma, Sean.L.; Abrikosov, Oleg; Lemoine, Jean-Michel; Marty, Jean Charles; Flechtner, Frank; Balmino, G.; Barthelmes, F.; Biancale, R. (2014): EIGEN-6C4 The latest combined global gravity field model including GOCE data up to degree and order 2190 of GFZ Potsdam and GRGS Toulouse. GFZ Data Services. doi:`10.5880/icgem.2015.1 <http://doi.org/10.5880/icgem.2015.1>`__
+.. [Grombein2013] Grombein, T., Seitz, K., Heck, B. (2013), Optimized formulas for the gravitational field of a tesseroid, Journal of Geodesy. doi:`10.1007/s00190-013-0636-1 <https://doi.org/10.1007/s00190-013-0636-1>`__
 .. [Hofmann-WellenhofMoritz2006] Hofmann-Wellenhof, B., & Moritz, H. (2006). Physical Geodesy (2nd, corr. ed. 2006 edition ed.). Wien ; New York: Springer.
 .. [LiGotze2001] Li, X. and H. J. Gotze, 2001, Tutorial: Ellipsoid, geoid, gravity, geodesy, and geophysics, Geophysics, 66(6), p. 1660-1668, doi:`10.1190/1.1487109 <https://doi.org/10.1190/1.1487109>`__
 .. [TurcotteSchubert2014] Turcotte, D. L., & Schubert, G. (2014). Geodynamics (3 edition). Cambridge, United Kingdom: Cambridge University Press.
diff --git a/examples/forward/point_mass.py b/examples/forward/point_mass.py
new file mode 100644
index 0000000..4c8f91a
--- /dev/null
+++ b/examples/forward/point_mass.py
@@ -0,0 +1,54 @@
+"""
+Point Mass in Cartesian Coordinates
+===================================
+
+The simplest geometry used to compute gravitational fields are point masses.
+Although modelling geologic structures with point masses can be challenging, they are
+very usefull for other purposes: creating synthetic models, solving inverse problems
+very quickly, generating equivalent sources for interpolation or gridding, etc.
+The gravitational fields generated by point masses can be quickly computed
+either in Cartesian or in geocentric spherical coordinate systems.
+We will compute the downward component of the gravitational acceleration generated by
+a set of point masses on a computation grid given in Cartesian coordinates. We will do
+it throught the :func:`harmonica.point_mass_gravity` function.
+"""
+import harmonica as hm
+import verde as vd
+import matplotlib.pyplot as plt
+import matplotlib.ticker
+
+
+# Define two point masses with oposite mass values of 10000000 kg
+easting = [5e3, 15e3]
+northing = [5e3, 15e3]
+down = [7e3, 2.5e3]
+points = [easting, northing, down]
+# We're using "negative" masses to represent a "mass deficit"
+masses = [10e6, -10e6]
+
+# Define computation points on a grid at 500m above the ground
+coordinates = vd.grid_coordinates(
+    region=[0, 20e3, 0, 20e3], shape=(80, 80), extra_coords=-500
+)
+
+# Compute the downward component of the acceleration
+gravity = hm.point_mass_gravity(
+    coordinates, points, masses, field="g_z", coordinate_system="cartesian"
+)
+print(gravity)
+
+# Plot the gravitational field
+fig, ax = plt.subplots(figsize=(8, 9))
+ax.set_aspect("equal")
+img = plt.pcolormesh(*coordinates[:2], gravity)
+# Add colorbar
+fmt = matplotlib.ticker.ScalarFormatter(useMathText=True)
+fmt.set_powerlimits((0, 0))
+plt.colorbar(img, ax=ax, format=fmt, pad=0.04, shrink=0.73, label="mGal")
+ax.set_title("Downward component of gravitational acceleration")
+# Convert axes units to km
+ax.set_xticklabels(ax.get_xticks() * 1e-3)
+ax.set_yticklabels(ax.get_yticks() * 1e-3)
+ax.set_xlabel("km")
+ax.set_ylabel("km")
+plt.show()
diff --git a/harmonica/forward/point_mass.py b/harmonica/forward/point_mass.py
index f6a200e..5b12660 100644
--- a/harmonica/forward/point_mass.py
+++ b/harmonica/forward/point_mass.py
@@ -7,28 +7,117 @@ from numba import jit
 from ..constants import GRAVITATIONAL_CONST
 
 
-def point_mass_gravity(coordinates, points, masses, field, dtype="float64"):
-    """
-    Compute gravitational fields of a point mass on spherical coordinates.
+def point_mass_gravity(
+    coordinates, points, masses, field, coordinate_system="cartesian", dtype="float64"
+):
+    r"""
+    Compute gravitational fields of point masses.
+
+    It can compute the gravitational fields of point masses on a set of computation
+    points defined either in Cartesian or geocentric spherical coordinates.
+
+    The potential gravity field generated by a point mass with mass :math:`m` located at
+    a point :math:`Q` on a computation point :math:`P` can be computed as:
+
+    .. math::
+
+        V(P) = \frac{G m}{l},
+
+    where :math:`G` is the gravitational constant and :math:`l` is the Euclidean
+    distance between :math:`P` and :math:`Q` [Blakely1995]_.
+
+    In Cartesian coordinates, the points :math:`P` and :math:`Q` are given by :math:`x`,
+    :math:`y` and :math:`z` coordinates, which can be translated into ``northing``,
+    ``easting`` and ``down``, respectively.
+    If :math:`P` is located at :math:`(x, y, z)`, and :math:`Q` at :math:`(x_p, y_p,
+    z_p)`, the distance :math:`l` can be computed as:
+
+    .. math::
+
+        l = \sqrt{ (x - x_p)^2 + (y - y_p)^2 + (z - z_p)^2 }.
+
+    The gradient of the potential, also known as the gravity acceleration vector
+    :math:`\vec{g}`, is defined as:
+
+    .. math::
+
+        \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 downward):
+
+    .. math::
+
+        g_z(P) = \frac{G m}{l^3} (z_p - z).
+
+    On a geocentric spherical coordinate system, the points :math:`P` and :math:`Q` are
+    given by the ``longitude``, ``latitude`` and ``radius`` coordinates, i.e.
+    :math:`\lambda`, :math:`\varphi` and :math:`r`, respectively. On this coordinate
+    system, the Euclidean distance between :math:`P(r, \varphi, \lambda)` and
+    :math:`Q(r_p, \varphi_p, \lambda_p)` can be calculated  as follows [Grombein2013]_:
+
+    .. math::
+
+        l = \sqrt{ r^2 + r_p^2 - 2 r r_p \cos \Psi },
+
+    where
+
+    .. math::
+
+        \cos \Psi = \sin \varphi \sin \varphi_p +
+        \cos \varphi \cos \varphi_p \cos(\lambda - \lambda_p).
+
+    The radial component of the acceleration vector on a local North-oriented
+    system whose origin is located on the point :math:`P(r, \varphi, \lambda)`
+    is given by [Grombein2013]_:
+
+    .. math::
+
+        g_r(P) = \frac{G m}{l^3} (r_p \cos \Psi - r).
+
+    .. warning::
+
+        When working in Cartesian coordinates, the **z direction points downwards**,
+        i.e. positive and negative values represent points below and above the surface,
+        respectively.
+
 
     Parameters
     ----------
     coordinates : list or array
-        List or array containing `longitude`, `latitude` and `radius` of computation
-        points defined on a spherical geocentric coordinate system.
-        Both `longitude` and `latitude` should be in degrees and `radius` in meters.
+        List or array containing the coordinates of computation points in the following
+        order: ``easting``, ``northing`` and ``down`` (if coordinates given in
+        Cartesian coordiantes), or ``longitude``, ``latitude`` and ``radius`` (if given
+        on a spherical geocentric coordinate system).
+        All ``easting``, ``northing`` and ``down`` should be in meters.
+        Both ``longitude`` and ``latitude`` should be in degrees and ``radius`` in
+        meters.
     points : list or array
-        List or array containing the coordinates of the point masses `longitude_p`,
-        `latitude_p`, `radius_p` defined on a spherical geocentric coordinate system.
-        Both `longitude_p` and `latitude_p` should be in degrees and `radius` in meters.
+        List or array containing the coordinates of the point masses in the following
+        order: ``easting``, ``northing`` and ``down`` (if coordinates given in
+        Cartesian coordiantes), or ``longitude``, ``latitude`` and ``radius`` (if given
+        on a spherical geocentric coordinate system).
+        All ``easting``, ``northing`` and ``down`` should be in meters.
+        Both ``longitude`` and ``latitude`` should be in degrees and ``radius`` in
+        meters.
     masses : list or array
         List or array containing the mass of each point mass in kg.
     field : str
         Gravitational field that wants to be computed.
-        The available fields are:
+        The available fields in Cartesian coordinates are:
+
+        - Gravitational potential: ``potential``
+        - Downward acceleration: ``g_z``
+
+        The available fields in spherical geocentric coordinates are:
 
         - Gravitational potential: ``potential``
         - Radial acceleration: ``g_r``
+
+    coordinate_system : str (optional)
+        Coordinate system of the coordinates of the computation points and the point
+        masses. Available coordinates systems: ``cartesian``, ``spherical``.
+        Default ``cartesian``.
     dtype : data-type (optional)
         Data type assigned to resulting gravitational field, and coordinates of point
         masses and computation points. Default to ``np.float64``.
@@ -37,50 +126,122 @@ def point_mass_gravity(coordinates, points, masses, field, dtype="float64"):
     Returns
     -------
     result : array
-        Gravitational field generated by the `point_mass` on the computation points
-        defined in `coordinates`.
+        Gravitational field generated by the ``point_mass`` on the computation points
+        defined in ``coordinates``.
         The potential is given in SI units, the accelerations in mGal and the Marussi
         tensor components in Eotvos.
     """
-    kernels = {"potential": kernel_potential, "g_r": kernel_g_r}
-    if field not in kernels:
+    # Organize dispatchers and kernel functions inside dictionaries
+    dispatchers = {
+        "cartesian": jit_point_mass_cartesian,
+        "spherical": jit_point_mass_spherical,
+    }
+    kernels = {
+        "cartesian": {"potential": kernel_potential_cartesian, "g_z": kernel_g_z},
+        "spherical": {"potential": kernel_potential_spherical, "g_r": kernel_g_r},
+    }
+    # Sanity checks for coordinate_system and field
+    if coordinate_system not in ("cartesian", "spherical"):
+        raise ValueError(
+            "Coordinate system {} not recognized".format(coordinate_system)
+        )
+    if field not in kernels[coordinate_system]:
         raise ValueError("Gravity field {} not recognized".format(field))
     # Figure out the shape and size of the output array
     cast = np.broadcast(*coordinates[:3])
     result = np.zeros(cast.size, dtype=dtype)
     # Prepare arrays to be passed to the jitted functions
-    longitude, latitude, radius = (
-        np.atleast_1d(i).ravel().astype(dtype) for i in coordinates[:3]
-    )
-    longitude_p, latitude_p, radius_p = (
-        np.atleast_1d(i).ravel().astype(dtype) for i in points[:3]
-    )
+    coordinates = (np.atleast_1d(i).ravel().astype(dtype) for i in coordinates[:3])
+    points = (np.atleast_1d(i).ravel().astype(dtype) for i in points[:3])
     masses = np.atleast_1d(masses).astype(dtype).ravel()
     # Compute gravitational field
-    jit_point_mass_gravity(
-        longitude,
-        latitude,
-        radius,
-        longitude_p,
-        latitude_p,
-        radius_p,
-        masses,
-        result,
-        kernels[field],
+    dispatchers[coordinate_system](
+        *coordinates, *points, masses, result, kernels[coordinate_system][field]
     )
     result *= GRAVITATIONAL_CONST
     # Convert to more convenient units
-    if field == "g_r":
+    if field in ("g_r", "g_z"):
         result *= 1e5  # SI to mGal
     return result.reshape(cast.shape)
 
 
 @jit(nopython=True)
-def jit_point_mass_gravity(
+def jit_point_mass_cartesian(
+    easting, northing, down, easting_p, northing_p, down_p, masses, out, kernel
+):  # pylint: disable=invalid-name
+    """
+    Compute gravity field of point masses on computation points in Cartesian coordinates
+
+    Parameters
+    ----------
+    easting, northing, down : 1d-arrays
+        Coordinates of computation points in Cartesian coordinate system.
+    easting_p, northing_p, down_p : 1d-arrays
+        Coordinates of point masses in Cartesian coordinate system.
+    masses : 1d-array
+        Mass of each point mass in SI units.
+    out : 1d-array
+        Array where the gravitational field on each computation point will be appended.
+        It must have the same size of ``easting``, ``northing`` and ``down``.
+    kernel : func
+        Kernel function that will be used to compute the gravity field on the
+        computation points.
+    """
+    for l in range(easting.size):
+        for m in range(easting_p.size):
+            out[l] += masses[m] * kernel(
+                easting[l], northing[l], down[l], easting_p[m], northing_p[m], down_p[m]
+            )
+
+
+@jit(nopython=True)
+def kernel_potential_cartesian(easting, northing, down, easting_p, northing_p, down_p):
+    """
+    Kernel function for potential gravity field in Cartesian coordinates
+    """
+    return 1 / _distance_cartesian(
+        [easting, northing, down], [easting_p, northing_p, down_p]
+    )
+
+
+@jit(nopython=True)
+def kernel_g_z(easting, northing, down, easting_p, northing_p, down_p):
+    """
+    Kernel function for downward component of gravity gradient in Cartesian coordinates
+    """
+    distance_sq = _distance_cartesian_sq(
+        [easting, northing, down], [easting_p, northing_p, down_p]
+    )
+    return (down_p - down) / distance_sq ** (3 / 2)
+
+
+@jit(nopython=True)
+def _distance_cartesian_sq(point_a, point_b):
+    """
+    Calculate the square distance between two points given in Cartesian coordinates
+    """
+    easting, northing, down = point_a[:]
+    easting_p, northing_p, down_p = point_b[:]
+    distance_sq = (
+        (easting - easting_p) ** 2 + (northing - northing_p) ** 2 + (down - down_p) ** 2
+    )
+    return distance_sq
+
+
+@jit(nopython=True)
+def _distance_cartesian(point_a, point_b):
+    """
+    Calculate the distance between two points given in Cartesian coordinates
+    """
+    return np.sqrt(_distance_cartesian_sq(point_a, point_b))
+
+
+@jit(nopython=True)
+def jit_point_mass_spherical(
     longitude, latitude, radius, longitude_p, latitude_p, radius_p, masses, out, kernel
 ):  # pylint: disable=invalid-name
     """
-    Compute gravity field of point masses on computation points.
+    Compute gravity field of point masses on computation points in spherical coordiantes
 
     Parameters
     ----------
@@ -123,11 +284,11 @@ def jit_point_mass_gravity(
 
 
 @jit(nopython=True)
-def kernel_potential(
+def kernel_potential_spherical(
     longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
 ):
     """
-    Kernel function for potential gravity field
+    Kernel function for potential gravity field in spherical coordinates
     """
     coslambda = np.cos(longitude_p - longitude)
     cospsi = sinphi_p * sinphi + cosphi_p * cosphi * coslambda
@@ -140,7 +301,7 @@ def kernel_g_r(
     longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
 ):
     """
-    Kernel function for radial component of gravity gradient
+    Kernel function for radial component of gravity gradient in spherical coordinates
     """
     coslambda = np.cos(longitude_p - longitude)
     cospsi = sinphi_p * sinphi + cosphi_p * cosphi * coslambda
diff --git a/harmonica/forward/tesseroid.py b/harmonica/forward/tesseroid.py
index 021d51f..d3d679e 100644
--- a/harmonica/forward/tesseroid.py
+++ b/harmonica/forward/tesseroid.py
@@ -6,7 +6,7 @@ import numpy as np
 from numpy.polynomial.legendre import leggauss
 
 from ..constants import GRAVITATIONAL_CONST
-from .point_mass import jit_point_mass_gravity, kernel_potential, kernel_g_r
+from .point_mass import jit_point_mass_spherical, kernel_potential_spherical, kernel_g_r
 
 STACK_SIZE = 100
 MAX_DISCRETIZATIONS = 100000
@@ -110,7 +110,7 @@ def tesseroid_gravity(
     array(-112.54539933)
 
     """
-    kernels = {"potential": kernel_potential, "g_r": kernel_g_r}
+    kernels = {"potential": kernel_potential_spherical, "g_r": kernel_g_r}
     if field not in kernels:
         raise ValueError("Gravity field {} not recognized".format(field))
     # Figure out the shape and size of the output array
@@ -258,7 +258,7 @@ def jit_tesseroid_gravity(
                 weights,
             )
             # Compute gravity fields
-            jit_point_mass_gravity(
+            jit_point_mass_spherical(
                 computation_point[0:1],
                 computation_point[1:2],
                 computation_point[2:3],


Test patch:
diff --git a/harmonica/tests/test_point_mass.py b/harmonica/tests/test_point_mass.py
index 3aba539..dfa1dd8 100644
--- a/harmonica/tests/test_point_mass.py
+++ b/harmonica/tests/test_point_mass.py
@@ -9,19 +9,106 @@ from ..constants import GRAVITATIONAL_CONST
 from ..forward.point_mass import point_mass_gravity
 
 
-def test_invalid_field():
-    "Check if an invalid gravitational field is passed as argument"
+def test_invalid_coordinate_system():
+    "Check if invalid coordinate system is passed"
+    coordinates = [0.0, 0.0, 0.0]
     point_mass = [0.0, 0.0, 0.0]
     mass = 1.0
-    longitude = np.array(0.0)
-    latitude = np.array(0.0)
-    height = np.array(0.0)
     with pytest.raises(ValueError):
         point_mass_gravity(
-            [longitude, latitude, height], point_mass, mass, "this-field-does-not-exist"
+            coordinates,
+            point_mass,
+            mass,
+            "potential",
+            "this-is-not-a-valid-coordinate-system",
         )
 
 
+def test_invalid_field():
+    "Check if an invalid gravitational field is passed as argument"
+    coordinates = [0.0, 0.0, 0.0]
+    point_mass = [0.0, 0.0, 0.0]
+    mass = 1.0
+    for coordinate_system in ("spherical", "cartesian"):
+        with pytest.raises(ValueError):
+            point_mass_gravity(
+                coordinates,
+                point_mass,
+                mass,
+                "this-field-does-not-exist",
+                coordinate_system,
+            )
+
+
+def test_invalid_field_for_coordinate_system():
+    "Check if an invalid field is passed for that coordinate system"
+    coordinates = [0.0, 0.0, 0.0]
+    point_mass = [0.0, 0.0, 0.0]
+    mass = 1.0
+    cartesian_exclusive_fields = ["g_z"]
+    spherical_exclusive_fields = ["g_r"]
+    for field in spherical_exclusive_fields:
+        with pytest.raises(ValueError):
+            point_mass_gravity(coordinates, point_mass, mass, field, "cartesian")
+    for field in cartesian_exclusive_fields:
+        with pytest.raises(ValueError):
+            point_mass_gravity(coordinates, point_mass, mass, field, "spherical")
+
+
+# ---------------------------
+# Cartesian coordinates tests
+# ---------------------------
+@pytest.mark.use_numba
+def test_potential_cartesian_symmetry():
+    """
+    Test if potential field of a point mass has symmetry in Cartesian coordinates
+    """
+    # Define a single point mass
+    point_mass = [1.1, 1.2, 1.3]
+    masses = [2670]
+    # Define a set of computation points at a fixed distance from the point mass
+    distance = 3.3
+    easting = point_mass[0] * np.ones(6)
+    northing = point_mass[1] * np.ones(6)
+    down = point_mass[2] * np.ones(6)
+    easting[0] += distance
+    easting[1] -= distance
+    northing[2] += distance
+    northing[3] -= distance
+    down[4] += distance
+    down[5] -= distance
+    coordinates = [easting, northing, down]
+    # Compute potential gravity field on each computation point
+    results = point_mass_gravity(
+        coordinates, point_mass, masses, "potential", "cartesian"
+    )
+    npt.assert_allclose(*results)
+
+
+@pytest.mark.use_numba
+def test_g_z_symmetry():
+    """
+    Test if g_z field of a point mass has symmetry in 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 bellow the point mass
+    distance = 3.3
+    easting = point_mass[0] * np.ones(2)
+    northing = point_mass[1] * np.ones(2)
+    down = point_mass[2] * np.ones(2)
+    down[0] += distance
+    down[1] -= distance
+    coordinates = [easting, northing, down]
+    # Compute g_z gravity field on each computation point
+    results = point_mass_gravity(coordinates, point_mass, masses, "g_z", "cartesian")
+    npt.assert_allclose(results[0], -results[1])
+
+
+# ---------------------------
+# Spherical coordinates tests
+# ---------------------------
 @pytest.mark.use_numba
 def test_point_mass_on_origin():
     "Check potential and g_r of point mass on origin"
@@ -39,7 +126,9 @@ def test_point_mass_on_origin():
     # Compare results with analytical solutions
     for field in analytical:
         npt.assert_allclose(
-            point_mass_gravity([longitude, latitude, radius], point_mass, mass, field),
+            point_mass_gravity(
+                [longitude, latitude, radius], point_mass, mass, field, "spherical"
+            ),
             analytical[field],
         )
 
@@ -67,7 +156,9 @@ def test_point_mass_same_radial_direction():
                 # Compare results with analytical solutions
                 for field in analytical:
                     npt.assert_allclose(
-                        point_mass_gravity(coordinates, point_mass, mass, field),
+                        point_mass_gravity(
+                            coordinates, point_mass, mass, field, "spherical"
+                        ),
                         analytical[field],
                     )
 
@@ -95,7 +186,9 @@ def test_point_mass_potential_on_equator():
                 analytical = {"potential": GRAVITATIONAL_CONST * mass / distance}
                 # Compare results with analytical solutions
                 npt.assert_allclose(
-                    point_mass_gravity(coordinates, point_mass, mass, "potential"),
+                    point_mass_gravity(
+                        coordinates, point_mass, mass, "potential", "spherical"
+                    ),
                     analytical["potential"],
                 )
 
@@ -123,6 +216,8 @@ def test_point_mass_potential_on_same_meridian():
                 analytical = {"potential": GRAVITATIONAL_CONST * mass / distance}
                 # Compare results with analytical solutions
                 npt.assert_allclose(
-                    point_mass_gravity(coordinates, point_mass, mass, "potential"),
+                    point_mass_gravity(
+                        coordinates, point_mass, mass, "potential", "spherical"
+                    ),
                     analytical["potential"],
                 )