# 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,
)
from ._utils import (
is_point_on_east_face,
is_point_on_easting_edge,
is_point_on_north_face,
is_point_on_northing_edge,
is_point_on_top_face,
is_point_on_upward_edge,
)

[docs]
@jit(nopython=True)
def gravity_pot(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
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_west : float
The West boundary of the prism. Must be in meters.
prism_east : float
The East boundary of the prism. Must be in meters.
prism_south : float
The South boundary of the prism. Must be in meters.
prism_north : float
The North boundary of the prism. Must be in meters.
prism_bottom : float
The bottom boundary of the prism. Must be in meters.
prism_top : float
The top boundary of the prism. Must be in meters.
density: float
Density of the rectangular prism in kilograms per cubic meter.

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 \, \operatorname{safe-ln} (z, r)
+ y z \, \operatorname{safe-ln} (x, r)
+ z x \, \operatorname{safe-ln} (y, r) \\
- \frac{x^2}{2} &\operatorname{safe-arctan} \left( yz, xr \right)
- \frac{y^2}{2} \operatorname{safe-arctan} \left( zx, yr \right)
- \frac{z^2}{2} \operatorname{safe-arctan} \left( 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:\operatorname{safe-ln} and :math:\operatorname{safe-arctan}
functions are defined as follows:

.. math::

\operatorname{safe-ln}(x, r) =
\begin{cases}
0 & x = 0, r = 0 \\
\ln(x + r) & x \ge 0 \\
\ln((r^2 - x^2) / (r - x)) & x < 0, r \ne -x \\
-\ln(-2 x) & x < 0, r = -x
\end{cases}

.. math::

\operatorname{safe-arctan} \left( 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_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
kernel_pot,
)
)

[docs]
@jit(nopython=True)
def gravity_e(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
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_west : float
The West boundary of the prism. Must be in meters.
prism_east : float
The East boundary of the prism. Must be in meters.
prism_south : float
The South boundary of the prism. Must be in meters.
prism_north : float
The North boundary of the prism. Must be in meters.
prism_bottom : float
The bottom boundary of the prism. Must be in meters.
prism_top : float
The top boundary of the prism. Must be in meters.
density: float
Density of the rectangular prism in kilograms per cubic meter.

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 \, \operatorname{safe-ln} (z, r)
+ z \, \operatorname{safe-ln} (y, r)
- x \, \operatorname{safe-arctan} \left( 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:\operatorname{safe-ln} and :math:\operatorname{safe-arctan}
functions are defined as follows:

.. math::

\operatorname{safe-ln}(x, r) =
\begin{cases}
0 & x = 0, r = 0 \\
\ln(x + r) & x \ge 0 \\
\ln((r^2 - x^2) / (r - x)) & x < 0, r \ne -x \\
-\ln(-2 x) & x < 0, r = -x
\end{cases}

.. math::

\operatorname{safe-arctan} \left( 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_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
kernel_e,
)
)

[docs]
@jit(nopython=True)
def gravity_n(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
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_west : float
The West boundary of the prism. Must be in meters.
prism_east : float
The East boundary of the prism. Must be in meters.
prism_south : float
The South boundary of the prism. Must be in meters.
prism_north : float
The North boundary of the prism. Must be in meters.
prism_bottom : float
The bottom boundary of the prism. Must be in meters.
prism_top : float
The top boundary of the prism. Must be in meters.
density: float
Density of the rectangular prism in kilograms per cubic meter.

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 \, \operatorname{safe-ln} (x, r)
+ x \, \operatorname{safe-ln} (z, r)
- y \, \operatorname{safe-arctan} \left( 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:\operatorname{safe-ln} and :math:\operatorname{safe-arctan}
functions are defined as follows:

.. math::

\operatorname{safe-ln}(x, r) =
\begin{cases}
0 & x = 0, r = 0 \\
\ln(x + r) & x \ge 0 \\
\ln((r^2 - x^2) / (r - x)) & x < 0, r \ne -x \\
-\ln(-2 x) & x < 0, r = -x
\end{cases}

.. math::

\operatorname{safe-arctan} \left( 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_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
kernel_n,
)
)

[docs]
@jit(nopython=True)
def gravity_u(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
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_west : float
The West boundary of the prism. Must be in meters.
prism_east : float
The East boundary of the prism. Must be in meters.
prism_south : float
The South boundary of the prism. Must be in meters.
prism_north : float
The North boundary of the prism. Must be in meters.
prism_bottom : float
The bottom boundary of the prism. Must be in meters.
prism_top : float
The top boundary of the prism. Must be in meters.
density: float
Density of the rectangular prism in kilograms per cubic meter.

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 \, \operatorname{safe-ln} (y, r)
+ y \, \operatorname{safe-ln} (x, r)
- z \, \operatorname{safe-arctan} \left( 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:\operatorname{safe-ln} and :math:\operatorname{safe-arctan}
functions are defined as follows:

.. math::

\operatorname{safe-ln}(x, r) =
\begin{cases}
0 & x = 0, r = 0 \\
\ln(x + r) & x \ge 0 \\
\ln((r^2 - x^2) / (r - x)) & x < 0, r \ne -x \\
-\ln(-2 x) & x < 0, r = -x
\end{cases}

.. math::

\operatorname{safe-arctan} \left( 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_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
kernel_u,
)
)

[docs]
@jit(nopython=True)
def gravity_ee(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
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_west : float
The West boundary of the prism. Must be in meters.
prism_east : float
The East boundary of the prism. Must be in meters.
prism_south : float
The South boundary of the prism. Must be in meters.
prism_north : float
The North boundary of the prism. Must be in meters.
prism_bottom : float
The bottom boundary of the prism. Must be in meters.
prism_top : float
The top boundary of the prism. Must be in meters.
density: float
Density of the rectangular prism in kilograms per cubic meter.

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.
Return numpy.nan if the observation point falls in a singular
point: prism vertices or prism edges perpendicular to the easting
direction.

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) = - \operatorname{safe-arctan} \left( 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:\operatorname{safe-arctan} function is defined as follows:

.. math::

\operatorname{safe-arctan} \left( 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 nan if the observation point falls on a singular point.
# For the g_ee this are edges perpendicular to the easting direction
# (parallel to northing and upward)
if is_point_on_northing_edge(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
) or is_point_on_upward_edge(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
):
return np.nan
# Evaluate kernel function on each vertex of the prism
result = _evaluate_kernel(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
kernel_ee,
)
# Add 4 pi if computing on the eastern face to return the limit approaching
# from outside (approaching from the east)
if is_point_on_east_face(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
):
result += 4 * np.pi
return GRAVITATIONAL_CONST * density * result

[docs]
@jit(nopython=True)
def gravity_nn(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
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_west : float
The West boundary of the prism. Must be in meters.
prism_east : float
The East boundary of the prism. Must be in meters.
prism_south : float
The South boundary of the prism. Must be in meters.
prism_north : float
The North boundary of the prism. Must be in meters.
prism_bottom : float
The bottom boundary of the prism. Must be in meters.
prism_top : float
The top boundary of the prism. Must be in meters.
density: float
Density of the rectangular prism in kilograms per cubic meter.

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.
Return numpy.nan if the observation point falls in a singular
point: prism vertices or prism edges perpendicular to the northing
direction.

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) = - \operatorname{safe-arctan} \left( 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:\operatorname{safe-arctan} function is defined as follows:

.. math::

\operatorname{safe-arctan} \left( 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 nan if the observation point falls on a singular point.
# For the g_nn this are edges perpendicular to the northing direction
# (parallel to easting and upward)
if is_point_on_easting_edge(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
) or is_point_on_upward_edge(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
):
return np.nan
# Evaluate kernel function on each vertex of the prism
result = _evaluate_kernel(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
kernel_nn,
)
# Add 4 pi if computing on the northern face to return the limit
# approaching from outside (approaching from the north)
if is_point_on_north_face(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
):
result += 4 * np.pi
return GRAVITATIONAL_CONST * density * result

[docs]
@jit(nopython=True)
def gravity_uu(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
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_west : float
The West boundary of the prism. Must be in meters.
prism_east : float
The East boundary of the prism. Must be in meters.
prism_south : float
The South boundary of the prism. Must be in meters.
prism_north : float
The North boundary of the prism. Must be in meters.
prism_bottom : float
The bottom boundary of the prism. Must be in meters.
prism_top : float
The top boundary of the prism. Must be in meters.
density: float
Density of the rectangular prism in kilograms per cubic meter.

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.
Return numpy.nan if the observation point falls in a singular
point: prism vertices or prism edges perpendicular to the upward
direction.

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) = - \operatorname{safe-arctan} \left( 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:\operatorname{safe-arctan} function is defined as follows:

.. math::

\operatorname{safe-arctan} \left( 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 nan if the observation point falls on a singular point.
# For the g_uu this are edges perpendicular to the upward direction
# (parallel to easting and northing)
if is_point_on_easting_edge(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
) or is_point_on_northing_edge(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
):
return np.nan
# Evaluate kernel function on each vertex of the prism
result = _evaluate_kernel(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
kernel_uu,
)
# Add 4 pi if computing on the top face to return the limit approaching
# from outside (approaching from the top)
if is_point_on_top_face(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
):
result += 4 * np.pi
return GRAVITATIONAL_CONST * density * result

[docs]
@jit(nopython=True)
def gravity_en(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
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_west : float
The West boundary of the prism. Must be in meters.
prism_east : float
The East boundary of the prism. Must be in meters.
prism_south : float
The South boundary of the prism. Must be in meters.
prism_north : float
The North boundary of the prism. Must be in meters.
prism_bottom : float
The bottom boundary of the prism. Must be in meters.
prism_top : float
The top boundary of the prism. Must be in meters.
density: float
Density of the rectangular prism in kilograms per cubic meter.

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.
Return numpy.nan if the observation point falls in a singular
point: prism vertices or prism edges parallel to the upward direction.

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) = \operatorname{safe-ln} \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:\operatorname{safe-ln} function is defined as follows:

.. math::

\operatorname{safe-ln}(x, r) =
\begin{cases}
0 & x = 0, r = 0 \\
\ln(x + r) & x \ge 0 \\
\ln((r^2 - x^2) / (r - x)) & x < 0, r \ne -x \\
-\ln(-2 x) & x < 0, r = -x
\end{cases}

It was defined after [Fukushima2020]_ and guarantee a good accuracy on any
observation point.

References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
# Return nan if the observation point falls on a singular point.
# For g_en this are edges parallel to to the upward direction
if is_point_on_upward_edge(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
):
return np.nan
return (
GRAVITATIONAL_CONST
* density
* _evaluate_kernel(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
kernel_en,
)
)

[docs]
@jit(nopython=True)
def gravity_eu(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
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_west : float
The West boundary of the prism. Must be in meters.
prism_east : float
The East boundary of the prism. Must be in meters.
prism_south : float
The South boundary of the prism. Must be in meters.
prism_north : float
The North boundary of the prism. Must be in meters.
prism_bottom : float
The bottom boundary of the prism. Must be in meters.
prism_top : float
The top boundary of the prism. Must be in meters.
density: float
Density of the rectangular prism in kilograms per cubic meter.

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.
Return numpy.nan if the observation point falls in a singular
point: prism vertices or prism edges parallel to the northing
direction.

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) = \operatorname{safe-ln} \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:\operatorname{safe-ln} function is defined as follows:

.. math::

\operatorname{safe-ln}(x, r) =
\begin{cases}
0 & x = 0, r = 0 \\
\ln(x + r) & x \ge 0 \\
\ln((r^2 - x^2) / (r - x)) & x < 0, r \ne -x \\
-\ln(-2 x) & x < 0, r = -x
\end{cases}

It was defined after [Fukushima2020]_ and guarantee a good accuracy on any
observation point.

References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
# Return nan if the observation point falls on a singular point.
# For g_eu this are edges parallel to to the northing direction
if is_point_on_northing_edge(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
):
return np.nan
return (
GRAVITATIONAL_CONST
* density
* _evaluate_kernel(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
kernel_eu,
)
)

[docs]
@jit(nopython=True)
def gravity_nu(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
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_west : float
The West boundary of the prism. Must be in meters.
prism_east : float
The East boundary of the prism. Must be in meters.
prism_south : float
The South boundary of the prism. Must be in meters.
prism_north : float
The North boundary of the prism. Must be in meters.
prism_bottom : float
The bottom boundary of the prism. Must be in meters.
prism_top : float
The top boundary of the prism. Must be in meters.
density: float
Density of the rectangular prism in kilograms per cubic meter.

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.
Return numpy.nan if the observation point falls in a singular
point: prism vertices or prism edges parallel to the easting
direction.

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) = \operatorname{safe-ln} \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:\operatorname{safe-ln} function is defined as follows:

.. math::

\operatorname{safe-ln}(x, r) =
\begin{cases}
0 & x = 0, r = 0 \\
\ln(x + r) & x \ge 0 \\
\ln((r^2 - x^2) / (r - x)) & x < 0, r \ne -x \\
-\ln(-2 x) & x < 0, r = -x
\end{cases}

It was defined after [Fukushima2020]_ and guarantee a good accuracy on any
observation point.

References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
# Return nan if the observation point falls on a singular point.
# For g_nu this are edges parallel to to the easting direction
if is_point_on_easting_edge(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
):
return np.nan
return (
GRAVITATIONAL_CONST
* density
* _evaluate_kernel(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
kernel_nu,
)
)

@jit(nopython=True)
def _evaluate_kernel(
easting,
northing,
upward,
prism_west,
prism_east,
prism_south,
prism_north,
prism_bottom,
prism_top,
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_west : float
The West boundary of the prism. Must be in meters.
prism_east : float
The East boundary of the prism. Must be in meters.
prism_south : float
The South boundary of the prism. Must be in meters.
prism_north : float
The North boundary of the prism. Must be in meters.
prism_bottom : float
The bottom boundary of the prism. Must be in meters.
prism_top : float
The top boundary of the prism. Must 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
if i == 0:
shift_east = prism_east - easting
else:
shift_east = prism_west - easting
shift_east_sq = shift_east**2
for j in range(2):
# Compute shifted northing coordinate
if j == 0:
shift_north = prism_north - northing
else:
shift_north = prism_south - northing
shift_north_sq = shift_north**2
for k in range(2):
# Compute shifted upward coordinate
if k == 0:
shift_upward = prism_top - upward
else:
shift_upward = prism_bottom - upward
shift_upward_sq = shift_upward**2