# Source code for choclo.prism._magnetic

# Copyright (c) 2022 The Choclo Developers.
#
# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
#
"""
Magnetic forward modelling functions for rectangular prisms
"""
import numpy as np
from numba import jit

from ..constants import VACUUM_MAGNETIC_PERMEABILITY
from ._kernels import kernel_ee, kernel_en, kernel_eu, kernel_nn, kernel_nu, kernel_uu

[docs]@jit(nopython=True)
def magnetic_field(easting, northing, upward, prism, magnetization):
r"""
Magnetic field due to a rectangular prism

Returns the three components of the magnetic field due to a single
rectangular prism on a single computation point.

.. note::

Use this function when all the three component of the magnetic fields
are needed. Running this function is faster than computing each
component separately. Use one of :func:magnetic_e,
:func:magnetic_n, :func:magnetic_u if you need only one of them.

Parameters
----------
easting : float
Easting coordinate of the observation point. Must be in meters.
northing : float
Northing coordinate of the observation point. Must be in meters.
upward : float
Upward coordinate of the observation point. Must be in meters.
prism : 1d-array
One dimensional array containing the coordinates of the prism in the
following order: west, east, south, north, bottom,
top in a Cartesian coordinate system.
All coordinates should be in meters.
magnetization : 1d-array
Magnetization vector of the prism. It should have three components in
the following order: magnetization_easting,
magnetization_northing, magnetization_upward.
Should be in :math:A m^{-1}.

Returns
-------
b_e, b_n, b_u : floats
The three components of the magnetic field generated by the prism on
the observation point in :math:\text{T}.
The components are returned in the following order: b_e, b_n,
b_u.

Notes
-----
Consider an observation point :math:\mathbf{p} and a prism :math:R with
a magnetization vector :math:\mathbf{M}. The magnetic field
:math:\mathbf{B} it generates on the observation point :math:\mathbf{p}
is defined as:

.. math::

\mathbf{B}(\mathbf{p}) =
- \frac{\mu_0}{4\pi} \nabla_\mathbf{p}
\left[
\int\limits_R \mathbf{M} \cdot \nabla_\mathbf{q}
\left( \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} \right)
dv
\right]

Since the magnetization vector is constant inside the boundaries of the
prism, we can write the easting component of :math:\mathbf{B} as:

.. math::

B_x(\mathbf{p}) =
- \frac{\mu_0}{4\pi}
\left[
M_x \int\limits_R
\frac{\partial}{\partial x_p}
\left[
\frac{\partial}{\partial x_q}
\left( \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} \right)
\right]
dv
+
M_y \int\limits_R
\frac{\partial}{\partial x_p}
\left[
\frac{\partial}{\partial y_q}
\left( \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} \right)
\right]
dv
+
M_z \int\limits_R
\frac{\partial}{\partial x_p}
\left[
\frac{\partial}{\partial z_q}
\left( \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} \right)
\right]
dv
\right]

where :math:M_x, :math:M_y and :math:M_z are the components of the
magnetization vector. The other components can be expressed in an analogous
way.

It can be proved that

.. math::

\frac{\partial}{\partial x_q}
\left( \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} \right)
=
- \frac{\partial}{\partial x_p}
\left( \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} \right)

and that it also holds for the two other directions. Therefore, we can
rewrite :math:B_x as:

.. math::

B_x(\mathbf{p}) =
+ \frac{\mu_0}{4\pi}
\left[
M_x
\frac{\partial^2}{\partial x_p^2}
\int\limits_R
\frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert}
dv
+
M_y
\frac{\partial^2}{\partial x_p \partial y_p}
\int\limits_R
\frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert}
dv
+
M_z
\frac{\partial^2}{\partial x_p \partial z_p}
\int\limits_R
\frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert}
dv
\right]

Solutions to each one of the integrals in the previous equation and their
second derivatives are given by [Nagy2000]_.

Following [Oliveira2015]_ we can define a symmetrical 3x3 matrix
:math:\mathbf{U} whose elements are the second derivatives of the
previous integrals, such as:

.. math::

u_{ij} =
\frac{\partial}{\partial i}
\frac{\partial}{\partial j}
\int\limits_R
\frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert}
dv

with :math:i, j \in \{x, y z\}.

We can then express the magnetic field :math:\mathbf{B}(\mathbf{p})
generated by the prism in a compact form:

.. math::

\mathbf{B}(\mathbf{p}) = \frac{\mu_0}{4\pi} \mathbf{U} \cdot \mathbf{M}

References
----------
- [Blakely1995]_
- [Oliveira2015]_
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_

--------
:func:choclo.prism.magnetic_e
:func:choclo.prism.magnetic_n
:func:choclo.prism.magnetic_u
"""
# Initialize magnetic field vector components
b_e, b_n, b_u = 0.0, 0.0, 0.0
# Iterate over the vertices of the prism
for i in range(2):
# Compute shifted easting coordinate
shift_east = prism[1 - i] - easting
shift_east_sq = shift_east**2
for j in range(2):
# Compute shifted northing coordinate
shift_north = prism[3 - j] - northing
shift_north_sq = shift_north**2
for k in range(2):
# Compute shifted upward coordinate
shift_upward = prism[5 - k] - upward
shift_upward_sq = shift_upward**2
radius = np.sqrt(shift_east_sq + shift_north_sq + shift_upward_sq)
# Compute all kernel tensor components for the current vertex
k_ee = kernel_ee(shift_east, shift_north, shift_upward, radius)
k_nn = kernel_nn(shift_east, shift_north, shift_upward, radius)
k_uu = kernel_uu(shift_east, shift_north, shift_upward, radius)
k_en = kernel_en(shift_east, shift_north, shift_upward, radius)
k_eu = kernel_eu(shift_east, shift_north, shift_upward, radius)
k_nu = kernel_nu(shift_east, shift_north, shift_upward, radius)
# Get the sign of this terms based on the current vertex
sign = (-1) ** (i + j + k)
# Compute the dot product between the kernel tensor and the
# magnetization vector of the prism
b_e += sign * (
magnetization[0] * k_ee
+ magnetization[1] * k_en
+ magnetization[2] * k_eu
)
b_n += sign * (
magnetization[0] * k_en
+ magnetization[1] * k_nn
+ magnetization[2] * k_nu
)
b_u += sign * (
magnetization[0] * k_eu
+ magnetization[1] * k_nu
+ magnetization[2] * k_uu
)
c_m = VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi
b_e *= c_m
b_n *= c_m
b_u *= c_m
return b_e, b_n, b_u

