# Source code for choclo.utils

# 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)
#
"""
Utility functions for the kernel and forward modelling functions
"""
import numpy as np
from numba import jit

[docs]
@jit(nopython=True)
def distance_cartesian(
easting_p, northing_p, upward_p, easting_q, northing_q, upward_q
):
r"""
Euclidean distance between two points given in Cartesian coordinates

.. warning::

All coordinates should be in the same units.

Parameters
----------
easting_p : float
Easting coordinate of point :math:\mathbf{p}.
northing_p : float
Northing coordinate of point :math:\mathbf{p}.
upward_p : float
Upward coordinate of point :math:\mathbf{p}.
easting_q : float
Easting coordinate of point :math:\mathbf{q}.
northing_q : float
Northing coordinate of point :math:\mathbf{q}.
upward_q : float
Upward coordinate of point :math:\mathbf{q}.

Returns
-------
distance : float
Euclidean distance between point_p and point_q.

Notes
-----
Given two points :math:\mathbf{p} = (x_p, y_p, z_p) and
:math:\mathbf{q} = (x_q, y_q, z_q) defined in a Cartesian coordinate
system :math:(x, y, z), return the Euclidean (L2) distance between
them:

.. math::

d = \sqrt{(x_p - x_q)^2 + (y_p - y_q)^2 + (z_p - z_q)^2}

"""
distance = np.sqrt(
(easting_p - easting_q) ** 2
+ (northing_p - northing_q) ** 2
+ (upward_p - upward_q) ** 2
)
return distance

[docs]
@jit(nopython=True)
def distance_spherical(
longitude_p, latitude_p, radius_p, longitude_q, latitude_q, radius_q
):
r"""
Euclidean distance between two points in spherical coordinates

.. important::

All angles must be in degrees and radii in meters.

Parameters
----------
longitude_p : float
Longitude coordinate of point :math:\mathbf{p} in degrees.
latitude_p : float
Latitude coordinate of point :math:\mathbf{p} in degrees.
radius_p : float
Radial coordinate of point :math:\mathbf{p} in meters.
longitude_q : float
Longitude coordinate of point :math:\mathbf{q} in degrees.
latitude_q : float
Latitude coordinate of point :math:\mathbf{q} in degrees.
radius_q : float
Radial coordinate of point :math:\mathbf{q} in meters.

Returns
-------
distance : float
Euclidean distance between point_p and point_q.

Notes
-----
Given two points :math:\mathbf{p} = (\lambda_p, \phi_p, r_p) and
:math:\mathbf{q} = (\lambda_q, \phi_q, r_q) defined in a spherical
coordinate system :math:(\lambda, \phi, r), return the Euclidean (L2)
distance between them:

.. math::

d = \sqrt{ (r_p - r_q) ^ 2 + 2 r_p r_q (1 - \cos\psi)}

where

.. math::

\cos\psi =
\sin\phi_p \sin\phi_q
+ \cos\phi_p \cos\phi_q \cos(\lambda_p - \lambda_q)

and :math:\lambda is the longitude angle, :math:\phi the spherical
latitude angle an :math:r is the radial coordinate.
"""
# Convert angles to radians
longitude_p, latitude_p = np.radians(longitude_p), np.radians(latitude_p)
longitude_q, latitude_q = np.radians(longitude_q), np.radians(latitude_q)
# Compute trigonometric quantities
cosphi_p = np.cos(latitude_p)
sinphi_p = np.sin(latitude_p)
cosphi_q = np.cos(latitude_q)
sinphi_q = np.sin(latitude_q)
distance, _, _ = distance_spherical_core(
longitude_p,
cosphi_p,
sinphi_p,
radius_p,
longitude_q,
cosphi_q,
sinphi_q,
radius_q,
)
return distance

[docs]
@jit(nopython=True)
def distance_spherical_core(
longitude_p, cosphi_p, sinphi_p, radius_p, longitude_q, cosphi_q, sinphi_q, radius_q
):
r"""
Core computation of distance between two points in spherical coordinates

.. important::

All longitudinal angles must be in degrees.

It computes the Euclidean distance between two points defined in spherical
coordinates given precomputed quantities related to the coordinates of both
points: the longitude in radians, the sine and cosine of the
latitude, and the radius in meters. Precomputing this quantities
may save computation time on some cases.

Parameters
----------
longitude_p : float
Longitude coordinate of the first point. Must be in radians.
cosphi_p : float
Cosine of the latitude coordinate of the first point.
sinphi_p : float
Sine of the latitude coordinate of the first point.
radius_p : float
Radial coordinate of the first point.
longitude_q : float
Longitude coordinate of the second point. Must be in radians.
cosphi_q : float
Cosine of the latitude coordinate of the second point.
sinphi_q : float
Sine of the latitude coordinate of the second point.
radius_q : float
Radial coordinate of the second point.

Returns
-------
distance : float
Distance between the two points.
cospsi : float
Cosine of the psi angle.
coslambda : float
Cosine of the diference between the longitudes of both points.

Notes
-----
Given two points :math:\mathbf{p} = (\lambda_p, \phi_p, r_p) and
:math:\mathbf{q} = (\lambda_q, \phi_q, r_q) defined in a spherical
coordinate system :math:(\lambda, \phi, r), return the Euclidean (L2)
distance between them:

.. math::

d = \sqrt{ (r_p - r_q) ^ 2 + 2 r_p r_q (1 - \cos\psi)},

the cosine of the :math:\psi angle:

.. math::

\cos\psi =
\sin\phi_p \sin\phi_q
+ \cos\phi_p \cos\phi_q \cos(\lambda_p - \lambda_q)

and the cosine of the difference between the longitude angles of the
observation and source points:

.. math::

\cos(\lambda_p - \lambda_q),

where :math:\lambda is the longitude angle, :math:\phi the spherical
latitude angle an :math:r is the radial coordinate.
"""
coslambda = np.cos(longitude_q - longitude_p)
cospsi = sinphi_q * sinphi_p + cosphi_q * cosphi_p * coslambda
distance = np.sqrt(
(radius_p - radius_q) ** 2 + 2 * radius_p * radius_q * (1 - cospsi)
)
return distance, cospsi, coslambda