# Copyright (c) 2022 The Choclo Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
#
# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
#
"""
Gravity forward modelling functions for rectangular prisms
"""
import numpy as np
from numba import jit
from ..constants import GRAVITATIONAL_CONST
from ._kernels import (
kernel_e,
kernel_ee,
kernel_en,
kernel_eu,
kernel_n,
kernel_nn,
kernel_nu,
kernel_pot,
kernel_u,
kernel_uu,
)
[docs]@jit(nopython=True)
def gravity_pot(easting, northing, upward, prism, density):
r"""
Gravitational potential field due to a rectangular prism
Returns the gravitational potential field produced by a single rectangular
prism on a single computation point.
Parameters
----------
easting : float
Easting coordinate of the observation point. Must be in meters.
northing : float
Northing coordinate of the observation point. Must be in meters.
upward : float
Upward coordinate of the observation point. Must be in meters.
prism : 1d-array
One dimensional array containing the coordinates of the prism in the
following order: ``west``, ``east``, ``south``, ``north``, ``bottom``,
``top`` in a Cartesian coordinate system.
All coordinates should be in meters.
Returns
-------
potential : float
Gravitational potential field generated by the rectangular prism on the
observation point in :math:`\text{J}/\text{kg}`.
Notes
-----
Returns the gravitational potential field :math:`V(\mathbf{p})` on the
observation point :math:`\mathbf{p} = (x_p, y_p, z_p)` generated by
a single rectangular prism defined by its boundaries :math:`x_1, x_2, y_1,
y_2, z_1, z_2` and with a density :math:`\rho`:
.. math::
V(\mathbf{p}) =
G
\rho \,\,
\Bigg\lvert \Bigg\lvert \Bigg\lvert
k_V(x, y, z)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
where
.. math::
k_V(x, y, z) &=
x y \, \text{ln2} (z + r)
+ y z \, \text{ln2} (x + r)
+ z x \, \text{ln2} (y + r) \\
- \frac{x^2}{2} &\text{arctan2} \left( \frac{yz}{xr} \right)
- \frac{y^2}{2} \text{arctan2} \left( \frac{zx}{yr} \right)
- \frac{z^2}{2} \text{arctan2} \left( \frac{xy}{zr} \right),
.. math::
r = \sqrt{x^2 + y^2 + z^2},
and
.. math::
X_1 = x_1 - x_p \\
X_2 = x_2 - x_p \\
Y_1 = y_1 - y_p \\
Y_2 = y_2 - y_p \\
Z_1 = z_1 - z_p \\
Z_2 = z_2 - z_p
are the shifted coordinates of the prism boundaries and :math:`G` is the
Universal Gravitational Constant.
The :math:`\text{ln2}` and :math:`\text{arctan2}` functions are defined as
follows:
.. math::
\text{ln2}(x) =
\begin{cases}
0 & |x| < 10^{-10} \\
\ln (x)
\end{cases}
.. math::
\text{arctan2} \left( \frac{y}{x} \right) =
\begin{cases}
\text{arctan}\left( \frac{y}{x} \right) & x \ne 0 \\
\frac{\pi}{2} & x = 0 \quad \text{and} \quad y > 0 \\
-\frac{\pi}{2} & x = 0 \quad \text{and} \quad y < 0 \\
0 & x = 0 \quad \text{and} \quad y = 0 \\
\end{cases}
These were defined after [Fukushima2020]_ and guarantee a good accuracy on
any observation point.
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
return (
GRAVITATIONAL_CONST
* density
* _evaluate_kernel(easting, northing, upward, prism, kernel_pot)
)
[docs]@jit(nopython=True)
def gravity_e(easting, northing, upward, prism, density):
r"""
Easting component of the gravitational acceleration due to a prism
Returns the easting component of the gravitational acceleration produced by
a single rectangular prism on a single computation point.
Parameters
----------
easting : float
Easting coordinate of the observation point. Must be in meters.
northing : float
Northing coordinate of the observation point. Must be in meters.
upward : float
Upward coordinate of the observation point. Must be in meters.
prism : 1d-array
One dimensional array containing the coordinates of the prism in the
following order: ``west``, ``east``, ``south``, ``north``, ``bottom``,
``top`` in a Cartesian coordinate system.
All coordinates should be in meters.
Returns
-------
g_e : float
Easting component of the gravitational acceleration generated by the
rectangular prism on the observation point in
:math:`\text{m}/\text{s}^2`.
Notes
-----
Returns the easting component :math:`g_x(\mathbf{p})` of the gravitational
acceleration :math:`\mathbf{g}` on the observation point :math:`\mathbf{p}
= (x_p, y_p, z_p)` generated by a single rectangular prism defined by its
boundaries :math:`x_1, x_2, y_1, y_2, z_1, z_2` and with a density
:math:`\rho`:
.. math::
g_x(\mathbf{p}) =
G
\rho \,\,
\Bigg\lvert \Bigg\lvert \Bigg\lvert
k_x(x, y, z)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
where
.. math::
k_x(x, y, z) =
\left[
y \, \text{ln2} (z + r)
+ z \, \text{ln2} (y + r)
- x \, \text{arctan2} \left( \frac{yz}{xr} \right)
\right]
.. math::
r = \sqrt{x^2 + y^2 + z^2},
and
.. math::
X_1 = x_1 - x_p \\
X_2 = x_2 - x_p \\
Y_1 = y_1 - y_p \\
Y_2 = y_2 - y_p \\
Z_1 = z_1 - z_p \\
Z_2 = z_2 - z_p
are the shifted coordinates of the prism boundaries and :math:`G` is the
Universal Gravitational Constant.
The :math:`\text{ln2}` and :math:`\text{arctan2}` functions are defined as
follows:
.. math::
\text{ln2}(x) =
\begin{cases}
0 & |x| < 10^{-10} \\
\ln (x)
\end{cases}
.. math::
\text{arctan2} \left( \frac{y}{x} \right) =
\begin{cases}
\text{arctan}\left( \frac{y}{x} \right) & x \ne 0 \\
\frac{\pi}{2} & x = 0 \quad \text{and} \quad y > 0 \\
-\frac{\pi}{2} & x = 0 \quad \text{and} \quad y < 0 \\
0 & x = 0 \quad \text{and} \quad y = 0 \\
\end{cases}
These were defined after [Fukushima2020]_ and guarantee a good accuracy on
any observation point.
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
return (
GRAVITATIONAL_CONST
* density
* _evaluate_kernel(easting, northing, upward, prism, kernel_e)
)
[docs]@jit(nopython=True)
def gravity_n(easting, northing, upward, prism, density):
r"""
Northing component of the gravitational acceleration due to a prism
Returns the northing component of the gravitational acceleration produced
by a single rectangular prism on a single computation point.
Parameters
----------
easting : float
Easting coordinate of the observation point. Must be in meters.
northing : float
Northing coordinate of the observation point. Must be in meters.
upward : float
Upward coordinate of the observation point. Must be in meters.
prism : 1d-array
One dimensional array containing the coordinates of the prism in the
following order: ``west``, ``east``, ``south``, ``north``, ``bottom``,
``top`` in a Cartesian coordinate system.
All coordinates should be in meters.
Returns
-------
g_n : float
Northing component of the gravitational acceleration generated by the
rectangular prism on the observation point in
:math:`\text{m}/\text{s}^2`.
Notes
-----
Returns the northing component :math:`g_y(\mathbf{p})` of the gravitational
acceleration :math:`\mathbf{g}` on the observation point :math:`\mathbf{p}
= (x_p, y_p, z_p)` generated by a single rectangular prism defined by its
boundaries :math:`x_1, x_2, y_1, y_2, z_1, z_2` and with a density
:math:`\rho`:
.. math::
g_y(\mathbf{p}) =
G
\rho \,\,
\Bigg\lvert \Bigg\lvert \Bigg\lvert
k_y(x, y, z)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
where
.. math::
k_y(x, y, z) =
\left[
z \, \text{ln2} (x + r)
+ x \, \text{ln2} (z + r)
- y \, \text{arctan2} \left( \frac{zx}{yr} \right)
\right]
.. math::
r = \sqrt{x^2 + y^2 + z^2},
and
.. math::
X_1 = x_1 - x_p \\
X_2 = x_2 - x_p \\
Y_1 = y_1 - y_p \\
Y_2 = y_2 - y_p \\
Z_1 = z_1 - z_p \\
Z_2 = z_2 - z_p
are the shifted coordinates of the prism boundaries and :math:`G` is the
Universal Gravitational Constant.
The :math:`\text{ln2}` and :math:`\text{arctan2}` functions are defined as
follows:
.. math::
\text{ln2}(x) =
\begin{cases}
0 & |x| < 10^{-10} \\
\ln (x)
\end{cases}
.. math::
\text{arctan2} \left( \frac{y}{x} \right) =
\begin{cases}
\text{arctan}\left( \frac{y}{x} \right) & x \ne 0 \\
\frac{\pi}{2} & x = 0 \quad \text{and} \quad y > 0 \\
-\frac{\pi}{2} & x = 0 \quad \text{and} \quad y < 0 \\
0 & x = 0 \quad \text{and} \quad y = 0 \\
\end{cases}
These were defined after [Fukushima2020]_ and guarantee a good accuracy on
any observation point.
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
return (
GRAVITATIONAL_CONST
* density
* _evaluate_kernel(easting, northing, upward, prism, kernel_n)
)
[docs]@jit(nopython=True)
def gravity_u(easting, northing, upward, prism, density):
r"""
Upward component of the gravitational acceleration due to a prism
Returns the upward component of the gravitational acceleration produced by
a single rectangular prism on a single computation point.
Parameters
----------
easting : float
Easting coordinate of the observation point. Must be in meters.
northing : float
Northing coordinate of the observation point. Must be in meters.
upward : float
Upward coordinate of the observation point. Must be in meters.
prism : 1d-array
One dimensional array containing the coordinates of the prism in the
following order: ``west``, ``east``, ``south``, ``north``, ``bottom``,
``top`` in a Cartesian coordinate system.
All coordinates should be in meters.
Returns
-------
g_u : float
Upward component of the gravitational acceleration generated by the
rectangular prism on the observation point in
:math:`\text{m}/\text{s}^2`.
Notes
-----
Returns the upward component :math:`g_z(\mathbf{p})` of the gravitational
acceleration :math:`\mathbf{g}` on the observation point :math:`\mathbf{p}
= (x_p, y_p, z_p)` generated by a single rectangular prism defined by its
boundaries :math:`x_1, x_2, y_1, y_2, z_1, z_2` and with a density
:math:`\rho`:
.. math::
g_z(\mathbf{p}) =
G
\rho \,\,
\Bigg\lvert \Bigg\lvert \Bigg\lvert
k_z(x, y, z)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
where
.. math::
k_z(x, y, z) =
- \left[
x \, \text{ln2} (y + r)
+ y \, \text{ln2} (x + r)
- z \, \text{arctan2} \left( \frac{xy}{zr} \right)
\right]
.. math::
r = \sqrt{x^2 + y^2 + z^2},
and
.. math::
X_1 = x_1 - x_p \\
X_2 = x_2 - x_p \\
Y_1 = y_1 - y_p \\
Y_2 = y_2 - y_p \\
Z_1 = z_1 - z_p \\
Z_2 = z_2 - z_p
are the shifted coordinates of the prism boundaries and :math:`G` is the
Universal Gravitational Constant.
.. important
The minus sign in the :math:`k_x(x, y, z)` function is needed to return
the **upward** component of the acceleration. [Nagy2000]_ and
[Nagy2002]_ equation corresponds to the *downward* coordinate.
The :math:`\text{ln2}` and :math:`\text{arctan2}` functions are defined as
follows:
.. math::
\text{ln2}(x) =
\begin{cases}
0 & |x| < 10^{-10} \\
\ln (x)
\end{cases}
.. math::
\text{arctan2} \left( \frac{y}{x} \right) =
\begin{cases}
\text{arctan}\left( \frac{y}{x} \right) & x \ne 0 \\
\frac{\pi}{2} & x = 0 \quad \text{and} \quad y > 0 \\
-\frac{\pi}{2} & x = 0 \quad \text{and} \quad y < 0 \\
0 & x = 0 \quad \text{and} \quad y = 0 \\
\end{cases}
These were defined after [Fukushima2020]_ and guarantee a good accuracy on
any observation point.
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
return (
GRAVITATIONAL_CONST
* density
* _evaluate_kernel(easting, northing, upward, prism, kernel_u)
)
[docs]@jit(nopython=True)
def gravity_ee(easting, northing, upward, prism, density):
r"""
Easting-easting component of the gravitational tensor due to a prism
Returns the easting-easting component of the gravitational tensor produced
by a single rectangular prism on a single computation point.
Parameters
----------
easting : float
Easting coordinate of the observation point. Must be in meters.
northing : float
Northing coordinate of the observation point. Must be in meters.
upward : float
Upward coordinate of the observation point. Must be in meters.
prism : 1d-array
One dimensional array containing the coordinates of the prism in the
following order: ``west``, ``east``, ``south``, ``north``, ``bottom``,
``top`` in a Cartesian coordinate system.
All coordinates should be in meters.
Returns
-------
g_ee : float
Easting-easting component of the gravitational tensor generated by the
rectangular prism on the observation point in
:math:`\text{m}/\text{s}^2`.
Notes
-----
Returns the easting-easting component :math:`g_{xx}(\mathbf{p})` of the
gravitational tensor :math:`\mathbf{T}` on the observation point
:math:`\mathbf{p} = (x_p, y_p, z_p)` generated by a single rectangular
prism defined by its boundaries :math:`x_1, x_2, y_1, y_2, z_1, z_2` and
with a density :math:`\rho`:
.. math::
g_{xx}(\mathbf{p}) =
G
\rho \,\,
\Bigg\lvert \Bigg\lvert \Bigg\lvert
k_{xx}(x, y, z)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
where
.. math::
k_{xx}(x, y, z) = - \text{arctan2} \left( \frac{yz}{xr} \right),
.. math::
r = \sqrt{x^2 + y^2 + z^2},
and
.. math::
X_1 = x_1 - x_p \\
X_2 = x_2 - x_p \\
Y_1 = y_1 - y_p \\
Y_2 = y_2 - y_p \\
Z_1 = z_1 - z_p \\
Z_2 = z_2 - z_p
are the shifted coordinates of the prism boundaries and :math:`G` is the
Universal Gravitational Constant.
The :math:`\text{arctan2}` function is defined as follows:
.. math::
\text{arctan2} \left( \frac{y}{x} \right) =
\begin{cases}
\text{arctan}\left( \frac{y}{x} \right) & x \ne 0 \\
\frac{\pi}{2} & x = 0 \quad \text{and} \quad y > 0 \\
-\frac{\pi}{2} & x = 0 \quad \text{and} \quad y < 0 \\
0 & x = 0 \quad \text{and} \quad y = 0 \\
\end{cases}
It was defined after [Fukushima2020]_ and guarantee a good accuracy on any
observation point.
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
return (
GRAVITATIONAL_CONST
* density
* _evaluate_kernel(easting, northing, upward, prism, kernel_ee)
)
[docs]@jit(nopython=True)
def gravity_nn(easting, northing, upward, prism, density):
r"""
Northing-northing component of the gravitational tensor due to a prism
Returns the northing-northing component of the gravitational tensor
produced by a single rectangular prism on a single computation point.
Parameters
----------
easting : float
Easting coordinate of the observation point. Must be in meters.
northing : float
Northing coordinate of the observation point. Must be in meters.
upward : float
Upward coordinate of the observation point. Must be in meters.
prism : 1d-array
One dimensional array containing the coordinates of the prism in the
following order: ``west``, ``east``, ``south``, ``north``, ``bottom``,
``top`` in a Cartesian coordinate system.
All coordinates should be in meters.
Returns
-------
g_nn : float
Northing-northing component of the gravitational tensor generated by
the rectangular prism on the observation point in
:math:`\text{m}/\text{s}^2`.
Notes
-----
Returns the northing-northing component :math:`g_{yy}(\mathbf{p})` of the
gravitational tensor :math:`\mathbf{T}` on the observation point
:math:`\mathbf{p} = (x_p, y_p, z_p)` generated by a single rectangular
prism defined by its boundaries :math:`x_1, x_2, y_1, y_2, z_1, z_2` and
with a density :math:`\rho`:
.. math::
g_{yy}(\mathbf{p}) =
G
\rho \,\,
\Bigg\lvert \Bigg\lvert \Bigg\lvert
k_{yy}(x, y, z)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
where
.. math::
k_{yy}(x, y, z) = - \text{arctan2} \left( \frac{zx}{yr} \right),
.. math::
r = \sqrt{x^2 + y^2 + z^2},
and
.. math::
X_1 = x_1 - x_p \\
X_2 = x_2 - x_p \\
Y_1 = y_1 - y_p \\
Y_2 = y_2 - y_p \\
Z_1 = z_1 - z_p \\
Z_2 = z_2 - z_p
are the shifted coordinates of the prism boundaries and :math:`G` is the
Universal Gravitational Constant.
The :math:`\text{arctan2}` function is defined as follows:
.. math::
\text{arctan2} \left( \frac{y}{x} \right) =
\begin{cases}
\text{arctan}\left( \frac{y}{x} \right) & x \ne 0 \\
\frac{\pi}{2} & x = 0 \quad \text{and} \quad y > 0 \\
-\frac{\pi}{2} & x = 0 \quad \text{and} \quad y < 0 \\
0 & x = 0 \quad \text{and} \quad y = 0 \\
\end{cases}
It was defined after [Fukushima2020]_ and guarantee a good accuracy on any
observation point.
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
return (
GRAVITATIONAL_CONST
* density
* _evaluate_kernel(easting, northing, upward, prism, kernel_nn)
)
[docs]@jit(nopython=True)
def gravity_uu(easting, northing, upward, prism, density):
r"""
Upward-upward component of the gravitational tensor due to a prism
Returns the northing-northing component of the gravitational tensor
produced by a single rectangular prism on a single computation point.
Parameters
----------
easting : float
Easting coordinate of the observation point. Must be in meters.
northing : float
Northing coordinate of the observation point. Must be in meters.
upward : float
Upward coordinate of the observation point. Must be in meters.
prism : 1d-array
One dimensional array containing the coordinates of the prism in the
following order: ``west``, ``east``, ``south``, ``north``, ``bottom``,
``top`` in a Cartesian coordinate system.
All coordinates should be in meters.
Returns
-------
g_uu : float
Upward-upward component of the gravitational tensor generated by
the rectangular prism on the observation point in
:math:`\text{m}/\text{s}^2`.
Notes
-----
Returns the upward-upward component :math:`g_{zz}(\mathbf{p})` of the
gravitational tensor :math:`\mathbf{T}` on the observation point
:math:`\mathbf{p} = (x_p, y_p, z_p)` generated by a single rectangular
prism defined by its boundaries :math:`x_1, x_2, y_1, y_2, z_1, z_2` and
with a density :math:`\rho`:
.. math::
g_{zz}(\mathbf{p}) =
G
\rho \,\,
\Bigg\lvert \Bigg\lvert \Bigg\lvert
k_{zz}(x, y, z)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
where
.. math::
k_{zz}(x, y, z) = - \text{arctan2} \left( \frac{xy}{zr} \right),
.. math::
r = \sqrt{x^2 + y^2 + z^2},
and
.. math::
X_1 = x_1 - x_p \\
X_2 = x_2 - x_p \\
Y_1 = y_1 - y_p \\
Y_2 = y_2 - y_p \\
Z_1 = z_1 - z_p \\
Z_2 = z_2 - z_p
are the shifted coordinates of the prism boundaries and :math:`G` is the
Universal Gravitational Constant.
The :math:`\text{arctan2}` function is defined as follows:
.. math::
\text{arctan2} \left( \frac{y}{x} \right) =
\begin{cases}
\text{arctan}\left( \frac{y}{x} \right) & x \ne 0 \\
\frac{\pi}{2} & x = 0 \quad \text{and} \quad y > 0 \\
-\frac{\pi}{2} & x = 0 \quad \text{and} \quad y < 0 \\
0 & x = 0 \quad \text{and} \quad y = 0 \\
\end{cases}
It was defined after [Fukushima2020]_ and guarantee a good accuracy on any
observation point.
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
return (
GRAVITATIONAL_CONST
* density
* _evaluate_kernel(easting, northing, upward, prism, kernel_uu)
)
[docs]@jit(nopython=True)
def gravity_en(easting, northing, upward, prism, density):
r"""
Easting-northing component of the gravitational tensor due to a prism
Returns the northing-northing component of the gravitational tensor
produced by a single rectangular prism on a single computation point.
Parameters
----------
easting : float
Easting coordinate of the observation point. Must be in meters.
northing : float
Northing coordinate of the observation point. Must be in meters.
upward : float
Upward coordinate of the observation point. Must be in meters.
prism : 1d-array
One dimensional array containing the coordinates of the prism in the
following order: ``west``, ``east``, ``south``, ``north``, ``bottom``,
``top`` in a Cartesian coordinate system.
All coordinates should be in meters.
Returns
-------
g_en : float
Easting-northing component of the gravitational tensor generated by the
rectangular prism on the observation point in
:math:`\text{m}/\text{s}^2`.
Notes
-----
Returns the easting-northing component :math:`g_{xy}(\mathbf{p})` of the
gravitational tensor :math:`\mathbf{T}` on the observation point
:math:`\mathbf{p} = (x_p, y_p, z_p)` generated by a single rectangular
prism defined by its boundaries :math:`x_1, x_2, y_1, y_2, z_1, z_2` and
with a density :math:`\rho`:
.. math::
g_{xy}(\mathbf{p}) =
G
\rho \,\,
\Bigg\lvert \Bigg\lvert \Bigg\lvert
k_{xy}(x, y, z)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
where
.. math::
k_{xy}(x, y, z) = \text{ln2} \left( z + r \right),
.. math::
r = \sqrt{x^2 + y^2 + z^2},
and
.. math::
X_1 = x_1 - x_p \\
X_2 = x_2 - x_p \\
Y_1 = y_1 - y_p \\
Y_2 = y_2 - y_p \\
Z_1 = z_1 - z_p \\
Z_2 = z_2 - z_p
are the shifted coordinates of the prism boundaries and :math:`G` is the
Universal Gravitational Constant.
The :math:`\text{ln2}` function is defined as follows:
.. math::
\text{ln2}(x) =
\begin{cases}
0 & |x| < 10^{-10} \\
\ln (x)
\end{cases}
It was defined after [Fukushima2020]_ and guarantee a good accuracy on any
observation point.
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
return (
GRAVITATIONAL_CONST
* density
* _evaluate_kernel(easting, northing, upward, prism, kernel_en)
)
[docs]@jit(nopython=True)
def gravity_eu(easting, northing, upward, prism, density):
r"""
Easting-upward component of the gravitational tensor due to a prism
Returns the easting-upward component of the gravitational tensor
produced by a single rectangular prism on a single computation point.
Parameters
----------
easting : float
Easting coordinate of the observation point. Must be in meters.
northing : float
Northing coordinate of the observation point. Must be in meters.
upward : float
Upward coordinate of the observation point. Must be in meters.
prism : 1d-array
One dimensional array containing the coordinates of the prism in the
following order: ``west``, ``east``, ``south``, ``north``, ``bottom``,
``top`` in a Cartesian coordinate system.
All coordinates should be in meters.
Returns
-------
g_eu : float
Easting-upward component of the gravitational tensor generated by the
rectangular prism on the observation point in
:math:`\text{m}/\text{s}^2`.
Notes
-----
Returns the easting-upward component :math:`g_{xz}(\mathbf{p})` of the
gravitational tensor :math:`\mathbf{T}` on the observation point
:math:`\mathbf{p} = (x_p, y_p, z_p)` generated by a single rectangular
prism defined by its boundaries :math:`x_1, x_2, y_1, y_2, z_1, z_2` and
with a density :math:`\rho`:
.. math::
g_{xz}(\mathbf{p}) =
G
\rho \,\,
\Bigg\lvert \Bigg\lvert \Bigg\lvert
k_{xz}(x, y, z)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
where
.. math::
k_{xz}(x, y, z) = \text{ln2} \left( y + r \right),
.. math::
r = \sqrt{x^2 + y^2 + z^2},
and
.. math::
X_1 = x_1 - x_p \\
X_2 = x_2 - x_p \\
Y_1 = y_1 - y_p \\
Y_2 = y_2 - y_p \\
Z_1 = z_1 - z_p \\
Z_2 = z_2 - z_p
are the shifted coordinates of the prism boundaries and :math:`G` is the
Universal Gravitational Constant.
The :math:`\text{ln2}` function is defined as follows:
.. math::
\text{ln2}(x) =
\begin{cases}
0 & |x| < 10^{-10} \\
\ln (x)
\end{cases}
It was defined after [Fukushima2020]_ and guarantee a good accuracy on any
observation point.
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
return (
GRAVITATIONAL_CONST
* density
* _evaluate_kernel(easting, northing, upward, prism, kernel_eu)
)
[docs]@jit(nopython=True)
def gravity_nu(easting, northing, upward, prism, density):
r"""
Northing-upward component of the gravitational tensor due to a prism
Returns the northing-upward component of the gravitational tensor
produced by a single rectangular prism on a single computation point.
Parameters
----------
easting : float
Easting coordinate of the observation point. Must be in meters.
northing : float
Northing coordinate of the observation point. Must be in meters.
upward : float
Upward coordinate of the observation point. Must be in meters.
prism : 1d-array
One dimensional array containing the coordinates of the prism in the
following order: ``west``, ``east``, ``south``, ``north``, ``bottom``,
``top`` in a Cartesian coordinate system.
All coordinates should be in meters.
Returns
-------
g_nu : float
Northing-upward component of the gravitational tensor generated by the
rectangular prism on the observation point in
:math:`\text{m}/\text{s}^2`.
Notes
-----
Returns the northing-upward component :math:`g_{yz}(\mathbf{p})` of the
gravitational tensor :math:`\mathbf{T}` on the observation point
:math:`\mathbf{p} = (x_p, y_p, z_p)` generated by a single rectangular
prism defined by its boundaries :math:`x_1, x_2, y_1, y_2, z_1, z_2` and
with a density :math:`\rho`:
.. math::
g_{yz}(\mathbf{p}) =
G
\rho \,\,
\Bigg\lvert \Bigg\lvert \Bigg\lvert
k_{yz}(x, y, z)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
where
.. math::
k_{yz}(x, y, z) = \text{ln2} \left( x + r \right),
.. math::
r = \sqrt{x^2 + y^2 + z^2},
and
.. math::
X_1 = x_1 - x_p \\
X_2 = x_2 - x_p \\
Y_1 = y_1 - y_p \\
Y_2 = y_2 - y_p \\
Z_1 = z_1 - z_p \\
Z_2 = z_2 - z_p
are the shifted coordinates of the prism boundaries and :math:`G` is the
Universal Gravitational Constant.
The :math:`\text{ln2}` function is defined as follows:
.. math::
\text{ln2}(x) =
\begin{cases}
0 & |x| < 10^{-10} \\
\ln (x)
\end{cases}
It was defined after [Fukushima2020]_ and guarantee a good accuracy on any
observation point.
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
return (
GRAVITATIONAL_CONST
* density
* _evaluate_kernel(easting, northing, upward, prism, kernel_nu)
)
@jit(nopython=True)
def _evaluate_kernel(easting, northing, upward, prism, kernel):
r"""
Evaluate a kernel function on every shifted vertex of a prism
Parameters
----------
easting : float
Easting coordinate of the observation point. Must be in meters.
northing : float
Northing coordinate of the observation point. Must be in meters.
upward : float
Upward coordinate of the observation point. Must be in meters.
prism : 1d-array
One dimensional array containing the coordinates of the prism in the
following order: ``west``, ``east``, ``south``, ``north``, ``bottom``,
``top`` in a Cartesian coordinate system.
All coordinates should be in meters.
kernel : callable
Kernel function that will be evaluated on each one of the shifted
vertices of the prism.
Returns
-------
result : float
Evaluation of the kernel function on each one of the vertices of the
prism.
Notes
-----
This function evaluates a numerical kernel :math:`k(x, y, z)` on each one
of the vertices of the prism:
.. math::
v(\mathbf{p}) =
\Bigg\lvert \Bigg\lvert \Bigg\lvert
k(x, y, z)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
where :math:`X_1`, :math:`X_2`, :math:`Y_1`, :math:`Y_2`, :math:`Z_1` and
:math:`Z_2` are boundaries of the rectangular prism in the *shifted
coordinates* defined by the Cartesian coordinate system with its origin
located on the observation point :math:`\mathbf{p}`.
References
----------
- [Nagy2000]_
- [Nagy2002]_
"""
# Initialize result float to zero
result = 0
# Iterate over the vertices of the prism
for i in range(2):
# Compute shifted easting coordinate
shift_east = prism[1 - i] - easting
shift_east_sq = shift_east**2
for j in range(2):
# Compute shifted northing coordinate
shift_north = prism[3 - j] - northing
shift_north_sq = shift_north**2
for k in range(2):
# Compute shifted upward coordinate
shift_upward = prism[5 - k] - upward
shift_upward_sq = shift_upward**2
# Compute the radius
radius = np.sqrt(shift_east_sq + shift_north_sq + shift_upward_sq)
# If i, j or k is 1, the corresponding shifted
# coordinate will refer to the lower boundary,
# meaning the corresponding term should have a minus
# sign.
result += (-1) ** (i + j + k) * kernel(
shift_east, shift_north, shift_upward, radius
)
return result