Repository: fatiando/harmonica

Issue:
Add variable density tesseroids forward modelling
**Description of the desired feature**

Would be nice to implement the forward modelling method to compute the gravitational fields of tesseroids with variable densities described on [Soler et al. (2019)](https://www.doi.org/10.1093/gji/ggz277).

The main thing to do is adding the distance-based discretization algorithm, which should be implemented inside Numba jitted functions without allocating memory.
Also, the `tesseroid_gravity` function should take a density function per tesseroid as argument, which should be jitted in order to speed things up.


**Are you willing to help implement and maintain this feature?** Yes


Discussion:
santisoler: Thanks for the review @leouieda! I've added the introductory text and also changed the `jit` decorator for `njit`.

Patch:
diff --git a/doc/conf.py b/doc/conf.py
index 2c6b136..a62ee87 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -42,6 +42,7 @@ extensions = [
 intersphinx_mapping = {
     "python": ("https://docs.python.org/3/", None),
     "numpy": ("https://docs.scipy.org/doc/numpy/", None),
+    "numba": ("https://numba.pydata.org/numba-doc/latest/", None),
     "scipy": ("https://docs.scipy.org/doc/scipy/reference", None),
     "pandas": ("http://pandas.pydata.org/pandas-docs/stable/", None),
     "xarray": ("http://xarray.pydata.org/en/stable/", None),
diff --git a/doc/install.rst b/doc/install.rst
index f4a3443..e814a77 100644
--- a/doc/install.rst
+++ b/doc/install.rst
@@ -22,6 +22,7 @@ Dependencies
 * `numpy <http://www.numpy.org/>`__
 * `pandas <http://pandas.pydata.org/>`__
 * `numba <https://numba.pydata.org/>`__
+* `scipy <https://www.scipy.org/>`__
 * `xarray <https://xarray.pydata.org/>`__
 * `scikit-learn <https://scikit-learn.org>`__
 * `pooch <http://www.fatiando.org/pooch/>`__
diff --git a/doc/references.rst b/doc/references.rst
index a5df54f..ce9b28f 100644
--- a/doc/references.rst
+++ b/doc/references.rst
@@ -12,6 +12,7 @@ References
 .. [Hofmann-WellenhofMoritz2006] Hofmann-Wellenhof, B., & Moritz, H. (2006). Physical Geodesy (2nd, corr. ed. 2006 edition ed.). Wien ; New York: Springer.
 .. [Nagy2000] Nagy, D., Papp, G. & Benedek, J.(2000). The gravitational potential and its derivatives for the prism. Journal of Geodesy 74: 552. doi:`10.1007/s001900000116 <https://doi.org/10.1007/s001900000116>`__
 .. [Nagy2002] Nagy, D., Papp, G. & Benedek, J.(2002). Corrections to "The gravitational potential and its derivatives for the prism". Journal of Geodesy. doi:`10.1007/s00190-002-0264-7 <https://doi.org/10.1007/s00190-002-0264-7>`__
+.. [Soler2019] Soler, S. R., Pesce, A., Gimenez, M. E., & Uieda, L. (2019). Gravitational field calculation in spherical coordinates using variable densities in depth, Geophysical Journal International. doi: `10.1093/gji/ggz277 <https://doi.org/10.1093/gji/ggz277>`__
 .. [Soler2021] Soler, S. R. and Uieda, L. (2021). Gradient-boosted equivalent sources, Geophysical Journal International. doi:`10.1093/gji/ggab297 <https://doi.org/10.1093/gji/ggab297>`__
 .. [Vajda2004] Vajda, P., Vaníček, P., Novák, P. and Meurers, B. (2004). On evaluation of Newton integrals in geodetic coordinates: Exact formulation and spherical approximation. Contributions to Geophysics and Geodesy, 34(4), 289-314.
 .. [TurcotteSchubert2014] Turcotte, D. L., & Schubert, G. (2014). Geodynamics (3 edition). Cambridge, United Kingdom: Cambridge University Press.
diff --git a/environment.yml b/environment.yml
index 805f8db..605121b 100644
--- a/environment.yml
+++ b/environment.yml
@@ -9,6 +9,7 @@ dependencies:
   - numpy
   - pandas
   - numba
+  - scipy
   - scikit-learn
   - pooch>=0.7.0
   - verde>=1.5.0
diff --git a/examples/forward/tesseroid_variable_density.py b/examples/forward/tesseroid_variable_density.py
new file mode 100644
index 0000000..990deb3
--- /dev/null
+++ b/examples/forward/tesseroid_variable_density.py
@@ -0,0 +1,82 @@
+# Copyright (c) 2018 The Harmonica 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)
+#
+"""
+Tesseroids with variable density
+================================
+
+The :func:`harmonica.tesseroid_gravity` is capable of computing the
+gravitational effects of tesseroids whose density is defined through
+a continuous function of the radial coordinate. This is achieved by the
+application of the method introduced in [Soler2021]_.
+
+To do so we need to define a regular Python function for the density, which
+should have a single argument (the ``radius`` coordinate) and return the
+density of the tesseroids at that radial coordinate.
+In addition, we need to decorate the density function with
+:func:`numba.jit(nopython=True)` or ``numba.njit`` for short.
+
+On this example we will show how we can compute the gravitational effect of
+four tesseroids whose densities are given by a custom linear ``density``
+function.
+"""
+import matplotlib.pyplot as plt
+import cartopy.crs as ccrs
+import verde as vd
+import boule as bl
+import harmonica as hm
+from numba import njit
+
+
+# Use the WGS84 ellipsoid to obtain the mean Earth radius which we'll use to
+# reference the tesseroid
+ellipsoid = bl.WGS84
+mean_radius = ellipsoid.mean_radius
+
+# Define tesseroid with top surface at the mean Earth radius, a thickness of
+# 10km and a linear density function
+tesseroids = (
+    [-70, -60, -40, -30, mean_radius - 3e3, mean_radius],
+    [-70, -60, -30, -20, mean_radius - 5e3, mean_radius],
+    [-60, -50, -40, -30, mean_radius - 7e3, mean_radius],
+    [-60, -50, -30, -20, mean_radius - 10e3, mean_radius],
+)
+
+# Define a linear density function. We should use the jit decorator so Numba
+# can run the forward model efficiently.
+
+
+@njit
+def density(radius):
+    """Linear density function"""
+    top = mean_radius
+    bottom = mean_radius - 10e3
+    density_top = 2670
+    density_bottom = 3000
+    slope = (density_top - density_bottom) / (top - bottom)
+    return slope * (radius - bottom) + density_bottom
+
+
+# Define computation points on a regular grid at 100km above the mean Earth
+# radius
+coordinates = vd.grid_coordinates(
+    region=[-80, -40, -50, -10],
+    shape=(80, 80),
+    extra_coords=100e3 + ellipsoid.mean_radius,
+)
+
+# Compute the radial component of the acceleration
+gravity = hm.tesseroid_gravity(coordinates, tesseroids, density, field="g_z")
+print(gravity)
+
+# Plot the gravitational field
+fig = plt.figure(figsize=(8, 9))
+ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-60))
+img = ax.pcolormesh(*coordinates[:2], gravity, transform=ccrs.PlateCarree())
+plt.colorbar(img, ax=ax, pad=0, aspect=50, orientation="horizontal", label="mGal")
+ax.coastlines()
+ax.set_title("Downward component of gravitational acceleration")
+plt.show()
diff --git a/harmonica/forward/_tesseroid_utils.py b/harmonica/forward/_tesseroid_utils.py
new file mode 100644
index 0000000..5f95256
--- /dev/null
+++ b/harmonica/forward/_tesseroid_utils.py
@@ -0,0 +1,499 @@
+# Copyright (c) 2018 The Harmonica 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)
+#
+"""
+Utils functions for tesseroid forward modelling
+"""
+from numba import jit
+import numpy as np
+from numpy.polynomial.legendre import leggauss
+
+from .utils import distance_spherical
+
+
+@jit(nopython=True)
+def gauss_legendre_quadrature(
+    longitude,
+    cosphi,
+    sinphi,
+    radius,
+    tesseroid,
+    density,
+    glq_nodes,
+    glq_weights,
+    kernel,
+):  # pylint: disable=too-many-locals
+    r"""
+    Compute the effect of a tesseroid on a single observation point through GLQ
+
+    The tesseroid is converted into a set of point masses located on the
+    scaled nodes of the Gauss-Legendre Quadrature. The number of point masses
+    created from each tesseroid is equal to the product of the GLQ degrees for
+    each direction (:math:`N_r`, :math:`N_\lambda`, :math:`N_\phi`). The mass
+    of each point mass is defined as the product of the tesseroid density
+    (:math:`\rho`), the GLQ weights for each direction (:math:`W_i^r`,
+    :math:`W_j^\phi`, :math:`W_k^\lambda`), the scale constant :math:`A` and
+    the :math:`\kappa` factor evaluated on the coordinates of the point mass.
+
+    Parameters
+    ----------
+    longitude : float
+        Longitudinal coordinate of the observation points in radians.
+    cosphi : float
+        Cosine of the latitudinal coordinate of the observation point in
+        radians.
+    sinphi : float
+        Sine of the latitudinal coordinate of the observation point in
+        radians.
+    radius : float
+        Radial coordinate of the observation point in meters.
+    tesseroids : 1d-array
+        Array containing the boundaries of the tesseroid:
+        ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``.
+        Horizontal boundaries should be in degrees and radial boundaries in
+        meters.
+    density : float
+        Density of the tesseroid in SI units.
+    glq_nodes : list
+        Unscaled location of GLQ nodes for each direction.
+    glq_weights : list
+        GLQ weigths for each node for each direction.
+    kernel : func
+        Kernel function for the gravitational field of point masses.
+
+    """
+    # Get tesseroid boundaries
+    w, e, s, n, bottom, top = tesseroid[:]
+    # Calculate the A factor for the tesseroid
+    a_factor = 1 / 8 * np.radians(e - w) * np.radians(n - s) * (top - bottom)
+    # Unpack nodes and weights
+    lon_nodes, lat_nodes, rad_nodes = glq_nodes[:]
+    lon_weights, lat_weights, rad_weights = glq_weights[:]
+    # Compute effect of the tesseroid on the observation point
+    # by iterating over the location of the point masses
+    # (move the iteration along the longitudinal nodes to the bottom for
+    # optimization: reduce the number of times that the trigonometric functions
+    # are evaluated)
+    result = 0.0
+    for j, lat_node in enumerate(lat_nodes):
+        # Get the latitude of the point mass
+        latitude_p = np.radians(0.5 * (n - s) * lat_node + 0.5 * (n + s))
+        cosphi_p = np.cos(latitude_p)
+        sinphi_p = np.sin(latitude_p)
+        for k, rad_node in enumerate(rad_nodes):
+            # Get the radius of the point mass
+            radius_p = 0.5 * (top - bottom) * rad_node + 0.5 * (top + bottom)
+            # Get kappa constant for the point mass
+            kappa = radius_p ** 2 * cosphi_p
+            for i, lon_node in enumerate(lon_nodes):
+                # Get the longitude of the point mass
+                longitude_p = np.radians(0.5 * (e - w) * lon_node + 0.5 * (e + w))
+                # Compute the mass of the point mass
+                mass = (
+                    density
+                    * a_factor
+                    * kappa
+                    * lon_weights[i]
+                    * lat_weights[j]
+                    * rad_weights[k]
+                )
+                # Add effect of the current point mass to the result
+                result += mass * kernel(
+                    longitude,
+                    cosphi,
+                    sinphi,
+                    radius,
+                    longitude_p,
+                    cosphi_p,
+                    sinphi_p,
+                    radius_p,
+                )
+    return result
+
+
+def glq_nodes_weights(glq_degrees):
+    """
+    Calculate GLQ unscaled nodes and weights
+
+    Parameters
+    ----------
+    glq_degrees : list
+        List of GLQ degrees for each direction: ``longitude``, ``latitude``,
+        ``radius``.
+
+    Returns
+    -------
+    glq_nodes : list
+        Unscaled GLQ nodes for each direction: ``longitude``, ``latitude``,
+        ``radius``.
+    glq_weights : list
+        GLQ weights for each node on each direction: ``longitude``,
+        ``latitude``, ``radius``.
+    """
+    # Unpack GLQ degrees
+    lon_degree, lat_degree, rad_degree = glq_degrees[:]
+    # Get nodes coordinates and weights
+    lon_node, lon_weights = leggauss(lon_degree)
+    lat_node, lat_weights = leggauss(lat_degree)
+    rad_node, rad_weights = leggauss(rad_degree)
+    # Reorder nodes and weights
+    glq_nodes = (lon_node, lat_node, rad_node)
+    glq_weights = (lon_weights, lat_weights, rad_weights)
+    return glq_nodes, glq_weights
+
+
+@jit(nopython=True)
+def _adaptive_discretization(
+    coordinates,
+    tesseroid,
+    distance_size_ratio,
+    stack,
+    small_tesseroids,
+    radial_discretization=False,
+):
+    """
+    Perform the adaptive discretization algorithm on a tesseroid
+
+    It apply the three or two dimensional adaptive discretization algorithm on
+    a tesseroid after a single computation point.
+
+    Parameters
+    ----------
+    coordinates : array
+        Array containing ``longitude``, ``latitude`` and ``radius`` of a single
+        computation point.
+    tesseroid : array
+        Array containing the boundaries of the tesseroid.
+    distance_size_ratio : float
+        Value for the distance-size ratio. A greater value will perform more
+        discretizations.
+    stack : 2d-array
+        Array with shape ``(6, stack_size)`` that will temporarly hold the
+        small tesseroids that are not yet processed.
+        If too many discretizations will take place, increase the
+        ``stack_size``.
+    small_tesseroids : 2d-array
+        Array with shape ``(6, MAX_DISCRETIZATIONS)`` that will contain every
+        small tesseroid produced by the adaptive discretization algorithm.
+        If too many discretizations will take place, increase the
+        ``MAX_DISCRETIZATIONS``.
+    radial_discretization : bool (optional)
+        If ``True`` the three dimensional adaptive discretization will be
+        applied.
+        If ``False`` the two dimensional adaptive discretization will be
+        applied, i.e. the tesseroid will only be split on the ``longitude`` and
+        ``latitude`` directions.
+        Default ``False``.
+
+    Returns
+    -------
+    n_splits : int
+        Total number of small tesseroids generated by the algorithm.
+    """
+    # Create stack of tesseroids
+    stack[0] = tesseroid
+    stack_top = 0
+    n_splits = 0
+    while stack_top >= 0:
+        # Pop the first tesseroid from the stack
+        tesseroid = stack[stack_top]
+        stack_top -= 1
+        # Get its dimensions
+        l_lon, l_lat, l_rad = _tesseroid_dimensions(tesseroid)
+        # Get distance between computation point and center of tesseroid
+        distance = _distance_tesseroid_point(coordinates, tesseroid)
+        # Check inequality
+        n_lon, n_lat, n_rad = 1, 1, 1
+        if distance / l_lon < distance_size_ratio:
+            n_lon = 2
+        if distance / l_lat < distance_size_ratio:
+            n_lat = 2
+        if distance / l_rad < distance_size_ratio and radial_discretization:
+            n_rad = 2
+        # Apply discretization
+        if n_lon * n_lat * n_rad > 1:
+            # Raise error if stack overflow
+            # Number of tesseroids in stack = stack_top + 1
+            if (stack_top + 1) + n_lon * n_lat * n_rad > stack.shape[0]:
+                raise OverflowError("Stack Overflow. Try to increase the stack size.")
+            stack_top = _split_tesseroid(
+                tesseroid, n_lon, n_lat, n_rad, stack, stack_top
+            )
+        else:
+            # Raise error if small_tesseroids overflow
+            if n_splits + 1 > small_tesseroids.shape[0]:
+                raise OverflowError(
+                    "Exceeded maximum discretizations."
+                    + " Please increase the MAX_DISCRETIZATIONS."
+                )
+            small_tesseroids[n_splits] = tesseroid
+            n_splits += 1
+    return n_splits
+
+
+@jit(nopython=True)
+def _split_tesseroid(
+    tesseroid, n_lon, n_lat, n_rad, stack, stack_top
+):  # pylint: disable=too-many-locals
+    """
+    Split tesseroid along each dimension
+    """
+    w, e, s, n, bottom, top = tesseroid[:]
+    # Compute differential distance
+    d_lon = (e - w) / n_lon
+    d_lat = (n - s) / n_lat
+    d_rad = (top - bottom) / n_rad
+    for i in range(n_lon):
+        for j in range(n_lat):
+            for k in range(n_rad):
+                stack_top += 1
+                stack[stack_top, 0] = w + d_lon * i
+                stack[stack_top, 1] = w + d_lon * (i + 1)
+                stack[stack_top, 2] = s + d_lat * j
+                stack[stack_top, 3] = s + d_lat * (j + 1)
+                stack[stack_top, 4] = bottom + d_rad * k
+                stack[stack_top, 5] = bottom + d_rad * (k + 1)
+    return stack_top
+
+
+@jit(nopython=True)
+def _tesseroid_dimensions(tesseroid):
+    """
+    Calculate the dimensions of the tesseroid.
+    """
+    w, e, s, n, bottom, top = tesseroid[:]
+    w, e, s, n = np.radians(w), np.radians(e), np.radians(s), np.radians(n)
+    latitude_center = (n + s) / 2
+    l_lat = top * np.arccos(np.sin(n) * np.sin(s) + np.cos(n) * np.cos(s))
+    l_lon = top * np.arccos(
+        np.sin(latitude_center) ** 2 + np.cos(latitude_center) ** 2 * np.cos(e - w)
+    )
+    l_rad = top - bottom
+    return l_lon, l_lat, l_rad
+
+
+@jit(nopython=True)
+def _distance_tesseroid_point(
+    coordinates, tesseroid
+):  # pylint: disable=too-many-locals
+    """
+    Distance between a computation point and the center of a tesseroid
+    """
+    # Get center of the tesseroid
+    w, e, s, n, bottom, top = tesseroid[:]
+    longitude_p = (w + e) / 2
+    latitude_p = (s + n) / 2
+    radius_p = (bottom + top) / 2
+    # Get distance between computation point and tesseroid center
+    distance = distance_spherical(coordinates, (longitude_p, latitude_p, radius_p))
+    return distance
+
+
+def _check_tesseroids(tesseroids):  # pylint: disable=too-many-branches
+    """
+    Check if tesseroids boundaries are well defined
+
+    A valid tesseroid should have:
+        - latitudinal boundaries within the [-90, 90] degrees interval,
+        - north boundaries greater or equal than the south boundaries,
+        - radial boundaries positive or zero,
+        - top boundaries greater or equal than the bottom boundaries,
+        - longitudinal boundaries within the [-180, 360] degrees interval,
+        - longitudinal interval must not be greater than one turn around the
+          globe.
+
+    Some valid tesseroids have its west boundary greater than the east one,
+    e.g. ``(350, 10, ...)``. On these cases the ``_longitude_continuity``
+    function is applied in order to move the longitudinal coordinates to the
+    [-180, 180) interval. Any valid tesseroid should have east boundaries
+    greater than the west boundaries before or after applying longitude
+    continuity.
+
+    Parameters
+    ----------
+    tesseroids : 2d-array
+        Array containing the boundaries of the tesseroids in the following
+        order: ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``.
+        Longitudinal and latitudinal boundaries must be in degrees.
+        The array must have the following shape: (``n_tesseroids``, 6), where
+        ``n_tesseroids`` is the total number of tesseroids.
+
+    Returns
+    -------
+    tesseroids :  2d-array
+        Array containing the boundaries of the tesseroids.
+        If no longitude continuity needs to be applied, the returned array is
+        the same one as the orignal.
+        Otherwise, it's copied and its longitudinal boundaries are modified.
+    """
+    west, east, south, north, bottom, top = tuple(tesseroids[:, i] for i in range(6))
+    err_msg = "Invalid tesseroid or tesseroids. "
+    # Check if latitudinal boundaries are inside the [-90, 90] interval
+    invalid = np.logical_or(
+        np.logical_or(south < -90, south > 90), np.logical_or(north < -90, north > 90)
+    )
+    if (invalid).any():
+        err_msg += (
+            "The latitudinal boundaries must be inside the [-90, 90] "
+            + "degrees interval.\n"
+        )
+        for tess in tesseroids[invalid]:
+            err_msg += "\tInvalid tesseroid: {}\n".format(tess)
+        raise ValueError(err_msg)
+    # Check if south boundary is not greater than the corresponding north
+    # boundary
+    invalid = south > north
+    if (invalid).any():
+        err_msg += "The south boundary can't be greater than the north one.\n"
+        for tess in tesseroids[invalid]:
+            err_msg += "\tInvalid tesseroid: {}\n".format(tess)
+        raise ValueError(err_msg)
+    # Check if radial boundaries are positive or zero
+    invalid = np.logical_or(bottom < 0, top < 0)
+    if (invalid).any():
+        err_msg += "The bottom and top radii should be positive or zero.\n"
+        for tess in tesseroids[invalid]:
+            err_msg += "\tInvalid tesseroid: {}\n".format(tess)
+        raise ValueError(err_msg)
+    # Check if top boundary is not greater than the corresponding bottom
+    # boundary
+    invalid = bottom > top
+    if (invalid).any():
+        err_msg += "The bottom radius boundary can't be greater than the top one.\n"
+        for tess in tesseroids[invalid]:
+            err_msg += "\tInvalid tesseroid: {}\n".format(tess)
+        raise ValueError(err_msg)
+    # Check if longitudinal boundaries are inside the [-180, 360] interval
+    invalid = np.logical_or(
+        np.logical_or(west < -180, west > 360), np.logical_or(east < -180, east > 360)
+    )
+    if (invalid).any():
+        err_msg += (
+            "The longitudinal boundaries must be inside the [-180, 360] "
+            + "degrees interval.\n"
+        )
+        for tess in tesseroids[invalid]:
+            err_msg += "\tInvalid tesseroid: {}\n".format(tess)
+        raise ValueError(err_msg)
+    # Apply longitude continuity if w > e
+    if (west > east).any():
+        tesseroids = _longitude_continuity(tesseroids)
+        west, east, south, north, bottom, top = tuple(
+            tesseroids[:, i] for i in range(6)
+        )
+    # Check if west boundary is not greater than the corresponding east
+    # boundary, even after applying the longitude continuity
+    invalid = west > east
+    if (invalid).any():
+        err_msg += "The west boundary can't be greater than the east one.\n"
+        for tess in tesseroids[invalid]:
+            err_msg += "\tInvalid tesseroid: {}\n".format(tess)
+        raise ValueError(err_msg)
+    # Check if the longitudinal interval is not grater than one turn around the
+    # globe
+    invalid = east - west > 360
+    if (invalid).any():
+        err_msg += (
+            "The difference between east and west boundaries cannot be greater than "
+            + "one turn around the globe.\n"
+        )
+        for tess in tesseroids[invalid]:
+            err_msg += "\tInvalid tesseroid: {}\n".format(tess)
+        raise ValueError(err_msg)
+    return tesseroids
+
+
+def _check_points_outside_tesseroids(
+    coordinates, tesseroids
+):  # pylint: disable=too-many-locals
+    """
+    Check if computation points are not inside the tesseroids
+
+    Parameters
+    ----------
+    coordinates : 2d-array
+        Array containing the coordinates of the computation points in the
+        following order: ``longitude``, ``latitude`` and ``radius``.
+        Both ``longitude`` and ``latitude`` must be in degrees.
+        The array must have the following shape: (3, ``n_points``), where
+        ``n_points`` is the total number of computation points.
+    tesseroids : 2d-array
+        Array containing the boundaries of the tesseroids in the following
+        order: ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``.
+        Longitudinal and latitudinal boundaries must be in degrees.
+        The array must have the following shape: (``n_tesseroids``, 6), where
+        ``n_tesseroids`` is the total number of tesseroids.
+        This array of tesseroids must have longitude continuity and valid
+        boundaries.
+        Run ``_check_tesseroids`` before.
+    """
+    longitude, latitude, radius = coordinates[:]
+    west, east, south, north, bottom, top = tuple(tesseroids[:, i] for i in range(6))
+    # Longitudinal boundaries of the tesseroid must be compared with
+    # longitudinal coordinates of computation points when moved to
+    # [0, 360) and [-180, 180).
+    longitude_360 = longitude % 360
+    longitude_180 = ((longitude + 180) % 360) - 180
+    inside_longitude = np.logical_or(
+        np.logical_and(
+            west < longitude_360[:, np.newaxis], longitude_360[:, np.newaxis] < east
+        ),
+        np.logical_and(
+            west < longitude_180[:, np.newaxis], longitude_180[:, np.newaxis] < east
+        ),
+    )
+    inside_latitude = np.logical_and(
+        south < latitude[:, np.newaxis], latitude[:, np.newaxis] < north
+    )
+    inside_radius = np.logical_and(
+        bottom < radius[:, np.newaxis], radius[:, np.newaxis] < top
+    )
+    # Build array of booleans.
+    # The (i, j) element is True if the computation point i is inside the
+    # tesseroid j.
+    inside = inside_longitude * inside_latitude * inside_radius
+    if inside.any():
+        err_msg = (
+            "Found computation point inside tesseroid. "
+            + "Computation points must be outside of tesseroids.\n"
+        )
+        for point_i, tess_i in np.argwhere(inside):
+            err_msg += "\tComputation point '{}' found inside tesseroid '{}'\n".format(
+                coordinates[:, point_i], tesseroids[tess_i, :]
+            )
+        raise ValueError(err_msg)
+
+
+def _longitude_continuity(tesseroids):
+    """
+    Modify longitudinal boundaries of tesseroids to ensure longitude continuity
+
+    Longitudinal boundaries of the tesseroids are moved to the ``[-180, 180)``
+    degrees interval in case the ``west`` boundary is numerically greater than
+    the ``east`` one.
+
+    Parameters
+    ----------
+    tesseroids : 2d-array
+        Longitudinal and latitudinal boundaries must be in degrees.
+        Array containing the boundaries of each tesseroid:
+        ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top`` under a geocentric
+        spherical coordinate system.
+        The array must have the following shape: (``n_tesseroids``, 6), where
+        ``n_tesseroids`` is the total number of tesseroids.
+
+    Returns
+    -------
+    tesseroids : 2d-array
+        Modified boundaries of the tesseroids.
+    """
+    # Copy the tesseroids to avoid modifying the original tesseroids array
+    tesseroids = tesseroids.copy()
+    west, east = tesseroids[:, 0], tesseroids[:, 1]
+    tess_to_be_changed = west > east
+    east[tess_to_be_changed] = ((east[tess_to_be_changed] + 180) % 360) - 180
+    west[tess_to_be_changed] = ((west[tess_to_be_changed] + 180) % 360) - 180
+    return tesseroids
diff --git a/harmonica/forward/_tesseroid_variable_density.py b/harmonica/forward/_tesseroid_variable_density.py
new file mode 100644
index 0000000..4b4c991
--- /dev/null
+++ b/harmonica/forward/_tesseroid_variable_density.py
@@ -0,0 +1,259 @@
+# Copyright (c) 2018 The Harmonica 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)
+#
+"""
+Utils functions for tesseroids with variable density
+"""
+from numba import jit
+import numpy as np
+from scipy.optimize import minimize_scalar
+
+DELTA_RATIO = 0.1
+
+
+@jit(nopython=True)
+def gauss_legendre_quadrature_variable_density(
+    longitude,
+    cosphi,
+    sinphi,
+    radius,
+    tesseroid,
+    density,
+    glq_nodes,
+    glq_weights,
+    kernel,
+):  # pylint: disable=too-many-locals
+    r"""
+    Compute the effect of a tesseroid on a single observation point through GLQ
+
+    The tesseroid is converted into a set of point masses located on the
+    scaled nodes of the Gauss-Legendre Quadrature. The number of point masses
+    created from each tesseroid is equal to the product of the GLQ degrees for
+    each direction (:math:`N_r`, :math:`N_\lambda`, :math:`N_\phi`). The mass
+    of each point mass is defined as the product of the tesseroid density
+    (:math:`\rho`), the GLQ weights for each direction (:math:`W_i^r`,
+    :math:`W_j^\phi`, :math:`W_k^\lambda`), the scale constant :math:`A` and
+    the :math:`\kappa` factor evaluated on the coordinates of the point mass.
+
+    Parameters
+    ----------
+    longitude : float
+        Longitudinal coordinate of the observation points in radians.
+    cosphi : float
+        Cosine of the latitudinal coordinate of the observation point in
+        radians.
+    sinphi : float
+        Sine of the latitudinal coordinate of the observation point in
+        radians.
+    radius : float
+        Radial coordinate of the observation point in meters.
+    tesseroids : 1d-array
+        Array containing the boundaries of the tesseroid:
+        ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``.
+        Horizontal boundaries should be in degrees and radial boundaries in
+        meters.
+    density : func
+        Density func of the tesseroid in SI units.
+    glq_nodes : list
+        Unscaled location of GLQ nodes for each direction.
+    glq_weights : list
+        GLQ weigths for each node for each direction.
+    kernel : func
+        Kernel function for the gravitational field of point masses.
+
+    """
+    # Get tesseroid boundaries
+    w, e, s, n, bottom, top = tesseroid[:]
+    # Calculate the A factor for the tesseroid
+    a_factor = 1 / 8 * np.radians(e - w) * np.radians(n - s) * (top - bottom)
+    # Unpack nodes and weights
+    lon_nodes, lat_nodes, rad_nodes = glq_nodes[:]
+    lon_weights, lat_weights, rad_weights = glq_weights[:]
+    # Compute effect of the tesseroid on the observation point
+    # by iterating over the location of the point masses
+    # (move the iteration along the longitudinal nodes to the bottom for
+    # optimization: reduce the number of times that the trigonometric functions
+    # are evaluated)
+    result = 0.0
+    for j, lat_node in enumerate(lat_nodes):
+        # Get the latitude of the point mass
+        latitude_p = np.radians(0.5 * (n - s) * lat_node + 0.5 * (n + s))
+        cosphi_p = np.cos(latitude_p)
+        sinphi_p = np.sin(latitude_p)
+        for k, rad_node in enumerate(rad_nodes):
+            # Get the radius of the point mass
+            radius_p = 0.5 * (top - bottom) * rad_node + 0.5 * (top + bottom)
+            density_p = density(radius_p)
+            # Get kappa constant for the point mass
+            kappa = radius_p ** 2 * cosphi_p
+            for i, lon_node in enumerate(lon_nodes):
+                # Get the longitude of the point mass
+                longitude_p = np.radians(0.5 * (e - w) * lon_node + 0.5 * (e + w))
+                # Compute the mass of the point mass
+                mass = (
+                    density_p
+                    * a_factor
+                    * kappa
+                    * lon_weights[i]
+                    * lat_weights[j]
+                    * rad_weights[k]
+                )
+                # Add effect of the current point mass to the result
+                result += mass * kernel(
+                    longitude,
+                    cosphi,
+                    sinphi,
+                    radius,
+                    longitude_p,
+                    cosphi_p,
+                    sinphi_p,
+                    radius_p,
+                )
+    return result
+
+
+# Density-based discretization functions
+# --------------------------------------
+def density_based_discretization(tesseroids, density):
+    """
+    Apply density_based discretization to a collection of tesseroids
+
+    Parameters
+    ----------
+    tesseroids : 2d-array
+        Array containing the coordinates of the tesseroid. Each row of the
+        array should contain the boundaries of each tesseroid in the following
+        order: ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``.
+        The longitudinal and latitudinal boundaries should be in degrees, while
+        the radial ones must be in meters.
+    density : func
+        Continuous density function of the tesseroid in SI units.
+
+    Returns
+    -------
+    discretized_tesseroids : 2d-array
+        Array containing the coordinates of radially discretized tesseriods.
+        Each row of the array will have the boundaries for each new tesseroid.
+    """
+    discretized_tesseroids = []
+    for tesseroid in tesseroids:
+        discretized_tesseroids.extend(_density_based_discretization(tesseroid, density))
+    return np.atleast_2d(discretized_tesseroids)
+
+
+def _density_based_discretization(
+    tesseroid, density
+):  # pylint: disable=too-many-locals
+    """
+    Applies density-based discretization to a single tesseroid
+
+    Splits the tesseroid on the points of maximum density variance
+
+    Parameters
+    ----------
+    tesseroid : tuple
+        Tuple containing the boundaries of the tesseroid:
+        ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``.
+        Horizontal boundaries should be in degrees and radial boundaries in
+        meters.
+    density : func
+        Density func of the tesseroid in SI units.
+
+    Returns
+    -------
+    tesseroids : list
+        List containing the boundaries of discretized tesseroids.
+    """
+    # Define normalized density
+    def normalized_density(radius):
+        return (density(radius) - density_min) / (density_max - density_min)
+
+    # Get boundaries of original tesseroid
+    w, e, s, n, bottom, top = tesseroid[:]
+    # Get minimum and maximum values of the density
+    density_min, density_max = density_minmax(density, bottom, top)
+    # Return the original tesseroid if max and min densities are equal
+    if np.isclose(density_min, density_max):
+        return [tesseroid]
+    # Store the size of the original tesseroid
+    size_original_tesseroid = top - bottom
+    # Initialize list of pending and output tesseroids
+    pending, tesseroids = [tesseroid], []
+    # Discretization of the tesseroid
+    while pending:
+        tesseroid = pending.pop(0)
+        bottom, top = tesseroid[-2:]
+        radius_split, max_diff = maximum_absolute_diff(normalized_density, bottom, top)
+        size_ratio = (top - bottom) / size_original_tesseroid
+        if max_diff * size_ratio > DELTA_RATIO:
+            pending.append([w, e, s, n, radius_split, top])
+            pending.append([w, e, s, n, bottom, radius_split])
+        else:
+            tesseroids.append([w, e, s, n, bottom, top])
+    return tesseroids
+
+
+def density_minmax(density, bottom, top):
+    """
+    Compute the minimum and maximum value of a bounded density
+    """
+    # Calculate min and max density values at the top and bottom boundaries
+    density_bounds_min, density_bounds_max = np.sort([density(bottom), density(top)])
+    # Estimate the minimum value of the density function withing bounds
+    kwargs = dict(bounds=[bottom, top], method="bounded")
+    minimum = np.min(
+        (
+            minimize_scalar(density, **kwargs).fun,
+            density_bounds_min,
+        )
+    )
+    # Estimate the maximum value of the density function withing bounds
+    maximum = np.max(
+        (
+            -minimize_scalar(lambda radius: -density(radius), **kwargs).fun,
+            density_bounds_max,
+        )
+    )
+    return minimum, maximum
+
+
+def maximum_absolute_diff(normalized_density, bottom, top):
+    """
+    Compute maximum abs difference between normalized density and straight line
+
+    The maximum difference is computed within the ``bottom`` and ``top``
+    boundaries.
+    """
+
+    def neg_absolute_difference(radius):
+        """
+        Define minus absolute diff between normalized density and straight line
+        """
+        return -np.abs(
+            normalized_density(radius)
+            - straight_line(radius, normalized_density, bottom, top)
+        )
+
+    # Use scipy.optimize.minimize_scalar for maximizing the absolute difference
+    result = minimize_scalar(
+        neg_absolute_difference,
+        bounds=[bottom, top],
+        method="bounded",
+    )
+    # Get maximum difference and the radius at which it takes place
+    radius_split = result.x
+    max_diff = -result.fun
+    return radius_split, max_diff
+
+
+def straight_line(radius, normalized_density, bottom, top):
+    """
+    Compute the reference straight line that joins points of normalized density
+    """
+    norm_density_bottom = normalized_density(bottom)
+    norm_density_top = normalized_density(top)
+    slope = (norm_density_top - norm_density_bottom) / (top - bottom)
+    return slope * (radius - bottom) + norm_density_bottom
diff --git a/harmonica/forward/tesseroid.py b/harmonica/forward/tesseroid.py
index 990ed61..9e14a52 100644
--- a/harmonica/forward/tesseroid.py
+++ b/harmonica/forward/tesseroid.py
@@ -9,14 +9,23 @@ Forward modelling for tesseroids
 """
 from numba import jit, prange
 import numpy as np
-from numpy.polynomial.legendre import leggauss
 
 from ..constants import GRAVITATIONAL_CONST
-from .utils import distance_spherical
 from .point_mass import (
     kernel_potential_spherical,
     kernel_g_z_spherical,
 )
+from ._tesseroid_utils import (
+    _check_tesseroids,
+    _check_points_outside_tesseroids,
+    _adaptive_discretization,
+    glq_nodes_weights,
+    gauss_legendre_quadrature,
+)
+from ._tesseroid_variable_density import (
+    density_based_discretization,
+    gauss_legendre_quadrature_variable_density,
+)
 
 STACK_SIZE = 100
 MAX_DISCRETIZATIONS = 100000
@@ -58,8 +67,13 @@ def tesseroid_gravity(
         spherical coordinate system.
         The longitudinal and latitudinal boundaries should be in degrees, while
         the radial ones must be in meters.
-    density : list or array
+    density : list, array or func decorated with :func:`numba.jit`
         List or array containing the density of each tesseroid in kg/m^3.
+        Alternatively, it can be a continuous function within the boundaries of
+        the tesseroids, in which case the variable density tesseroids method
+        introduced in [Soler2019]_ will be used.
+        If ``density`` is a function, it should be decorated with
+        :func:`numba.jit`.
     field : str
         Gravitational field that wants to be computed.
         The available fields are:
@@ -93,6 +107,10 @@ def tesseroid_gravity(
         Gravitational field generated by the tesseroids on the computation
         points.
 
+    References
+    ----------
+    [Soler2019]_
+
     Examples
     --------
 
@@ -114,6 +132,19 @@ def tesseroid_gravity(
     >>> tesseroid_gravity(coordinates, tesseroid, density, field="g_z")
     array(112.54539933)
 
+    >>> # Define a linear density function for the same tesseroid.
+    >>> # It should be decorated with numba.njit
+    >>> from numba import jit
+    >>> @jit
+    ... def linear_density(radius):
+    ...     density_top = 2670.
+    ...     density_bottom = 3300.
+    ...     slope = (density_top - density_bottom) / (top - bottom)
+    ...     return slope * (radius - bottom) + density_bottom
+    >>> # Compute the downward acceleration it generates
+    >>> tesseroid_gravity(coordinates, tesseroid, linear_density, field="g_z")
+    array(125.80498632)
+
     """
     kernels = {"potential": kernel_potential_spherical, "g_z": kernel_g_z_spherical}
     if field not in kernels:
@@ -124,21 +155,26 @@ def tesseroid_gravity(
     # Convert coordinates, tesseroids and density to arrays
     coordinates = tuple(np.atleast_1d(i).ravel() for i in coordinates[:3])
     tesseroids = np.atleast_2d(tesseroids)
-    density = np.atleast_1d(density).ravel()
     # Sanity checks for tesseroids and computation points
     if not disable_checks:
-        if density.size != tesseroids.shape[0]:
+        tesseroids = _check_tesseroids(tesseroids)
+        _check_points_outside_tesseroids(coordinates, tesseroids)
+    # Check if density is a function or constant values
+    if callable(density):
+        # Run density-based discretization on each tesseroid
+        tesseroids = density_based_discretization(tesseroids, density)
+    else:
+        density = np.atleast_1d(density).ravel()
+        if not disable_checks and density.size != tesseroids.shape[0]:
             raise ValueError(
                 "Number of elements in density ({}) ".format(density.size)
                 + "mismatch the number of tesseroids ({})".format(tesseroids.shape[0])
             )
-        tesseroids = _check_tesseroids(tesseroids)
-        _check_points_outside_tesseroids(coordinates, tesseroids)
     # Get GLQ unscaled nodes, weights and number of nodes for each small
     # tesseroid
     glq_nodes, glq_weights = glq_nodes_weights(GLQ_DEGREES)
     # Compute gravitational field
-    dispatcher(parallel)(
+    dispatcher(parallel, density)(
         coordinates,
         tesseroids,
         density,
@@ -157,14 +193,23 @@ def tesseroid_gravity(
     return result.reshape(cast.shape)
 
 
-def dispatcher(parallel):
+def dispatcher(parallel, density):
     """
-    Return the parallelized or serialized forward modelling function
+    Return the jitted compiled forward modelling function
+
+    The choice of the forward modelling function is based on whether the
+    density is a function and if the model should be run in parallel.
     """
-    dispatchers = {
-        True: jit_tesseroid_gravity_parallel,
-        False: jit_tesseroid_gravity_serial,
-    }
+    if callable(density):
+        dispatchers = {
+            True: jit_tesseroid_gravity_variable_density_parallel,
+            False: jit_tesseroid_gravity_variable_density_serial,
+        }
+    else:
+        dispatchers = {
+            True: jit_tesseroid_gravity_parallel,
+            False: jit_tesseroid_gravity_serial,
+        }
     return dispatchers[parallel]
 
 
@@ -207,12 +252,6 @@ def jit_tesseroid_gravity(
         should be in meters.
     density : 1d-array
         Density of each tesseroid in SI units.
-    stack : 2d-array
-        Empty array where tesseroids created by adaptive discretization
-        algorithm will be processed.
-    small_tesseroids : 2d-array
-        Empty array where smaller tesseroids created by adaptive discretization
-        algorithm will be stored.
     result : 1d-array
         Array where the gravitational effect of each tesseroid will be added.
     distance_size_ratio : float
@@ -269,489 +308,99 @@ def jit_tesseroid_gravity(
                 )
 
 
-@jit(nopython=True)
-def gauss_legendre_quadrature(
-    longitude,
-    cosphi,
-    sinphi,
-    radius,
-    tesseroid,
+def jit_tesseroid_gravity_variable_density(
+    coordinates,
+    tesseroids,
     density,
+    result,
+    distance_size_ratio,
+    radial_adaptive_discretization,
     glq_nodes,
     glq_weights,
     kernel,
-):  # pylint: disable=too-many-locals
-    r"""
-    Compute the effect of a tesseroid on a single observation point through GLQ
-
-    The tesseroid is converted into a set of point masses located on the
-    scaled nodes of the Gauss-Legendre Quadrature. The number of point masses
-    created from each tesseroid is equal to the product of the GLQ degrees for
-    each direction (:math:`N_r`, :math:`N_\lambda`, :math:`N_\phi`). The mass
-    of each point mass is defined as the product of the tesseroid density
-    (:math:`\rho`), the GLQ weights for each direction (:math:`W_i^r`,
-    :math:`W_j^\phi`, :math:`W_k^\lambda`), the scale constant :math:`A` and
-    the :math:`\kappa` factor evaluated on the coordinates of the point mass.
-
-    Parameters
-    ----------
-    longitude : float
-        Longitudinal coordinate of the observation points in radians.
-    cosphi : float
-        Cosine of the latitudinal coordinate of the observation point in
-        radians.
-    sinphi : float
-        Sine of the latitudinal coordinate of the observation point in
-        radians.
-    radius : float
-        Radial coordinate of the observation point in meters.
-    tesseroids : 1d-array
-        Array containing the boundaries of the tesseroid:
-        ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``.
-        Horizontal boundaries should be in degrees and radial boundaries in
-        meters.
-    density : float
-        Density of the tesseroid in SI units.
-    glq_nodes : list
-        Unscaled location of GLQ nodes for each direction.
-    glq_weights : list
-        GLQ weigths for each node for each direction.
-    kernel : func
-        Kernel function for the gravitational field of point masses.
-
-    """
-    # Get tesseroid boundaries
-    w, e, s, n, bottom, top = tesseroid[:]
-    # Calculate the A factor for the tesseroid
-    a_factor = 1 / 8 * np.radians(e - w) * np.radians(n - s) * (top - bottom)
-    # Unpack nodes and weights
-    lon_nodes, lat_nodes, rad_nodes = glq_nodes[:]
-    lon_weights, lat_weights, rad_weights = glq_weights[:]
-    # Compute effect of the tesseroid on the observation point
-    # by iterating over the location of the point masses
-    # (move the iteration along the longitudinal nodes to the bottom for
-    # optimization: reduce the number of times that the trigonometric functions
-    # are evaluated)
-    result = 0.0
-    for j, lat_node in enumerate(lat_nodes):
-        # Get the latitude of the point mass
-        latitude_p = np.radians(0.5 * (n - s) * lat_node + 0.5 * (n + s))
-        cosphi_p = np.cos(latitude_p)
-        sinphi_p = np.sin(latitude_p)
-        for k, rad_node in enumerate(rad_nodes):
-            # Get the radius of the point mass
-            radius_p = 0.5 * (top - bottom) * rad_node + 0.5 * (top + bottom)
-            # Get kappa constant for the point mass
-            kappa = radius_p ** 2 * cosphi_p
-            for i, lon_node in enumerate(lon_nodes):
-                # Get the longitude of the point mass
-                longitude_p = np.radians(0.5 * (e - w) * lon_node + 0.5 * (e + w))
-                # Compute the mass of the point mass
-                mass = (
-                    density
-                    * a_factor
-                    * kappa
-                    * lon_weights[i]
-                    * lat_weights[j]
-                    * rad_weights[k]
-                )
-                # Add effect of the current point mass to the result
-                result += mass * kernel(
-                    longitude,
-                    cosphi,
-                    sinphi,
-                    radius,
-                    longitude_p,
-                    cosphi_p,
-                    sinphi_p,
-                    radius_p,
-                )
-    return result
-
-
-def glq_nodes_weights(glq_degrees):
-    """
-    Calculate GLQ unscaled nodes and weights
-
-    Parameters
-    ----------
-    glq_degrees : list
-        List of GLQ degrees for each direction: ``longitude``, ``latitude``,
-        ``radius``.
-
-    Returns
-    -------
-    glq_nodes : list
-        Unscaled GLQ nodes for each direction: ``longitude``, ``latitude``,
-        ``radius``.
-    glq_weights : list
-        GLQ weights for each node on each direction: ``longitude``,
-        ``latitude``, ``radius``.
-    """
-    # Unpack GLQ degrees
-    lon_degree, lat_degree, rad_degree = glq_degrees[:]
-    # Get nodes coordinates and weights
-    lon_node, lon_weights = leggauss(lon_degree)
-    lat_node, lat_weights = leggauss(lat_degree)
-    rad_node, rad_weights = leggauss(rad_degree)
-    # Reorder nodes and weights
-    glq_nodes = (lon_node, lat_node, rad_node)
-    glq_weights = (lon_weights, lat_weights, rad_weights)
-    return glq_nodes, glq_weights
-
-
-@jit(nopython=True)
-def _adaptive_discretization(
-    coordinates,
-    tesseroid,
-    distance_size_ratio,
-    stack,
-    small_tesseroids,
-    radial_discretization=False,
-):
-    """
-    Perform the adaptive discretization algorithm on a tesseroid
-
-    It apply the three or two dimensional adaptive discretization algorithm on
-    a tesseroid after a single computation point.
-
-    Parameters
-    ----------
-    coordinates : array
-        Array containing ``longitude``, ``latitude`` and ``radius`` of a single
-        computation point.
-    tesseroid : array
-        Array containing the boundaries of the tesseroid.
-    distance_size_ratio : float
-        Value for the distance-size ratio. A greater value will perform more
-        discretizations.
-    stack : 2d-array
-        Array with shape ``(6, stack_size)`` that will temporarly hold the
-        small tesseroids that are not yet processed.
-        If too many discretizations will take place, increase the
-        ``stack_size``.
-    small_tesseroids : 2d-array
-        Array with shape ``(6, MAX_DISCRETIZATIONS)`` that will contain every
-        small tesseroid produced by the adaptive discretization algorithm.
-        If too many discretizations will take place, increase the
-        ``MAX_DISCRETIZATIONS``.
-    radial_discretization : bool (optional)
-        If ``True`` the three dimensional adaptive discretization will be
-        applied.
-        If ``False`` the two dimensional adaptive discretization will be
-        applied, i.e. the tesseroid will only be split on the ``longitude`` and
-        ``latitude`` directions.
-        Default ``False``.
-
-    Returns
-    -------
-    n_splits : int
-        Total number of small tesseroids generated by the algorithm.
-    """
-    # Create stack of tesseroids
-    stack[0] = tesseroid
-    stack_top = 0
-    n_splits = 0
-    while stack_top >= 0:
-        # Pop the first tesseroid from the stack
-        tesseroid = stack[stack_top]
-        stack_top -= 1
-        # Get its dimensions
-        l_lon, l_lat, l_rad = _tesseroid_dimensions(tesseroid)
-        # Get distance between computation point and center of tesseroid
-        distance = _distance_tesseroid_point(coordinates, tesseroid)
-        # Check inequality
-        n_lon, n_lat, n_rad = 1, 1, 1
-        if distance / l_lon < distance_size_ratio:
-            n_lon = 2
-        if distance / l_lat < distance_size_ratio:
-            n_lat = 2
-        if distance / l_rad < distance_size_ratio and radial_discretization:
-            n_rad = 2
-        # Apply discretization
-        if n_lon * n_lat * n_rad > 1:
-            # Raise error if stack overflow
-            # Number of tesseroids in stack = stack_top + 1
-            if (stack_top + 1) + n_lon * n_lat * n_rad > stack.shape[0]:
-                raise OverflowError("Stack Overflow. Try to increase the stack size.")
-            stack_top = _split_tesseroid(
-                tesseroid, n_lon, n_lat, n_rad, stack, stack_top
-            )
-        else:
-            # Raise error if small_tesseroids overflow
-            if n_splits + 1 > small_tesseroids.shape[0]:
-                raise OverflowError(
-                    "Exceeded maximum discretizations."
-                    + " Please increase the MAX_DISCRETIZATIONS."
-                )
-            small_tesseroids[n_splits] = tesseroid
-            n_splits += 1
-    return n_splits
-
-
-@jit(nopython=True)
-def _split_tesseroid(
-    tesseroid, n_lon, n_lat, n_rad, stack, stack_top
-):  # pylint: disable=too-many-locals
-    """
-    Split tesseroid along each dimension
-    """
-    w, e, s, n, bottom, top = tesseroid[:]
-    # Compute differential distance
-    d_lon = (e - w) / n_lon
-    d_lat = (n - s) / n_lat
-    d_rad = (top - bottom) / n_rad
-    for i in range(n_lon):
-        for j in range(n_lat):
-            for k in range(n_rad):
-                stack_top += 1
-                stack[stack_top, 0] = w + d_lon * i
-                stack[stack_top, 1] = w + d_lon * (i + 1)
-                stack[stack_top, 2] = s + d_lat * j
-                stack[stack_top, 3] = s + d_lat * (j + 1)
-                stack[stack_top, 4] = bottom + d_rad * k
-                stack[stack_top, 5] = bottom + d_rad * (k + 1)
-    return stack_top
-
-
-@jit(nopython=True)
-def _tesseroid_dimensions(tesseroid):
-    """
-    Calculate the dimensions of the tesseroid.
-    """
-    w, e, s, n, bottom, top = tesseroid[:]
-    w, e, s, n = np.radians(w), np.radians(e), np.radians(s), np.radians(n)
-    latitude_center = (n + s) / 2
-    l_lat = top * np.arccos(np.sin(n) * np.sin(s) + np.cos(n) * np.cos(s))
-    l_lon = top * np.arccos(
-        np.sin(latitude_center) ** 2 + np.cos(latitude_center) ** 2 * np.cos(e - w)
-    )
-    l_rad = top - bottom
-    return l_lon, l_lat, l_rad
-
-
-@jit(nopython=True)
-def _distance_tesseroid_point(
-    coordinates, tesseroid
-):  # pylint: disable=too-many-locals
-    """
-    Distance between a computation point and the center of a tesseroid
-    """
-    # Get center of the tesseroid
-    w, e, s, n, bottom, top = tesseroid[:]
-    longitude_p = (w + e) / 2
-    latitude_p = (s + n) / 2
-    radius_p = (bottom + top) / 2
-    # Get distance between computation point and tesseroid center
-    distance = distance_spherical(coordinates, (longitude_p, latitude_p, radius_p))
-    return distance
-
-
-def _check_tesseroids(tesseroids):  # pylint: disable=too-many-branches
-    """
-    Check if tesseroids boundaries are well defined
-
-    A valid tesseroid should have:
-        - latitudinal boundaries within the [-90, 90] degrees interval,
-        - north boundaries greater or equal than the south boundaries,
-        - radial boundaries positive or zero,
-        - top boundaries greater or equal than the bottom boundaries,
-        - longitudinal boundaries within the [-180, 360] degrees interval,
-        - longitudinal interval must not be greater than one turn around the
-          globe.
-
-    Some valid tesseroids have its west boundary greater than the east one,
-    e.g. ``(350, 10, ...)``. On these cases the ``_longitude_continuity``
-    function is applied in order to move the longitudinal coordinates to the
-    [-180, 180) interval. Any valid tesseroid should have east boundaries
-    greater than the west boundaries before or after applying longitude
-    continuity.
-
-    Parameters
-    ----------
-    tesseroids : 2d-array
-        Array containing the boundaries of the tesseroids in the following
-        order: ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``.
-        Longitudinal and latitudinal boundaries must be in degrees.
-        The array must have the following shape: (``n_tesseroids``, 6), where
-        ``n_tesseroids`` is the total number of tesseroids.
-
-    Returns
-    -------
-    tesseroids :  2d-array
-        Array containing the boundaries of the tesseroids.
-        If no longitude continuity needs to be applied, the returned array is
-        the same one as the orignal.
-        Otherwise, it's copied and its longitudinal boundaries are modified.
-    """
-    west, east, south, north, bottom, top = tuple(tesseroids[:, i] for i in range(6))
-    err_msg = "Invalid tesseroid or tesseroids. "
-    # Check if latitudinal boundaries are inside the [-90, 90] interval
-    invalid = np.logical_or(
-        np.logical_or(south < -90, south > 90), np.logical_or(north < -90, north > 90)
-    )
-    if (invalid).any():
-        err_msg += (
-            "The latitudinal boundaries must be inside the [-90, 90] "
-            + "degrees interval.\n"
-        )
-        for tess in tesseroids[invalid]:
-            err_msg += "\tInvalid tesseroid: {}\n".format(tess)
-        raise ValueError(err_msg)
-    # Check if south boundary is not greater than the corresponding north
-    # boundary
-    invalid = south > north
-    if (invalid).any():
-        err_msg += "The south boundary can't be greater than the north one.\n"
-        for tess in tesseroids[invalid]:
-            err_msg += "\tInvalid tesseroid: {}\n".format(tess)
-        raise ValueError(err_msg)
-    # Check if radial boundaries are positive or zero
-    invalid = np.logical_or(bottom < 0, top < 0)
-    if (invalid).any():
-        err_msg += "The bottom and top radii should be positive or zero.\n"
-        for tess in tesseroids[invalid]:
-            err_msg += "\tInvalid tesseroid: {}\n".format(tess)
-        raise ValueError(err_msg)
-    # Check if top boundary is not greater than the corresponding bottom
-    # boundary
-    invalid = bottom > top
-    if (invalid).any():
-        err_msg += "The bottom radius boundary can't be greater than the top one.\n"
-        for tess in tesseroids[invalid]:
-            err_msg += "\tInvalid tesseroid: {}\n".format(tess)
-        raise ValueError(err_msg)
-    # Check if longitudinal boundaries are inside the [-180, 360] interval
-    invalid = np.logical_or(
-        np.logical_or(west < -180, west > 360), np.logical_or(east < -180, east > 360)
-    )
-    if (invalid).any():
-        err_msg += (
-            "The longitudinal boundaries must be inside the [-180, 360] "
-            + "degrees interval.\n"
-        )
-        for tess in tesseroids[invalid]:
-            err_msg += "\tInvalid tesseroid: {}\n".format(tess)
-        raise ValueError(err_msg)
-    # Apply longitude continuity if w > e
-    if (west > east).any():
-        tesseroids = _longitude_continuity(tesseroids)
-        west, east, south, north, bottom, top = tuple(
-            tesseroids[:, i] for i in range(6)
-        )
-    # Check if west boundary is not greater than the corresponding east
-    # boundary, even after applying the longitude continuity
-    invalid = west > east
-    if (invalid).any():
-        err_msg += "The west boundary can't be greater than the east one.\n"
-        for tess in tesseroids[invalid]:
-            err_msg += "\tInvalid tesseroid: {}\n".format(tess)
-        raise ValueError(err_msg)
-    # Check if the longitudinal interval is not grater than one turn around the
-    # globe
-    invalid = east - west > 360
-    if (invalid).any():
-        err_msg += (
-            "The difference between east and west boundaries cannot be greater than "
-            + "one turn around the globe.\n"
-        )
-        for tess in tesseroids[invalid]:
-            err_msg += "\tInvalid tesseroid: {}\n".format(tess)
-        raise ValueError(err_msg)
-    return tesseroids
-
-
-def _check_points_outside_tesseroids(
-    coordinates, tesseroids
-):  # pylint: disable=too-many-locals
-    """
-    Check if computation points are not inside the tesseroids
-
-    Parameters
-    ----------
-    coordinates : 2d-array
-        Array containing the coordinates of the computation points in the
-        following order: ``longitude``, ``latitude`` and ``radius``.
-        Both ``longitude`` and ``latitude`` must be in degrees.
-        The array must have the following shape: (3, ``n_points``), where
-        ``n_points`` is the total number of computation points.
-    tesseroids : 2d-array
-        Array containing the boundaries of the tesseroids in the following
-        order: ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``.
-        Longitudinal and latitudinal boundaries must be in degrees.
-        The array must have the following shape: (``n_tesseroids``, 6), where
-        ``n_tesseroids`` is the total number of tesseroids.
-        This array of tesseroids must have longitude continuity and valid
-        boundaries.
-        Run ``_check_tesseroids`` before.
-    """
-    longitude, latitude, radius = coordinates[:]
-    west, east, south, north, bottom, top = tuple(tesseroids[:, i] for i in range(6))
-    # Longitudinal boundaries of the tesseroid must be compared with
-    # longitudinal coordinates of computation points when moved to
-    # [0, 360) and [-180, 180).
-    longitude_360 = longitude % 360
-    longitude_180 = ((longitude + 180) % 360) - 180
-    inside_longitude = np.logical_or(
-        np.logical_and(
-            west < longitude_360[:, np.newaxis], longitude_360[:, np.newaxis] < east
-        ),
-        np.logical_and(
-            west < longitude_180[:, np.newaxis], longitude_180[:, np.newaxis] < east
-        ),
-    )
-    inside_latitude = np.logical_and(
-        south < latitude[:, np.newaxis], latitude[:, np.newaxis] < north
-    )
-    inside_radius = np.logical_and(
-        bottom < radius[:, np.newaxis], radius[:, np.newaxis] < top
-    )
-    # Build array of booleans.
-    # The (i, j) element is True if the computation point i is inside the
-    # tesseroid j.
-    inside = inside_longitude * inside_latitude * inside_radius
-    if inside.any():
-        err_msg = (
-            "Found computation point inside tesseroid. "
-            + "Computation points must be outside of tesseroids.\n"
-        )
-        for point_i, tess_i in np.argwhere(inside):
-            err_msg += "\tComputation point '{}' found inside tesseroid '{}'\n".format(
-                coordinates[:, point_i], tesseroids[tess_i, :]
-            )
-        raise ValueError(err_msg)
-
-
-def _longitude_continuity(tesseroids):
+    dtype,
+):  # pylint: disable=too-many-locals,invalid-name,not-an-iterable
     """
-    Modify longitudinal boundaries of tesseroids to ensure longitude continuity
+    Compute gravitational field of tesseroids on computations points
 
-    Longitudinal boundaries of the tesseroids are moved to the ``[-180, 180)``
-    degrees interval in case the ``west`` boundary is numerically greater than
-    the ``east`` one.
+    Perform adaptive discretization, convert each small tesseroid to equivalent
+    point masses through GLQ and use point masses kernel functions to compute
+    the gravitational field.
 
     Parameters
     ----------
+    coordinates : tuple
+        Tuple containing the coordinates of the computation points in spherical
+        geocentric coordinate system in the following order:
+        ``longitude``, ``latitude``, ``radius``.
+        Each element of the tuple must be a 1d array.
+        Both ``longitude`` and ``latitude`` should be in degrees and ``radius``
+        in meters.
     tesseroids : 2d-array
-        Longitudinal and latitudinal boundaries must be in degrees.
         Array containing the boundaries of each tesseroid:
         ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top`` under a geocentric
         spherical coordinate system.
         The array must have the following shape: (``n_tesseroids``, 6), where
         ``n_tesseroids`` is the total number of tesseroids.
-
-    Returns
-    -------
-    tesseroids : 2d-array
-        Modified boundaries of the tesseroids.
+        All tesseroids must have valid boundary coordinates.
+        Horizontal boundaries should be in degrees while radial boundaries
+        should be in meters.
+    density : func
+        Density function of every tesseroid in SI units.
+    result : 1d-array
+        Array where the gravitational effect of each tesseroid will be added.
+    distance_size_ratio : float
+        Value of the distance size ratio.
+    radial_adaptive_discretization : bool
+        If ``False``, the adaptive discretization algorithm will split the
+        tesseroid only on the horizontal direction.
+        If ``True``, it will perform a three dimensional adaptive
+        discretization, splitting the tesseroids on every direction.
+    glq_nodes : list
+        List containing unscaled GLQ nodes.
+    glq_weights : list
+        List containing GLQ weights of the nodes.
+    kernel : func
+        Kernel function for the gravitational field of point masses.
+    dtype : data-type
+        Data type assigned to the resulting gravitational field.
     """
-    # Copy the tesseroids to avoid modifying the original tesseroids array
-    tesseroids = tesseroids.copy()
-    west, east = tesseroids[:, 0], tesseroids[:, 1]
-    tess_to_be_changed = west > east
-    east[tess_to_be_changed] = ((east[tess_to_be_changed] + 180) % 360) - 180
-    west[tess_to_be_changed] = ((west[tess_to_be_changed] + 180) % 360) - 180
-    return tesseroids
+    # Get coordinates of the observation points
+    # and precompute trigonometric functions
+    longitude, latitude, radius = coordinates[:]
+    longitude_rad = np.radians(longitude)
+    cosphi = np.cos(np.radians(latitude))
+    sinphi = np.sin(np.radians(latitude))
+    # Loop over computation points
+    for l in prange(longitude.size):
+        # Initialize arrays to perform memory allocation only once
+        stack = np.empty((STACK_SIZE, 6), dtype=dtype)
+        small_tesseroids = np.empty((MAX_DISCRETIZATIONS, 6), dtype=dtype)
+        # Loop over tesseroids
+        for m in range(tesseroids.shape[0]):
+            # Apply adaptive discretization on tesseroid
+            n_splits = _adaptive_discretization(
+                (longitude[l], latitude[l], radius[l]),
+                tesseroids[m, :],
+                distance_size_ratio,
+                stack,
+                small_tesseroids,
+                radial_adaptive_discretization,
+            )
+            # Compute effect of the tesseroid through GLQ
+            for tess_index in range(n_splits):
+                tesseroid = small_tesseroids[tess_index, :]
+                result[l] += gauss_legendre_quadrature_variable_density(
+                    longitude_rad[l],
+                    cosphi[l],
+                    sinphi[l],
+                    radius[l],
+                    tesseroid,
+                    density,
+                    glq_nodes,
+                    glq_weights,
+                    kernel,
+                )
 
 
 # Define jitted versions of the forward modelling function
@@ -760,3 +409,9 @@ jit_tesseroid_gravity_serial = jit(nopython=True)(jit_tesseroid_gravity)
 jit_tesseroid_gravity_parallel = jit(nopython=True, parallel=True)(
     jit_tesseroid_gravity
 )
+jit_tesseroid_gravity_variable_density_serial = jit(nopython=True)(
+    jit_tesseroid_gravity_variable_density
+)
+jit_tesseroid_gravity_variable_density_parallel = jit(nopython=True, parallel=True)(
+    jit_tesseroid_gravity_variable_density
+)
diff --git a/requirements.txt b/requirements.txt
index 5efe7d2..1abc0db 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,6 @@
 numpy
 pandas
+scipy
 scikit-learn
 numba
 pooch>=0.7.0


Test patch:
diff --git a/harmonica/tests/test_tesseroid.py b/harmonica/tests/test_tesseroid.py
index 740b68f..45473e6 100644
--- a/harmonica/tests/test_tesseroid.py
+++ b/harmonica/tests/test_tesseroid.py
@@ -19,15 +19,18 @@ from ..forward.tesseroid import (
     tesseroid_gravity,
     _check_tesseroids,
     _check_points_outside_tesseroids,
+    _adaptive_discretization,
+    STACK_SIZE,
+    MAX_DISCRETIZATIONS,
+)
+from ..forward._tesseroid_utils import (
     _distance_tesseroid_point,
     _tesseroid_dimensions,
     _split_tesseroid,
-    _adaptive_discretization,
     _longitude_continuity,
-    STACK_SIZE,
-    MAX_DISCRETIZATIONS,
 )
 
+
 # Define the accuracy threshold for tesseroids (0.1%) as a
 # relative error (0.001)
 ACCURACY_THRESHOLD = 1e-3
diff --git a/harmonica/tests/test_tesseroid_variable_density.py b/harmonica/tests/test_tesseroid_variable_density.py
new file mode 100644
index 0000000..650911d
--- /dev/null
+++ b/harmonica/tests/test_tesseroid_variable_density.py
@@ -0,0 +1,502 @@
+# Copyright (c) 2018 The Harmonica 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 forward modelling for tesseroids with variable density
+"""
+import pytest
+import numpy as np
+import numpy.testing as npt
+from verde import grid_coordinates
+from numba import jit
+import harmonica
+
+from .utils import run_only_with_numba
+from ..constants import GRAVITATIONAL_CONST
+from .. import tesseroid_gravity
+from ..forward._tesseroid_variable_density import (
+    straight_line,
+    maximum_absolute_diff,
+    density_minmax,
+    _density_based_discretization,
+)
+
+
+# Define the accuracy threshold for tesseroids (0.1%) as a
+# relative error (0.001)
+ACCURACY_THRESHOLD = 1e-3
+
+
+# ---------------------------------
+# Test density-based discretization
+# ---------------------------------
+
+
+@pytest.fixture(name="bottom")
+def fixture_bottom():
+    """Return a bottom boundary"""
+    bottom = 2e3
+    return bottom
+
+
+@pytest.fixture(name="top")
+def fixture_top():
+    """Return a top boundary"""
+    top = 5e3
+    return top
+
+
+@pytest.fixture(name="quadratic_params")
+def fixture_quadratic_params():
+    """Return parameters for building a quadratic density function"""
+    factor = 1e-3
+    vertex_radius = 3e3
+    vertex_density = 1900.0
+    return factor, vertex_radius, vertex_density
+
+
+@pytest.fixture(name="quadratic_density")
+def fixture_quadratic_density(quadratic_params):
+    """
+    Return a quadratic density function
+
+    .. math::
+
+        f(x) = a (x - h)^2 + k
+
+    Where :math:`a` is the ``factor``, :math:`h` is the ``vertex_radius``,
+    :math:`k` is the ``vertex_density`` and :math:`x` is the ``radius``.
+
+    ::
+
+       +------------------------------------------------+
+       |                                               *|
+       |                                              * |
+       |                                             *  |
+       |                                           **   |
+       |                                          *     |
+       |                                        **      |
+       |                                       *        |
+       |                                     **         |
+       |                                   **           |
+       |                                 **             |
+       |*                              **               |
+       | **                         ***                 |
+       |   ****                  ***                    |
+       |       ******     *******                       |
+       |             *****                              |
+       +------------------------------------------------+
+        |              |                               |
+      bottom     vertex_radius                        top
+
+    """
+    factor, vertex_radius, vertex_density = quadratic_params
+
+    @jit
+    def density(radius):
+        """Quadratic density function"""
+        return factor * (radius - vertex_radius) ** 2 + vertex_density
+
+    return density
+
+
+@pytest.fixture(name="straight_line_analytic")
+def fixture_straight_line_analytic(bottom, top, quadratic_density):
+    """
+    Return the analytic solution for the straigh line of the quadratic density
+    """
+    density_bottom, density_top = quadratic_density(bottom), quadratic_density(top)
+    slope = (density_top - density_bottom) / (top - bottom)
+
+    @jit
+    def line(radius):
+        return slope * (radius - bottom) + density_bottom
+
+    return line
+
+
+@pytest.fixture(name="max_abs_diff_analytic")
+def fixture_max_abs_diff_analytic(
+    bottom, top, quadratic_params, quadratic_density, straight_line_analytic
+):
+    r"""
+    Return analytic solution for maximum absolute difference between quadratic
+    density and the straight line
+
+    .. math:
+
+        x_m = \frac{m}{2a} + h
+
+    """
+    factor, vertex_radius, _ = quadratic_params
+    density_bottom, density_top = quadratic_density(bottom), quadratic_density(top)
+    slope = (density_top - density_bottom) / (top - bottom)
+    radius_split = 0.5 * slope / factor + vertex_radius
+    max_diff = np.abs(
+        quadratic_density(radius_split) - straight_line_analytic(radius_split)
+    )
+    return radius_split, max_diff
+
+
+@pytest.fixture(name="quadratic_density_minmax")
+def fixture_quadratic_density_minmax(top, quadratic_params, quadratic_density):
+    """
+    Return the analytic maximum and minimum value of the quadratic density
+    between top and bottom
+    """
+    _, _, vertex_density = quadratic_params
+    minimum = vertex_density
+    maximum = quadratic_density(top)
+    return minimum, maximum
+
+
+def test_straight_line(bottom, top, quadratic_density, straight_line_analytic):
+    """
+    Test the straight_line function
+    """
+    radii = np.linspace(bottom, top, 51)
+    npt.assert_allclose(
+        straight_line(radii, quadratic_density, bottom, top),
+        straight_line_analytic(radii),
+    )
+
+
+def test_max_abs_diff(bottom, top, quadratic_density, max_abs_diff_analytic):
+    """
+    Test the maximum absolute difference
+
+    Test against the sine density defined in density_sine_portion fixture.
+    The solution is analytic.
+    """
+    radius_split_expected, max_diff_expected = max_abs_diff_analytic
+    radius_split, max_diff = maximum_absolute_diff(quadratic_density, bottom, top)
+    npt.assert_allclose(radius_split_expected, radius_split)
+    npt.assert_allclose(max_diff_expected, max_diff)
+
+
+def test_density_minmax(bottom, top, quadratic_density, quadratic_density_minmax):
+    """
+    Test the density_minmax function
+    """
+    density_min, density_max = density_minmax(quadratic_density, bottom, top)
+    expected_min, expected_max = quadratic_density_minmax
+    npt.assert_allclose(density_min, expected_min)
+    npt.assert_allclose(density_max, expected_max)
+
+
+def test_density_minmax_exponential_function(bottom, top):
+    """
+    Test density_minmax with an extreme function
+
+    A previous implementation of density_minmax failed this test
+    """
+    # Define an expontential density that has a strong L shape
+    thickness = top - bottom
+    b_factor = 50
+
+    @jit
+    def exponential_density(radius):
+        """
+        Create a dummy exponential density
+        """
+        density_outer, density_inner = 2500.0, 3300.0
+        a_factor = (density_inner - density_outer) / (1 - np.exp(-b_factor))
+        constant_term = density_inner - a_factor
+        return (
+            a_factor * np.exp(-b_factor * (radius - bottom) / thickness) + constant_term
+        )
+
+    density_min, density_max = density_minmax(exponential_density, bottom, top)
+    npt.assert_allclose(density_min, exponential_density(top))
+    npt.assert_allclose(density_max, exponential_density(bottom))
+
+
+def test_single_density_based_discretization(
+    bottom, top, quadratic_density, max_abs_diff_analytic
+):
+    """
+    Test the density-based discretization algorithm
+    """
+    # Define some dummy horizontal coordinates for the tesseroid
+    w, e, s, n = -3, 2, -4, 5
+    tesseroid = w, e, s, n, bottom, top
+    tesseroids = _density_based_discretization(tesseroid, quadratic_density)
+    # With the default value of DELTA=0.1, it should produce only 2 tesseroids
+    assert len(tesseroids) == 2
+    # Check the horizontal coordinates of the tesseroids
+    for tess in tesseroids:
+        for coord, original_coord in zip(tess[:4], tesseroid):
+            npt.assert_allclose(coord, original_coord)
+
+    # Check the radial coordinates
+    # ----------------------------
+    # Take the analytical solution for the radius_split
+    radius_split, _ = max_abs_diff_analytic
+    # The first tesseroid in the list should be the one in the top
+    npt.assert_allclose(tesseroids[0][-2], radius_split)
+    npt.assert_allclose(tesseroids[0][-1], top)
+    # The second tesseroid in the list should be the one in the bottom
+    npt.assert_allclose(tesseroids[1][-2], bottom)
+    npt.assert_allclose(tesseroids[1][-1], radius_split)
+
+
+def test_density_based_discret_with_delta(
+    bottom,
+    top,
+    quadratic_density,
+):  # pylint: disable=protected-access
+    """
+    Test the density-based discretization algorithm against values of DELTA
+    """
+    # Define some dummy horizontal coordinates for the tesseroid
+    w, e, s, n = -3, 2, -4, 5
+    tesseroid = w, e, s, n, bottom, top
+    # Define a collection of values for the delta ratio
+    deltas = [1e-3, 1e-2, 1e-1, 1e0]
+    # Define an empty list to count the number of splits for each delta
+    splits = []
+    for delta in deltas:
+        # Override the DELTA_RATIO global variable
+        harmonica.forward._tesseroid_variable_density.DELTA_RATIO = delta
+        # Count the splits generated by density based discretization
+        splits.append(len(_density_based_discretization(tesseroid, quadratic_density)))
+    splits = np.array(splits)
+    # Check if numbers of splits increases for lower deltas
+    assert (splits[1:] < splits[:-1]).all()
+
+
+def test_density_based_discret_linear_density():
+    """
+    Test if density-based discretization generates no splits when linear
+    density is passed
+    """
+    w, e, s, n, bottom, top = -3, 2, -4, 5, 30, 50
+    tesseroid = [w, e, s, n, bottom, top]
+
+    @jit
+    def linear_density(radius):
+        """Define a dummy linear density"""
+        return 3 * radius + 2
+
+    tesseroids = _density_based_discretization(tesseroid, linear_density)
+    assert len(tesseroids) == 1
+    npt.assert_allclose(tesseroids[0], tesseroid)
+
+
+def test_density_based_discret_constant_density():
+    """
+    Test if density-based discretization generates no splits when a constant
+    density function is passed (this should not be done IRL, pass an array of
+    floats as density instead)
+    """
+    w, e, s, n, bottom, top = -3, 2, -4, 5, 30, 50
+    tesseroid = [w, e, s, n, bottom, top]
+
+    def stupid_constant_density(radius):  # pylint: disable=unused-argument
+        """Define a dummy constant density function"""
+        return 3
+
+    tesseroids = _density_based_discretization(tesseroid, stupid_constant_density)
+    assert len(tesseroids) == 1
+    npt.assert_allclose(tesseroids[0], tesseroid)
+
+
+# ----------------------------------------------
+# Test single tesseroid against constant density
+# ----------------------------------------------
+
+
+@pytest.mark.use_numba
+@pytest.mark.parametrize("field", ("potential", "g_z"))
+def test_single_tesseroid_against_constant_density(field):
+    """
+    Test the output of a single tesseroid against one with constant density
+    """
+    # Define a single tesseroid with constant density
+    bottom, top = 5400e3, 6300e3
+    tesseroid = [-3, 3, -2, 2, bottom, top]
+    density = 2900.0
+
+    # Define a constant density
+    @jit
+    def constant_density(radius):  # pylint: disable=unused-argument
+        return density
+
+    # Define a set of observation points
+    coordinates = grid_coordinates(region=(-5, 5, -5, 5), spacing=1, extra_coords=top)
+
+    # Compare effects against the constant density implementation
+    npt.assert_allclose(
+        tesseroid_gravity(coordinates, [tesseroid], density=[density], field=field),
+        tesseroid_gravity(
+            coordinates, [tesseroid], density=constant_density, field=field
+        ),
+    )
+
+
+# ----------------------------
+# Test against spherical shell
+# ----------------------------
+
+
+def analytical_spherical_shell_linear(radius, bottom, top, slope, constant_term):
+    """
+    Analytical solutions of a spherical shell with linear density
+    """
+    constant = np.pi * GRAVITATIONAL_CONST * slope * (
+        top ** 4 - bottom ** 4
+    ) + 4 / 3.0 * np.pi * GRAVITATIONAL_CONST * constant_term * (top ** 3 - bottom ** 3)
+    potential = constant / radius
+    data = {
+        "potential": potential,
+        "g_z": 1e5 * (potential / radius),
+    }
+    return data
+
+
+def analytical_spherical_shell_exponential(
+    radius, bottom, top, a_factor, b_factor, constant_term
+):  # pylint: disable=too-many-locals
+    r"""
+    Analytical solutions of a spherical shell with exponential density
+
+    .. math :
+        \rho(r') = a e^{- b * (r' - R_1) / T} + c
+    """
+    thickness = top - bottom
+    k = b_factor / thickness
+    potential = (
+        4
+        * np.pi
+        * GRAVITATIONAL_CONST
+        * a_factor
+        / k ** 3
+        / radius
+        * (
+            ((bottom * k) ** 2 + 2 * bottom * k + 2)
+            - ((top * k) ** 2 + 2 * top * k + 2) * np.exp(-k * thickness)
+        )
+        + 4
+        / 3
+        * np.pi
+        * GRAVITATIONAL_CONST
+        * constant_term
+        * (top ** 3 - bottom ** 3)
+        / radius
+    )
+    data = {
+        "potential": potential,
+        "g_z": 1e5 * (potential / radius),
+    }
+    return data
+
+
+def build_spherical_shell(
+    bottom, top, shape=(6, 12)
+):  # pylint: disable=too-many-locals
+    """
+    Return a set of tesseroids modelling a spherical shell
+
+    Parameters
+    ----------
+    bottom : float
+        Inner radius of the spherical shell
+    top : float
+        Outer radius of the spherical shell
+    shape : tuple (n_latitude, n_longitude)
+    """
+    region = (-180, 180, -90, 90)
+    n_lat, n_lon = shape[:]
+    longitude = np.linspace(*region[:2], n_lon + 1, dtype="float64")
+    latitude = np.linspace(*region[2:], n_lat + 1, dtype="float64")
+    west, east = longitude[:-1], longitude[1:]
+    south, north = latitude[:-1], latitude[1:]
+    tesseroids = []
+    for w, e in zip(west, east):
+        for s, n in zip(south, north):
+            tesseroids.append([w, e, s, n, bottom, top])
+    return tesseroids
+
+
+@run_only_with_numba
+@pytest.mark.parametrize("field", ("potential", "g_z"))
+@pytest.mark.parametrize("thickness", (1e2, 1e3, 1e6))
+def test_spherical_shell_linear_density(field, thickness):
+    """
+    Test numerical results against analytical solution for linear density
+    """
+    # Define a spherical shell made out of tesseroids
+    top = 6371e3
+    bottom = top - thickness
+    tesseroids = build_spherical_shell(bottom, top)
+    # Define a linear density
+    density_outer, density_inner = 2500.0, 3300.0
+    slope = (density_outer - density_inner) / (top - bottom)
+    constant_term = density_outer - slope * top
+
+    @jit
+    def linear_density(radius):
+        """
+        Create a dummy linear density
+        """
+        return slope * radius + constant_term
+
+    # Create a set of observation points located on the shell outer radius
+    region = (-180, 180, -90, 90)
+    shape = (19, 13)
+    coordinates = grid_coordinates(region=region, shape=shape, extra_coords=top)
+    # Get analytic solution
+    analytic = analytical_spherical_shell_linear(
+        coordinates[-1], bottom, top, slope, constant_term
+    )
+    # Check if numerical results are accurate enough
+    npt.assert_allclose(
+        tesseroid_gravity(coordinates, tesseroids, linear_density, field=field),
+        analytic[field],
+        rtol=ACCURACY_THRESHOLD,
+    )
+
+
+@run_only_with_numba
+@pytest.mark.parametrize("field", ("potential", "g_z"))
+@pytest.mark.parametrize("thickness", (1e2, 1e3, 1e6))
+@pytest.mark.parametrize("b_factor", (5, 100))
+def test_spherical_shell_exponential_density(field, thickness, b_factor):
+    """
+    Test numerical results against analytical solution for exponential density
+    """
+    # Define a spherical shell made out of tesseroids
+    top = 6371e3
+    bottom = top - thickness
+    tesseroids = build_spherical_shell(bottom, top)
+    # Define a linear density
+    density_outer, density_inner = 2670.0, 3300.0
+    a_factor = (density_inner - density_outer) / (1 - np.exp(-b_factor))
+    constant_term = density_inner - a_factor
+
+    @jit
+    def exponential_density(radius):
+        """
+        Create a dummy exponential density
+        """
+        return (
+            a_factor * np.exp(-b_factor * (radius - bottom) / thickness) + constant_term
+        )
+
+    # Create a set of observation points located on the shell outer radius
+    region = (-180, 180, -90, 90)
+    shape = (19, 13)
+    coordinates = grid_coordinates(region=region, shape=shape, extra_coords=top)
+    # Get analytic solution
+    analytic = analytical_spherical_shell_exponential(
+        coordinates[-1], bottom, top, a_factor, b_factor, constant_term
+    )
+    # Check if numerical results are accurate enough
+    npt.assert_allclose(
+        tesseroid_gravity(coordinates, tesseroids, exponential_density, field=field),
+        analytic[field],
+        rtol=ACCURACY_THRESHOLD,
+    )