Repository: fatiando/choclo (language: python)

Issue:
Clarify what is meant by "kernel"
<!-- Please describe your suggestion/problem/question in as much detail as possible. -->

It would be great to clarify this in the Project Goals. For example, I noticed that the functions don't include universal constants or physical properties. For gravity, this is ok since getting the gravity field is a matter of simple multiplication. For magnetics, the combination is a bit more difficult since it involves a dot product with the vector magnetisation. So it would be worth having that in a function already instead of requiring users to code it. It also may be necessary since doing it out of the jit compiled function would mean losing parallel and low memory execution. 

It would be worth clarifying the scope of the package with this in mind. A good way to go about it would be to try to add a magnetic dipole implementation to see how it fits with the current paradigm. 

Or maybe the "kernel" could mean the version without the for loops? It may be worth making that available in the API for use in external jit-compiled code (like tesseroids). 


Discussion / hints:
I agree that the use of "kernel" is a too loose at the moment. I still struggle to give a closed definition to it. I've been using this name for grouping a set of low-level functions used in the forward models for gravity and magnetics.

The name "kernel" is inherited from the integral kernels involved in the forward modelling functions of potential fields.
For example, for any three dimensional body $\Gamma$, a potential field $\phi$ on an observation point $\mathbf{p}$ can be computed as:

$$ \phi(\mathbf{p}) = \int\limits_\Gamma g(\mathbf{p}) k(\mathbf{p}, \mathbf{q}) \text{d}\mathbf{q} $$

where $k(\mathbf{p}, \mathbf{q})$ is the integral kernel and $\mathbf{g(\mathbf{p})}$ is a function that doesn't depend on the integration variable and it's usually constant for gravity and magnetic forward modelling.

For the gravity potential it would be:

$$ 
V_\text{grav}(\mathbf{p}) = 
    G 
    \int\limits_\Gamma 
    \frac{
        \rho(\mathbf{q})
    }{
        \lVert \mathbf{p} - \mathbf{q} \rVert
    }
\text{d}\mathbf{q}
$$

And for magnetic potential (in absence of currents):

$$ 
V_\text{mag}(\mathbf{p}) = 
    \frac{\mu_0}{4\pi}
    \int\limits_\Gamma 
    \mathbf{M(\mathbf{q})}
    \cdot
    \nabla_\mathbf{q} \left( \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} \right)
\text{d}\mathbf{q}
$$

If we want to continue with this idea of a "kernel" function, then the physical properties (density and magnetization) must be included in them.

> It would be worth clarifying the scope of the package with this in mind. A good way to go about it would be to try to add a magnetic dipole implementation to see how it fits with the current paradigm.

I agree. In my head the implementation of the magnetic dipole would ask the coordinates of the observation point, the coordinates of the dipole and its magnetization vector.

> Or maybe the "kernel" could mean the version without the for loops? It may be worth making that available in the API for use in external jit-compiled code (like tesseroids).

For now I'm only adding functions without the for loops, so: single computation point and single source. There is a broad interest in having those functions without the for loops.

---

Nevertheless, there are some caveats. We perform the prism forward model by evaluating the analytical solution to these integrals, so we are not actually computing the integral kernel. Should we use a different name for that function then? Which one? Could we still include it in the "kernel" group?

After a discussion we had during a meeting this week we decided to stick with the following design for Choclo.

# Everything is JIT compiled

We could offer non-jit compiled versions of our functions, although there are some caveats. Some functions are written to support floats (like the `safe_log` and `safe_atan` ones). If the would have to support comparisons for arrays, we would need to introduce some checks there. We want these functions to run quick and dirty, no extra steps.

Besides, the forward modelling functions we will offer in Choclo (single source, single computation point) are meant to be run using Numba for loops, so it makes sense for everything to be JIT compiled.

# No external loops

Choclo will not offer (at least for now) functions that compute forward models over sets of computations points and/or sets of sources. It will only offer forward modelling functions for **single computation point** and **single source**.
The reason behind this is that the outer for loops can be optimized and written in many different combinations depending on the specific problem the user want to tackle. It's easier to leave that part to them rather than offering a collection of different functions that perform almost the same, but through different ways. Besides, with the functions that Choclo will offer, these outer loops only cost two to three extra lines, e.g.:

## Full forward model

```python
import choclo as ch

@numba.jit(parallel=True)
def gravity_potential(...):
    result = np.zeros(N)
    for point in prange(N):
        for source in range(M):
            result[point] += ch.point.gravity_pot(...)
    return result
```

## Build Jacobian matrix

```python
import choclo as ch

@numba.jit(parallel=True)
def jacobian(...):
    A = np.empty(shape)
    for i in prange(shape[0]):
        for j in range(shape[1]):
            A[i, j] = ch.point.gravity_pot(..., density=1)
```

# Architecture

Choclo will be split in submodules, one for each source type plus a `utils` one (will include the distance functions).

Each submodule will have a two level system of functions.
1. **Level 1:** Low-level kernel functions that only perform the math. They don't call any other function within Choclo. For example, for point sources, they compute `1/r` and its derivatives. They will use the following name convention: `kernel_pot` (for the kernel function used in the potential field), `kernel_e`, `kernel_n`, `kernel_u` (for the kernels used in the gradient of the potential) and `kernel_ee`, `kernel_nn`, `kernel_uu`, `kernel_en`, `kernel_ez`, `kernel_nz` (for the kernels used in the tensor components and magnetic field).
2. **Level 2:** Higher-level functions that perform the forward modelling of a single source on a single computation point. These will depend on the Level 1 functions and on the utils functions. These functions will take the physical property of the source (density, magnetization vector) and will make use of universal constants. They will return every result in SI units. Their naming conventions will be:
```
gravity_pot(coordinates, source, density)
gravity_e
gravity_n
gravity_u
gravity_ee
gravity_en
gravity_eu
gravity_nn
gravity_nu
gravity_uu
magnetic_e(coordinates, source, magnetization_vector)
magnetic_n
magnetic_u
magnetic_tf_anomaly(coordinates, source, magnetization_vector, igrf_inc, igrf_dec)
```

# Docs

The docs will include instructions on how to use Choclo functions properly and in an optimized way, showing examples on how to build a full forward model and a Jacobian matrix.

Fix patch:
diff --git a/choclo/__init__.py b/choclo/__init__.py
index a9734453..fb21a163 100644
--- a/choclo/__init__.py
+++ b/choclo/__init__.py
@@ -5,17 +5,4 @@
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #
 # Import functions/classes to make the public API
-from ._distance import distance_cartesian, distance_spherical, distance_spherical_core
-from ._point import (
-    kernel_point_g_easting,
-    kernel_point_g_ee,
-    kernel_point_g_en,
-    kernel_point_g_ez,
-    kernel_point_g_nn,
-    kernel_point_g_northing,
-    kernel_point_g_nz,
-    kernel_point_g_upward,
-    kernel_point_g_zz,
-    kernel_point_potential,
-)
 from ._version import __version__
