From 36b36e16cb346545f01f8a95aae8df01a5d64fee Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 7 Oct 2024 14:02:10 -0700 Subject: [PATCH] Add magnetic gradiometry forwards for prisms (#97) Add forward modelling functions for the magnetic gradiometry components due to rectangular prisms: spatial derivatives of the magnetic field components. Include third-order kernel functions that are needed to compute the derivatives of the magnetic field components. Create private function to evaluate the kernels on the vertices of the prism provided its magnetization so we can reduce the amount of repeated code in the magnetic forward modelling functions. Add tests for the new functions. --- choclo/prism/__init__.py | 23 +- choclo/prism/_kernels.py | 690 +++++++++++++++ choclo/prism/_magnetic.py | 1234 ++++++++++++++++++++++++--- choclo/tests/test_prism_magnetic.py | 214 ++++- doc/api/index.rst | 16 + 5 files changed, 2047 insertions(+), 130 deletions(-) diff --git a/choclo/prism/__init__.py b/choclo/prism/__init__.py index 04f37b8a..4288e4d4 100644 --- a/choclo/prism/__init__.py +++ b/choclo/prism/__init__.py @@ -22,13 +22,34 @@ from ._kernels import ( kernel_e, kernel_ee, + kernel_eee, + kernel_een, + kernel_eeu, kernel_en, + kernel_enn, + kernel_enu, kernel_eu, + kernel_euu, kernel_n, kernel_nn, + kernel_nnn, + kernel_nnu, kernel_nu, + kernel_nuu, kernel_pot, kernel_u, kernel_uu, + kernel_uuu, +) +from ._magnetic import ( + magnetic_e, + magnetic_ee, + magnetic_en, + magnetic_eu, + magnetic_field, + magnetic_n, + magnetic_nn, + magnetic_nu, + magnetic_u, + magnetic_uu, ) -from ._magnetic import magnetic_e, magnetic_field, magnetic_n, magnetic_u diff --git a/choclo/prism/_kernels.py b/choclo/prism/_kernels.py index 4b1eea65..b3f8f845 100644 --- a/choclo/prism/_kernels.py +++ b/choclo/prism/_kernels.py @@ -899,3 +899,693 @@ def _safe_log(x, y, z, r): result = np.log((y**2 + z**2) / (r - x)) return result return np.log(x + r) + + +@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 : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate 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. + + 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) + + +@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 : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate 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. + + 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) + + +@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 : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate 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. + + 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) + + +@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 : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate 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. + + 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) + + +@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 : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate 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. + + 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) + + +@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 : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate 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. + + 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) + + +@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 : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate 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. + + 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) + + +@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 : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate 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. + + 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) + + +@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 : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate 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. + + 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) + + +@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 : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate 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. + + 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 + """ + 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 : floats + 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. + + 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 (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 : floats + 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. + + 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 x_i == 0 and x_j == 0: + return 0.0 + return -x_i / ((x_k + radius) * radius) diff --git a/choclo/prism/_magnetic.py b/choclo/prism/_magnetic.py index 3d358d7a..d8ff39c9 100644 --- a/choclo/prism/_magnetic.py +++ b/choclo/prism/_magnetic.py @@ -11,7 +11,24 @@ 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 +from ._kernels import ( + kernel_ee, + kernel_eee, + kernel_een, + kernel_eeu, + kernel_en, + kernel_enn, + kernel_enu, + kernel_eu, + kernel_euu, + kernel_nn, + kernel_nnn, + kernel_nnu, + kernel_nu, + kernel_nuu, + kernel_uu, + kernel_uuu, +) from ._utils import ( is_interior_point, is_point_on_east_face, @@ -198,8 +215,7 @@ def magnetic_field( .. math:: u_{ij} = - \frac{\partial}{\partial i} - \frac{\partial}{\partial j} + \frac{\partial^2}{\partial i \partial j} \int\limits_R \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} dv @@ -226,6 +242,12 @@ def magnetic_field( :func:`choclo.prism.magnetic_e` :func:`choclo.prism.magnetic_n` :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` """ # Check if observation point falls in a singular point if is_point_on_edge( @@ -428,8 +450,7 @@ def magnetic_e( .. math:: u_{ij} = - \frac{\partial}{\partial i} - \frac{\partial}{\partial j} + \frac{\partial^2}{\partial i \partial j} \int\limits_R \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} dv @@ -474,6 +495,12 @@ def magnetic_e( :func:`choclo.prism.magnetic_field` :func:`choclo.prism.magnetic_n` :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` """ # Check if observation point falls in a singular point if is_point_on_edge( @@ -498,43 +525,24 @@ def magnetic_e( prism_top, ): return np.nan - # 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 - 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) - # 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_east * k_ee - + magnetization_north * k_en - + magnetization_up * k_eu - ) + # Compute magnetic field vector component + b_e = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_ee, + kernel_en, + kernel_eu, + ) # Add 4 pi to Be if computing on the eastmost face, to correctly evaluate # the limit approaching from outside (approaching from the east) if is_point_on_east_face( @@ -629,8 +637,7 @@ def magnetic_n( .. math:: u_{ij} = - \frac{\partial}{\partial i} - \frac{\partial}{\partial j} + \frac{\partial^2}{\partial i \partial j} \int\limits_R \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} dv @@ -675,6 +682,12 @@ def magnetic_n( :func:`choclo.prism.magnetic_field` :func:`choclo.prism.magnetic_e` :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` """ # Check if observation point falls in a singular point if is_point_on_edge( @@ -699,43 +712,24 @@ def magnetic_n( prism_top, ): return np.nan - # 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 - 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) - # 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_east * k_en - + magnetization_north * k_nn - + magnetization_up * k_nu - ) + # Compute magnetic field vector component + b_n = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_en, + kernel_nn, + kernel_nu, + ) # Add 4 pi to Bn if computing on the northmost face, to correctly evaluate # the limit approaching from outside (approaching from the north) if is_point_on_north_face( @@ -830,8 +824,7 @@ def magnetic_u( .. math:: u_{ij} = - \frac{\partial}{\partial i} - \frac{\partial}{\partial j} + \frac{\partial^2}{\partial i \partial j} \int\limits_R \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} dv @@ -876,6 +869,12 @@ def magnetic_u( :func:`choclo.prism.magnetic_field` :func:`choclo.prism.magnetic_e` :func:`choclo.prism.magnetic_n` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` """ # Check if observation point falls in a singular point if is_point_on_edge( @@ -900,43 +899,24 @@ def magnetic_u( prism_top, ): return np.nan - # 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 - 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) - # 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_east * k_eu - + magnetization_north * k_nu - + magnetization_up * k_uu - ) + # Compute magnetic field vector component + b_u = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_eu, + kernel_nu, + kernel_uu, + ) # Add 4 pi to Bu if computing on the north face, to correctly evaluate the # limit approaching from outside (approaching from the top) if is_point_on_top_face( @@ -952,3 +932,1017 @@ def magnetic_u( ): b_u += 4 * np.pi * magnetization_up return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_u + + +@jit(nopython=True) +def magnetic_ee( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, +): + r""" + Easting derivative of the easting component of the magnetic field. + + Returns the easting derivative of 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_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. + magnetization_east : float + The East component of the magnetization vector of the prism. Must be in + :math:`A m^{-1}`. + magnetization_north : float + The North component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + magnetization_up : float + The upward component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + + Returns + ------- + b_ee : float + Easting derivative of the easting component of the magnetic field + generated by the prism on the observation point in :math:`\text{T}`. + Return ``numpy.nan`` if the observation point falls in + a singular point: prism vertices, prism edges or interior points. + + Notes + ----- + Computes the easting derivative of 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_{xx}(\mathbf{p}) = + \frac{\mu_0}{4\pi} + \left( M_x u_{xxx} + M_y u_{xxy} + M_z u_{xxz} \right) + + Where :math:`u_{ijk}` are: + + .. math:: + + u_{ijk} = + \frac{\partial^3}{\partial i \partial j \partial k} + \int\limits_R + \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} + dv + + with :math:`i,j,k \in \{x, y, z\}`. + + References + ---------- + - [Blakely1995]_ + - [Oliveira2015]_ + - [Nagy2000]_ + - [Nagy2002]_ + - [Fukushima2020]_ + + See Also + -------- + :func:`choclo.prism.magnetic_field` + :func:`choclo.prism.magnetic_e` + :func:`choclo.prism.magnetic_n` + :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` + """ + # Check if observation point falls in a singular point + if is_point_on_edge( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ) or is_interior_point( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ): + return np.nan + # Compute magnetic gradiometry component + b_ee = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_eee, + kernel_een, + kernel_eeu, + ) + return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_ee + + +@jit(nopython=True) +def magnetic_nn( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, +): + r""" + Northing derivative of the northing component of the magnetic field. + + Returns the northing derivative of 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_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. + magnetization_east : float + The East component of the magnetization vector of the prism. Must be in + :math:`A m^{-1}`. + magnetization_north : float + The North component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + magnetization_up : float + The upward component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + + Returns + ------- + b_nn : float + Northing derivative of the northing component of the magnetic field + generated by the prism on the observation point in :math:`\text{T}`. + Return ``numpy.nan`` if the observation point falls in + a singular point: prism vertices, prism edges or interior points. + + Notes + ----- + Computes the northing derivative of 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_{yy}(\mathbf{p}) = + \frac{\mu_0}{4\pi} + \left( M_x u_{xyy} + M_y u_{yyy} + M_z u_{yyz} \right) + + Where :math:`u_{ijk}` are: + + .. math:: + + u_{ijk} = + \frac{\partial^3}{\partial i \partial j \partial k} + \int\limits_R + \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} + dv + + with :math:`i,j,k \in \{x, y, z\}`. + + References + ---------- + - [Blakely1995]_ + - [Oliveira2015]_ + - [Nagy2000]_ + - [Nagy2002]_ + - [Fukushima2020]_ + + See Also + -------- + :func:`choclo.prism.magnetic_field` + :func:`choclo.prism.magnetic_e` + :func:`choclo.prism.magnetic_n` + :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` + """ + # Check if observation point falls in a singular point + if is_point_on_edge( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ) or is_interior_point( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ): + return np.nan + # Compute magnetic gradiometry component + b_nn = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_enn, + kernel_nnn, + kernel_nnu, + ) + return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_nn + + +@jit(nopython=True) +def magnetic_uu( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, +): + r""" + Upward derivative of the upward component of the magnetic field. + + Returns the upward derivative of 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_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. + magnetization_east : float + The East component of the magnetization vector of the prism. Must be in + :math:`A m^{-1}`. + magnetization_north : float + The North component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + magnetization_up : float + The upward component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + + Returns + ------- + b_uu : float + Upward derivative of the upward component of the magnetic field + generated by the prism on the observation point in :math:`\text{T}`. + Return ``numpy.nan`` if the observation point falls in + a singular point: prism vertices, prism edges or interior points. + + Notes + ----- + Computes the upward derivative of 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_{zz}(\mathbf{p}) = + \frac{\mu_0}{4\pi} + \left( M_x u_{xzz} + M_y u_{yzz} + M_z u_{zzz} \right) + + Where :math:`u_{ijk}` are: + + .. math:: + + u_{ijk} = + \frac{\partial^3}{\partial i \partial j \partial k} + \int\limits_R + \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} + dv + + with :math:`i,j,k \in \{x, y, z\}`. + + References + ---------- + - [Blakely1995]_ + - [Oliveira2015]_ + - [Nagy2000]_ + - [Nagy2002]_ + - [Fukushima2020]_ + + See Also + -------- + :func:`choclo.prism.magnetic_field` + :func:`choclo.prism.magnetic_e` + :func:`choclo.prism.magnetic_n` + :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` + """ + # Check if observation point falls in a singular point + if is_point_on_edge( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ) or is_interior_point( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ): + return np.nan + # Compute magnetic gradiometry component + b_uu = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_euu, + kernel_nuu, + kernel_uuu, + ) + return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_uu + + +@jit(nopython=True) +def magnetic_en( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, +): + r""" + Northing derivative of the easting component of the magnetic field. + + Returns the northing derivative of 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_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. + magnetization_east : float + The East component of the magnetization vector of the prism. Must be in + :math:`A m^{-1}`. + magnetization_north : float + The North component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + magnetization_up : float + The upward component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + + Returns + ------- + b_en : float + Northing derivative of the easting component of the magnetic field + generated by the prism on the observation point in :math:`\text{T}`. + Return ``numpy.nan`` if the observation point falls in a singular + point: prism vertices, prism edges or interior points. + + Notes + ----- + Computes the northing derivative of 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_{xy}(\mathbf{p}) = + \frac{\mu_0}{4\pi} + \left( M_x u_{xxy} + M_y u_{xyy} + M_z u_{xyz} \right) + + Where :math:`u_{ijk}` are: + + .. math:: + + u_{ijk} = + \frac{\partial^3}{\partial i \partial j \partial k} + \int\limits_R + \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} + dv + + with :math:`i,j,k \in \{x, y, z\}`. + + References + ---------- + - [Blakely1995]_ + - [Oliveira2015]_ + - [Nagy2000]_ + - [Nagy2002]_ + - [Fukushima2020]_ + + See Also + -------- + :func:`choclo.prism.magnetic_field` + :func:`choclo.prism.magnetic_e` + :func:`choclo.prism.magnetic_n` + :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` + """ + # Check if observation point falls in a singular point + if is_point_on_edge( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ) or is_interior_point( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ): + return np.nan + # Compute magnetic gradiometry component + b_en = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_een, + kernel_enn, + kernel_enu, + ) + return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_en + + +@jit(nopython=True) +def magnetic_eu( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, +): + r""" + Upward derivative of the easting component of the magnetic field. + + Returns the upward derivative of 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_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. + magnetization_east : float + The East component of the magnetization vector of the prism. Must be in + :math:`A m^{-1}`. + magnetization_north : float + The North component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + magnetization_up : float + The upward component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + + Returns + ------- + b_eu : float + Upward derivative of the easting component of the magnetic field + generated by the prism on the observation point in :math:`\text{T}`. + Return ``numpy.nan`` if the observation point falls in a singular + point: prism vertices, prism edges or interior points. + + Notes + ----- + Computes the northing derivative of 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_{xz}(\mathbf{p}) = + \frac{\mu_0}{4\pi} + \left( M_x u_{xxz} + M_y u_{xyz} + M_z u_{xzz} \right) + + Where :math:`u_{ijk}` are: + + .. math:: + + u_{ijk} = + \frac{\partial^3}{\partial i \partial j \partial k} + \int\limits_R + \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} + dv + + with :math:`i,j,k \in \{x, y, z\}`. + + References + ---------- + - [Blakely1995]_ + - [Oliveira2015]_ + - [Nagy2000]_ + - [Nagy2002]_ + - [Fukushima2020]_ + + See Also + -------- + :func:`choclo.prism.magnetic_field` + :func:`choclo.prism.magnetic_e` + :func:`choclo.prism.magnetic_n` + :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_nu` + """ + # Check if observation point falls in a singular point + if is_point_on_edge( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ) or is_interior_point( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ): + return np.nan + # Compute magnetic gradiometry component + b_eu = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_eeu, + kernel_enu, + kernel_euu, + ) + return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_eu + + +@jit(nopython=True) +def magnetic_nu( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, +): + r""" + Upward derivative of the northing component of the magnetic field. + + Returns the upward derivative of 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_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. + magnetization_east : float + The East component of the magnetization vector of the prism. Must be in + :math:`A m^{-1}`. + magnetization_north : float + The North component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + magnetization_up : float + The upward component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + + Returns + ------- + b_nu : float + Upward derivative of the northing component of the magnetic field + generated by the prism on the observation point in :math:`\text{T}`. + Return ``numpy.nan`` if the observation point falls in a singular + point: prism vertices, prism edges or interior points. + + Notes + ----- + Computes the northing derivative of 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_{yz}(\mathbf{p}) = + \frac{\mu_0}{4\pi} + \left( M_x u_{xyz} + M_y u_{yyz} + M_z u_{yzz} \right) + + Where :math:`u_{ijk}` are: + + .. math:: + + u_{ijk} = + \frac{\partial^3}{\partial i \partial j \partial k} + \int\limits_R + \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} + dv + + with :math:`i,j,k \in \{x, y, z\}`. + + References + ---------- + - [Blakely1995]_ + - [Oliveira2015]_ + - [Nagy2000]_ + - [Nagy2002]_ + - [Fukushima2020]_ + + See Also + -------- + :func:`choclo.prism.magnetic_field` + :func:`choclo.prism.magnetic_e` + :func:`choclo.prism.magnetic_n` + :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + """ + # Check if observation point falls in a singular point + if is_point_on_edge( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ) or is_interior_point( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ): + return np.nan + # Compute magnetic gradiometry component + b_nu = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_enu, + kernel_nnu, + kernel_nuu, + ) + return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_nu + + +@jit(nopython=True) +def _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_i, + kernel_j, + kernel_k, +): + r""" + Calculate field component for a single prism and observation point. + + Evaluate the provided kernels on the shifted coordinates of prism vertices + to compute the magnetic field component given a magnetization vector of the + prism. + + Parameters + ---------- + easting, northing, upward : floats + Coordinates of the observation point. Must be in meters. + prism_west, prism_east : floats + The easting boundaries of the prism. Must be in meters. + prism_south, prism_north : floats + The northing boundaries of the prism. Must be in meters. + prism_bottom, prism_top : floats + The upward boundaries of the prism. Must be in meters. + magnetization_east, magnetization_north, magnetization_up : floats + The components of the magnetization vector of the prism. Must be in + :math:`A m^{-1}`. + kernel_i, kernel_j, kernel_k : callables + Kernel functions to be evaluated on each vertex of the prism. + + Returns + ------- + float + + Notes + ----- + Given the kernels :math:`k_i(\hat{x}, \hat{y}, \hat{z})`, + :math:`k_j(\hat{x}, \hat{y}, \hat{z})`, and :math:`k_k(\hat{x}, \hat{y}, + \hat{z})`; a prism provided by its boundaries :math:`x_1`, :math:`x_2`, + :math:`y_1`, :math:`y_2`, :math:`z_1`, and :math:`z_2`; a magnetization + vector :math:`\mathbf{M}=(M_x, M_y, M_z)`; and an observation point + :math:`\mathbf{p}=(x, y, z)`, this function returns: + + .. math:: + + M_x u_x(x, y, z) + M_y u_y(x, y, z) + M_z u_z(x, y, z), + + where + + .. math:: + + u_x(x, y, z) = + \Bigg\lvert\Bigg\lvert\Bigg\lvert + k_i(\hat{x}, \hat{y}, \hat{z}) + \Bigg\rvert_{x_1 - x}^{x_2 - x} + \Bigg\rvert_{y_1 - y}^{y_2 - y} + \Bigg\rvert_{z_1 - z}^{z_2 - z} + + .. math:: + + u_y(x, y, z) = + \Bigg\lvert\Bigg\lvert\Bigg\lvert + k_j(\hat{x}, \hat{y}, \hat{z}) + \Bigg\rvert_{x_1 - x}^{x_2 - x} + \Bigg\rvert_{y_1 - y}^{y_2 - y} + \Bigg\rvert_{z_1 - z}^{z_2 - z} + + .. math:: + + u_z(x, y, z) = + \Bigg\lvert\Bigg\lvert\Bigg\lvert + k_k(\hat{x}, \hat{y}, \hat{z}) + \Bigg\rvert_{x_1 - x}^{x_2 - x} + \Bigg\rvert_{y_1 - y}^{y_2 - y} + \Bigg\rvert_{z_1 - z}^{z_2 - z} + """ + result = 0.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) + # Compute kernel tensor components for the current vertex + k_e = kernel_i(shift_east, shift_north, shift_upward, radius) + k_n = kernel_j(shift_east, shift_north, shift_upward, radius) + k_u = kernel_k(shift_east, shift_north, shift_upward, radius) + # Compute b_en using the dot product between the kernel tensor + # and the magnetization vector of the prism + result += (-1) ** (i + j + k) * ( + magnetization_east * k_e + + magnetization_north * k_n + + magnetization_up * k_u + ) + return result diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index 224cb9c1..0b46c3d6 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -7,6 +7,8 @@ """ Test magnetic forward modelling functions for rectangular prisms """ +import itertools + import numpy as np import numpy.testing as npt import pytest @@ -16,11 +18,35 @@ kernel_eu, kernel_nu, magnetic_e, + magnetic_ee, + magnetic_en, + magnetic_eu, magnetic_field, magnetic_n, + magnetic_nn, + magnetic_nu, magnetic_u, + magnetic_uu, ) +COMPONENTS = { + "e": magnetic_e, + "n": magnetic_n, + "u": magnetic_u, + "ee": magnetic_ee, + "nn": magnetic_nn, + "uu": magnetic_uu, + "en": magnetic_en, + "eu": magnetic_eu, + "nu": magnetic_nu, +} + + +def get_forward_function(component: str): + """Return desired magnetic forward function.""" + component = "".join(sorted(component)) + return COMPONENTS[component] + @pytest.fixture(name="sample_prism") def fixture_sample_prism(): @@ -599,6 +625,21 @@ class TestDivergenceOfB: Test if the divergence of the magnetic field is equal to zero. """ + def test_divergence_of_b(self, sample_3d_grid, sample_prism, sample_magnetization): + """ + Test div of B through magnetic gradiometry analytical functions. + """ + kwargs = dict( + coordinates=sample_3d_grid, + prism=sample_prism, + magnetization=sample_magnetization, + ) + b_ee = evaluate(magnetic_ee, **kwargs) + b_nn = evaluate(magnetic_nn, **kwargs) + b_uu = evaluate(magnetic_uu, **kwargs) + # Check if the divergence of B is zero + npt.assert_allclose(-b_uu, b_ee + b_nn, atol=1e-11) + def test_divergence_of_b_finite_differences( self, sample_3d_grid, sample_prism, sample_magnetization ): @@ -635,6 +676,19 @@ class TestMagneticFieldSingularities: we approach from outside of the prism. """ + COMPONENTS = ( + magnetic_field, + magnetic_e, + magnetic_n, + magnetic_u, + magnetic_ee, + magnetic_en, + magnetic_eu, + magnetic_nn, + magnetic_nu, + magnetic_uu, + ) + @pytest.fixture() def sample_prism(self): """ @@ -711,9 +765,7 @@ def get_interior_points(self, prism): coordinates = tuple(c.ravel() for c in np.meshgrid(easting, northing, upward)) return coordinates - @pytest.mark.parametrize( - "forward_func", (magnetic_field, magnetic_e, magnetic_n, magnetic_u) - ) + @pytest.mark.parametrize("forward_func", COMPONENTS) def test_on_vertices(self, sample_prism, forward_func): """ Test if magnetic field components on vertices are equal to NaN @@ -726,9 +778,7 @@ def test_on_vertices(self, sample_prism, forward_func): ) assert np.isnan(results).all() - @pytest.mark.parametrize( - "forward_func", (magnetic_field, magnetic_e, magnetic_n, magnetic_u) - ) + @pytest.mark.parametrize("forward_func", COMPONENTS) @pytest.mark.parametrize("direction", ("easting", "northing", "upward")) def test_on_edges(self, sample_prism, direction, forward_func): """ @@ -744,9 +794,7 @@ def test_on_edges(self, sample_prism, direction, forward_func): ) assert np.isnan(results).all() - @pytest.mark.parametrize( - "forward_func", (magnetic_field, magnetic_e, magnetic_n, magnetic_u) - ) + @pytest.mark.parametrize("forward_func", COMPONENTS) def test_on_interior_points(self, sample_prism, forward_func): """ Test if magnetic field components are NaN on internal points @@ -787,6 +835,154 @@ def test_symmetry_on_faces(self, sample_prism, direction, forward_func): ) npt.assert_allclose(results, results[0]) + @pytest.mark.parametrize( + "forward_func", + ( + magnetic_ee, + magnetic_nn, + magnetic_uu, + magnetic_en, + magnetic_eu, + magnetic_nu, + ), + ) + @pytest.mark.parametrize("direction", ("easting", "northing", "upward")) + def test_gradiometry_symmetry_on_faces(self, sample_prism, direction, forward_func): + """ + Tests symmetry of magnetic gradiometry components on the center of + faces normal to the component direction + + For example, check if ``magnetic_ee`` has the opposite value on points + in the top face and points of the bottom face. + """ + easting, northing, upward = self.get_faces_centers(sample_prism, direction) + magnetization = np.array([1.0, 1.0, 1.0]) + results = list( + forward_func(e, n, u, *sample_prism, *magnetization) + for (e, n, u) in zip(easting, northing, upward) + ) + npt.assert_allclose(-results[0], results[1:], atol=1e-23) + + +class TestMagGradiometryFiniteDifferences: + """ + Test the magnetic gradiometry components by comparing them to numerical + derivatives computed through finite differences of the magnetic components. + """ + + delta = 1e-6 # displacement used in the finite difference calculations + rtol, atol = 5e-4, 5e-12 # tolerances used in the comparisons + + @pytest.mark.parametrize("i", ["e", "n", "u"]) + @pytest.mark.parametrize("j", ["e", "n", "u"]) + def test_derivatives_of_magnetic_components( + self, sample_3d_grid, sample_prism, sample_magnetization, i, j + ): + """ + Compare the magnetic gradiometry functions with a finite difference + approximation of them computed from the magnetic field components. + + With this function we can compare the result of ``magnetic_{i}{j}`` + against a finite difference approximation using the ``magnetic_{i}`` + forward function, where ``i`` and ``j`` can be ``e``, ``n`` or ``u``. + """ + # Get forward functions + mag_func = get_forward_function(i) + mag_grad_func = get_forward_function(i + j) + direction = j # direction along which to apply the finite differences + # Evaluate the mag grad function on the sample grid + mag_grad = evaluate( + mag_grad_func, sample_3d_grid, sample_prism, sample_magnetization + ) + # Compute the mag grad function using finite differences + mag_grad_finite_diff = finite_differences( + sample_3d_grid, + sample_prism, + sample_magnetization, + direction=direction, + forward_func=mag_func, + delta=self.delta, + ) + # Compare results + npt.assert_allclose( + mag_grad, mag_grad_finite_diff, rtol=self.rtol, atol=self.atol + ) + + +class TestMagGradiometryLimits: + """ + Test if limits are well evaluated in magnetic gradiometry kernels. + + Most of the third-order kernels have singular points along the lines that + matches with the edges of the prism, but the magnetic gradiometry functions + are defined (and finite) in those regions. These tests ensure that these + limits are well resolved. + """ + + delta = 1e-6 # displacement used in the finite difference calculations + rtol, atol = 5e-4, 5e-12 # tolerances used in the comparisons + + @pytest.fixture + def coordinates(self, sample_prism): + """ + Return a set of coordinates located along the lines of the prism edges. + """ + west, east, south, north, bottom, top = sample_prism + distance = 13.4 # define distance from the nearest vertex for obs points + points = [] + # Define points on lines parallel to upward edges + for easting, northing in itertools.product((west, east), (south, north)): + points.append([easting, northing, bottom - distance]) + points.append([easting, northing, top + distance]) + # Define points on lines parallel to easting edges + for northing, upward in itertools.product((south, north), (bottom, top)): + points.append([west - distance, northing, upward]) + points.append([east + distance, northing, upward]) + # Define points on lines parallel to northing edges + for easting, upward in itertools.product((west, east), (bottom, top)): + points.append([easting, south - distance, upward]) + points.append([easting, north + distance, upward]) + points = np.array(points).T + easting, northing, upward = points[:] + return easting, northing, upward + + @pytest.mark.parametrize("i", ["e", "n", "u"]) + @pytest.mark.parametrize("j", ["e", "n", "u"]) + def test_derivatives_of_magnetic_components_on_edges( + self, coordinates, sample_prism, sample_magnetization, i, j + ): + """ + Compare the magnetic gradiometry functions with a finite difference + approximation of it computed from the magnetic field components on + observation points that are located along the same lines of the prism + edges. + + With this function we can compare the result of ``magnetic_{i}{j}`` + against a finite difference approximation using the ``magnetic_{i}`` + forward function, where ``i`` and ``j`` can be ``e``, ``n`` or ``u``. + """ + # Get forward functions + mag_func = get_forward_function(i) + mag_grad_func = get_forward_function(f"{i}{j}") + direction = j # direction along which to apply the finite differences + # Evaluate the mag grad function on the sample grid + mag_grad = evaluate( + mag_grad_func, coordinates, sample_prism, sample_magnetization + ) + # Compute the mag grad function using finite differences + mag_grad_finite_diff = finite_differences( + coordinates, + sample_prism, + sample_magnetization, + direction=direction, + forward_func=mag_func, + delta=self.delta, + ) + # Compare results + npt.assert_allclose( + mag_grad, mag_grad_finite_diff, rtol=self.rtol, atol=self.atol + ) + class TestBugfixKernelEvaluation: r""" diff --git a/doc/api/index.rst b/doc/api/index.rst index c7ae9bbf..6107e2bd 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -95,6 +95,12 @@ Magnetic prism.magnetic_e prism.magnetic_n prism.magnetic_u + prism.magnetic_ee + prism.magnetic_nn + prism.magnetic_uu + prism.magnetic_en + prism.magnetic_eu + prism.magnetic_nu Kernels ^^^^^^^ @@ -112,6 +118,16 @@ Kernels prism.kernel_en prism.kernel_eu prism.kernel_nu + prism.kernel_eee + prism.kernel_nnn + prism.kernel_uuu + prism.kernel_een + prism.kernel_eeu + prism.kernel_enn + prism.kernel_nnu + prism.kernel_euu + prism.kernel_nuu + prism.kernel_enu Utilities