# 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)
#
"""
Kernel functions for rectangular prisms
"""
import numpy as np
from numba import jit
[docs]
@jit(nopython=True)
def kernel_pot(easting, northing, upward, radius):
r"""
Kernel for the potential field due to a rectangular prism
Evaluates the integration kernel for the potential field generated by
a prism [Nagy2000] on a single vertex of the prism. The coordinates that
must be passed are shifted coordinates: the coordinates of the vertex from
a Cartesian coordinate system whose origin is located in the observation
point.
This function makes use of a safe natural logarithmic function and a safe
arctangent function [Fukushima2020]_ that guarantee a good accuracy on
every observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the numerical kernel function for the potential field due to
a rectangular prism evaluated on a single vertex.
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. 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)
where
.. math::
\operatorname{safe\_ln}(x, r) =
\begin{cases}
0 & r = 0 \\
\ln(x + r) & x \ge 0 \\
\ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\
-\ln(-2 x) & x < 0, r = |x|
\end{cases}
and
.. 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}
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
kernel = (
easting * northing * _safe_log(upward, easting, northing, radius)
+ northing * upward * _safe_log(easting, northing, upward, radius)
+ easting * upward * _safe_log(northing, upward, easting, radius)
- 0.5 * easting**2 * _safe_atan2(upward * northing, easting * radius)
- 0.5 * northing**2 * _safe_atan2(upward * easting, northing * radius)
- 0.5 * upward**2 * _safe_atan2(easting * northing, upward * radius)
)
return kernel
[docs]
@jit(nopython=True)
def kernel_e(easting, northing, upward, radius):
r"""
Kernel for easting component of the gradient due to a rectangular prism
Evaluates the integration kernel for the easting component of the gradient
of the potential field generated by a prism [Nagy2000]_ on a single vertex
of the prism. The coordinates that must be passed are shifted coordinates:
the coordinates of the vertex from a Cartesian coordinate system whose
origin is located in the observation point.
This function makes use of a safe natural logarithmic function and a safe
arctangent function [Fukushima2020]_ that guarantee a good accuracy on
every observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the numerical kernel function for the easting component of the
gradient of the potential field due to a rectangular prism evaluated on
a single vertex.
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. 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]
where
.. math::
\operatorname{safe\_ln}(x, r) =
\begin{cases}
0 & r = 0 \\
\ln(x + r) & x \ge 0 \\
\ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\
-\ln(-2 x) & x < 0, r = |x|
\end{cases}
and
.. 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}
.. important::
In the first equation a minus sign has been added to the one obtained
by [Nagy2000]_ in order to compute the numerical kernel for the
**eastward** component instead for the westward one.
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
kernel = -(
northing * _safe_log(upward, easting, northing, radius)
+ upward * _safe_log(northing, upward, easting, radius)
- easting * _safe_atan2(northing * upward, easting * radius)
)
return kernel
[docs]
@jit(nopython=True)
def kernel_n(easting, northing, upward, radius):
r"""
Kernel for northing component of the gradient due to a rectangular prism
Evaluates the integration kernel for the northing component of the gradient
of the potential field generated by a prism [Nagy2000]_ on a single vertex
of the prism. The coordinates that must be passed are shifted coordinates:
the coordinates of the vertex from a Cartesian coordinate system whose
origin is located in the observation point.
This function makes use of a safe natural logarithmic function and a safe
arctangent function [Fukushima2020]_ that guarantee a good accuracy on
every observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the numerical kernel function for the northing component of
the gradient of the potential field due to a rectangular prism
evaluated on a single vertex.
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. 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]
where
.. math::
\operatorname{safe\_ln}(x, r) =
\begin{cases}
0 & r = 0 \\
\ln(x + r) & x \ge 0 \\
\ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\
-\ln(-2 x) & x < 0, r = |x|
\end{cases}
and
.. 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}
.. important::
In the first equation a minus sign has been added to the one obtained
by [Nagy2000]_ in order to compute the numerical kernel for the
**northward** component instead for the southward one.
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
kernel = -(
upward * _safe_log(easting, northing, upward, radius)
+ easting * _safe_log(upward, easting, northing, radius)
- northing * _safe_atan2(easting * upward, northing * radius)
)
return kernel
[docs]
@jit(nopython=True)
def kernel_u(easting, northing, upward, radius):
r"""
Kernel for upward component of the gradient due to a rectangular prism
Evaluates the integration kernel for the upward component of the gradient
of the potential field generated by a prism [Nagy2000]_ on a single vertex
of the prism. The coordinates that must be passed are shifted coordinates:
the coordinates of the vertex from a Cartesian coordinate system whose
origin is located in the observation point.
This function makes use of a safe natural logarithmic function and a safe
arctangent function [Fukushima2020]_ that guarantee a good accuracy on
every observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the numerical kernel function for the upward component of the
gradient of the potential field due to a rectangular prism evaluated on
a single vertex.
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. 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]
where
.. math::
\operatorname{safe\_ln}(x, r) =
\begin{cases}
0 & r = 0 \\
\ln(x + r) & x \ge 0 \\
\ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\
-\ln(-2 x) & x < 0, r = |x|
\end{cases}
and
.. 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}
.. important::
In the first equation a minus sign has been added to the one obtained
by [Nagy2000]_ in order to compute the numerical kernel for the
**upward** component instead for the downward one.
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
# The minus sign is to return the kernel for the upward component instead
# of the downward one.
kernel = -(
easting * _safe_log(northing, upward, easting, radius)
+ northing * _safe_log(easting, northing, upward, radius)
- upward * _safe_atan2(easting * northing, upward * radius)
)
return kernel
[docs]
@jit(nopython=True)
def kernel_ee(easting, northing, upward, radius):
r"""
Kernel for easting-easting component of the tensor due to a prism
Evaluates the integration kernel for the easting-easting component of the
tensor generated by a prism [Nagy2000]_ on a single vertex of the prism.
The coordinates that must be passed are shifted coordinates: the
coordinates of the vertex from a Cartesian coordinate system whose origin
is located in the observation point.
This function makes use of a safe arctangent function [Fukushima2020]_ that
guarantee a good accuracy on every observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the easting-easting component of the
tensor due to a rectangular prism evaluated on a single vertex.
Return ``np.nan`` if ``radius`` is zero.
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{xx}(x, y, z) = - \operatorname{safe-arctan} \left( yz, xr \right)
where
.. 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}
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
if radius == 0.0:
return np.nan
return -_safe_atan2(northing * upward, easting * radius)
[docs]
@jit(nopython=True)
def kernel_nn(easting, northing, upward, radius):
r"""
Kernel for northing-northing component of the tensor due to a prism
Evaluates the integration kernel for the northing-northing component of the
tensor generated by a prism [Nagy2000]_ on a single vertex of the prism.
The coordinates that must be passed are shifted coordinates: the
coordinates of the vertex from a Cartesian coordinate system whose origin
is located in the observation point.
This function makes use of a safe arctangent function [Fukushima2020]_ that
guarantee a good accuracy on every observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the northing-northing component of the
tensor due to a rectangular prism evaluated on a single vertex.
Return ``np.nan`` if ``radius`` is zero.
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{yy}(x, y, z) = - \operatorname{safe-arctan} \left( zx, yr \right)
where
.. 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}
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
if radius == 0.0:
return np.nan
return -_safe_atan2(easting * upward, northing * radius)
[docs]
@jit(nopython=True)
def kernel_uu(easting, northing, upward, radius):
r"""
Kernel for upward-upward component of the tensor due to a prism
Evaluates the integration kernel for the upward-upward component of the
tensor generated by a prism [Nagy2000]_ on a single vertex of the prism.
The coordinates that must be passed are shifted coordinates: the
coordinates of the vertex from a Cartesian coordinate system whose origin
is located in the observation point.
This function makes use of a safe arctangent function [Fukushima2020]_ that
guarantee a good accuracy on every observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the upward-upward component of the
tensor due to a rectangular prism evaluated on a single vertex.
Return ``np.nan`` if ``radius`` is zero.
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{zz}(x, y, z) = - \operatorname{safe-arctan} \left( xy, zr \right)
where
.. 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}
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
if radius == 0.0:
return np.nan
return -_safe_atan2(easting * northing, upward * radius)
[docs]
@jit(nopython=True)
def kernel_en(easting, northing, upward, radius):
r"""
Kernel for easting-northing component of the tensor due to a prism
Evaluates the integration kernel for the easting-northing component of the
tensor generated by a prism [Nagy2000]_ on a single vertex of the prism.
The coordinates that must be passed are shifted coordinates: the
coordinates of the vertex from a Cartesian coordinate system whose origin
is located in the observation point.
This function makes use of a safe natural logarithmic function
[Fukushima2020]_ that guarantee a good accuracy on every observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the easting-northing component of the
tensor due to a rectangular prism evaluated on a single vertex.
Return ``np.nan`` if ``radius`` is zero.
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{xy}(x, y, z) = \operatorname{safe\_ln} \left(z, r \right)
where
.. math::
\operatorname{safe\_ln}(x, r) =
\begin{cases}
0 & r = 0 \\
\ln(x + r) & x \ge 0 \\
\ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\
-\ln(-2 x) & x < 0, r = |x|
\end{cases}
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
if radius == 0.0:
return np.nan
return _safe_log(upward, easting, northing, radius)
[docs]
@jit(nopython=True)
def kernel_eu(easting, northing, upward, radius):
r"""
Kernel for easting-upward component of the tensor due to a prism
Evaluates the integration kernel for the easting-upward component of the
tensor generated by a prism [Nagy2000]_ on a single vertex of the prism.
The coordinates that must be passed are shifted coordinates: the
coordinates of the vertex from a Cartesian coordinate system whose origin
is located in the observation point.
This function makes use of a safe natural logarithmic function
[Fukushima2020]_ that guarantee a good accuracy on every observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the easting-upward component of the
tensor due to a rectangular prism evaluated on a single vertex.
Return ``np.nan`` if ``radius`` is zero.
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{xz}(x, y, z) = \operatorname{safe\_ln} \left(y, r \right)
where
.. math::
\operatorname{safe\_ln}(x, r) =
\begin{cases}
0 & r = 0 \\
\ln(x + r) & x \ge 0 \\
\ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\
-\ln(-2 x) & x < 0, r = |x|
\end{cases}
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
if radius == 0.0:
return np.nan
return _safe_log(northing, easting, upward, radius)
[docs]
@jit(nopython=True)
def kernel_nu(easting, northing, upward, radius):
r"""
Kernel for northing-upward component of the tensor due to a prism
Evaluates the integration kernel for the northing-upward component of the
tensor generated by a prism [Nagy2000]_ on a single vertex of the prism.
The coordinates that must be passed are shifted coordinates: the
coordinates of the vertex from a Cartesian coordinate system whose origin
is located in the observation point.
This function makes use of a safe natural logarithmic function
[Fukushima2020]_ that guarantee a good accuracy on every observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the northing-upward component of the
tensor due to a rectangular prism evaluated on a single vertex.
Return ``np.nan`` if ``radius`` is zero.
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{yz}(x, y, z) = \operatorname{safe\_ln} \left(x, r \right)
where
.. math::
\operatorname{safe\_ln}(x, r) =
\begin{cases}
0 & r = 0 \\
\ln(x + r) & x \ge 0 \\
\ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\
-\ln(-2 x) & x < 0, r = |x|
\end{cases}
References
----------
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_
"""
if radius == 0.0:
return np.nan
return _safe_log(easting, northing, upward, radius)
@jit(nopython=True)
def _safe_atan2(y, x):
r"""
Principal value of the arctangent expressed as a two variable function
This modification has to be made to the arctangent function so the
gravitational field of the prism satisfies the Poisson's equation.
Therefore, it guarantees that the fields satisfies the symmetry properties
of the prism. This modified function has been defined according to
[Fukushima2020]_.
Notes
-----
.. 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}
References
----------
- [Fukushima2020]_
"""
if x == 0:
if y > 0:
result = np.pi / 2
elif y < 0:
result = -np.pi / 2
else:
result = 0
return result
return np.arctan(y / x)
@jit(nopython=True)
def _safe_log(x, y, z, r):
r"""
Safe log function to use in the prism kernels.
Evaluates the :math:`\ln(x + r)` where :math:`x` is one of the shifted
coordinate of the prism vertex and :math:`r` is the Euclidean distance
(always non-negative) from the vertex to the observation point.
Parameters
----------
x : float
Shifted coordinate of the prism vertex that will be used for evaluating
the log.
y, z : float
The other two shifted coordinates of the prism vertex. They are used to
determine if ``abs(x) == r`` with more accuracy, and to evaluate the
log function to avoid floating point errors when :math:`x < 0` and
:math:`r \ne |x|`.
r : float
Euclidean distance from the vertex to the observation point. We require
this argument because this quantity is usually precomputed before
evaluating the kernel, thus saving computation time.
Returns
-------
float
Evaluation of the :math:`\ln{x + r}` to be used on kernels.
Notes
-----
The :math:`\text{safe\_ln}(x, r)` function is defined as:
.. math::
\text{safe\_ln}(x, r) =
\begin{cases}
0 & r = 0 \\
\ln(x + r) & x \ge 0 \\
\ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\
-\ln(-2 x) & x < 0, r = |x|
\end{cases}
The function evaluates :math:`\ln(x + r)` when :math:`x \ge 0`. Otherwise,
the returned value takes into account the limit cases:
* If :math:`r = 0` (observation point is on one of the vertices of the
prism), it returns zero. This accounts for the limits of the
zeroth-order and first-order kernels when evaluated on the vertices of
the prism. The second-order kernels are not defined on the vertices.
* If :math:`x < 0, r \ne |x|`, then it evaluates
:math:`\ln((y^2 + z^2) / (r - x))` to reduce floating point errors
([Fukushima2020]_).
* If :math:`x < 0, r = |x|`, then it returns one of the terms of the limit
of the second-order kernels (when the observation point is inline with
two nodes): :math:`-\ln(-2x)`.
When checking if :math:`r = |x|`, we will evaluate if ``y == 0.0 and z ==
0.0``, to avoid issues with floating point errors in cases in which
:math:`|x| \gg |y|` and/or :math:`|x| \gg |z|`. In such cases, the ``r``
value could be exactly the same as ``x`` (up to machine precision) and
not reflect that ``y`` and/or ``z`` are not-null.
This modified version of the log function was inspired by [Nagy2000]_ and
[Fukushima2020]_.
References
----------
- [Nagy2000]_
- [Fukushima2020]_
"""
if r == 0:
return 0
if x < 0:
if y == 0.0 and z == 0.0:
result = -np.log(-2 * x)
else:
result = np.log((y**2 + z**2) / (r - x))
return result
return np.log(x + r)
[docs]
@jit(nopython=True)
def kernel_eee(easting, northing, upward, radius):
r"""
Easting-easting-easting kernel due to a prism.
Evaluates the integration kernel for the easting-easting-easting component
of the grad tensor generated by a prism [Nagy2000]_ on a single vertex of
the prism. The coordinates that must be passed are shifted coordinates: the
coordinates of the vertex from a Cartesian coordinate system whose origin
is located in the observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the easting-easting-easting component
of the grad tensor due to a rectangular prism evaluated on a single
vertex. Return ``np.nan`` if ``radius`` is zero (i.e. observation point
is on the vertex).
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{xxx}(x, y, z) =
\frac{
- y z (2 x^2 + y^2 + z^2)
}{
(x^2 + y^2)(x^2 + z^2)
\sqrt{x^2 + y^2 + z^2}
}
Examples
--------
>>> x, y, z = 3.1, 5.2, -3.0
>>> r = np.sqrt(x**2 + y**2 + z**2)
>>> float(kernel_eee(x, y, z, r)) # doctest: +NUMBER
0.18706595
"""
return _kernel_iii(easting, northing, upward, radius)
[docs]
@jit(nopython=True)
def kernel_nnn(easting, northing, upward, radius):
r"""
Northing-northing-northing kernel due to a prism.
Evaluates the integration kernel for the northing-northing-northing
component of the grad tensor generated by a prism [Nagy2000]_ on a single
vertex of the prism. The coordinates that must be passed are shifted
coordinates: the coordinates of the vertex from a Cartesian coordinate
system whose origin is located in the observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the northing-northing-northing
component of the grad tensor due to a rectangular prism evaluated on
a single vertex. Return ``np.nan`` if ``radius`` is zero (i.e.
observation point is on the vertex).
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{yyy}(x, y, z) =
\frac{
- x z (x^2 + 2 y^2 + z^2)
}{
(x^2 + y^2)(y^2 + z^2)
\sqrt{x^2 + y^2 + z^2}
}
Examples
--------
>>> x, y, z = 3.1, 5.2, -3.0
>>> r = np.sqrt(x**2 + y**2 + z**2)
>>> float(kernel_nnn(x, y, z, r)) # doctest: +NUMBER
0.07574927
"""
return _kernel_iii(northing, upward, easting, radius)
[docs]
@jit(nopython=True)
def kernel_uuu(easting, northing, upward, radius):
r"""
Upward-upward-upward kernel due to a prism.
Evaluates the integration kernel for the upward-upward-upward
component of the grad tensor generated by a prism [Nagy2000]_ on a single
vertex of the prism. The coordinates that must be passed are shifted
coordinates: the coordinates of the vertex from a Cartesian coordinate
system whose origin is located in the observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the upward-upward-upward component of
the grad tensor due to a rectangular prism evaluated on a single
vertex. Return ``np.nan`` if ``radius`` is zero (i.e. observation point
is on the vertex).
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{zzz}(x, y, z) =
\frac{
- x y (x^2 + y^2 + 2 z^2)
}{
(x^2 + z^2)(y^2 + z^2)
\sqrt{x^2 + y^2 + z^2}
}
Examples
--------
>>> x, y, z = 3.1, 5.2, -3.0
>>> r = np.sqrt(x**2 + y**2 + z**2)
>>> float(kernel_uuu(x, y, z, r)) # doctest: +NUMBER
-0.19440331
"""
return _kernel_iii(upward, northing, easting, radius)
[docs]
@jit(nopython=True)
def kernel_een(easting, northing, upward, radius):
r"""
Easting-easting-northing kernel due to a prism.
Evaluates the integration kernel for the easting-easting-northing
component of the grad tensor generated by a prism [Nagy2000]_ on a single
vertex of the prism. The coordinates that must be passed are shifted
coordinates: the coordinates of the vertex from a Cartesian coordinate
system whose origin is located in the observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the easting-easting-northing component
of the grad tensor due to a rectangular prism evaluated on a single
vertex. Return ``np.nan`` if ``radius`` is zero (i.e. observation point
is on the vertex).
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{xxy}(x, y, z) =
\frac{
- x
}{
(z + \sqrt{x^2 + y^2 + z^2})
\sqrt{x^2 + y^2 + z^2}
}
Examples
--------
>>> x, y, z = 3.1, 5.2, -3.0
>>> r = np.sqrt(x**2 + y**2 + z**2)
>>> float(kernel_een(x, y, z, r)) # doctest: +NUMBER
-0.12214070
"""
return _kernel_iij(easting, northing, upward, radius)
[docs]
@jit(nopython=True)
def kernel_eeu(easting, northing, upward, radius):
r"""
Easting-easting-upward kernel due to a prism.
Evaluates the integration kernel for the easting-easting-upward
component of the grad tensor generated by a prism [Nagy2000]_ on a single
vertex of the prism. The coordinates that must be passed are shifted
coordinates: the coordinates of the vertex from a Cartesian coordinate
system whose origin is located in the observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the easting-easting-upward component
of the grad tensor due to a rectangular prism evaluated on a single
vertex. Return ``np.nan`` if ``radius`` is zero (i.e. observation point
is on the vertex).
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{xxz}(x, y, z) =
\frac{
- x
}{
(y + \sqrt{x^2 + y^2 + z^2})
\sqrt{x^2 + y^2 + z^2}
}
Examples
--------
>>> x, y, z = 3.1, 5.2, -3.0
>>> r = np.sqrt(x**2 + y**2 + z**2)
>>> float(kernel_eeu(x, y, z, r)) # doctest: +NUMBER
-0.03837408
"""
return _kernel_iij(easting, upward, northing, radius)
[docs]
@jit(nopython=True)
def kernel_enn(easting, northing, upward, radius):
r"""
Easting-northing-northing kernel due to a prism.
Evaluates the integration kernel for the easting-northing-northing
component of the grad tensor generated by a prism [Nagy2000]_ on a single
vertex of the prism. The coordinates that must be passed are shifted
coordinates: the coordinates of the vertex from a Cartesian coordinate
system whose origin is located in the observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the easting-northing-northing
component of the grad tensor due to a rectangular prism evaluated on
a single vertex. Return ``np.nan`` if ``radius`` is zero (i.e.
observation point is on the vertex).
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{xyy}(x, y, z) =
\frac{
- y
}{
(z + \sqrt{x^2 + y^2 + z^2})
\sqrt{x^2 + y^2 + z^2}
}
Examples
--------
>>> x, y, z = 3.1, 5.2, -3.0
>>> r = np.sqrt(x**2 + y**2 + z**2)
>>> float(kernel_enn(x, y, z, r)) # doctest: +NUMBER
-0.20488118
"""
return _kernel_iij(northing, easting, upward, radius)
[docs]
@jit(nopython=True)
def kernel_nnu(easting, northing, upward, radius):
r"""
Northing-northing-upward kernel due to a prism.
Evaluates the integration kernel for the northing-northing-upward
component of the grad tensor generated by a prism [Nagy2000]_ on a single
vertex of the prism. The coordinates that must be passed are shifted
coordinates: the coordinates of the vertex from a Cartesian coordinate
system whose origin is located in the observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the northing-northing-upward
component of the grad tensor due to a rectangular prism evaluated on
a single vertex. Return ``np.nan`` if ``radius`` is zero (i.e.
observation point is on the vertex).
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{yyz}(x, y, z) =
\frac{
- y
}{
(x + \sqrt{x^2 + y^2 + z^2})
\sqrt{x^2 + y^2 + z^2}
}
Examples
--------
>>> x, y, z = 3.1, 5.2, -3.0
>>> r = np.sqrt(x**2 + y**2 + z**2)
>>> float(kernel_nnu(x, y, z, r)) # doctest: +NUMBER
-0.07808384
"""
return _kernel_iij(northing, upward, easting, radius)
[docs]
@jit(nopython=True)
def kernel_euu(easting, northing, upward, radius):
r"""
Easting-upward-upward kernel due to a prism.
Evaluates the integration kernel for the easting-upward-upward
component of the grad tensor generated by a prism [Nagy2000]_ on a single
vertex of the prism. The coordinates that must be passed are shifted
coordinates: the coordinates of the vertex from a Cartesian coordinate
system whose origin is located in the observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the easting-upward-upward
component of the grad tensor due to a rectangular prism evaluated on
a single vertex. Return ``np.nan`` if ``radius`` is zero (i.e.
observation point is on the vertex).
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{xzz}(x, y, z) =
\frac{
- z
}{
(y + \sqrt{x^2 + y^2 + z^2})
\sqrt{x^2 + y^2 + z^2}
}
Examples
--------
>>> x, y, z = 3.1, 5.2, -3.0
>>> r = np.sqrt(x**2 + y**2 + z**2)
>>> float(kernel_euu(x, y, z, r)) # doctest: +NUMBER
0.03713621
"""
return _kernel_iij(upward, easting, northing, radius)
[docs]
@jit(nopython=True)
def kernel_nuu(easting, northing, upward, radius):
r"""
Northing-upward-upward kernel due to a prism.
Evaluates the integration kernel for the northing-upward-upward
component of the grad tensor generated by a prism [Nagy2000]_ on a single
vertex of the prism. The coordinates that must be passed are shifted
coordinates: the coordinates of the vertex from a Cartesian coordinate
system whose origin is located in the observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the northing-upward-upward
component of the grad tensor due to a rectangular prism evaluated on
a single vertex. Return ``np.nan`` if ``radius`` is zero (i.e.
observation point is on the vertex).
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{yzz}(x, y, z) =
\frac{
- z
}{
(x + \sqrt{x^2 + y^2 + z^2})
\sqrt{x^2 + y^2 + z^2}
}
Examples
--------
>>> x, y, z = 3.1, 5.2, -3.0
>>> r = np.sqrt(x**2 + y**2 + z**2)
>>> float(kernel_nuu(x, y, z, r)) # doctest: +NUMBER
0.04504837
"""
return _kernel_iij(upward, northing, easting, radius)
[docs]
@jit(nopython=True)
def kernel_enu(easting, northing, upward, radius):
r"""
Easting-northing-upward kernel due to a prism.
Evaluates the integration kernel for the easting-northing-upward
component of the grad tensor generated by a prism [Nagy2000]_ on a single
vertex of the prism. The coordinates that must be passed are shifted
coordinates: the coordinates of the vertex from a Cartesian coordinate
system whose origin is located in the observation point.
Parameters
----------
easting, northing, upward : float
Shifted easting, northing and upward coordinates of the vertex of the
prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``easting``, ``northing``
and ``upward`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function for the easting-northing-upward component
of the grad tensor due to a rectangular prism evaluated on a single
vertex. Return ``np.nan`` if ``radius`` is zero (i.e. observation point
is on the vertex).
Notes
-----
Computes the following numerical kernel on the passed *shifted
coordinates*:
.. math::
k_{xyz}(x, y, z) = \frac{-1}{\sqrt{x^2 + y^2 + z^2}}
Examples
--------
>>> x, y, z = 3.1, 5.2, -3.0
>>> r = np.sqrt(x**2 + y**2 + z**2)
>>> float(kernel_enu(x, y, z, r)) # doctest: +NUMBER
-0.1480061
"""
if radius == 0.0:
return np.nan
return -1 / radius
@jit(nopython=True)
def _kernel_iii(x_i, x_j, x_k, radius):
r"""
Diagonal 3rd order kernels of a prism.
Evaluates the following integration kernel:
.. math::
k_{iii}(x_i, x_j, x_k) =
\frac{
- x_j x_k (2 x_i^2 + x_j^2 + x_k^2)
}{
(x_i^2 + x_j^2)(x_i^2 + x_k^2)
\sqrt{x_i^2 + x_j^2 + x_k^2}
}
This function allows us to evaluate the third order kernels
:math:`k_{xxx}(x, y, z)`, :math:`k_{yyy}(x, y, z)`, and :math:`k_{zzz}(x,
y, z)` by cycling through the coordinates of the nodes of the prism.
Parameters
----------
x_i, x_j, x_k : float
Shifted coordinates of the vertex of the prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``x_i``, ``x_j``, and
``x_k`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function. Return ``np.nan`` if ``radius`` is zero.
Examples
--------
Given the shifted coordinates of a prism vertex ``easting``, ``northing``
and ``upward``, we can use this function to compute the diagonal third
order kernels by cycling the order of the shifted coordinates of the
vertex:
>>> import numpy as np
>>> easting, northing, upward = 1.7, 2.8, 3.9
>>> radius = np.sqrt(easting**2 + northing**2 + upward**2)
>>> k_eee = _kernel_iii(easting, northing, upward, radius)
>>> k_nnn = _kernel_iii(northing, upward, easting, radius)
>>> k_uuu = _kernel_iii(upward, easting, northing, radius)
"""
if radius == 0.0:
return np.nan
if (x_i == 0 and x_j == 0) or (x_i == 0 and x_k == 0):
return 0.0
x_i_sq, x_j_sq, x_k_sq = x_i**2, x_j**2, x_k**2
numerator = -x_j * x_k * (2 * x_i_sq + x_j_sq + x_k_sq)
denominator = (x_i_sq + x_j_sq) * (x_i_sq + x_k_sq) * radius
return numerator / denominator
@jit(nopython=True)
def _kernel_iij(x_i, x_j, x_k, radius):
r"""
Non-diagonal 3rd order kernels of a prism.
Evaluates the following integration kernel:
.. math::
easting / ((upward + radius) * radius)
k_{iij}(x_i, x_j, x_k) =
\frac{
- x_i
}{
(x_k + \sqrt{x_i^2 + x_j^2 + x_k^2})
\sqrt{x_i^2 + x_j^2 + x_k^2}
}
This function allows us to evaluate the non-diagonal third order kernels
:math:`k_{xxy}(x, y, z)`, :math:`k_{xxz}(x, y, z)`,
:math:`k_{xyy}(x, y, z)`, :math:`k_{yyz}(x, y, z)`,
:math:`k_{xzz}(x, y, z)`, and :math:`k_{yzz}(x, y, z)` by cycling
through the coordinates of the nodes of the prism.
Parameters
----------
x_i, x_j, x_k : float
Shifted coordinates of the vertex of the prism. Must be in meters.
radius : float
Square root of the sum of the squares of the ``x_i``, ``x_j``, and
``x_k`` shifted coordinates.
Returns
-------
kernel : float
Value of the kernel function. Return ``np.nan`` if ``radius`` is zero.
Examples
--------
Given the shifted coordinates of a prism vertex ``easting``, ``northing``
and ``upward``, we can use this function to compute the non-diagonal third
order kernels by changing the order of the shifted coordinates of the
vertex:
>>> import numpy as np
>>> easting, northing, upward = 1.7, 2.8, 3.9
>>> radius = np.sqrt(easting**2 + northing**2 + upward**2)
>>> k_een = _kernel_iij(easting, northing, upward, radius)
>>> k_eeu = _kernel_iij(easting, upward, northing, radius)
>>> k_enn = _kernel_iij(northing, easting, upward, radius)
>>> k_nnu = _kernel_iij(northing, upward, easting, radius)
>>> k_euu = _kernel_iij(upward, easting, northing, radius)
>>> k_nuu = _kernel_iij(upward, northing, upward, radius)
"""
if radius == 0.0:
return np.nan
if x_i == 0 and x_j == 0:
return 0.0
return -x_i / ((x_k + radius) * radius)