diff --git a/choclo/_point.py b/choclo/_point.py
deleted file mode 100644
index 78677347..00000000
--- a/choclo/_point.py
+++ /dev/null
@@ -1,654 +0,0 @@
-# Copyright (c) 2022 The Choclo Developers.
-# Distributed under the terms of the BSD 3-Clause License.
-# SPDX-License-Identifier: BSD-3-Clause
-#
-# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
-#
-"""
-Kernel functions for point sources
-"""
-from numba import jit
-
-from ._distance import distance_cartesian
-
-
-@jit(nopython=True)
-def kernel_point_potential(
-    easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
-):
-    r"""
-    Kernel for the potential field due to a point source
-
-    .. important ::
-
-        The observation point and the location of the point source must be in
-        Cartesian coordinates and have the same units.
-
-    Parameters
-    ----------
-    easting_p : float
-        Easting coordinate of the observation point.
-    northing_p : float
-        Northing coordinate of the observation point.
-    upward_p : float
-        Upward coordinate of the observation point.
-    easting_q : float
-        Easting coordinate of the point source location.
-    northing_q : float
-        Northing coordinate of the point source location.
-    upward_q : float
-        Upward coordinate of the point source location.
-
-    Returns
-    -------
-    kernel : float
-        Value of the kernel function for the potential field of a point source.
-
-    Notes
-    -----
-    Given an observation point located in :math:`\mathbf{p} = (x_p, y_p, z_p)`
-    and a point source located in :math:`\mathbf{q} = (x_q, y_q, z_q)` defined
-    in a Cartesian coordinate system, compute the kernel function for the
-    potential field that the point source generates on the observation point:
-
-    .. math::
-
-        k_V(\mathbf{p}, \mathbf{q}) =
-        \frac{
-            1
-        }{
-            \lVert \mathbf{p} - \mathbf{q} \rVert_2
-        }
-
-    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
-    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`).
-    """
-    distance = distance_cartesian(
-        (easting_p, northing_p, upward_p), (easting_q, northing_q, upward_q)
-    )
-    return 1 / distance
-
-
-@jit(nopython=True)
-def kernel_point_g_easting(
-    easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
-):
-    r"""
-    Kernel for easting component of the gradient due to a point source
-
-    .. important ::
-
-        The observation point and the location of the point source must be in
-        Cartesian coordinates and have the same units.
-
-    Parameters
-    ----------
-    easting_p : float
-        Easting coordinate of the observation point.
-    northing_p : float
-        Northing coordinate of the observation point.
-    upward_p : float
-        Upward coordinate of the observation point.
-    easting_q : float
-        Easting coordinate of the point source location.
-    northing_q : float
-        Northing coordinate of the point source location.
-    upward_q : float
-        Upward coordinate of the point source location.
-
-    Returns
-    -------
-    kernel : float
-        Value of the kernel function for the easting component of the
-        gradient of the potential field of a point source.
-
-    Notes
-    -----
-    Given an observation point located in :math:`\mathbf{p} = (x_p, y_p, z_p)`
-    and a point source located in :math:`\mathbf{q} = (x_q, y_q, z_q)` defined
-    in a Cartesian coordinate system, compute the kernel function for the
-    easting component of the gradient of the potential field that the point
-    source generates on the observation point:
-
-    .. math::
-
-        k_{g_x}(\mathbf{p}, \mathbf{q}) =
-        - \frac{
-            x_p - x_q
-        }{
-            \lVert \mathbf{p} - \mathbf{q} \rVert_2
-        }
-
-    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
-    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`).
-    """
-    distance = distance_cartesian(
-        (easting_p, northing_p, upward_p), (easting_q, northing_q, upward_q)
-    )
-    return -(easting_p - easting_q) / distance**3
-
-
-@jit(nopython=True)
-def kernel_point_g_northing(
-    easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
-):
-    r"""
-    Kernel for northing component of the gradient due to a point source
-
-    .. important ::
-
-        The observation point and the location of the point source must be in
-        Cartesian coordinates
-
-    Parameters
-    ----------
-    easting_p : float
-        Easting coordinate of the observation point.
-    northing_p : float
-        Northing coordinate of the observation point.
-    upward_p : float
-        Upward coordinate of the observation point.
-    easting_q : float
-        Easting coordinate of the point source location.
-    northing_q : float
-        Northing coordinate of the point source location.
-    upward_q : float
-        Upward coordinate of the point source location.
-
-    Returns
-    -------
-    kernel : float
-        Value of the kernel function for the northing component of the
-        gradient of the potential field of a point source.
-
-    Notes
-    -----
-    Given an observation point located in :math:`\mathbf{p} = (x_p, y_p, z_p)`
-    and a point source located in :math:`\mathbf{q} = (x_q, y_q, z_q)` defined
-    in a Cartesian coordinate system, compute the kernel function for the
-    northing component of the gradient of the potential field that the point
-    source generates on the observation point:
-
-    .. math::
-
-        k_{g_y}(\mathbf{p}, \mathbf{q}) =
-        - \frac{
-            y_p - y_q
-        }{
-            \lVert \mathbf{p} - \mathbf{q} \rVert_2
-        }
-
-    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
-    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`).
-    """
-    distance = distance_cartesian(
-        (easting_p, northing_p, upward_p), (easting_q, northing_q, upward_q)
-    )
-    return -(northing_p - northing_q) / distance**3
-
-
-@jit(nopython=True)
-def kernel_point_g_upward(
-    easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
-):
-    r"""
-    Kernel for upward component of the gradient due to a point source
-
-    .. important ::
-
-        The observation point and the location of the point source must be in
-        Cartesian coordinates and have the same units.
-
-    .. warning ::
-
-        This kernel corresponds to the **upward** component of the
-        gradient of the potential field. In need of the _downward_ component,
-        multiply this result by -1.
-
-    Parameters
-    ----------
-    easting_p : float
-        Easting coordinate of the observation point.
-    northing_p : float
-        Northing coordinate of the observation point.
-    upward_p : float
-        Upward coordinate of the observation point.
-    easting_q : float
-        Easting coordinate of the point source location.
-    northing_q : float
-        Northing coordinate of the point source location.
-    upward_q : float
-        Upward coordinate of the point source location.
-
-    Returns
-    -------
-    kernel : float
-        Value of the kernel function for the upward component of the gradient
-        of the potential field of a point source.
-
-    Notes
-    -----
-    Given an observation point located in :math:`\mathbf{p} = (x_p, y_p, z_p)`
-    and a point source located in :math:`\mathbf{q} = (x_q, y_q, z_q)` defined
-    in a Cartesian coordinate system, compute the kernel function for the
-    upward component of the gradient of the potential field that the point
-    source generates on the observation point:
-
-    .. math::
-
-        k_{g_z}(\mathbf{p}, \mathbf{q}) =
-        - \frac{
-            z_p - z_q
-        }{
-            \lVert \mathbf{p} - \mathbf{q} \rVert_2
-        }
-
-    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
-    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`).
-    """
-    distance = distance_cartesian(
-        (easting_p, northing_p, upward_p), (easting_q, northing_q, upward_q)
-    )
-    return -(upward_p - upward_q) / distance**3
-
-
-@jit(nopython=True)
-def kernel_point_g_ee(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q):
-    r"""
-    Kernel for the :math:`G_\text{ee}` tensor component due to a point source
-
-    .. important ::
-
-        The observation point and the location of the point source must be in
-        Cartesian coordinates and have the same units.
-
-    Parameters
-    ----------
-    easting_p : float
-        Easting coordinate of the observation point.
-    northing_p : float
-        Northing coordinate of the observation point.
-    upward_p : float
-        Upward coordinate of the observation point.
-    easting_q : float
-        Easting coordinate of the point source location.
-    northing_q : float
-        Northing coordinate of the point source location.
-    upward_q : float
-        Upward coordinate of the point source location.
-
-    Returns
-    -------
-    kernel : float
-        Value of the kernel function for the ``G_ee`` component of the
-        potential field tensor due to a point source.
-
-    Notes
-    -----
-    Given an observation point located in :math:`\mathbf{p} = (x_p, y_p, z_p)`
-    and a point source located in :math:`\mathbf{q} = (x_q, y_q, z_q)` defined
-    in a Cartesian coordinate system, compute the kernel function for the
-    diagonal :math:`G_\text{ee}` component of the potential field tensor that
-    the point source generates on the observation point:
-
-    .. math::
-
-        k_{G_\text{ee}}(\mathbf{p}, \mathbf{q}) =
-        \frac{
-            3 (x_p - x_q)^2
-        }{
-            \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
-        }
-        - \frac{
-            1
-        }{
-            \lVert \mathbf{p} - \mathbf{q} \rVert_2^3
-        }
-
-    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
-    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`) and
-    :math:`G_\text{ee}` is defined as:
-
-    .. math::
-
-        G_\text{ee}(\mathbf{p}) = \frac{\partial^2 V(\mathbf{p})}{\partial x^2}
-    """
-    distance = distance_cartesian(
-        (easting_p, northing_p, upward_p), (easting_q, northing_q, upward_q)
-    )
-    return 3 * (easting_p - easting_q) ** 2 / distance**5 - 1 / distance**3
-
-
-@jit(nopython=True)
-def kernel_point_g_nn(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q):
-    r"""
-    Kernel for the :math:`G_\text{nn}` tensor component due to a point source
-
-    .. important ::
-
-        The observation point and the location of the point source must be in
-        Cartesian coordinates and have the same units.
-
-    Parameters
-    ----------
-    easting_p : float
-        Easting coordinate of the observation point.
-    northing_p : float
-        Northing coordinate of the observation point.
-    upward_p : float
-        Upward coordinate of the observation point.
-    easting_q : float
-        Easting coordinate of the point source location.
-    northing_q : float
-        Northing coordinate of the point source location.
-    upward_q : float
-        Upward coordinate of the point source location.
-
-    Returns
-    -------
-    kernel : float
-        Value of the kernel function for the ``G_nn`` component of the
-        potential field tensor due to a point source.
-
-    Notes
-    -----
-    Given an observation point located in :math:`\mathbf{p} = (x_p, y_p, z_p)`
-    and a point source located in :math:`\mathbf{q} = (x_q, y_q, z_q)` defined
-    in a Cartesian coordinate system, compute the kernel function for the
-    diagonal :math:`G_\text{nn}` component of the potential field tensor that
-    the point source generates on the observation point:
-
-    .. math::
-
-        k_{G_\text{nn}}(\mathbf{p}, \mathbf{q}) =
-        \frac{
-            3 (y_p - y_q)^2
-        }{
-            \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
-        }
-        - \frac{
-            1
-        }{
-            \lVert \mathbf{p} - \mathbf{q} \rVert_2^3
-        }
-
-    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
-    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`) and
-    :math:`G_\text{nn}` is defined as:
-
-    .. math::
-
-        G_\text{ee}(\mathbf{p}) = \frac{\partial^2 V(\mathbf{p})}{\partial x^2}
-    """
-    distance = distance_cartesian(
-        (easting_p, northing_p, upward_p), (easting_q, northing_q, upward_q)
-    )
-    return 3 * (northing_p - northing_q) ** 2 / distance**5 - 1 / distance**3
-
-
-@jit(nopython=True)
-def kernel_point_g_zz(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q):
-    r"""
-    Kernel for the :math:`G_\text{zz}` tensor component due to a point source
-
-    .. important ::
-
-        The observation point and the location of the point source must be in
-        Cartesian coordinates and have the same units.
-
-    Parameters
-    ----------
-    easting_p : float
-        Easting coordinate of the observation point.
-    northing_p : float
-        Northing coordinate of the observation point.
-    upward_p : float
-        Upward coordinate of the observation point.
-    easting_q : float
-        Easting coordinate of the point source location.
-    northing_q : float
-        Northing coordinate of the point source location.
-    upward_q : float
-        Upward coordinate of the point source location.
-
-    Returns
-    -------
-    kernel : float
-        Value of the kernel function for the ``G_zz`` component of the
-        potential field tensor due to a point source.
-
-    Notes
-    -----
-    Given an observation point located in :math:`\mathbf{p} = (x_p, y_p, z_p)`
-    and a point source located in :math:`\mathbf{q} = (x_q, y_q, z_q)` defined
-    in a Cartesian coordinate system, compute the kernel function for the
-    diagonal :math:`G_\text{zz}` component of the potential field tensor that
-    the point source generates on the observation point:
-
-    .. math::
-
-        k_{G_\text{zz}}(\mathbf{p}, \mathbf{q}) =
-        \frac{
-            3 (z_p - z_q)^2
-        }{
-            \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
-        }
-        - \frac{
-            1
-        }{
-            \lVert \mathbf{p} - \mathbf{q} \rVert_2^3
-        }
-
-    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
-    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`) and
-    :math:`G_\text{zz}` is defined as:
-
-    .. math::
-
-        G_\text{zz}(\mathbf{p}) = \frac{\partial^2 V(\mathbf{p})}{\partial z^2}
-    """
-    distance = distance_cartesian(
-        (easting_p, northing_p, upward_p), (easting_q, northing_q, upward_q)
-    )
-    return 3 * (upward_p - upward_q) ** 2 / distance**5 - 1 / distance**3
-
-
-@jit(nopython=True)
-def kernel_point_g_en(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q):
-    r"""
-    Kernel for the :math:`G_\text{en}` tensor component due to a point source
-
-    .. important ::
-
-        The observation point and the location of the point source must be in
-        Cartesian coordinates and have the same units.
-
-    Parameters
-    ----------
-    easting_p : float
-        Easting coordinate of the observation point.
-    northing_p : float
-        Northing coordinate of the observation point.
-    upward_p : float
-        Upward coordinate of the observation point.
-    easting_q : float
-        Easting coordinate of the point source location.
-    northing_q : float
-        Northing coordinate of the point source location.
-    upward_q : float
-        Upward coordinate of the point source location.
-
-    Returns
-    -------
-    kernel : float
-        Value of the kernel function for the ``G_en`` component of the
-        potential field tensor due to a point source.
-
-    Notes
-    -----
-    Given an observation point located in :math:`\mathbf{p} = (x_p, y_p, z_p)`
-    and a point source located in :math:`\mathbf{q} = (x_q, y_q, z_q)` defined
-    in a Cartesian coordinate system, compute the kernel function for the
-    non-diagonal :math:`G_\text{en}` component of the potential field tensor
-    that the point source generates on the observation point:
-
-    .. math::
-
-        k_{G_\text{en}}(\mathbf{p}, \mathbf{q}) =
-        \frac{
-            3 (x_p - x_q) (y_p - y_q)
-        }{
-            \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
-        }
-
-    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
-    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`) and
-    :math:`G_\text{en}` is defined as:
-
-    .. math::
-
-        G_\text{en}(\mathbf{p}) =
-        \frac{
-            \partial^2 V(\mathbf{p})
-        }{
-            \partial x \partial y
-        }
-    """
-    distance = distance_cartesian(
-        (easting_p, northing_p, upward_p), (easting_q, northing_q, upward_q)
-    )
-    return 3 * (easting_p - easting_q) * (northing_p - northing_q) / distance**5
-
-
-@jit(nopython=True)
-def kernel_point_g_ez(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q):
-    r"""
-    Kernel for the :math:`G_\text{ez}` tensor component due to a point source
-
-    .. important ::
-
-        The observation point and the location of the point source must be in
-        Cartesian coordinates and have the same units.
-
-    Parameters
-    ----------
-    easting_p : float
-        Easting coordinate of the observation point.
-    northing_p : float
-        Northing coordinate of the observation point.
-    upward_p : float
-        Upward coordinate of the observation point.
-    easting_q : float
-        Easting coordinate of the point source location.
-    northing_q : float
-        Northing coordinate of the point source location.
-    upward_q : float
-        Upward coordinate of the point source location.
-
-    Returns
-    -------
-    kernel : float
-        Value of the kernel function for the ``G_ez`` component of the
-        potential field tensor due to a point source.
-
-    Notes
-    -----
-    Given an observation point located in :math:`\mathbf{p} = (x_p, y_p, z_p)`
-    and a point source located in :math:`\mathbf{q} = (x_q, y_q, z_q)` defined
-    in a Cartesian coordinate system, compute the kernel function for the
-    non-diagonal :math:`G_\text{ez}` component of the potential field tensor
-    that the point source generates on the observation point:
-
-    .. math::
-
-        k_{G_\text{ez}}(\mathbf{p}, \mathbf{q}) =
-        \frac{
-            3 (x_p - x_q) (z_p - z_q)
-        }{
-            \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
-        }
-
-    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
-    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`) and
-    :math:`G_\text{en}` is defined as:
-
-    .. math::
-
-        G_\text{ez}(\mathbf{p}) =
-        \frac{
-            \partial^2 V(\mathbf{p})
-        }{
-            \partial x \partial z
-        }
-    """
-    distance = distance_cartesian(
-        (easting_p, northing_p, upward_p), (easting_q, northing_q, upward_q)
-    )
-    return 3 * (easting_p - easting_q) * (upward_p - upward_q) / distance**5
-
-
-@jit(nopython=True)
-def kernel_point_g_nz(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q):
-    r"""
-    Kernel for the :math:`G_\text{nz}` tensor component due to a point source
-
-    .. important ::
-
-        The observation point and the location of the point source must be in
-        Cartesian coordinates and have the same units.
-
-    Parameters
-    ----------
-    easting_p : float
-        Easting coordinate of the observation point.
-    northing_p : float
-        Northing coordinate of the observation point.
-    upward_p : float
-        Upward coordinate of the observation point.
-    easting_q : float
-        Easting coordinate of the point source location.
-    northing_q : float
-        Northing coordinate of the point source location.
-    upward_q : float
-        Upward coordinate of the point source location.
-
-    Returns
-    -------
-    kernel : float
-        Value of the kernel function for the ``G_nz`` component of the
-        potential field tensor due to a point source.
-
-    Notes
-    -----
-    Given an observation point located in :math:`\mathbf{p} = (x_p, y_p, z_p)`
-    and a point source located in :math:`\mathbf{q} = (x_q, y_q, z_q)` defined
-    in a Cartesian coordinate system, compute the kernel function for the
-    non-diagonal :math:`G_\text{nz}` component of the potential field tensor
-    that the point source generates on the observation point:
-
-    .. math::
-
-        k_{G_\text{nz}}(\mathbf{p}, \mathbf{q}) =
-        \frac{
-            3 (y_p - y_q) (z_p - z_q)
-        }{
-            \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
-        }
-
-    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
-    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`) and
-    :math:`G_\text{en}` is defined as:
-
-    .. math::
-
-        G_\text{nz}(\mathbf{p}) =
-        \frac{
-            \partial^2 V(\mathbf{p})
-        }{
-            \partial y \partial z
-        }
-    """
-    distance = distance_cartesian(
-        (easting_p, northing_p, upward_p), (easting_q, northing_q, upward_q)
-    )
-    return 3 * (northing_p - northing_q) * (upward_p - upward_q) / distance**5
diff --git a/choclo/constants.py b/choclo/constants.py
new file mode 100644
index 00000000..3bff83c2
--- /dev/null
+++ b/choclo/constants.py
@@ -0,0 +1,11 @@
+# Copyright (c) 2022 The Choclo Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Define universal constants
+"""
+#: The gravitational constant in SI units :math:`m^3 kg^{-1} s^{-1}`
+GRAVITATIONAL_CONST = 0.00000000006673
diff --git a/choclo/point/__init__.py b/choclo/point/__init__.py
new file mode 100644
index 00000000..a80eb695
--- /dev/null
+++ b/choclo/point/__init__.py
@@ -0,0 +1,33 @@
+# Copyright (c) 2022 The Choclo Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Kernels and forward modelling functions for point sources
+"""
+from ._forward import (
+    gravity_e,
+    gravity_ee,
+    gravity_en,
+    gravity_eu,
+    gravity_n,
+    gravity_nn,
+    gravity_nu,
+    gravity_pot,
+    gravity_u,
+    gravity_uu,
+)
+from ._kernels import (
+    kernel_e,
+    kernel_ee,
+    kernel_en,
+    kernel_eu,
+    kernel_n,
+    kernel_nn,
+    kernel_nu,
+    kernel_pot,
+    kernel_u,
+    kernel_uu,
+)
diff --git a/choclo/point/_forward.py b/choclo/point/_forward.py
new file mode 100644
index 00000000..f7393f69
--- /dev/null
+++ b/choclo/point/_forward.py
@@ -0,0 +1,651 @@
+# Copyright (c) 2022 The Choclo Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Forward modelling function for point sources
+"""
+from numba import jit
+
+from ..constants import GRAVITATIONAL_CONST
+from ..utils import distance_cartesian
+from ._kernels import (
+    kernel_e,
+    kernel_ee,
+    kernel_en,
+    kernel_eu,
+    kernel_n,
+    kernel_nn,
+    kernel_nu,
+    kernel_pot,
+    kernel_u,
+    kernel_uu,
+)
+
+
+@jit(nopython=True)
+def gravity_pot(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, mass):
+    r"""
+    Gravitational potential field due to a point source
+
+    Returns the gravitational potential field produced by a single point source
+    on a single computation point
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of the observation point in meters.
+    northing_p : float
+        Northing coordinate of the observation point in meters.
+    upward_p : float
+        Upward coordinate of the observation point in meters.
+    easting_q : float
+        Easting coordinate of the point source in meters.
+    northing_q : float
+        Northing coordinate of the point source in meters.
+    upward_q : float
+        Upward coordinate of the point source in meters.
+    mass : float
+        Mass of the point source in kilograms.
+
+    Returns
+    -------
+    potential : float
+        Gravitational potential field generated by the point source on the
+        observation point in :math:`\text{J}/\text{kg}`.
+
+    Notes
+    -----
+    Returns the gravitational potential field :math:`V(\mathbf{p})` on the
+    observation point :math:`\mathbf{p} = (x_p, y_p, z_p)` generated by
+    a single point source located in :math:`\mathbf{q} = (x_q, y_q, z_q)` and
+    mass :math:`m`.
+
+    .. math::
+
+        V(\mathbf{p}) = G m \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert_2^2}
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`)
+    and :math:`G` is the Universal Gravitational Constant.
+    """
+    distance = distance_cartesian(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
+    )
+    kernel = kernel_pot(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
+    return GRAVITATIONAL_CONST * mass * kernel
+
+
+@jit(nopython=True)
+def gravity_e(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, mass):
+    r"""
+    Easting component of the gravitational acceleration due to a point source
+
+    Returns the easting component of the gravitational acceleration produced by
+    a single point source on a single computation point
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of the observation point in meters.
+    northing_p : float
+        Northing coordinate of the observation point in meters.
+    upward_p : float
+        Upward coordinate of the observation point in meters.
+    easting_q : float
+        Easting coordinate of the point source in meters.
+    northing_q : float
+        Northing coordinate of the point source in meters.
+    upward_q : float
+        Upward coordinate of the point source in meters.
+    mass : float
+        Mass of the point source in kilograms.
+
+    Returns
+    -------
+    g_e : float
+        Easting component of the gravitational acceleration generated by the
+        point source on the observation point in :math:`\text{m}/\text{s}^2`.
+
+    Notes
+    -----
+    Returns the easting component :math:`g_x(\mathbf{p})` of the gravitational
+    acceleration :math:`\mathbf{g}` on the observation point :math:`\mathbf{p}
+    = (x_p, y_p, z_p)` generated by a single point source located in
+    :math:`\mathbf{q} = (x_q, y_q, z_q)` and mass :math:`m`.
+
+    .. math::
+
+        g_x(\mathbf{p}) = - G m
+            \frac{
+                x_p - x_q
+            }{
+                \lVert \mathbf{p} - \mathbf{q} \rVert_2^2
+            }
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`)
+    and :math:`G` is the Universal Gravitational Constant.
+    """
+    distance = distance_cartesian(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
+    )
+    kernel = kernel_e(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
+    return GRAVITATIONAL_CONST * mass * kernel
+
+
+@jit(nopython=True)
+def gravity_n(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, mass):
+    r"""
+    Northing component of the gravitational acceleration due to a point source
+
+    Returns the northing component of the gravitational acceleration produced
+    by a single point source on a single computation point
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of the observation point in meters.
+    northing_p : float
+        Northing coordinate of the observation point in meters.
+    upward_p : float
+        Upward coordinate of the observation point in meters.
+    easting_q : float
+        Easting coordinate of the point source in meters.
+    northing_q : float
+        Northing coordinate of the point source in meters.
+    upward_q : float
+        Upward coordinate of the point source in meters.
+    mass : float
+        Mass of the point source in kilograms.
+
+    Returns
+    -------
+    g_e : float
+        Northing component of the gravitational acceleration generated by the
+        point source on the observation point in :math:`\text{m}/\text{s}^2`.
+
+    Notes
+    -----
+    Returns the northing component :math:`g_y(\mathbf{p})` of the gravitational
+    acceleration :math:`\mathbf{g}` on the observation point :math:`\mathbf{p}
+    = (x_p, y_p, z_p)` generated by a single point source located in
+    :math:`\mathbf{q} = (x_q, y_q, z_q)` and mass :math:`m`.
+
+    .. math::
+
+        g_y(\mathbf{p}) = - G m
+            \frac{
+                y_p - y_q
+            }{
+                \lVert \mathbf{p} - \mathbf{q} \rVert_2^2
+            }
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`)
+    and :math:`G` is the Universal Gravitational Constant.
+    """
+    distance = distance_cartesian(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
+    )
+    kernel = kernel_n(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
+    return GRAVITATIONAL_CONST * mass * kernel
+
+
+@jit(nopython=True)
+def gravity_u(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, mass):
+    r"""
+    Upward component of the gravitational acceleration due to a point source
+
+    Returns the upward component of the gravitational acceleration produced by
+    a single point source on a single computation point
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of the observation point in meters.
+    northing_p : float
+        Northing coordinate of the observation point in meters.
+    upward_p : float
+        Upward coordinate of the observation point in meters.
+    easting_q : float
+        Easting coordinate of the point source in meters.
+    northing_q : float
+        Northing coordinate of the point source in meters.
+    upward_q : float
+        Upward coordinate of the point source in meters.
+    mass : float
+        Mass of the point source in kilograms.
+
+    Returns
+    -------
+    g_u : float
+        Upward component of the gravitational acceleration generated by the
+        point source on the observation point in :math:`\text{m}/\text{s}^2`.
+
+    Notes
+    -----
+    Returns the upward component :math:`g_y(\mathbf{p})` of the gravitational
+    acceleration :math:`\mathbf{g}` on the observation point :math:`\mathbf{p}
+    = (x_p, y_p, z_p)` generated by a single point source located in
+    :math:`\mathbf{q} = (x_q, y_q, z_q)` and mass :math:`m`.
+
+    .. math::
+
+        g_z(\mathbf{p}) = - G m
+            \frac{
+                z_p - z_q
+            }{
+                \lVert \mathbf{p} - \mathbf{q} \rVert_2^2
+            }
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`)
+    and :math:`G` is the Universal Gravitational Constant.
+    """
+    distance = distance_cartesian(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
+    )
+    kernel = kernel_u(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
+    return GRAVITATIONAL_CONST * mass * kernel
+
+
+@jit(nopython=True)
+def gravity_ee(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, mass):
+    r"""
+    Easting-easting component of the gravitational tensor due to a point source
+
+    Returns the easting-easting component of the gravitational tensor produced
+    by a single point source on a single computation point
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of the observation point in meters.
+    northing_p : float
+        Northing coordinate of the observation point in meters.
+    upward_p : float
+        Upward coordinate of the observation point in meters.
+    easting_q : float
+        Easting coordinate of the point source in meters.
+    northing_q : float
+        Northing coordinate of the point source in meters.
+    upward_q : float
+        Upward coordinate of the point source in meters.
+    mass : float
+        Mass of the point source in kilograms.
+
+    Returns
+    -------
+    g_ee : float
+        Easting-easting component of the gravitational tensor generated by the
+        point source on the observation point in :math:`\text{m}/\text{s}^2`.
+
+    Notes
+    -----
+    Returns the easting-easting component :math:`g_{xx}(\mathbf{p})` of the
+    gravitational tensor :math:`\mathbf{T}` on the observation point
+    :math:`\mathbf{p} = (x_p, y_p, z_p)` generated by a single point source
+    located in :math:`\mathbf{q} = (x_q, y_q, z_q)` and mass :math:`m`.
+
+    .. math::
+
+        g_{xx}(\mathbf{p}) = G m
+            \left[
+            \frac{
+                3 (x_p - x_q)^2
+            }{
+                \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
+            }
+            - \frac{
+                1
+            }{
+                \lVert \mathbf{p} - \mathbf{q} \rVert_2^3
+            }
+            \right]
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`)
+    and :math:`G` is the Universal Gravitational Constant.
+    """
+    distance = distance_cartesian(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
+    )
+    kernel = kernel_ee(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
+    return GRAVITATIONAL_CONST * mass * kernel
+
+
+@jit(nopython=True)
+def gravity_nn(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, mass):
+    r"""
+    Northing-northing component of the gravitational tensor due to point source
+
+    Returns the northing-northing component of the gravitational tensor
+    produced by a single point source on a single computation point
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of the observation point in meters.
+    northing_p : float
+        Northing coordinate of the observation point in meters.
+    upward_p : float
+        Upward coordinate of the observation point in meters.
+    easting_q : float
+        Easting coordinate of the point source in meters.
+    northing_q : float
+        Northing coordinate of the point source in meters.
+    upward_q : float
+        Upward coordinate of the point source in meters.
+    mass : float
+        Mass of the point source in kilograms.
+
+    Returns
+    -------
+    g_nn : float
+        Northing-northing component of the gravitational tensor generated by
+        the point source on the observation point in
+        :math:`\text{m}/\text{s}^2`.
+
+    Notes
+    -----
+    Returns the northing-northing component :math:`g_{yy}(\mathbf{p})` of the
+    gravitational tensor :math:`\mathbf{T}` on the observation point
+    :math:`\mathbf{p} = (x_p, y_p, z_p)` generated by a single point source
+    located in :math:`\mathbf{q} = (x_q, y_q, z_q)` and mass :math:`m`.
+
+    .. math::
+
+        g_{yy}(\mathbf{p}) = G m
+            \left[
+            \frac{
+                3 (y_p - y_q)^2
+            }{
+                \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
+            }
+            - \frac{
+                1
+            }{
+                \lVert \mathbf{p} - \mathbf{q} \rVert_2^3
+            }
+            \right]
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`)
+    and :math:`G` is the Universal Gravitational Constant.
+    """
+    distance = distance_cartesian(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
+    )
+    kernel = kernel_nn(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
+    return GRAVITATIONAL_CONST * mass * kernel
+
+
+@jit(nopython=True)
+def gravity_uu(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, mass):
+    r"""
+    Upward-upward component of the gravitational tensor due to a point source
+
+    Returns the upward-upward component of the gravitational tensor
+    produced by a single point source on a single computation point
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of the observation point in meters.
+    northing_p : float
+        Northing coordinate of the observation point in meters.
+    upward_p : float
+        Upward coordinate of the observation point in meters.
+    easting_q : float
+        Easting coordinate of the point source in meters.
+    northing_q : float
+        Northing coordinate of the point source in meters.
+    upward_q : float
+        Upward coordinate of the point source in meters.
+    mass : float
+        Mass of the point source in kilograms.
+
+    Returns
+    -------
+    g_uu : float
+        Upward-upward component of the gravitational tensor generated by
+        the point source on the observation point in
+        :math:`\text{m}/\text{s}^2`.
+
+    Notes
+    -----
+    Returns the upward-upward component :math:`g_{zz}(\mathbf{p})` of the
+    gravitational tensor :math:`\mathbf{T}` on the observation point
+    :math:`\mathbf{p} = (x_p, y_p, z_p)` generated by a single point source
+    located in :math:`\mathbf{q} = (x_q, y_q, z_q)` and mass :math:`m`.
+
+    .. math::
+
+        g_{zz}(\mathbf{p}) = G m
+            \left[
+            \frac{
+                3 (z_p - z_q)^2
+            }{
+                \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
+            }
+            - \frac{
+                1
+            }{
+                \lVert \mathbf{p} - \mathbf{q} \rVert_2^3
+            }
+            \right]
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`)
+    and :math:`G` is the Universal Gravitational Constant.
+    """
+    distance = distance_cartesian(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
+    )
+    kernel = kernel_uu(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
+    return GRAVITATIONAL_CONST * mass * kernel
+
+
+@jit(nopython=True)
+def gravity_en(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, mass):
+    r"""
+    Easting-northing component of the gravitational tensor due to point source
+
+    Returns the easting-northing component of the gravitational tensor
+    produced by a single point source on a single computation point
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of the observation point in meters.
+    northing_p : float
+        Northing coordinate of the observation point in meters.
+    upward_p : float
+        Upward coordinate of the observation point in meters.
+    easting_q : float
+        Easting coordinate of the point source in meters.
+    northing_q : float
+        Northing coordinate of the point source in meters.
+    upward_q : float
+        Upward coordinate of the point source in meters.
+    mass : float
+        Mass of the point source in kilograms.
+
+    Returns
+    -------
+    g_en : float
+        Easting-northing component of the gravitational tensor generated by
+        the point source on the observation point in
+        :math:`\text{m}/\text{s}^2`.
+
+    Notes
+    -----
+    Returns the easting-northing component :math:`g_{xy}(\mathbf{p})` of the
+    gravitational tensor :math:`\mathbf{T}` on the observation point
+    :math:`\mathbf{p} = (x_p, y_p, z_p)` generated by a single point source
+    located in :math:`\mathbf{q} = (x_q, y_q, z_q)` and mass :math:`m`.
+
+    .. math::
+
+        g_{xy}(\mathbf{p}) = G m
+            \left[
+            \frac{
+                3 (x_p - x_q) (y_p - y_q)
+            }{
+                \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
+            }
+            \right]
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`)
+    and :math:`G` is the Universal Gravitational Constant.
+    """
+    distance = distance_cartesian(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
+    )
+    kernel = kernel_en(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
+    return GRAVITATIONAL_CONST * mass * kernel
+
+
+@jit(nopython=True)
+def gravity_eu(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, mass):
+    r"""
+    Easting-upward component of the gravitational tensor due to point source
+
+    Returns the easting-upward component of the gravitational tensor
+    produced by a single point source on a single computation point
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of the observation point in meters.
+    northing_p : float
+        Northing coordinate of the observation point in meters.
+    upward_p : float
+        Upward coordinate of the observation point in meters.
+    easting_q : float
+        Easting coordinate of the point source in meters.
+    northing_q : float
+        Northing coordinate of the point source in meters.
+    upward_q : float
+        Upward coordinate of the point source in meters.
+    mass : float
+        Mass of the point source in kilograms.
+
+    Returns
+    -------
+    g_eu : float
+        Easting-upward component of the gravitational tensor generated by
+        the point source on the observation point in
+        :math:`\text{m}/\text{s}^2`.
+
+    Notes
+    -----
+    Returns the easting-upward component :math:`g_{xz}(\mathbf{p})` of the
+    gravitational tensor :math:`\mathbf{T}` on the observation point
+    :math:`\mathbf{p} = (x_p, y_p, z_p)` generated by a single point source
+    located in :math:`\mathbf{q} = (x_q, y_q, z_q)` and mass :math:`m`.
+
+    .. math::
+
+        g_{xz}(\mathbf{p}) = G m
+            \left[
+            \frac{
+                3 (x_p - x_q) (z_p - z_q)
+            }{
+                \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
+            }
+            \right]
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`)
+    and :math:`G` is the Universal Gravitational Constant.
+    """
+    distance = distance_cartesian(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
+    )
+    kernel = kernel_eu(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
+    return GRAVITATIONAL_CONST * mass * kernel
+
+
+@jit(nopython=True)
+def gravity_nu(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, mass):
+    r"""
+    Northing-upward component of the gravitational tensor due to point source
+
+    Returns the northing-upward component of the gravitational tensor
+    produced by a single point source on a single computation point
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of the observation point in meters.
+    northing_p : float
+        Northing coordinate of the observation point in meters.
+    upward_p : float
+        Upward coordinate of the observation point in meters.
+    easting_q : float
+        Easting coordinate of the point source in meters.
+    northing_q : float
+        Northing coordinate of the point source in meters.
+    upward_q : float
+        Upward coordinate of the point source in meters.
+    mass : float
+        Mass of the point source in kilograms.
+
+    Returns
+    -------
+    g_nu : float
+        Northing-upward component of the gravitational tensor generated by the
+        point source on the observation point in :math:`\text{m}/\text{s}^2`.
+
+    Notes
+    -----
+    Returns the northing-upward component :math:`g_{xz}(\mathbf{p})` of the
+    gravitational tensor :math:`\mathbf{T}` on the observation point
+    :math:`\mathbf{p} = (x_p, y_p, z_p)` generated by a single point source
+    located in :math:`\mathbf{q} = (x_q, y_q, z_q)` and mass :math:`m`.
+
+    .. math::
+
+        g_{yz}(\mathbf{p}) = G m
+            \left[
+            \frac{
+                3 (y_p - y_q) (z_p - z_q)
+            }{
+                \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
+            }
+            \right]
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`)
+    and :math:`G` is the Universal Gravitational Constant.
+    """
+    distance = distance_cartesian(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
+    )
+    kernel = kernel_nu(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
+    return GRAVITATIONAL_CONST * mass * kernel
diff --git a/choclo/point/_kernels.py b/choclo/point/_kernels.py
new file mode 100644
index 00000000..dc061d5a
--- /dev/null
+++ b/choclo/point/_kernels.py
@@ -0,0 +1,630 @@
+# Copyright (c) 2022 The Choclo Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Kernel functions for point sources
+"""
+from numba import jit
+
+
+@jit(nopython=True)
+def kernel_pot(
+    easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+):
+    r"""
+    The inverse of the distance between the two points
+
+    .. important ::
+
+        The coordinates of the two points must be in Cartesian coordinates and
+        have the same units.
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of point :math:`\mathbf{p}`.
+    northing_p : float
+        Northing coordinate of point :math:`\mathbf{p}`.
+    upward_p : float
+        Upward coordinate of point :math:`\mathbf{p}`.
+    easting_q : float
+        Easting coordinate of point :math:`\mathbf{q}`.
+    northing_q : float
+        Northing coordinate of point :math:`\mathbf{q}`.
+    upward_q : float
+        Upward coordinate of point :math:`\mathbf{q}`.
+    distance : float
+        Euclidean distance between points :math:`\mathbf{p}` and
+        :math:`\mathbf{q}`.
+
+    Returns
+    -------
+    kernel : float
+        Value of the kernel function.
+
+    Notes
+    -----
+    Given two points :math:`\mathbf{p} = (x_p, y_p, z_p)` and :math:`\mathbf{q}
+    = (x_q, y_q, z_q)` defined in a Cartesian coordinate system, compute the
+    following kernel function:
+
+    .. math::
+
+        k_V(\mathbf{p}, \mathbf{q}) =
+        \frac{
+            1
+        }{
+            \lVert \mathbf{p} - \mathbf{q} \rVert_2
+        }
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`).
+    """
+    return 1 / distance
+
+
+@jit(nopython=True)
+def kernel_e(
+    easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+):
+    r"""
+    Easting component of the gradient of the inverse of the distance
+
+    .. important ::
+
+        The coordinates of the two points must be in Cartesian coordinates and
+        have the same units.
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of point :math:`\mathbf{p}`.
+    northing_p : float
+        Northing coordinate of point :math:`\mathbf{p}`.
+    upward_p : float
+        Upward coordinate of point :math:`\mathbf{p}`.
+    easting_q : float
+        Easting coordinate of point :math:`\mathbf{q}`.
+    northing_q : float
+        Northing coordinate of point :math:`\mathbf{q}`.
+    upward_q : float
+        Upward coordinate of point :math:`\mathbf{q}`.
+    distance : float
+        Euclidean distance between points :math:`\mathbf{p}` and
+        :math:`\mathbf{q}`.
+
+    Returns
+    -------
+    kernel : float
+        Value of the kernel function.
+
+    Notes
+    -----
+    Given two points :math:`\mathbf{p} = (x_p, y_p, z_p)` and :math:`\mathbf{q}
+    = (x_q, y_q, z_q)` defined in a Cartesian coordinate system, compute the
+    following kernel function:
+
+    .. math::
+
+        k_x(\mathbf{p}, \mathbf{q}) =
+        \frac{\partial}{\partial x}
+        \left(
+            \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert_2}
+        \right)
+        =
+        - \frac{
+            x_p - x_q
+        }{
+            \lVert \mathbf{p} - \mathbf{q} \rVert_2
+        }
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`).
+    """
+    return -(easting_p - easting_q) / distance**3
+
+
+@jit(nopython=True)
+def kernel_n(
+    easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+):
+    r"""
+    Northing component of the gradient of the inverse of the distance
+
+    .. important ::
+
+        The coordinates of the two points must be in Cartesian coordinates and
+        have the same units.
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of point :math:`\mathbf{p}`.
+    northing_p : float
+        Northing coordinate of point :math:`\mathbf{p}`.
+    upward_p : float
+        Upward coordinate of point :math:`\mathbf{p}`.
+    easting_q : float
+        Easting coordinate of point :math:`\mathbf{q}`.
+    northing_q : float
+        Northing coordinate of point :math:`\mathbf{q}`.
+    upward_q : float
+        Upward coordinate of point :math:`\mathbf{q}`.
+    distance : float
+        Euclidean distance between points :math:`\mathbf{p}` and
+        :math:`\mathbf{q}`.
+
+    Returns
+    -------
+    kernel : float
+        Value of the kernel function.
+
+    Notes
+    -----
+    Given two points :math:`\mathbf{p} = (x_p, y_p, z_p)` and :math:`\mathbf{q}
+    = (x_q, y_q, z_q)` defined in a Cartesian coordinate system, compute the
+    following kernel function:
+
+    .. math::
+
+        k_y(\mathbf{p}, \mathbf{q}) =
+        \frac{\partial}{\partial y}
+        \left(
+            \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert_2}
+        \right)
+        =
+        - \frac{
+            y_p - y_q
+        }{
+            \lVert \mathbf{p} - \mathbf{q} \rVert_2
+        }
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`).
+    """
+    return -(northing_p - northing_q) / distance**3
+
+
+@jit(nopython=True)
+def kernel_u(
+    easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+):
+    r"""
+    Upward component of the gradient of the inverse of the distance
+
+    .. important ::
+
+        The coordinates of the two points must be in Cartesian coordinates and
+        have the same units.
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of point :math:`\mathbf{p}`.
+    northing_p : float
+        Northing coordinate of point :math:`\mathbf{p}`.
+    upward_p : float
+        Upward coordinate of point :math:`\mathbf{p}`.
+    easting_q : float
+        Easting coordinate of point :math:`\mathbf{q}`.
+    northing_q : float
+        Northing coordinate of point :math:`\mathbf{q}`.
+    upward_q : float
+        Upward coordinate of point :math:`\mathbf{q}`.
+    distance : float
+        Euclidean distance between points :math:`\mathbf{p}` and
+        :math:`\mathbf{q}`.
+
+    Returns
+    -------
+    kernel : float
+        Value of the kernel function.
+
+    Notes
+    -----
+    Given two points :math:`\mathbf{p} = (x_p, y_p, z_p)` and :math:`\mathbf{q}
+    = (x_q, y_q, z_q)` defined in a Cartesian coordinate system, compute the
+    following kernel function:
+
+    .. math::
+
+        k_z(\mathbf{p}, \mathbf{q}) =
+        \frac{\partial}{\partial z}
+        \left(
+            \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert_2}
+        \right)
+        =
+        - \frac{
+            z_p - z_q
+        }{
+            \lVert \mathbf{p} - \mathbf{q} \rVert_2
+        }
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`).
+    """
+    return -(upward_p - upward_q) / distance**3
+
+
+@jit(nopython=True)
+def kernel_ee(
+    easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+):
+    r"""
+    Second derivative of the inverse of the distance along easting-easting
+
+    .. important ::
+
+        The coordinates of the two points must be in Cartesian coordinates and
+        have the same units.
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of point :math:`\mathbf{p}`.
+    northing_p : float
+        Northing coordinate of point :math:`\mathbf{p}`.
+    upward_p : float
+        Upward coordinate of point :math:`\mathbf{p}`.
+    easting_q : float
+        Easting coordinate of point :math:`\mathbf{q}`.
+    northing_q : float
+        Northing coordinate of point :math:`\mathbf{q}`.
+    upward_q : float
+        Upward coordinate of point :math:`\mathbf{q}`.
+    distance : float
+        Euclidean distance between points :math:`\mathbf{p}` and
+        :math:`\mathbf{q}`.
+
+    Returns
+    -------
+    kernel : float
+        Value of the kernel function.
+
+    Notes
+    -----
+    Given two points :math:`\mathbf{p} = (x_p, y_p, z_p)` and :math:`\mathbf{q}
+    = (x_q, y_q, z_q)` defined in a Cartesian coordinate system, compute the
+    following kernel function:
+
+    .. math::
+
+        k_{xx}(\mathbf{p}, \mathbf{q}) =
+        \frac{\partial^2}{\partial x^2}
+        \left(
+            \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert_2}
+        \right)
+        =
+        \frac{
+            3 (x_p - x_q)^2
+        }{
+            \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
+        }
+        - \frac{
+            1
+        }{
+            \lVert \mathbf{p} - \mathbf{q} \rVert_2^3
+        }
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`).
+    """
+    return 3 * (easting_p - easting_q) ** 2 / distance**5 - 1 / distance**3
+
+
+@jit(nopython=True)
+def kernel_nn(
+    easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+):
+    r"""
+    Second derivative of the inverse of the distance along northing-northing
+
+    .. important ::
+
+        The coordinates of the two points must be in Cartesian coordinates and
+        have the same units.
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of point :math:`\mathbf{p}`.
+    northing_p : float
+        Northing coordinate of point :math:`\mathbf{p}`.
+    upward_p : float
+        Upward coordinate of point :math:`\mathbf{p}`.
+    easting_q : float
+        Easting coordinate of point :math:`\mathbf{q}`.
+    northing_q : float
+        Northing coordinate of point :math:`\mathbf{q}`.
+    upward_q : float
+        Upward coordinate of point :math:`\mathbf{q}`.
+    distance : float
+        Euclidean distance between points :math:`\mathbf{p}` and
+        :math:`\mathbf{q}`.
+
+    Returns
+    -------
+    kernel : float
+        Value of the kernel function.
+
+    Notes
+    -----
+    Given two points :math:`\mathbf{p} = (x_p, y_p, z_p)` and :math:`\mathbf{q}
+    = (x_q, y_q, z_q)` defined in a Cartesian coordinate system, compute the
+    following kernel function:
+
+    .. math::
+
+        k_{yy}(\mathbf{p}, \mathbf{q}) =
+        \frac{\partial^2}{\partial y^2}
+        \left(
+            \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert_2}
+        \right)
+        =
+        \frac{
+            3 (y_p - y_q)^2
+        }{
+            \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
+        }
+        - \frac{
+            1
+        }{
+            \lVert \mathbf{p} - \mathbf{q} \rVert_2^3
+        }
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`).
+    """
+    return 3 * (northing_p - northing_q) ** 2 / distance**5 - 1 / distance**3
+
+
+@jit(nopython=True)
+def kernel_uu(
+    easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+):
+    r"""
+    Second derivative of the inverse of the distance along upward-upward
+
+    .. important ::
+
+        The coordinates of the two points must be in Cartesian coordinates and
+        have the same units.
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of point :math:`\mathbf{p}`.
+    northing_p : float
+        Northing coordinate of point :math:`\mathbf{p}`.
+    upward_p : float
+        Upward coordinate of point :math:`\mathbf{p}`.
+    easting_q : float
+        Easting coordinate of point :math:`\mathbf{q}`.
+    northing_q : float
+        Northing coordinate of point :math:`\mathbf{q}`.
+    upward_q : float
+        Upward coordinate of point :math:`\mathbf{q}`.
+    distance : float
+        Euclidean distance between points :math:`\mathbf{p}` and
+        :math:`\mathbf{q}`.
+
+    Returns
+    -------
+    kernel : float
+        Value of the kernel function.
+
+    Notes
+    -----
+    Given two points :math:`\mathbf{p} = (x_p, y_p, z_p)` and :math:`\mathbf{q}
+    = (x_q, y_q, z_q)` defined in a Cartesian coordinate system, compute the
+    following kernel function:
+
+    .. math::
+
+        k_{zz}(\mathbf{p}, \mathbf{q}) =
+        \frac{\partial^2}{\partial z^2}
+        \left(
+            \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert_2}
+        \right)
+        =
+        \frac{
+            3 (z_p - z_q)^2
+        }{
+            \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
+        }
+        - \frac{
+            1
+        }{
+            \lVert \mathbf{p} - \mathbf{q} \rVert_2^3
+        }
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`).
+    """
+    return 3 * (upward_p - upward_q) ** 2 / distance**5 - 1 / distance**3
+
+
+@jit(nopython=True)
+def kernel_en(
+    easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+):
+    r"""
+    Second derivative of the inverse of the distance along easting-northing
+
+    .. important ::
+
+        The coordinates of the two points must be in Cartesian coordinates and
+        have the same units.
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of point :math:`\mathbf{p}`.
+    northing_p : float
+        Northing coordinate of point :math:`\mathbf{p}`.
+    upward_p : float
+        Upward coordinate of point :math:`\mathbf{p}`.
+    easting_q : float
+        Easting coordinate of point :math:`\mathbf{q}`.
+    northing_q : float
+        Northing coordinate of point :math:`\mathbf{q}`.
+    upward_q : float
+        Upward coordinate of point :math:`\mathbf{q}`.
+    distance : float
+        Euclidean distance between points :math:`\mathbf{p}` and
+        :math:`\mathbf{q}`.
+
+    Returns
+    -------
+    kernel : float
+        Value of the kernel function.
+
+    Notes
+    -----
+    Given two points :math:`\mathbf{p} = (x_p, y_p, z_p)` and :math:`\mathbf{q}
+    = (x_q, y_q, z_q)` defined in a Cartesian coordinate system, compute the
+    following kernel function:
+
+    .. math::
+
+        k_{xy}(\mathbf{p}, \mathbf{q}) =
+        \frac{\partial}{\partial x \partial y}
+        \left(
+            \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert_2}
+        \right)
+        =
+        \frac{
+            3 (x_p - x_q) (y_p - y_q)
+        }{
+            \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
+        }
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`).
+    """
+    return 3 * (easting_p - easting_q) * (northing_p - northing_q) / distance**5
+
+
+@jit(nopython=True)
+def kernel_eu(
+    easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+):
+    r"""
+    Second derivative of the inverse of the distance along easting-upward
+
+    .. important ::
+
+        The coordinates of the two points must be in Cartesian coordinates and
+        have the same units.
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of point :math:`\mathbf{p}`.
+    northing_p : float
+        Northing coordinate of point :math:`\mathbf{p}`.
+    upward_p : float
+        Upward coordinate of point :math:`\mathbf{p}`.
+    easting_q : float
+        Easting coordinate of point :math:`\mathbf{q}`.
+    northing_q : float
+        Northing coordinate of point :math:`\mathbf{q}`.
+    upward_q : float
+        Upward coordinate of point :math:`\mathbf{q}`.
+    distance : float
+        Euclidean distance between points :math:`\mathbf{p}` and
+        :math:`\mathbf{q}`.
+
+    Returns
+    -------
+    kernel : float
+        Value of the kernel function.
+
+    Notes
+    -----
+    Given two points :math:`\mathbf{p} = (x_p, y_p, z_p)` and :math:`\mathbf{q}
+    = (x_q, y_q, z_q)` defined in a Cartesian coordinate system, compute the
+    following kernel function:
+
+    .. math::
+
+        k_{xz}(\mathbf{p}, \mathbf{q}) =
+        \frac{\partial}{\partial x \partial z}
+        \left(
+            \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert_2}
+        \right)
+        =
+        \frac{
+            3 (x_p - x_q) (z_p - z_q)
+        }{
+            \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
+        }
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`).
+    """
+    return 3 * (easting_p - easting_q) * (upward_p - upward_q) / distance**5
+
+
+@jit(nopython=True)
+def kernel_nu(
+    easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+):
+    r"""
+    Second derivative of the inverse of the distance along northing-upward
+
+    .. important ::
+
+        The coordinates of the two points must be in Cartesian coordinates and
+        have the same units.
+
+    Parameters
+    ----------
+    easting_p : float
+        Easting coordinate of point :math:`\mathbf{p}`.
+    northing_p : float
+        Northing coordinate of point :math:`\mathbf{p}`.
+    upward_p : float
+        Upward coordinate of point :math:`\mathbf{p}`.
+    easting_q : float
+        Easting coordinate of point :math:`\mathbf{q}`.
+    northing_q : float
+        Northing coordinate of point :math:`\mathbf{q}`.
+    upward_q : float
+        Upward coordinate of point :math:`\mathbf{q}`.
+    distance : float
+        Euclidean distance between points :math:`\mathbf{p}` and
+        :math:`\mathbf{q}`.
+
+    Returns
+    -------
+    kernel : float
+        Value of the kernel function.
+
+    Notes
+    -----
+    Given two points :math:`\mathbf{p} = (x_p, y_p, z_p)` and :math:`\mathbf{q}
+    = (x_q, y_q, z_q)` defined in a Cartesian coordinate system, compute the
+    following kernel function:
+
+    .. math::
+
+        k_{yz}(\mathbf{p}, \mathbf{q}) =
+        \frac{\partial}{\partial y \partial z}
+        \left(
+            \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert_2}
+        \right)
+        =
+        \frac{
+            3 (y_p - y_q) (z_p - z_q)
+        }{
+            \lVert \mathbf{p} - \mathbf{q} \rVert_2^5
+        }
+
+    where :math:`\lVert \cdot \rVert_2` refer to the :math:`L_2` norm (the
+    Euclidean distance between :math:`\mathbf{p}` and :math:`\mathbf{q}`).
+    """
+    return 3 * (northing_p - northing_q) * (upward_p - upward_q) / distance**5
diff --git a/choclo/_distance.py b/choclo/utils.py
similarity index 78%
rename from choclo/_distance.py
rename to choclo/utils.py
index 81ef26fb..1a07325f 100644
--- a/choclo/_distance.py
+++ b/choclo/utils.py
@@ -5,27 +5,37 @@
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #
 """
-Compute the distance between two points in different coordinate systems
+Utility functions for the kernel and forward modelling functions
 """
 import numpy as np
 from numba import jit
 
 
 @jit(nopython=True)