[docs]@jit(nopython=True)
def magnetic_e(easting, northing, upward, prism, magnetization):
r"""
Easting component of the magnetic field due to a prism

Returns the easting component of the magnetic field due to a single
rectangular prism on a single computation point.

Parameters
----------
easting : float
Easting coordinate of the observation point. Must be in meters.
northing : float
Northing coordinate of the observation point. Must be in meters.
upward : float
Upward coordinate of the observation point. Must be in meters.
prism : 1d-array
One dimensional array containing the coordinates of the prism in the
following order: west, east, south, north, bottom,
top in a Cartesian coordinate system.
All coordinates should be in meters.
magnetization : 1d-array
Magnetization vector of the prism. It should have three components in
the following order: magnetization_easting,
magnetization_northing, magnetization_upward.
Should be in :math:A m^{-1}.

Returns
-------
b_e : float
Easting component of the magnetic field generated by the prism
on the observation point in :math:\text{T}.

Notes
-----
Computes the easting component of the magnetic field
:math:\mathbf{B}(\mathbf{p}) generated by a rectangular prism :math:R
with a magnetization vector :math:M on the observation point
:math:\mathbf{p} as follows:

.. math::

B_x(\mathbf{p}) =
\frac{\mu_0}{4\pi}
\left( M_x u_{xx} + M_y u_{xy} + M_z u_{xz} \right)

Where :math:u_{ij} are:

.. math::

u_{ij} =
\frac{\partial}{\partial i}
\frac{\partial}{\partial j}
\int\limits_R
\frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert}
dv

with :math:i,j \in \{x, y, z\}.
Solutions of the second derivatives of these integrals are given by
[Nagy2000]_:

.. math::

u_{xx} &=
\Bigg\lvert\Bigg\lvert\Bigg\lvert
- \arctan \left( \frac{yz}{xr} \right)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
\\
u_{xy} &=
\Bigg\lvert\Bigg\lvert\Bigg\lvert
\ln (z + r)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
\\
u_{xz} &=
\Bigg\lvert\Bigg\lvert\Bigg\lvert
\ln (y + r)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}

References
----------
- [Blakely1995]_
- [Oliveira2015]_
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_

--------
:func:choclo.prism.magnetic_field
:func:choclo.prism.magnetic_n
:func:choclo.prism.magnetic_u
"""
# Initialize magnetic field vector component
b_e = 0.0
# Iterate over the vertices of the prism
for i in range(2):
# Compute shifted easting coordinate
shift_east = prism[1 - i] - easting
shift_east_sq = shift_east**2
for j in range(2):
# Compute shifted northing coordinate
shift_north = prism[3 - j] - northing
shift_north_sq = shift_north**2
for k in range(2):
# Compute shifted upward coordinate
shift_upward = prism[5 - k] - upward
shift_upward_sq = shift_upward**2
radius = np.sqrt(shift_east_sq + shift_north_sq + shift_upward_sq)
# Compute kernel tensor components for the current vertex
k_ee = kernel_ee(shift_east, shift_north, shift_upward, radius)
k_en = kernel_en(shift_east, shift_north, shift_upward, radius)
k_eu = kernel_eu(shift_east, shift_north, shift_upward, radius)
# Compute b_e using the dot product between the kernel tensor
# and the magnetization vector of the prism
b_e += (-1) ** (i + j + k) * (
magnetization[0] * k_ee
+ magnetization[1] * k_en
+ magnetization[2] * k_eu
)
return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_e

