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