.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "gallery/forward/tesseroid_variable_density.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        Click :ref:`here <sphx_glr_download_gallery_forward_tesseroid_variable_density.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_gallery_forward_tesseroid_variable_density.py:


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.

.. GENERATED FROM PYTHON SOURCE LINES 26-96



.. image-sg:: /gallery/forward/images/sphx_glr_tesseroid_variable_density_001.png
   :alt: tesseroid variable density
   :srcset: /gallery/forward/images/sphx_glr_tesseroid_variable_density_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none

    [[15.64931026 15.84632554 16.04629019 ... 17.73106017 17.50445005
      17.28024265]
     [15.93321711 16.14345748 16.35722954 ... 18.16928493 17.9230556
      17.68001566]
     [16.2246854  16.4491229  16.6777669  ... 18.62722145 18.3593706
      18.09566506]
     ...
     [15.52398623 15.81724318 16.12003228 ... 18.50055111 18.12075672
      17.75302149]
     [15.24656009 15.52448978 15.8108793  ... 18.06017674 17.705425
      17.36104098]
     [14.97689905 15.24046063 15.51153133 ... 17.6376896  17.30600261
      16.98321625]]
    <IPython.core.display.Image object>






|

.. code-block:: default

    import boule as bl
    import pygmt
    import verde as vd
    from numba import njit

    import harmonica as hm

    # 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)
    grid = vd.make_xarray_grid(
        coordinates, gravity, data_names="gravity", extra_coords_names="extra"
    )

    # Plot the gravitational field
    fig = pygmt.Figure()

    title = "Downward component of gravitational acceleration"

    with pygmt.config(FONT_TITLE="16p"):
        fig.grdimage(
            region=[-80, -40, -50, -10],
            projection="M-60/-30/10c",
            grid=grid.gravity,
            frame=["a", f"+t{title}"],
            cmap="viridis",
        )

    fig.colorbar(cmap=True, frame=["a200f50", "x+lmGal"])

    fig.coast(shorelines="1p,black")

    fig.show()


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 0 minutes  4.574 seconds)


.. _sphx_glr_download_gallery_forward_tesseroid_variable_density.py:


.. only :: html

 .. container:: sphx-glr-footer
    :class: sphx-glr-footer-example



  .. container:: sphx-glr-download sphx-glr-download-python

     :download:`Download Python source code: tesseroid_variable_density.py <tesseroid_variable_density.py>`



  .. container:: sphx-glr-download sphx-glr-download-jupyter

     :download:`Download Jupyter notebook: tesseroid_variable_density.ipynb <tesseroid_variable_density.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_