-def distance_cartesian(point_p, point_q):
+def distance_cartesian(
+    easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
+):
     r"""
     Euclidean distance between two points given in Cartesian coordinates
 
+    .. warning::
+
+        All coordinates should be in the same units.
+
     Parameters
     ----------
-    point_p : tuple or 1d-array
-        Tuple or array containing the coordinates of the first point in the
-        following order: (``easting``, ``northing`` and ``upward``)
-        All coordinates must be in the same units.
-    point_q : tuple or 1d-array
-        Tuple or array containing the coordinates of the second point in the
-        following order: (``easting``, ``northing`` and ``upward``)
-        All coordinates must be in the same units.
+    easting_p : float
+        Easting coordinate of point :math:`\mathbf{p}`.
+    northing_p : float
+        Northing coordinate of point :math:`\mathbf{p}`.
+    upward_p : float
+        Upward coordinate of point :math:`\mathbf{p}`.
+    easting_q : float
+        Easting coordinate of point :math:`\mathbf{q}`.
+    northing_q : float
+        Northing coordinate of point :math:`\mathbf{q}`.
+    upward_q : float
+        Upward coordinate of point :math:`\mathbf{q}`.
 
     Returns
     -------
@@ -44,8 +54,6 @@ def distance_cartesian(point_p, point_q):
         d = \sqrt{(x_p - x_q)^2 + (y_p - y_q)^2 + (z_p - z_q)^2}
 
     """
