# 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,
)
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
# 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