[docs]@jit(nopython=True)
def magnetic_n(easting, northing, upward, prism, magnetization):
r"""
Northing component of the magnetic field due to a prism

Returns the northing component of the magnetic field due to a single
rectangular prism on a single computation point.

Parameters
----------
easting : float
Easting coordinate of the observation point. Must be in meters.
northing : float
Northing coordinate of the observation point. Must be in meters.
upward : float
Upward coordinate of the observation point. Must be in meters.
prism : 1d-array
One dimensional array containing the coordinates of the prism in the
following order: west, east, south, north, bottom,
top in a Cartesian coordinate system.
All coordinates should be in meters.
magnetization : 1d-array
Magnetization vector of the prism. It should have three components in
the following order: magnetization_easting,
magnetization_northing, magnetization_upward.
Should be in :math:A m^{-1}.

Returns
-------
b_n : float
Northing component of the magnetic field generated by the prism
on the observation point in :math:\text{T}.

Notes
-----
Computes the northing component of the magnetic field
:math:\mathbf{B}(\mathbf{p}) generated by a rectangular prism :math:R
with a magnetization vector :math:M on the observation point
:math:\mathbf{p} as follows:

.. math::

B_y(\mathbf{p}) =
\frac{\mu_0}{4\pi}
\left( M_x u_{xy} + M_y u_{yy} + M_z u_{yz} \right)

Where :math:u_{ij} are:

.. math::

u_{ij} =
\frac{\partial}{\partial i}
\frac{\partial}{\partial j}
\int\limits_R
\frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert}
dv

with :math:i,j \in \{x, y, z\}.
Solutions of the second derivatives of these integrals are given by
[Nagy2000]_:

.. math::

u_{xy} &=
\Bigg\lvert\Bigg\lvert\Bigg\lvert
\ln (z + r)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
\\
u_{yy} &=
\Bigg\lvert\Bigg\lvert\Bigg\lvert
- \arctan \left( \frac{xz}{yr} \right)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
\\
u_{yz} &=
\Bigg\lvert\Bigg\lvert\Bigg\lvert
\ln (x + r)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}

References
----------
- [Blakely1995]_
- [Oliveira2015]_
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_

--------
:func:choclo.prism.magnetic_field
:func:choclo.prism.magnetic_e
:func:choclo.prism.magnetic_u
"""
# Initialize magnetic field vector component
b_n = 0.0
# Iterate over the vertices of the prism
for i in range(2):
# Compute shifted easting coordinate
shift_east = prism[1 - i] - easting
shift_east_sq = shift_east**2
for j in range(2):
# Compute shifted northing coordinate
shift_north = prism[3 - j] - northing
shift_north_sq = shift_north**2
for k in range(2):
# Compute shifted upward coordinate
shift_upward = prism[5 - k] - upward
shift_upward_sq = shift_upward**2
radius = np.sqrt(shift_east_sq + shift_north_sq + shift_upward_sq)
# Compute kernel tensor components for the current vertex
k_nn = kernel_nn(shift_east, shift_north, shift_upward, radius)
k_en = kernel_en(shift_east, shift_north, shift_upward, radius)
k_nu = kernel_nu(shift_east, shift_north, shift_upward, radius)
# Compute b_n using the dot product between the kernel tensor
# and the magnetization vector of the prism
b_n += (-1) ** (i + j + k) * (
magnetization[0] * k_en
+ magnetization[1] * k_nn
+ magnetization[2] * k_nu
)
return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_n

