harmonica.tesseroid_gravity

harmonica.tesseroid_gravity#

harmonica.tesseroid_gravity(coordinates, tesseroids, density, field, parallel=True, radial_adaptive_discretization=False, dtype=<class 'numpy.float64'>, progressbar=False, disable_checks=False)[source]#

Compute gravitational field of tesseroids on computation points.

Warning

The g_z field returns the downward component of the gravitational acceleration on the local North oriented coordinate system. It is equivalent to the opposite of the radial component, therefore it’s positive if the acceleration vector points inside the spheroid.

Important

  • The gravitational potential is returned in \(\text{J}/\text{kg}\).

  • The gravity acceleration components are returned in mGal (\(\text{m}/\text{s}^2\)).

Parameters:
coordinateslist of arrays

List of arrays containing the longitude, latitude and radius coordinates of the computation points, defined on a spherical geocentric coordinate system. Both longitude and latitude should be in degrees and radius in meters.

tesseroidslist or 2d-array

List or array containing the coordinates of the tesseroid: w, e, s, n, bottom, top under a geocentric spherical coordinate system. The longitudinal and latitudinal boundaries should be in degrees, while the radial ones must be in meters.

densitylist, array or func decorated with 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 numba.jit.

fieldstr

Gravitational field that wants to be computed. The available fields are:

  • Gravitational potential: potential

  • Downward acceleration: g_z

parallelbool (optional)

If True the computations will run in parallel using Numba built-in parallelization. If False, the forward model will run on a single core. Might be useful to disable parallelization if the forward model is run by an already parallelized workflow. Default to True.

radial_adaptive_discretizationbool (optional)

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. Default False.

dtypedata-type (optional)

Data type assigned to the resulting gravitational field. Default to np.float64.

progressbarbool (optional)

If True, a progress bar of the computation will be printed to standard error (stderr). Requires numba_progress to be installed. Default to False.

disable_checksbool (optional)

Flag that controls whether to perform a sanity check on the model. Should be set to True only when it is certain that the input model is valid and it does not need to be checked. Default to False.

Returns:
resultarray

Gravitational field generated by the tesseroids on the computation points. Gravitational potential is returned in \(\text{J}/\text{kg}\) and acceleration components in mGal.

References

[Soler2019]

Examples

>>> # Get WGS84 ellipsoid from the Boule package
>>> import boule
>>> ellipsoid = boule.WGS84
>>> # Define tesseroid of 1km of thickness with top surface on the mean
>>> # Earth radius
>>> thickness = 1000
>>> top = ellipsoid.mean_radius
>>> bottom = top - thickness
>>> w, e, s, n = -1.0, 1.0, -1.0, 1.0
>>> tesseroid = [w, e, s, n, bottom, top]
>>> # Set a density of 2670 kg/m^3
>>> density = 2670.0
>>> # Define computation point located on the top surface of the tesseroid
>>> coordinates = [0, 0, ellipsoid.mean_radius]
>>> # Compute radial component of the gravitational gradient in mGal
>>> g_z = tesseroid_gravity(coordinates, tesseroid, density, field="g_z")

>>> # Define a linear density function for the same tesseroid.
>>> # It should be decorated with numba.njit
>>> from numba import jit
>>> @jit(nopython=True)
... 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
>>> g_z = tesseroid_gravity(
...     coordinates, tesseroid, linear_density, field="g_z"
... )

Examples using harmonica.tesseroid_gravity#

Tesseroids (spherical prisms)

Tesseroids (spherical prisms)

Tesseroids with variable density

Tesseroids with variable density
Contents