diff --git a/choclo/dipole/_forward.py b/choclo/dipole/_forward.py
index 9676bf1a..90527257 100644
--- a/choclo/dipole/_forward.py
+++ b/choclo/dipole/_forward.py
@@ -113,9 +113,7 @@ def magnetic_field(
     b_n = c_m * (
         3 * dotproduct * r_n / distance**5 - magnetic_moment_north / distance**3
     )
-    b_u = c_m * (
-        3 * dotproduct * r_u / distance**5 - magnetic_moment_up / distance**3
-    )
+    b_u = c_m * (3 * dotproduct * r_u / distance**5 - magnetic_moment_up / distance**3)
     return b_e, b_n, b_u
 
 
@@ -297,9 +295,7 @@ def magnetic_n(
         + magnetic_moment_north * r_n
         + magnetic_moment_up * r_u
     )
-    result = (
-        3 * dotproduct * r_n / distance**5 - magnetic_moment_north / distance**3
-    )
+    result = 3 * dotproduct * r_n / distance**5 - magnetic_moment_north / distance**3
     return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * result
 
 
diff --git a/choclo/point/__init__.py b/choclo/point/__init__.py
index a72144c3..a80eb695 100644
--- a/choclo/point/__init__.py
+++ b/choclo/point/__init__.py
@@ -19,3 +19,15 @@
     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
index bd2023b8..4332e8f2 100644
--- a/choclo/point/_forward.py
+++ b/choclo/point/_forward.py
@@ -11,6 +11,18 @@
 
 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)
@@ -62,7 +74,10 @@ def gravity_pot(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
     distance = distance_cartesian(
         easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
     )
-    return GRAVITATIONAL_CONST * mass / distance
+    kernel = kernel_pot(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
+    return GRAVITATIONAL_CONST * mass * kernel
 
 
 @jit(nopython=True)
@@ -119,7 +134,9 @@ def gravity_e(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q,
     distance = distance_cartesian(
         easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
     )
-    kernel = -(easting_p - easting_q) / distance**3
+    kernel = kernel_e(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
     return GRAVITATIONAL_CONST * mass * kernel
 
 
@@ -177,7 +194,9 @@ def gravity_n(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q,
     distance = distance_cartesian(
         easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
     )
-    kernel = -(northing_p - northing_q) / distance**3
+    kernel = kernel_n(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
     return GRAVITATIONAL_CONST * mass * kernel
 
 
@@ -235,7 +254,9 @@ def gravity_u(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q,
     distance = distance_cartesian(
         easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
     )
-    kernel = -(upward_p - upward_q) / distance**3
+    kernel = kernel_u(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
     return GRAVITATIONAL_CONST * mass * kernel
 
 
@@ -300,7 +321,9 @@ def gravity_ee(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q,
     distance = distance_cartesian(
         easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
     )
-    kernel = 3 * (easting_p - easting_q) ** 2 / distance**5 - 1 / distance**3
+    kernel = kernel_ee(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
     return GRAVITATIONAL_CONST * mass * kernel
 
 
@@ -365,7 +388,9 @@ def gravity_nn(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q,
     distance = distance_cartesian(
         easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
     )
-    kernel = 3 * (northing_p - northing_q) ** 2 / distance**5 - 1 / distance**3
+    kernel = kernel_nn(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
     return GRAVITATIONAL_CONST * mass * kernel
 
 
@@ -430,7 +455,9 @@ def gravity_uu(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q,
     distance = distance_cartesian(
         easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
     )
-    kernel = 3 * (upward_p - upward_q) ** 2 / distance**5 - 1 / distance**3
+    kernel = kernel_uu(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
     return GRAVITATIONAL_CONST * mass * kernel
 
 
@@ -490,7 +517,9 @@ def gravity_en(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q,
     distance = distance_cartesian(
         easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
     )
-    kernel = 3 * (easting_p - easting_q) * (northing_p - northing_q) / distance**5
+    kernel = kernel_en(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
     return GRAVITATIONAL_CONST * mass * kernel
 
 
@@ -550,7 +579,9 @@ def gravity_eu(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q,
     distance = distance_cartesian(
         easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
     )
-    kernel = 3 * (easting_p - easting_q) * (upward_p - upward_q) / distance**5
+    kernel = kernel_eu(
+        easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, distance
+    )
     return GRAVITATIONAL_CONST * mass * kernel
 
 
@@ -610,5 +641,7 @@ def gravity_nu(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q,
     distance = distance_cartesian(
         easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
     )
-    kernel = 3 * (northing_p - northing_q) * (upward_p - upward_q) / distance**5
+    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/doc/api/index.rst b/doc/api/index.rst
index 721f53cb..c7ae9bbf 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -47,6 +47,24 @@ Magnetic
     dipole.magnetic_u
     dipole.magnetic_field
 
+Kernels
+^^^^^^^
+
+.. autosummary::
+   :toctree: generated/
+
+    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
+
+
 Rectangular Prisms
 ~~~~~~~~~~~~~~~~~~