[docs]@jit(nopython=True)
def magnetic_u(easting, northing, upward, prism, magnetization):
r"""
Upward component of the magnetic field due to a prism

Returns the upward component of the magnetic field due to a single
rectangular prism on a single computation point.

Parameters
----------
easting : float
Easting coordinate of the observation point. Must be in meters.
northing : float
Northing coordinate of the observation point. Must be in meters.
upward : float
Upward coordinate of the observation point. Must be in meters.
prism : 1d-array
One dimensional array containing the coordinates of the prism in the
following order: west, east, south, north, bottom,
top in a Cartesian coordinate system.
All coordinates should be in meters.
magnetization : 1d-array
Magnetization vector of the prism. It should have three components in
the following order: magnetization_easting,
magnetization_northing, magnetization_upward.
Should be in :math:A m^{-1}.

Returns
-------
b_u : float
Upward component of the magnetic field generated by the prism
on the observation point in :math:\text{T}.

Notes
-----
Computes the upward component of the magnetic field
:math:\mathbf{B}(\mathbf{p}) generated by a rectangular prism :math:R
with a magnetization vector :math:M on the observation point
:math:\mathbf{p} as follows:

.. math::

B_z(\mathbf{p}) =
\frac{\mu_0}{4\pi}
\left( M_x u_{xz} + M_y u_{yz} + M_z u_{zz} \right)

Where :math:u_{ij} are:

.. math::

u_{ij} =
\frac{\partial}{\partial i}
\frac{\partial}{\partial j}
\int\limits_R
\frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert}
dv

with :math:i,j \in \{x, y, z\}.
Solutions of the second derivatives of these integrals are given by
[Nagy2000]_:

.. math::

u_{xz} &=
\Bigg\lvert\Bigg\lvert\Bigg\lvert
\ln (y + r)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
\\
u_{yz} &=
\Bigg\lvert\Bigg\lvert\Bigg\lvert
\ln (x + r)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}
\\
u_{zz} &=
\Bigg\lvert\Bigg\lvert\Bigg\lvert
- \arctan \left( \frac{xy}{zr} \right)
\Bigg\rvert_{X_1}^{X_2}
\Bigg\rvert_{Y_1}^{Y_2}
\Bigg\rvert_{Z_1}^{Z_2}

References
----------
- [Blakely1995]_
- [Oliveira2015]_
- [Nagy2000]_
- [Nagy2002]_
- [Fukushima2020]_

--------
:func:choclo.prism.magnetic_field
:func:choclo.prism.magnetic_e
:func:choclo.prism.magnetic_n
"""
# Initialize magnetic field vector component
b_u = 0.0
# Iterate over the vertices of the prism
for i in range(2):
# Compute shifted easting coordinate
shift_east = prism[1 - i] - easting
shift_east_sq = shift_east**2
for j in range(2):
# Compute shifted northing coordinate
shift_north = prism[3 - j] - northing
shift_north_sq = shift_north**2
for k in range(2):
# Compute shifted upward coordinate
shift_upward = prism[5 - k] - upward
shift_upward_sq = shift_upward**2
radius = np.sqrt(shift_east_sq + shift_north_sq + shift_upward_sq)
# Compute kernel tensor components for the current vertex
k_uu = kernel_uu(shift_east, shift_north, shift_upward, radius)
k_eu = kernel_eu(shift_east, shift_north, shift_upward, radius)
k_nu = kernel_nu(shift_east, shift_north, shift_upward, radius)
# Compute b_n using the dot product between the kernel tensor
# and the magnetization vector of the prism
b_u += (-1) ** (i + j + k) * (
magnetization[0] * k_eu
+ magnetization[1] * k_nu
+ magnetization[2] * k_uu
)
return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_u