# Source code for choclo.prism._gravity

# Copyright (c) 2022 The Choclo Developers.
#
# 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