-    easting_p, northing_p, upward_p = point_p[:]
-    easting_q, northing_q, upward_q = point_q[:]
     distance = np.sqrt(
         (easting_p - easting_q) ** 2
         + (northing_p - northing_q) ** 2
@@ -55,7 +63,9 @@ def distance_cartesian(point_p, point_q):
 
 
 @jit(nopython=True)
-def distance_spherical(point_p, point_q):
+def distance_spherical(
+    longitude_p, latitude_p, radius_p, longitude_q, latitude_q, radius_q
+):
     r"""
     Euclidean distance between two points in spherical coordinates
 
@@ -65,16 +75,18 @@ def distance_spherical(point_p, point_q):
 
     Parameters
     ----------
-    point_p : tuple or 1d-array
-        Tuple or array containing the coordinates of the first point in the
-        following order: (``longitude``, ``latitude`` and ``radius``).
-        Both ``longitude`` and ``latitude`` must be in degrees, while
-        ``radius`` in meters.
-    point_q : tuple or 1d-array
-        Tuple or array containing the coordinates of the second point in the
-        following order: (``longitude``, ``latitude`` and ``radius``).
-        Both ``longitude`` and ``latitude`` must be in degrees, while
-        ``radius`` in meters.
+    longitude_p : float
+        Longitude coordinate of point :math:`\mathbf{p}` in degrees.
+    latitude_p : float
+        Latitude coordinate of point :math:`\mathbf{p}` in degrees.
+    radius_p : float
+        Radial coordinate of point :math:`\mathbf{p}` in meters.
+    longitude_q : float
+        Longitude coordinate of point :math:`\mathbf{q}` in degrees.
+    latitude_q : float
+        Latitude coordinate of point :math:`\mathbf{q}` in degrees.
+    radius_q : float
+        Radial coordinate of point :math:`\mathbf{q}` in meters.
 
     Returns
     -------
@@ -103,9 +115,6 @@ def distance_spherical(point_p, point_q):
     and :math:`\lambda` is the longitude angle, :math:`\phi` the spherical
     latitude angle an :math:`r` is the radial coordinate.
     """
-    # Get coordinates of the two points
-    longitude_p, latitude_p, radius_p = point_p[:]
-    longitude_q, latitude_q, radius_q = point_q[:]
     # Convert angles to radians
     longitude_p, latitude_p = np.radians(longitude_p), np.radians(latitude_p)
     longitude_q, latitude_q = np.radians(longitude_q), np.radians(latitude_q)
diff --git a/doc/api/index.rst b/doc/api/index.rst
index 64a4d3c7..42da3b5d 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -7,51 +7,67 @@ List of functions and classes (API)
 
 .. currentmodule:: choclo
 
-Kernel functions
-----------------
+Kernel and forward modelling functions
+--------------------------------------
 
-Here you will find the list of available kernel functions for gravity and
-magnetic forward modellings.
+Here you will find the list of available kernel and forward modelling functions
+for gravity and magnetic fields generated by different type of sources.
 
-Point sources
-~~~~~~~~~~~~~
+Point sources and dipoles
+~~~~~~~~~~~~~~~~~~~~~~~~~
 
 For point sources and observation points defined in **Cartesian coordinates**.
 
-Potential field
-^^^^^^^^^^^^^^^
+Gravity
+^^^^^^^
 
 .. autosummary::
    :toctree: generated/
 
-    kernel_point_potential
+    point.gravity_pot
+    point.gravity_e
+    point.gravity_n
+    point.gravity_u
+    point.gravity_ee
+    point.gravity_nn
+    point.gravity_uu
+    point.gravity_en
+    point.gravity_eu
+    point.gravity_nu
 
-Gradient components
-^^^^^^^^^^^^^^^^^^^
+Magnetic
+^^^^^^^^
 
 .. autosummary::
    :toctree: generated/
 
-    kernel_point_g_easting
-    kernel_point_g_northing
-    kernel_point_g_upward
+    point.magnetic_e
+    point.magnetic_n
+    point.magnetic_u
 
-Tensor components
-^^^^^^^^^^^^^^^^^
+Kernels
+^^^^^^^
 
 .. autosummary::
    :toctree: generated/
 
-    kernel_point_g_ee
-    kernel_point_g_nn
-    kernel_point_g_zz
-    kernel_point_g_en
-    kernel_point_g_ez
-    kernel_point_g_nz
+    point.kernel_pot
+    point.kernel_e
+    point.kernel_n
+    point.kernel_u
+    point.kernel_ee
+    point.kernel_nn
+    point.kernel_uu
+    point.kernel_en
+    point.kernel_eu
+    point.kernel_nu
 
 
+Utilities
+---------
+
 Euclidean distances
--------------------
+~~~~~~~~~~~~~~~~~~~
 
 Use these functions to compute Euclidean distance in Cartesian and spherical
 coordinates:
@@ -59,7 +75,17 @@ coordinates:
 .. autosummary::
    :toctree: generated/
 
-    distance_cartesian
-    distance_spherical
-    distance_spherical_core
+    utils.distance_cartesian
+    utils.distance_spherical
+    utils.distance_spherical_core
+
+
+Universal Constants
+~~~~~~~~~~~~~~~~~~~
+
+Access universal physical constants given in SI units:
+
+.. autosummary::
+   :toctree: generated/
 
+    constants.GRAVITATIONAL_CONST
diff --git a/setup.cfg b/setup.cfg
index 5381b782..2d476ad9 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -63,6 +63,8 @@ exclude =
 per-file-ignores =
     # disable unused-imports errors on __init__.py
     __init__.py: F401
+    # disable unused arguments in kernel files
+    _kernels.py: U100
 
 # Configure flake8-rst-docstrings
 # -------------------------------


Test patch:
diff --git a/choclo/tests/test_distance.py b/choclo/tests/test_distance.py
index 0424fd3d..369233ea 100644
--- a/choclo/tests/test_distance.py
+++ b/choclo/tests/test_distance.py
@@ -7,7 +7,7 @@
 import numpy.testing as npt
 import pytest
 
-from .._distance import distance_cartesian, distance_spherical
+from ..utils import distance_cartesian, distance_spherical
 from .utils import dumb_spherical_distance
 
 
@@ -24,7 +24,7 @@ def test_distance_cartesian(point_a, point_b, expected_distance):
     """
     Test if distance_cartesian works as expected
     """
-    npt.assert_allclose(distance_cartesian(point_a, point_b), expected_distance)
+    npt.assert_allclose(distance_cartesian(*point_a, *point_b), expected_distance)
 
 
 class TestDistanceSpherical:
@@ -67,7 +67,7 @@ def test_pole_and_equator(self, point_equator, point_pole):
         """
         expected_distance = 2 ** (3 / 2)
         npt.assert_allclose(
-            distance_spherical(point_equator, point_pole), expected_distance
+            distance_spherical(*point_equator, *point_pole), expected_distance
         )
 
     def test_sample_points(self, point_a, point_b):
@@ -75,6 +75,6 @@ def test_sample_points(self, point_a, point_b):
         Test spherical distance using two sample points
         """
         npt.assert_allclose(
-            distance_spherical(point_a, point_b),
+            distance_spherical(*point_a, *point_b),
             dumb_spherical_distance(point_a, point_b),
         )
diff --git a/choclo/tests/test_point_gravity.py b/choclo/tests/test_point_gravity.py
new file mode 100644
index 00000000..ffdd5067
--- /dev/null
+++ b/choclo/tests/test_point_gravity.py
@@ -0,0 +1,701 @@
+# Copyright (c) 2022 The Choclo Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Test gravity forward modelling functions for point sources
+"""
+import numpy as np
+import numpy.testing as npt
+import pytest
+
+from ..point import (
+    gravity_e,
+    gravity_ee,
+    gravity_en,
+    gravity_eu,
+    gravity_n,
+    gravity_nn,
+    gravity_nu,
+    gravity_pot,
+    gravity_u,
+    gravity_uu,
+)
+from .utils import NUMBA_IS_DISABLED
+
+
+@pytest.fixture(name="sample_point_source")
+def fixture_sample_point_source():
+    """
+    Return a sample point source
+    """
+    return (3.7, 4.3, -5.8)
+
+
+@pytest.fixture(name="sample_mass")
+def fixture_sample_mass():
+    """
+    Return the mass for the sample point source
+    """
+    return 9e4
+
+
+@pytest.fixture(name="sample_coordinate")
+def fixture_sample_coordinate():
+    """
+    Define a sample observation point
+    """
+    return (16.7, 13.2, 7.8)
+
+
+class TestSymmetryPotential:
+    """
+    Test the symmetry of gravity potential due to a point source
+    """
+
+    def build_points_in_sphere(self, radius, point):
+        """
+        Generate a set of observation points located in a sphere of radius
+        ``radius`` and centered around the ``point``.
+        """
+        longitude = np.arange(0, 2 * np.pi, np.pi / 8)
+        latitude = np.arange(-np.pi / 2, np.pi / 2, np.pi / 8)
+        longitude, latitude = np.meshgrid(longitude, latitude)
+        x = radius * np.cos(longitude) * np.cos(latitude)
+        y = radius * np.sin(longitude) * np.cos(latitude)
+        z = radius * np.sin(latitude)
+        x += point[0]
+        y += point[1]
+        z += point[2]
+        return tuple(coord.ravel() for coord in (x, y, z))
+
+    def test_symmetry_on_sphere(self, sample_point_source, sample_mass):
+        """
+        Test the symmetry of the gravity potential in points of a sphere
+        """
+        radius = 3.5
+        observation_points = self.build_points_in_sphere(radius, sample_point_source)
+        # Evaluate gravity_pot on every observation point
+        potential = tuple(
+            gravity_pot(*coords, *sample_point_source, sample_mass)
+            for coords in zip(*observation_points)
+        )
+        # Check if all values are equal
+        npt.assert_allclose(potential[0], potential)
+
+    def test_potential_between_two_spheres(self, sample_point_source, sample_mass):
+        """
+        Test the gravity potential between observation points in two spheres
+        """
+        radius_1, radius_2 = 3.5, 8.7
+        sphere_1 = self.build_points_in_sphere(radius_1, sample_point_source)
+        sphere_2 = self.build_points_in_sphere(radius_2, sample_point_source)
+        # Evaluate kernel_pot on every observation point
+        potential_1 = gravity_pot(*sphere_1, *sample_point_source, sample_mass)
+        potential_2 = gravity_pot(*sphere_2, *sample_point_source, sample_mass)
+        # Check if all values are equal
+        npt.assert_allclose(potential_1, radius_2 / radius_1 * potential_2)
+
+    @pytest.mark.skipif(not NUMBA_IS_DISABLED, reason="Numba is not disabled")
+    def test_infinite_potential(self, sample_mass):
+        """
+        Test if we get an infinite kernel if the computation point is in the
+        same location as the observation point.
+
+        This test should be run only with Numba disabled.
+        """
+        point = (4.6, -8.9, -50.3)
+        assert np.isinf(gravity_pot(*point, *point, sample_mass))
+
+    @pytest.mark.skipif(NUMBA_IS_DISABLED, reason="Numba is disabled")
+    def test_division_by_zero(self, sample_mass):
+        """
+        Test if we get a division by zero if the computation point is in the
+        same location as the observation point.
+
+        This test should be run only with Numba enabled.
+        """
+        point = (4.6, -8.9, -50.3)
+        with pytest.raises(ZeroDivisionError):
+            gravity_pot(*point, *point, sample_mass)
+
+
+class TestSymmetryGravityE:
+    """
+    Test symmetry of the gravity_e
+    """
+
+    @pytest.fixture
+    def coords_northing_upward_plane(self, sample_point_source):
+        """
+        Define set of observation points in the northing-upward plane
+        """
+        # Define northing and upward coordinates (avoid the zero, that's where
+        # the sample_point_source is located)
+        northing = np.linspace(-11, 11, 12)
+        upward = np.linspace(-7, 7, 8)
+        # Generate meshgrid
+        northing, upward = np.meshgrid(northing, upward)
+        # Compute easting
+        easting = np.zeros_like(northing)
+        # Add the coordinates of the sample point source
+        easting += sample_point_source[0]
+        northing += sample_point_source[1]
+        upward += sample_point_source[2]
+        return easting, northing, upward
+
+    @pytest.fixture
+    def mirrored_points(self, sample_point_source):
+        """
+        Define two set of mirrored points across the northing-upward plane
+        """
+        # Define the northing and upward coordinates of the points
+        northing = np.linspace(-11, 11, 27)
+        upward = np.linspace(-38, 38, 23)
+        # Define the easting coordinates for the first set
+        easting_1 = np.linspace(1, 15, 15)
+        # And mirror it for the second set
+        easting_2 = -easting_1
+        # Shift the coordinates to the sample_point_source
+        easting_1 += sample_point_source[0]
+        easting_2 += sample_point_source[0]
+        northing += sample_point_source[1]
+        upward += sample_point_source[2]
+        return (easting_1, northing, upward), (easting_2, northing, upward)
+
+    def test_zero_in_northing_upward(
+        self, sample_point_source, coords_northing_upward_plane, sample_mass
+    ):
+        """
+        Test if kernel_e is zero in the northing-upward plane
+        """
+        g_e = np.array(
+            list(
+                gravity_e(*coords, *sample_point_source, sample_mass)
+                for coords in zip(*coords_northing_upward_plane)
+            )
+        )
+        assert (g_e <= 1e-30).all()
+
+    def test_mirror_symmetry(self, sample_point_source, mirrored_points, sample_mass):
+        """
+        Test points with opposite easting coordinate
+        """
+        coords_1, coords_2 = mirrored_points
+        g_e_1 = np.array(
+            list(
+                gravity_e(*coords, *sample_point_source, sample_mass)
+                for coords in zip(*coords_1)
+            )
+        )
+        g_e_2 = np.array(
+            list(
+                gravity_e(*coords, *sample_point_source, sample_mass)
+                for coords in zip(*coords_2)
+            )
+        )
+        npt.assert_allclose(g_e_1, -g_e_2)
+
+
+class TestSymmetryGravityN:
+    """
+    Test symmetry of the gravity_n
+    """
+
+    @pytest.fixture
+    def coords_easting_upward_plane(self, sample_point_source):
+        """
+        Define set of observation points in the easting-upward plane
+        """
+        # Define easting and upward coordinates (avoid the zero, that's where
+        # the sample_point_source is located)
+        easting = np.linspace(-11, 11, 12)
+        upward = np.linspace(-7, 7, 8)
+        # Generate meshgrid
+        easting, upward = np.meshgrid(easting, upward)
+        # Compute northing
+        northing = np.zeros_like(easting)
+        # Add the coordinates of the sample point source
+        easting += sample_point_source[0]
+        northing += sample_point_source[1]
+        upward += sample_point_source[2]
+        return easting, northing, upward
+
+    @pytest.fixture
+    def mirrored_points(self, sample_point_source):
+        """
+        Define two set of mirrored points across the easting-upward plane
+        """
+        # Define the easting and upward coordinates of the points
+        easting = np.linspace(-11, 11, 27)
+        upward = np.linspace(-38, 38, 23)
+        # Define the northing coordinates for the first set
+        northing_1 = np.linspace(1, 15, 15)
+        # And mirror it for the second set
+        northing_2 = -northing_1
+        # Shift the coordinates to the sample_point_source
+        easting += sample_point_source[0]
+        northing_1 += sample_point_source[1]
+        northing_2 += sample_point_source[1]
+        upward += sample_point_source[2]
+        return (easting, northing_1, upward), (easting, northing_2, upward)
+
+    def test_zero_in_easting_upward(
+        self, sample_point_source, coords_easting_upward_plane, sample_mass
+    ):
+        """
+        Test if kernel_e is zero in the northing-upward plane
+        """
+        g_n = np.array(
+            list(
+                gravity_n(*coords, *sample_point_source, sample_mass)
+                for coords in zip(*coords_easting_upward_plane)
+            )
+        )
+        assert (g_n <= 1e-30).all()
+
+    def test_mirror_symmetry(self, sample_point_source, mirrored_points, sample_mass):
+        """
+        Test points with opposite easting coordinate
+        """
+        coords_1, coords_2 = mirrored_points
+        g_n_1 = np.array(
+            list(
+                gravity_n(*coords, *sample_point_source, sample_mass)
+                for coords in zip(*coords_1)
+            )
+        )
+        g_n_2 = np.array(
+            list(
+                gravity_n(*coords, *sample_point_source, sample_mass)
+                for coords in zip(*coords_2)
+            )
+        )
+        npt.assert_allclose(g_n_1, -g_n_2)
+
+
+class TestSymmetryGravityU:
+    """
+    Test symmetry of the gravity_u
+    """
+
+    @pytest.fixture
+    def coords_easting_northing_plane(self, sample_point_source):
+        """
+        Define set of observation points in the easting-northing plane
+        """
+        # Define easting and northing coordinates (avoid the zero, that's where
+        # the sample_point_source is located)
+        easting = np.linspace(-11, 11, 12)
+        northing = np.linspace(-7, 7, 8)
+        # Generate meshgrid
+        easting, northing = np.meshgrid(easting, northing)
+        # Compute upward
+        upward = np.zeros_like(easting)
+        # Add the coordinates of the sample point source
+        easting += sample_point_source[0]
+        northing += sample_point_source[1]
+        upward += sample_point_source[2]
+        return easting, northing, upward
+
+    @pytest.fixture
+    def mirrored_points(self, sample_point_source):
+        """
+        Define two set of mirrored points across the easting-northing plane
+        """
+        # Define the easting and upward coordinates of the points
+        easting = np.linspace(-11, 11, 27)
+        northing = np.linspace(-38, 38, 23)
+        # Define the northing coordinates for the first set
+        upward_1 = np.linspace(1, 15, 15)
+        # And mirror it for the second set
+        upward_2 = -upward_1
+        # Shift the coordinates to the sample_point_source
+        easting += sample_point_source[0]
+        northing += sample_point_source[1]
+        upward_1 += sample_point_source[2]
+        upward_2 += sample_point_source[2]
+        return (easting, northing, upward_1), (easting, northing, upward_2)
+
+    def test_zero_in_easting_northing(
+        self, sample_point_source, coords_easting_northing_plane, sample_mass
+    ):
+        """
+        Test if gravity_u is zero in the easting-upward plane
+        """
+        g_u = np.array(
+            list(
+                gravity_u(*coords, *sample_point_source, sample_mass)
+                for coords in zip(*coords_easting_northing_plane)
+            )
+        )
+        assert (g_u <= 1e-30).all()
+
+    def test_mirror_symmetry(self, sample_point_source, mirrored_points, sample_mass):
+        """
+        Test points with opposite upward coordinate
+        """
+        coords_1, coords_2 = mirrored_points
+        g_u_1 = np.array(
+            list(
+                gravity_u(*coords, *sample_point_source, sample_mass)
+                for coords in zip(*coords_1)
+            )
+        )
+        g_u_2 = np.array(
+            list(
+                gravity_u(*coords, *sample_point_source, sample_mass)
+                for coords in zip(*coords_2)
+            )
+        )
+        npt.assert_allclose(g_u_1, -g_u_2)
+
+
+class TestAccelerationFiniteDifferences:
+    """
+    Test acceleration components against finite-differences approximations of
+    the gravity potential
+    """
+
+    # Define percentage for the finite difference displacement
+    delta_percentage = 1e-8
+
+    # Define expected relative error tolerance in the comparisons
+    rtol = 1e-5
+
+    @pytest.fixture
+    def finite_diff_gravity_e(
+        self, sample_coordinate, sample_point_source, sample_mass
+    ):
+        """
+        Compute gravity_e through finite differences of the gravity_pot
+        """
+        easting_p, northing_p, upward_p = sample_coordinate
+        easting_q, _, _ = sample_point_source
+        # Compute a small increment in the easting coordinate
+        d_easting = self.delta_percentage * (easting_p - easting_q)
+        # Compute shifted coordinate
+        shifted_coordinate = (easting_p + d_easting, northing_p, upward_p)
+        # Calculate g_e through finite differences
+        g_e = (
+            gravity_pot(*shifted_coordinate, *sample_point_source, sample_mass)
+            - gravity_pot(*sample_coordinate, *sample_point_source, sample_mass)
+        ) / d_easting
+        return g_e
+
+    @pytest.fixture
+    def finite_diff_gravity_n(
+        self, sample_coordinate, sample_point_source, sample_mass
+    ):
+        """
+        Compute gravity_n through finite differences of the gravity_pot
+        """
+        easting_p, northing_p, upward_p = sample_coordinate
+        _, northing_q, _ = sample_point_source
+        # Compute a small increment in the easting coordinate
+        d_northing = self.delta_percentage * (northing_p - northing_q)
+        # Compute shifted coordinate
+        shifted_coordinate = (easting_p, northing_p + d_northing, upward_p)
+        # Calculate g_n through finite differences
+        g_n = (
+            gravity_pot(*shifted_coordinate, *sample_point_source, sample_mass)
+            - gravity_pot(*sample_coordinate, *sample_point_source, sample_mass)
+        ) / d_northing
+        return g_n
+
+    @pytest.fixture
+    def finite_diff_gravity_u(
+        self, sample_coordinate, sample_point_source, sample_mass
+    ):
+        """
+        Compute gravity_u through finite differences of the gravity_pot
+        """
+        easting_p, northing_p, upward_p = sample_coordinate
+        _, _, upward_q = sample_point_source
+        # Compute a small increment in the easting coordinate
+        d_upward = self.delta_percentage * (upward_p - upward_q)
+        # Compute shifted coordinate
+        shifted_coordinate = (easting_p, northing_p, upward_p + d_upward)
+        # Calculate g_u through finite differences
+        g_u = (
+            gravity_pot(*shifted_coordinate, *sample_point_source, sample_mass)
+            - gravity_pot(*sample_coordinate, *sample_point_source, sample_mass)
+        ) / d_upward
+        return g_u
+
+    def test_gravity_e(
+        self, sample_coordinate, sample_point_source, finite_diff_gravity_e, sample_mass
+    ):
+        """
+        Test gravity_e against finite differences of the gravity_pot
+        """
+        npt.assert_allclose(
+            finite_diff_gravity_e,
+            gravity_e(*sample_coordinate, *sample_point_source, sample_mass),
+            rtol=self.rtol,
+        )
+
+    def test_gravity_n(
+        self, sample_coordinate, sample_point_source, finite_diff_gravity_n, sample_mass
+    ):
+        """
+        Test gravity_n against finite differences of the gravity_pot
+        """
+        npt.assert_allclose(
+            finite_diff_gravity_n,
+            gravity_n(*sample_coordinate, *sample_point_source, sample_mass),
+            rtol=self.rtol,
+        )
+
+    def test_gravity_u(
+        self, sample_coordinate, sample_point_source, finite_diff_gravity_u, sample_mass
+    ):
+        """
+        Test gravity_u against finite differences of the gravity_pot
+        """
+        npt.assert_allclose(
+            finite_diff_gravity_u,
+            gravity_u(*sample_coordinate, *sample_point_source, sample_mass),
+            rtol=self.rtol,
+        )
+
+
+class TestTensorFiniteDifferences:
+    """
+    Test tensor gravity functions against finite-differences approximations of
+    the gradient components
+    """
+
+    # Define percentage for the finite difference displacement
+    delta_percentage = 1e-8
+
+    # Define expected relative error tolerance in the comparisons
+    rtol = 1e-5
+
+    @pytest.fixture
+    def finite_diff_gravity_ee(
+        self, sample_coordinate, sample_point_source, sample_mass
+    ):
+        """
+        Compute gravity_ee through finite differences of the gravity_e
+        """
+        easting_p, northing_p, upward_p = sample_coordinate
+        easting_q, _, _ = sample_point_source
+        # Compute a small increment in the easting coordinate
+        d_easting = self.delta_percentage * (easting_p - easting_q)
+        # Compute shifted coordinate
+        shifted_coordinate = (easting_p + d_easting, northing_p, upward_p)
+        # Calculate g_ee through finite differences
+        g_ee = (
+            gravity_e(*shifted_coordinate, *sample_point_source, sample_mass)
+            - gravity_e(*sample_coordinate, *sample_point_source, sample_mass)
+        ) / d_easting
+        return g_ee
+
+    @pytest.fixture
+    def finite_diff_gravity_nn(
+        self, sample_coordinate, sample_point_source, sample_mass
+    ):
+        """
+        Compute gravity_nn through finite differences of the gravity_n
+        """
+        easting_p, northing_p, upward_p = sample_coordinate
+        _, northing_q, _ = sample_point_source
+        # Compute a small increment in the easting coordinate
+        d_northing = self.delta_percentage * (northing_p - northing_q)
+        # Compute shifted coordinate
+        shifted_coordinate = (easting_p, northing_p + d_northing, upward_p)
+        # Calculate g_nn through finite differences
+        g_nn = (
+            gravity_n(*shifted_coordinate, *sample_point_source, sample_mass)
+            - gravity_n(*sample_coordinate, *sample_point_source, sample_mass)
+        ) / d_northing
+        return g_nn
+
+    @pytest.fixture
+    def finite_diff_gravity_uu(
+        self, sample_coordinate, sample_point_source, sample_mass
+    ):
+        """
+        Compute gravity_uu through finite differences of the gravity_u
+        """
+        easting_p, northing_p, upward_p = sample_coordinate
+        _, _, upward_q = sample_point_source
+        # Compute a small increment in the easting coordinate
+        d_upward = self.delta_percentage * (upward_p - upward_q)
+        # Compute shifted coordinate
+        shifted_coordinate = (easting_p, northing_p, upward_p + d_upward)
+        # Calculate g_uu through finite differences
+        g_uu = (
+            gravity_u(*shifted_coordinate, *sample_point_source, sample_mass)
+            - gravity_u(*sample_coordinate, *sample_point_source, sample_mass)
+        ) / d_upward
+        return g_uu
+
+    @pytest.fixture
+    def finite_diff_gravity_en(
+        self, sample_coordinate, sample_point_source, sample_mass
+    ):
+        """
+        Compute gravity_en through finite differences of the gravity_e
+        """
+        easting_p, northing_p, upward_p = sample_coordinate
+        _, northing_q, _ = sample_point_source
+        # Compute a small increment in the easting coordinate
+        d_northing = self.delta_percentage * (northing_p - northing_q)
+        # Compute shifted coordinate
+        shifted_coordinate = (easting_p, northing_p + d_northing, upward_p)
+        # Calculate g_en through finite differences
+        g_en = (
+            gravity_e(*shifted_coordinate, *sample_point_source, sample_mass)
+            - gravity_e(*sample_coordinate, *sample_point_source, sample_mass)
+        ) / d_northing
+        return g_en
+
+    @pytest.fixture
+    def finite_diff_gravity_eu(
+        self, sample_coordinate, sample_point_source, sample_mass
+    ):
+        """
+        Compute gravity_eu through finite differences of the gravity_e
+        """
+        easting_p, northing_p, upward_p = sample_coordinate
+        _, _, upward_q = sample_point_source
+        # Compute a small increment in the easting coordinate
+        d_upward = self.delta_percentage * (upward_p - upward_q)
+        # Compute shifted coordinate
+        shifted_coordinate = (easting_p, northing_p, upward_p + d_upward)
+        # Calculate g_eu through finite differences
+        g_eu = (
+            gravity_e(*shifted_coordinate, *sample_point_source, sample_mass)
+            - gravity_e(*sample_coordinate, *sample_point_source, sample_mass)
+        ) / d_upward
+        return g_eu
+
+    @pytest.fixture
+    def finite_diff_gravity_nu(
+        self, sample_coordinate, sample_point_source, sample_mass
+    ):
+        """
+        Test gravity_nu against finite differences of the gravity_n
+        """
+        easting_p, northing_p, upward_p = sample_coordinate
+        _, _, upward_q = sample_point_source
+        # Compute a small increment in the easting coordinate
+        d_upward = self.delta_percentage * (upward_p - upward_q)
+        # Compute shifted coordinate
+        shifted_coordinate = (easting_p, northing_p, upward_p + d_upward)
+        # Calculate g_nu through finite differences
+        g_nu = (
+            gravity_n(*shifted_coordinate, *sample_point_source, sample_mass)
+            - gravity_n(*sample_coordinate, *sample_point_source, sample_mass)
+        ) / d_upward
+        return g_nu
+
+    def test_gravity_ee(
+        self,
+        sample_coordinate,
+        sample_point_source,
+        finite_diff_gravity_ee,
+        sample_mass,
+    ):
+        """
+        Test gravity_ee against finite differences of the gravity_e
+        """
+        npt.assert_allclose(
+            finite_diff_gravity_ee,
+            gravity_ee(*sample_coordinate, *sample_point_source, sample_mass),
+            rtol=self.rtol,
+        )
+
+    def test_gravity_nn(
+        self,
+        sample_coordinate,
+        sample_point_source,
+        finite_diff_gravity_nn,
+        sample_mass,
+    ):
+        """
+        Test gravity_nn against finite differences of the gravity_n
+        """
+        npt.assert_allclose(
+            finite_diff_gravity_nn,
+            gravity_nn(*sample_coordinate, *sample_point_source, sample_mass),
+            rtol=self.rtol,
+        )
+
+    def test_gravity_uu(
+        self,
+        sample_coordinate,
+        sample_point_source,
+        finite_diff_gravity_uu,
+        sample_mass,
+    ):
+        """
+        Test gravity_uu against finite differences of the gravity_u
+        """
+        npt.assert_allclose(
+            finite_diff_gravity_uu,
+            gravity_uu(*sample_coordinate, *sample_point_source, sample_mass),
+            rtol=self.rtol,
+        )
+
+    def test_gravity_en(
+        self,
+        sample_coordinate,
+        sample_point_source,
+        finite_diff_gravity_en,
+        sample_mass,
+    ):
+        """
+        Test gravity_en against finite differences of the gravity_e
+        """
+        npt.assert_allclose(
+            finite_diff_gravity_en,
+            gravity_en(*sample_coordinate, *sample_point_source, sample_mass),
+            rtol=self.rtol,
+        )
+
+    def test_gravity_eu(
+        self,
+        sample_coordinate,
+        sample_point_source,
+        finite_diff_gravity_eu,
+        sample_mass,
+    ):
+        """
+        Test gravity_eu against finite differences of the gravity_e
+        """
+        npt.assert_allclose(
+            finite_diff_gravity_eu,
+            gravity_eu(*sample_coordinate, *sample_point_source, sample_mass),
+            rtol=self.rtol,
+        )
+
+    def test_gravity_nu(
+        self,
+        sample_coordinate,
+        sample_point_source,
+        finite_diff_gravity_nu,
+        sample_mass,
+    ):
+        """
+        Test gravity_nu against finite differences of the gravity_n
+        """
+        npt.assert_allclose(
+            finite_diff_gravity_nu,
+            gravity_nu(*sample_coordinate, *sample_point_source, sample_mass),
+            rtol=self.rtol,
+        )
+
+
+def test_laplacian(sample_coordinate, sample_point_source, sample_mass):
+    """
+    Test if diagonal tensor functions satisfy Laplace's equation
+    """
+    g_ee = gravity_ee(*sample_coordinate, *sample_point_source, sample_mass)
+    g_nn = gravity_nn(*sample_coordinate, *sample_point_source, sample_mass)
+    g_zz = gravity_uu(*sample_coordinate, *sample_point_source, sample_mass)
+    npt.assert_allclose(-g_zz, g_ee + g_nn)
diff --git a/choclo/tests/test_point.py b/choclo/tests/test_point_kernels.py
similarity index 55%
rename from choclo/tests/test_point.py
rename to choclo/tests/test_point_kernels.py
index 994e0fa2..1554a77d 100644
--- a/choclo/tests/test_point.py
+++ b/choclo/tests/test_point_kernels.py
@@ -11,21 +11,43 @@
 import numpy.testing as npt
 import pytest
 
-from .._point import (
-    kernel_point_g_easting,
-    kernel_point_g_ee,
-    kernel_point_g_en,
-    kernel_point_g_ez,
-    kernel_point_g_nn,
-    kernel_point_g_northing,
-    kernel_point_g_nz,
-    kernel_point_g_upward,
-    kernel_point_g_zz,
-    kernel_point_potential,
+from ..point import (
+    kernel_e,
+    kernel_ee,
+    kernel_en,
+    kernel_eu,
+    kernel_n,
+    kernel_nn,
+    kernel_nu,
+    kernel_pot,
+    kernel_u,
+    kernel_uu,
 )
+from ..utils import distance_cartesian
 from .utils import NUMBA_IS_DISABLED
 
 
+def evaluate_kernel(coordinates, source, kernel):
+    """
+    Evaluate kernel on a set of observation points using a single source
+
+    Parameters
+    ----------
+    coordinates : list
+        Coordinates of the observation points in the following order:
+            ``easting``, ``northing``, ``upward``
+    source : tuple
+        Tuple containing the coordinates of the point source
+    kernel : func
+        Kernel function that wants to be evaluated
+    """
+    result = []
+    for coords in zip(*coordinates):
+        distance = distance_cartesian(*coords, *source)
+        result.append(kernel(*coords, *source, distance))
+    return np.array(result)
+
+
 @pytest.fixture(name="sample_point_source")
 def fixture_sample_point_source():
     """
@@ -44,7 +66,7 @@ def fixture_sample_coordinate():
 
 class TestSymmetryPotential:
     """
-    Test the symmetry of the kernel for the potential of a point source
+    Test the symmetry of kernel_pot
     """
 
     def build_points_in_sphere(self, radius, point):
@@ -65,49 +87,39 @@ def build_points_in_sphere(self, radius, point):
 
     def test_symmetry_on_sphere(self, sample_point_source):
         """
-        Test the symmetry of the potential in points of a sphere
+        Test the symmetry of kernel_pot in points of a sphere
         """
         radius = 3.5
         observation_points = self.build_points_in_sphere(radius, sample_point_source)
-        kernel_potential = [
-            kernel_point_potential(*coords, *sample_point_source)
-            for coords in zip(*observation_points)
-        ]
-        npt.assert_allclose(kernel_potential[0], kernel_potential)
+        # Evaluate kernel_pot on every observation point
+        kernel = evaluate_kernel(observation_points, sample_point_source, kernel_pot)
+        # Check if all values are equal
+        npt.assert_allclose(kernel[0], kernel)
 
     def test_potential_between_two_spheres(self, sample_point_source):
         """
-        Test the potential between observation points in two spheres
+        Test kernel_pot between observation points in two spheres
         """
         radius_1, radius_2 = 3.5, 8.7
         sphere_1 = self.build_points_in_sphere(radius_1, sample_point_source)
         sphere_2 = self.build_points_in_sphere(radius_2, sample_point_source)
-        kernel_potential_1 = np.array(
-            [
-                kernel_point_potential(*coords, *sample_point_source)
-                for coords in zip(*sphere_1)
-            ]
-        )
-        kernel_potential_2 = np.array(
-            [
-                kernel_point_potential(*coords, *sample_point_source)
-                for coords in zip(*sphere_2)
-            ]
-        )
-        npt.assert_allclose(
-            kernel_potential_1, radius_2 / radius_1 * kernel_potential_2
-        )
+        # Evaluate kernel_pot on every observation point
+        kernel_1 = evaluate_kernel(sphere_1, sample_point_source, kernel_pot)
+        kernel_2 = evaluate_kernel(sphere_2, sample_point_source, kernel_pot)
+        # Check if all values are equal
+        npt.assert_allclose(kernel_1, radius_2 / radius_1 * kernel_2)
 
     @pytest.mark.skipif(not NUMBA_IS_DISABLED, reason="Numba is not disabled")
     def test_infinite_potential(self):
         """
-        Test if we get an infinite kernel if the computation point is in the
+        Test if we get an infinite kernel_pot the computation point is in the
         same location as the observation point.
 
         This test should be run only with Numba disabled.
         """
         point = (4.6, -8.9, -50.3)
-        assert np.isinf(kernel_point_potential(*point, *point))
+        distance = distance_cartesian(*point, *point)
+        assert np.isinf(kernel_pot(*point, *point, distance))
 
     @pytest.mark.skipif(NUMBA_IS_DISABLED, reason="Numba is disabled")
     def test_division_by_zero(self):
@@ -118,13 +130,14 @@ def test_division_by_zero(self):
         This test should be run only with Numba enabled.
         """
         point = (4.6, -8.9, -50.3)
+        distance = distance_cartesian(*point, *point)
         with pytest.raises(ZeroDivisionError):
-            kernel_point_potential(*point, *point)
+            kernel_pot(*point, *point, distance)
 
 
-class TestSymmetryGradientEasting:
+class TestSymmetryKernelE:
     """
-    Test symmetry of the kernel for the g_easting component of a point source
+    Test symmetry of the kernel_e
     """
 
     @pytest.fixture
@@ -169,39 +182,26 @@ def test_zero_in_northing_upward(
         self, sample_point_source, coords_northing_upward_plane
     ):
         """
-        Test if g_easting is zero in the northing-upward plane
+        Test if kernel_e is zero in the northing-upward plane
         """
-        g_easting = np.array(
-            [
-                kernel_point_g_easting(easting, northing, upward, *sample_point_source)
-                for easting, northing, upward in zip(*coords_northing_upward_plane)
-            ]
+        kernel = evaluate_kernel(
+            coords_northing_upward_plane, sample_point_source, kernel_e
         )
-        assert (g_easting <= 1e-30).all()
+        assert (kernel <= 1e-30).all()
 
     def test_mirror_symmetry(self, sample_point_source, mirrored_points):
         """
         Test points with opposite easting coordinate
         """
         coords_1, coords_2 = mirrored_points
-        g_easting_1 = np.array(
-            [
-                kernel_point_g_easting(easting, northing, upward, *sample_point_source)
-                for easting, northing, upward in zip(*coords_1)
-            ]
-        )
-        g_easting_2 = np.array(
-            [
-                kernel_point_g_easting(easting, northing, upward, *sample_point_source)
-                for easting, northing, upward in zip(*coords_2)
-            ]
-        )
-        npt.assert_allclose(g_easting_1, -g_easting_2)
+        kernel_1 = evaluate_kernel(coords_1, sample_point_source, kernel_e)
+        kernel_2 = evaluate_kernel(coords_2, sample_point_source, kernel_e)
+        npt.assert_allclose(kernel_1, -kernel_2)
 
 
-class TestSymmetryGradientNorthing:
+class TestSymmetryKernelN:
     """
-    Test symmetry of the kernel for the g_northing component of a point source
+    Test symmetry of the kernel_n
     """
 
     @pytest.fixture
@@ -246,39 +246,26 @@ def test_zero_in_easting_upward(
         self, sample_point_source, coords_easting_upward_plane
     ):
         """
-        Test if g_northing is zero in the easting-upward plane
+        Test if kernel_n is zero in the easting-upward plane
         """
-        g_northing = np.array(
-            [
-                kernel_point_g_northing(easting, northing, upward, *sample_point_source)
-                for easting, northing, upward in zip(*coords_easting_upward_plane)
-            ]
+        kernel = evaluate_kernel(
+            coords_easting_upward_plane, sample_point_source, kernel_n
         )
-        assert (g_northing <= 1e-30).all()
+        assert (kernel <= 1e-30).all()
 
     def test_mirror_symmetry(self, sample_point_source, mirrored_points):
         """
         Test points with opposite northing coordinate
         """
         coords_1, coords_2 = mirrored_points
-        g_northing_1 = np.array(
-            [
-                kernel_point_g_northing(easting, northing, upward, *sample_point_source)
-                for easting, northing, upward in zip(*coords_1)
-            ]
-        )
-        g_northing_2 = np.array(
-            [
-                kernel_point_g_northing(easting, northing, upward, *sample_point_source)
-                for easting, northing, upward in zip(*coords_2)
-            ]
-        )
-        npt.assert_allclose(g_northing_1, -g_northing_2)
+        kernel_1 = evaluate_kernel(coords_1, sample_point_source, kernel_n)
+        kernel_2 = evaluate_kernel(coords_2, sample_point_source, kernel_n)
+        npt.assert_allclose(kernel_1, -kernel_2)
 
 
-class TestSymmetryGradientUpward:
+class TestSymmetryKernelU:
     """
-    Test symmetry of the kernel for the g_upward component of a point source
+    Test symmetry of the kernel_u
     """
 
     @pytest.fixture
@@ -323,34 +310,21 @@ def test_zero_in_easting_northing(
         self, sample_point_source, coords_easting_northing_plane
     ):
         """
-        Test if g_upward is zero in the easting-upward plane
+        Test if kernel_u is zero in the easting-upward plane
         """
-        g_upward = np.array(
-            [
-                kernel_point_g_upward(easting, northing, upward, *sample_point_source)
-                for easting, northing, upward in zip(*coords_easting_northing_plane)
-            ]
+        kernel = evaluate_kernel(
+            coords_easting_northing_plane, sample_point_source, kernel_u
         )
-        assert (g_upward <= 1e-30).all()
+        assert (kernel <= 1e-30).all()
 
     def test_mirror_symmetry(self, sample_point_source, mirrored_points):
         """
         Test points with opposite upward coordinate
         """
         coords_1, coords_2 = mirrored_points
-        g_upward_1 = np.array(
-            [
-                kernel_point_g_upward(easting, northing, upward, *sample_point_source)
-                for easting, northing, upward in zip(*coords_1)
-            ]
-        )
-        g_upward_2 = np.array(
-            [
-                kernel_point_g_upward(easting, northing, upward, *sample_point_source)
-                for easting, northing, upward in zip(*coords_2)
-            ]
-        )
-        npt.assert_allclose(g_upward_1, -g_upward_2)
+        kernel_1 = evaluate_kernel(coords_1, sample_point_source, kernel_u)
+        kernel_2 = evaluate_kernel(coords_2, sample_point_source, kernel_u)
+        npt.assert_allclose(kernel_1, -kernel_2)
 
 
 class TestGradientFiniteDifferences:
@@ -366,9 +340,9 @@ class TestGradientFiniteDifferences:
     rtol = 1e-5
 
     @pytest.fixture
-    def finite_diff_g_easting(self, sample_coordinate, sample_point_source):
+    def finite_diff_kernel_e(self, sample_coordinate, sample_point_source):
         """
-        Compute g_easting through finite differences of the potential
+        Compute kernel_e through finite differences of the kernel_pot
         """
         easting_p, northing_p, upward_p = sample_coordinate
         easting_q, _, _ = sample_point_source
@@ -376,17 +350,19 @@ def finite_diff_g_easting(self, sample_coordinate, sample_point_source):
         d_easting = self.delta_percentage * (easting_p - easting_q)
         # Compute shifted coordinate
         shifted_coordinate = (easting_p + d_easting, northing_p, upward_p)
-        # Calculate g_easting through finite differences
-        g_easting = (
-            kernel_point_potential(*shifted_coordinate, *sample_point_source)
-            - kernel_point_potential(*sample_coordinate, *sample_point_source)
+        # Calculate kernel_e through finite differences
+        distance_shifted = distance_cartesian(*shifted_coordinate, *sample_point_source)
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
+        k_easting = (
+            kernel_pot(*shifted_coordinate, *sample_point_source, distance_shifted)
+            - kernel_pot(*sample_coordinate, *sample_point_source, distance)
         ) / d_easting
-        return g_easting
+        return k_easting
 
     @pytest.fixture
-    def finite_diff_g_northing(self, sample_coordinate, sample_point_source):
+    def finite_diff_kernel_n(self, sample_coordinate, sample_point_source):
         """
-        Compute g_northing through finite differences of the potential
+        Compute kernel_n through finite differences of the kernel_pot
         """
         easting_p, northing_p, upward_p = sample_coordinate
         _, northing_q, _ = sample_point_source
@@ -394,17 +370,19 @@ def finite_diff_g_northing(self, sample_coordinate, sample_point_source):
         d_northing = self.delta_percentage * (northing_p - northing_q)
         # Compute shifted coordinate
         shifted_coordinate = (easting_p, northing_p + d_northing, upward_p)
-        # Calculate g_easting through finite differences
-        g_northing = (
-            kernel_point_potential(*shifted_coordinate, *sample_point_source)
-            - kernel_point_potential(*sample_coordinate, *sample_point_source)
+        # Calculate kernel_n through finite differences
+        distance_shifted = distance_cartesian(*shifted_coordinate, *sample_point_source)
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
+        k_northing = (
+            kernel_pot(*shifted_coordinate, *sample_point_source, distance_shifted)
+            - kernel_pot(*sample_coordinate, *sample_point_source, distance)
         ) / d_northing
-        return g_northing
+        return k_northing
 
     @pytest.fixture
-    def finite_diff_g_upward(self, sample_coordinate, sample_point_source):
+    def finite_diff_kernel_u(self, sample_coordinate, sample_point_source):
         """
-        Compute g_upward through finite differences of the potential
+        Compute kernel_u through finite differences of the kernel_pot
         """
         easting_p, northing_p, upward_p = sample_coordinate
         _, _, upward_q = sample_point_source
@@ -412,46 +390,51 @@ def finite_diff_g_upward(self, sample_coordinate, sample_point_source):
         d_upward = self.delta_percentage * (upward_p - upward_q)
         # Compute shifted coordinate
         shifted_coordinate = (easting_p, northing_p, upward_p + d_upward)
-        # Calculate g_easting through finite differences
-        g_upward = (
-            kernel_point_potential(*shifted_coordinate, *sample_point_source)
-            - kernel_point_potential(*sample_coordinate, *sample_point_source)
+        # Calculate kernel_u through finite differences
+        distance_shifted = distance_cartesian(*shifted_coordinate, *sample_point_source)
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
+        k_upward = (
+            kernel_pot(*shifted_coordinate, *sample_point_source, distance_shifted)
+            - kernel_pot(*sample_coordinate, *sample_point_source, distance)
         ) / d_upward
-        return g_upward
+        return k_upward
 
-    def test_g_easting(
-        self, sample_coordinate, sample_point_source, finite_diff_g_easting
+    def test_kernel_e(
+        self, sample_coordinate, sample_point_source, finite_diff_kernel_e
     ):
         """
-        Test kernel of g_easting against finite differences of the potential
+        Test kernel_e against finite differences of the kernel_pot
         """
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
         npt.assert_allclose(
-            finite_diff_g_easting,
-            kernel_point_g_easting(*sample_coordinate, *sample_point_source),
+            finite_diff_kernel_e,
+            kernel_e(*sample_coordinate, *sample_point_source, distance),
             rtol=self.rtol,
         )
 
-    def test_g_northing(
-        self, sample_coordinate, sample_point_source, finite_diff_g_northing
+    def test_kernel_n(
+        self, sample_coordinate, sample_point_source, finite_diff_kernel_n
     ):
         """
-        Test kernel of g_northing against finite differences of the potential
+        Test kernel_n against finite differences of the kernel_pot
         """
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
         npt.assert_allclose(
-            finite_diff_g_northing,
-            kernel_point_g_northing(*sample_coordinate, *sample_point_source),
+            finite_diff_kernel_n,
+            kernel_n(*sample_coordinate, *sample_point_source, distance),
             rtol=self.rtol,
         )
 
-    def test_g_upward(
-        self, sample_coordinate, sample_point_source, finite_diff_g_upward
+    def test_kernel_u(
+        self, sample_coordinate, sample_point_source, finite_diff_kernel_u
     ):
         """
-        Test kernel of g_upward against finite differences of the potential
+        Test kernel_u against finite differences of the kernel_pot
         """
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
         npt.assert_allclose(
-            finite_diff_g_upward,
-            kernel_point_g_upward(*sample_coordinate, *sample_point_source),
+            finite_diff_kernel_u,
+            kernel_u(*sample_coordinate, *sample_point_source, distance),
             rtol=self.rtol,
         )
 
@@ -469,9 +452,9 @@ class TestTensorFiniteDifferences:
     rtol = 1e-5
 
     @pytest.fixture
-    def finite_diff_g_ee(self, sample_coordinate, sample_point_source):
+    def finite_diff_kernel_ee(self, sample_coordinate, sample_point_source):
         """
-        Compute g_ee through finite differences of the g_easting
+        Compute kernel_ee through finite differences of the kernel_e
         """
         easting_p, northing_p, upward_p = sample_coordinate
         easting_q, _, _ = sample_point_source
@@ -479,17 +462,19 @@ def finite_diff_g_ee(self, sample_coordinate, sample_point_source):
         d_easting = self.delta_percentage * (easting_p - easting_q)
         # Compute shifted coordinate
         shifted_coordinate = (easting_p + d_easting, northing_p, upward_p)
-        # Calculate g_easting through finite differences
-        g_ee = (
-            kernel_point_g_easting(*shifted_coordinate, *sample_point_source)
-            - kernel_point_g_easting(*sample_coordinate, *sample_point_source)
+        # Calculate kernel_ee through finite differences
+        distance_shifted = distance_cartesian(*shifted_coordinate, *sample_point_source)
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
+        kernel_ee = (
+            kernel_e(*shifted_coordinate, *sample_point_source, distance_shifted)
+            - kernel_e(*sample_coordinate, *sample_point_source, distance)
         ) / d_easting
-        return g_ee
+        return kernel_ee
 
     @pytest.fixture
-    def finite_diff_g_nn(self, sample_coordinate, sample_point_source):
+    def finite_diff_kernel_nn(self, sample_coordinate, sample_point_source):
         """
-        Compute g_nn through finite differences of the g_northing
+        Compute kernel_nn through finite differences of the kernel_n
         """
         easting_p, northing_p, upward_p = sample_coordinate
         _, northing_q, _ = sample_point_source
@@ -497,17 +482,19 @@ def finite_diff_g_nn(self, sample_coordinate, sample_point_source):
         d_northing = self.delta_percentage * (northing_p - northing_q)
         # Compute shifted coordinate
         shifted_coordinate = (easting_p, northing_p + d_northing, upward_p)
-        # Calculate g_easting through finite differences
-        g_nn = (
-            kernel_point_g_northing(*shifted_coordinate, *sample_point_source)
-            - kernel_point_g_northing(*sample_coordinate, *sample_point_source)
+        # Calculate kernel_nn through finite differences
+        distance_shifted = distance_cartesian(*shifted_coordinate, *sample_point_source)
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
+        kernel_nn = (
+            kernel_n(*shifted_coordinate, *sample_point_source, distance_shifted)
+            - kernel_n(*sample_coordinate, *sample_point_source, distance)
         ) / d_northing
-        return g_nn
+        return kernel_nn
 
     @pytest.fixture
-    def finite_diff_g_zz(self, sample_coordinate, sample_point_source):
+    def finite_diff_kernel_uu(self, sample_coordinate, sample_point_source):
         """
-        Compute g_zz through finite differences of the g_upward
+        Compute kernel_uu through finite differences of the kernel_u
         """
         easting_p, northing_p, upward_p = sample_coordinate
         _, _, upward_q = sample_point_source
@@ -515,17 +502,19 @@ def finite_diff_g_zz(self, sample_coordinate, sample_point_source):
         d_upward = self.delta_percentage * (upward_p - upward_q)
         # Compute shifted coordinate
         shifted_coordinate = (easting_p, northing_p, upward_p + d_upward)
-        # Calculate g_easting through finite differences
-        g_zz = (
-            kernel_point_g_upward(*shifted_coordinate, *sample_point_source)
-            - kernel_point_g_upward(*sample_coordinate, *sample_point_source)
+        # Calculate kernel_uu through finite differences
+        distance_shifted = distance_cartesian(*shifted_coordinate, *sample_point_source)
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
+        k_uu = (
+            kernel_u(*shifted_coordinate, *sample_point_source, distance_shifted)
+            - kernel_u(*sample_coordinate, *sample_point_source, distance)
         ) / d_upward
-        return g_zz
+        return k_uu
 
     @pytest.fixture
-    def finite_diff_g_en(self, sample_coordinate, sample_point_source):
+    def finite_diff_kernel_en(self, sample_coordinate, sample_point_source):
         """
-        Compute g_en through finite differences of the g_easting
+        Compute kernel_en through finite differences of the kernel_e
         """
         easting_p, northing_p, upward_p = sample_coordinate
         _, northing_q, _ = sample_point_source
@@ -533,17 +522,19 @@ def finite_diff_g_en(self, sample_coordinate, sample_point_source):
         d_northing = self.delta_percentage * (northing_p - northing_q)
         # Compute shifted coordinate
         shifted_coordinate = (easting_p, northing_p + d_northing, upward_p)
-        # Calculate g_easting through finite differences
-        g_en = (
-            kernel_point_g_easting(*shifted_coordinate, *sample_point_source)
-            - kernel_point_g_easting(*sample_coordinate, *sample_point_source)
+        # Calculate kernel_en through finite differences
+        distance_shifted = distance_cartesian(*shifted_coordinate, *sample_point_source)
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
+        k_en = (
+            kernel_e(*shifted_coordinate, *sample_point_source, distance_shifted)
+            - kernel_e(*sample_coordinate, *sample_point_source, distance)
         ) / d_northing
-        return g_en
+        return k_en
 
     @pytest.fixture
-    def finite_diff_g_ez(self, sample_coordinate, sample_point_source):
+    def finite_diff_kernel_eu(self, sample_coordinate, sample_point_source):
         """
-        Compute g_ez through finite differences of the g_easting
+        Compute kernel_eu through finite differences of the kernel_e
         """
         easting_p, northing_p, upward_p = sample_coordinate
         _, _, upward_q = sample_point_source
@@ -551,17 +542,19 @@ def finite_diff_g_ez(self, sample_coordinate, sample_point_source):
         d_upward = self.delta_percentage * (upward_p - upward_q)
         # Compute shifted coordinate
         shifted_coordinate = (easting_p, northing_p, upward_p + d_upward)
-        # Calculate g_easting through finite differences
-        g_ez = (
-            kernel_point_g_easting(*shifted_coordinate, *sample_point_source)
-            - kernel_point_g_easting(*sample_coordinate, *sample_point_source)
+        # Calculate kernel_eu through finite differences
+        distance_shifted = distance_cartesian(*shifted_coordinate, *sample_point_source)
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
+        k_eu = (
+            kernel_e(*shifted_coordinate, *sample_point_source, distance_shifted)
+            - kernel_e(*sample_coordinate, *sample_point_source, distance)
         ) / d_upward
-        return g_ez
+        return k_eu
 
     @pytest.fixture
-    def finite_diff_g_nz(self, sample_coordinate, sample_point_source):
+    def finite_diff_kernel_nu(self, sample_coordinate, sample_point_source):
         """
-        Compute g_nz through finite differences of the g_northing
+        Test kernel_nu against finite differences of the kernel_n
         """
         easting_p, northing_p, upward_p = sample_coordinate
         _, _, upward_q = sample_point_source
@@ -569,79 +562,88 @@ def finite_diff_g_nz(self, sample_coordinate, sample_point_source):
         d_upward = self.delta_percentage * (upward_p - upward_q)
         # Compute shifted coordinate
         shifted_coordinate = (easting_p, northing_p, upward_p + d_upward)
-        # Calculate g_easting through finite differences
-        g_nz = (
-            kernel_point_g_northing(*shifted_coordinate, *sample_point_source)
-            - kernel_point_g_northing(*sample_coordinate, *sample_point_source)
+        # Calculate k_nu through finite differences
+        distance_shifted = distance_cartesian(*shifted_coordinate, *sample_point_source)
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
+        k_nu = (
+            kernel_n(*shifted_coordinate, *sample_point_source, distance_shifted)
+            - kernel_n(*sample_coordinate, *sample_point_source, distance)
         ) / d_upward
-        return g_nz
+        return k_nu
 
-    def test_g_ee(self, sample_coordinate, sample_point_source, finite_diff_g_ee):
+    def test_g_ee(self, sample_coordinate, sample_point_source, finite_diff_kernel_ee):
         """
-        Test kernel of g_ee against finite differences of the g_easting
+        Test kernel_ee against finite differences of the kernel_e
         """
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
         npt.assert_allclose(
-            finite_diff_g_ee,
-            kernel_point_g_ee(*sample_coordinate, *sample_point_source),
+            finite_diff_kernel_ee,
+            kernel_ee(*sample_coordinate, *sample_point_source, distance),
             rtol=self.rtol,
         )
 
-    def test_g_nn(self, sample_coordinate, sample_point_source, finite_diff_g_nn):
+    def test_g_nn(self, sample_coordinate, sample_point_source, finite_diff_kernel_nn):
         """
-        Test kernel of g_nn against finite differences of the g_northing
+        Test kernel_nn against finite differences of the kernel_n
         """
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
         npt.assert_allclose(
-            finite_diff_g_nn,
-            kernel_point_g_nn(*sample_coordinate, *sample_point_source),
+            finite_diff_kernel_nn,
+            kernel_nn(*sample_coordinate, *sample_point_source, distance),
             rtol=self.rtol,
         )
 
-    def test_g_zz(self, sample_coordinate, sample_point_source, finite_diff_g_zz):
+    def test_g_zz(self, sample_coordinate, sample_point_source, finite_diff_kernel_uu):
         """
-        Test kernel of g_zz against finite differences of the g_upward
+        Test kernel_uu against finite differences of the kernel_u
         """
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
         npt.assert_allclose(
-            finite_diff_g_zz,
-            kernel_point_g_zz(*sample_coordinate, *sample_point_source),
+            finite_diff_kernel_uu,
+            kernel_uu(*sample_coordinate, *sample_point_source, distance),
             rtol=self.rtol,
         )
 
-    def test_g_en(self, sample_coordinate, sample_point_source, finite_diff_g_en):
+    def test_g_en(self, sample_coordinate, sample_point_source, finite_diff_kernel_en):
         """
-        Test kernel of g_en against finite differences of the g_easting
+        Test kernel_en against finite differences of the kernel_e
         """
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
         npt.assert_allclose(
-            finite_diff_g_en,
-            kernel_point_g_en(*sample_coordinate, *sample_point_source),
+            finite_diff_kernel_en,
+            kernel_en(*sample_coordinate, *sample_point_source, distance),
             rtol=self.rtol,
         )
 
-    def test_g_ez(self, sample_coordinate, sample_point_source, finite_diff_g_ez):
+    def test_g_ez(self, sample_coordinate, sample_point_source, finite_diff_kernel_eu):
         """
-        Test kernel of g_ez against finite differences of the g_easting
+        Test kernel_eu against finite differences of the kernel_e
         """
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
         npt.assert_allclose(
-            finite_diff_g_ez,
-            kernel_point_g_ez(*sample_coordinate, *sample_point_source),
+            finite_diff_kernel_eu,
+            kernel_eu(*sample_coordinate, *sample_point_source, distance),
             rtol=self.rtol,
         )
 
-    def test_g_nz(self, sample_coordinate, sample_point_source, finite_diff_g_nz):
+    def test_g_nz(self, sample_coordinate, sample_point_source, finite_diff_kernel_nu):
         """
-        Test kernel of g_nz against finite differences of the g_northing
+        Test kernel_nu against finite differences of the kernel_n
         """
+        distance = distance_cartesian(*sample_coordinate, *sample_point_source)
         npt.assert_allclose(
-            finite_diff_g_nz,
-            kernel_point_g_nz(*sample_coordinate, *sample_point_source),
+            finite_diff_kernel_nu,
+            kernel_nu(*sample_coordinate, *sample_point_source, distance),
             rtol=self.rtol,
         )
 
 
 def test_laplacian(sample_coordinate, sample_point_source):
     """
-    Test if diagonal tensor components satisfy Laplace's equation
+    Test if diagonal tensor kernels satisfy Laplace's equation
     """
-    g_ee = kernel_point_g_ee(*sample_coordinate, *sample_point_source)
-    g_nn = kernel_point_g_nn(*sample_coordinate, *sample_point_source)
-    g_zz = kernel_point_g_zz(*sample_coordinate, *sample_point_source)
+    distance = distance_cartesian(*sample_coordinate, *sample_point_source)
+    g_ee = kernel_ee(*sample_coordinate, *sample_point_source, distance)
+    g_nn = kernel_nn(*sample_coordinate, *sample_point_source, distance)
+    g_zz = kernel_uu(*sample_coordinate, *sample_point_source, distance)
     npt.assert_allclose(-g_zz, g_ee + g_nn)