From c845e9023769039a3dc1cdc460fe27a242ad05ed Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 20 Sep 2024 11:34:38 -0700 Subject: [PATCH 01/34] Add diagonal third order kernels --- choclo/prism/_kernels.py | 208 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 208 insertions(+) diff --git a/choclo/prism/_kernels.py b/choclo/prism/_kernels.py index 9646c72..e8274c6 100644 --- a/choclo/prism/_kernels.py +++ b/choclo/prism/_kernels.py @@ -845,3 +845,211 @@ def _safe_log(x, r): result = np.log((r**2 - x**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} + } + + """ + 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} + } + + """ + 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} + } + + """ + return _kernel_iii(upward, northing, easting, 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 + >>> 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) + """ + 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 From cc9f30dbe74559e6824df46d5f8c7da76b1f6dac Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 20 Sep 2024 14:10:56 -0700 Subject: [PATCH 02/34] Add non-diagonal kernel functions --- choclo/prism/_kernels.py | 374 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 373 insertions(+), 1 deletion(-) diff --git a/choclo/prism/_kernels.py b/choclo/prism/_kernels.py index e8274c6..31a1cc6 100644 --- a/choclo/prism/_kernels.py +++ b/choclo/prism/_kernels.py @@ -882,7 +882,8 @@ def kernel_eee(easting, northing, upward, radius): Notes ----- - Computes the following numerical kernel on the passed *shifted coordinates*: + Computes the following numerical kernel on the passed *shifted + coordinates*: .. math:: @@ -1002,6 +1003,318 @@ def kernel_uuu(easting, northing, upward, radius): 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} + } + + """ + 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} + } + + """ + 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} + } + + """ + 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} + } + + """ + 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} + } + + """ + 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} + } + + """ + return _kernel_iij(upward, northing, easting, radius) + + @jit(nopython=True) def _kernel_iii(x_i, x_j, x_k, radius): r""" @@ -1044,6 +1357,7 @@ def _kernel_iii(x_i, x_j, x_k, radius): 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) @@ -1053,3 +1367,61 @@ def _kernel_iii(x_i, x_j, x_k, radius): 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) + """ + return x_i / ((x_k + radius) * radius) From af40c19f659d3416fc5f203138b3c4e3ff15b911 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 20 Sep 2024 14:13:19 -0700 Subject: [PATCH 03/34] Add function for k_xyz --- choclo/prism/_kernels.py | 46 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/choclo/prism/_kernels.py b/choclo/prism/_kernels.py index 31a1cc6..96deb62 100644 --- a/choclo/prism/_kernels.py +++ b/choclo/prism/_kernels.py @@ -1315,6 +1315,52 @@ def kernel_nuu(easting, northing, upward, radius): 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}} + + """ + return 1 / radius + + @jit(nopython=True) def _kernel_iii(x_i, x_j, x_k, radius): r""" From 6c9c2162c37c1ca20d47dec9ae91dae9b30a5ae1 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 23 Sep 2024 16:02:35 -0700 Subject: [PATCH 04/34] Start adding forward functions for grad components --- choclo/prism/_magnetic.py | 160 +++++++++++++++++++++++++++++++++++++- 1 file changed, 159 insertions(+), 1 deletion(-) diff --git a/choclo/prism/_magnetic.py b/choclo/prism/_magnetic.py index 3d358d7..7819eb5 100644 --- a/choclo/prism/_magnetic.py +++ b/choclo/prism/_magnetic.py @@ -11,7 +11,23 @@ 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_en, + kernel_eu, + kernel_nn, + kernel_nu, + kernel_uu, + kernel_eee, + kernel_nnn, + kernel_uuu, + kernel_een, + kernel_eeu, + kernel_enn, + kernel_nnu, + kernel_euu, + kernel_nuu, +) from ._utils import ( is_interior_point, is_point_on_east_face, @@ -952,3 +968,145 @@ 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 due to a prism + + 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 (TODO: specify which are singular 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}{\partial i} + \frac{\partial}{\partial j} + \frac{\partial}{\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` + """ + # TODO: Check if observation point falls in a singular point + # Initialize magnetic gradiometry component + b_ee = 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_eee = kernel_eee(shift_east, shift_north, shift_upward, radius) + k_een = kernel_een(shift_east, shift_north, shift_upward, radius) + k_eeu = kernel_eeu(shift_east, shift_north, shift_upward, radius) + # Compute b_ee using the dot product between the kernel tensor + # and the magnetization vector of the prism + b_ee += (-1) ** (i + j + k) * ( + magnetization_east * k_eee + + magnetization_north * k_een + + magnetization_up * k_eeu + ) + return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_ee From f90dd306708b81a62525b0807ba6475f48714306 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 24 Sep 2024 10:57:15 -0700 Subject: [PATCH 05/34] Start adding the other diagonal components --- choclo/prism/_magnetic.py | 292 +++++++++++++++++++++++++++++++++++++- 1 file changed, 291 insertions(+), 1 deletion(-) diff --git a/choclo/prism/_magnetic.py b/choclo/prism/_magnetic.py index 7819eb5..86b1162 100644 --- a/choclo/prism/_magnetic.py +++ b/choclo/prism/_magnetic.py @@ -986,7 +986,7 @@ def magnetic_ee( magnetization_up, ): r""" - Easting derivative of the easting component of the magnetic field due to a prism + 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. @@ -1070,6 +1070,8 @@ def magnetic_ee( :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` """ # TODO: Check if observation point falls in a singular point # Initialize magnetic gradiometry component @@ -1110,3 +1112,291 @@ def magnetic_ee( + magnetization_up * k_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 (TODO: specify which are singular 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}{\partial i} + \frac{\partial}{\partial j} + \frac{\partial}{\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` + """ + # TODO: Check if observation point falls in a singular point + # Initialize magnetic gradiometry component + b_nn = 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_enn = kernel_enn(shift_east, shift_north, shift_upward, radius) + k_nnn = kernel_nnn(shift_east, shift_north, shift_upward, radius) + k_nnu = kernel_nnu(shift_east, shift_north, shift_upward, radius) + # Compute b_nn using the dot product between the kernel tensor + # and the magnetization vector of the prism + b_nn += (-1) ** (i + j + k) * ( + magnetization_east * k_enn + + magnetization_north * k_nnn + + magnetization_up * k_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 (TODO: specify which are singular 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}{\partial i} + \frac{\partial}{\partial j} + \frac{\partial}{\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` + """ + # TODO: Check if observation point falls in a singular point + # Initialize magnetic gradiometry component + b_uu = 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_euu = kernel_euu(shift_east, shift_north, shift_upward, radius) + k_nuu = kernel_nuu(shift_east, shift_north, shift_upward, radius) + k_uuu = kernel_uuu(shift_east, shift_north, shift_upward, radius) + # Compute b_uu using the dot product between the kernel tensor + # and the magnetization vector of the prism + b_uu += (-1) ** (i + j + k) * ( + magnetization_east * k_euu + + magnetization_north * k_nuu + + magnetization_up * k_uuu + ) + return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_uu From 1fcfdebe20ebe656ef71160a532d55c6ca0263d6 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 24 Sep 2024 11:00:56 -0700 Subject: [PATCH 06/34] Add new functions to __init__.py --- choclo/prism/__init__.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/choclo/prism/__init__.py b/choclo/prism/__init__.py index 04f37b8..289bbea 100644 --- a/choclo/prism/__init__.py +++ b/choclo/prism/__init__.py @@ -31,4 +31,12 @@ kernel_u, kernel_uu, ) -from ._magnetic import magnetic_e, magnetic_field, magnetic_n, magnetic_u +from ._magnetic import ( + magnetic_e, + magnetic_field, + magnetic_n, + magnetic_u, + magnetic_ee, + magnetic_nn, + magnetic_uu, +) From 0bbd8c37eead410b18a071ac0de7f88ca619c8c4 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 24 Sep 2024 11:01:07 -0700 Subject: [PATCH 07/34] Add divergence of B test with analytic solutions --- choclo/tests/test_prism_magnetic.py | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index fcecddd..ff6b1bd 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -11,7 +11,15 @@ import numpy.testing as npt import pytest -from ..prism import magnetic_e, magnetic_field, magnetic_n, magnetic_u +from ..prism import ( + magnetic_e, + magnetic_field, + magnetic_n, + magnetic_u, + magnetic_ee, + magnetic_nn, + magnetic_uu, +) @pytest.fixture(name="sample_prism") @@ -591,6 +599,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 ): From 06b51769a6e970787e9407c02b415987a2354570 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 24 Sep 2024 16:21:58 -0700 Subject: [PATCH 08/34] Add kernels to __init__.py --- choclo/prism/__init__.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/choclo/prism/__init__.py b/choclo/prism/__init__.py index 289bbea..a740f03 100644 --- a/choclo/prism/__init__.py +++ b/choclo/prism/__init__.py @@ -30,6 +30,15 @@ kernel_pot, kernel_u, kernel_uu, + kernel_eee, + kernel_nnn, + kernel_uuu, + kernel_een, + kernel_eeu, + kernel_enn, + kernel_nnu, + kernel_euu, + kernel_nuu, ) from ._magnetic import ( magnetic_e, From e03721cf57fdfe522086cc2d389e6cfdcc859cae Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 24 Sep 2024 16:22:07 -0700 Subject: [PATCH 09/34] Start adding some tests --- choclo/tests/test_prism_magnetic.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index ff6b1bd..27235bc 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -801,3 +801,24 @@ def test_symmetry_on_faces(self, sample_prism, direction, forward_func): for (e, n, u) in zip(easting, northing, upward) ) npt.assert_allclose(results, results[0]) + + +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 + + def test(self, sample_3d_grid, sample_prism, sample_magnetization): + b_ee = evaluate(magnetic_ee, sample_3d_grid, sample_prism, sample_magnetization) + b_ee_finite_diff = finite_differences( + sample_3d_grid, + sample_prism, + sample_magnetization, + direction="e", + forward_func=magnetic_e, + delta=self.delta, + ) + npt.assert_allclose(b_ee, b_ee_finite_diff, atol=1e-14) From 4dd6d3d5ad0a41f4a398bc82177f338cb7b73f65 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 26 Sep 2024 16:34:48 -0700 Subject: [PATCH 10/34] Add missing kernel_enu to `__init__.py` --- choclo/prism/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/choclo/prism/__init__.py b/choclo/prism/__init__.py index a740f03..8f24af6 100644 --- a/choclo/prism/__init__.py +++ b/choclo/prism/__init__.py @@ -39,6 +39,7 @@ kernel_nnu, kernel_euu, kernel_nuu, + kernel_enu, ) from ._magnetic import ( magnetic_e, From c36099f8f687f446cb83898183ec94e56f75778b Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 26 Sep 2024 16:40:36 -0700 Subject: [PATCH 11/34] Fix sign of third-order kernels Update the signs in the documentation and in the implementation of each kernel. Also make the kernels to return zero on singular points (not on vertices) so when computing on lines shared with the edges are not falling in division by zero errors and the limits for those cases are well resolved. --- choclo/prism/_kernels.py | 34 +++++++++++++++++++--------------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/choclo/prism/_kernels.py b/choclo/prism/_kernels.py index 96deb62..148d254 100644 --- a/choclo/prism/_kernels.py +++ b/choclo/prism/_kernels.py @@ -889,7 +889,7 @@ def kernel_eee(easting, northing, upward, radius): k_{xxx}(x, y, z) = \frac{ - y z (2 x^2 + y^2 + z^2) + - y z (2 x^2 + y^2 + z^2) }{ (x^2 + y^2)(x^2 + z^2) \sqrt{x^2 + y^2 + z^2} @@ -941,7 +941,7 @@ def kernel_nnn(easting, northing, upward, radius): k_{yyy}(x, y, z) = \frac{ - x z (x^2 + 2 y^2 + z^2) + - x z (x^2 + 2 y^2 + z^2) }{ (x^2 + y^2)(y^2 + z^2) \sqrt{x^2 + y^2 + z^2} @@ -993,7 +993,7 @@ def kernel_uuu(easting, northing, upward, radius): k_{zzz}(x, y, z) = \frac{ - x y (x^2 + y^2 + 2 z^2) + - x y (x^2 + y^2 + 2 z^2) }{ (x^2 + z^2)(y^2 + z^2) \sqrt{x^2 + y^2 + z^2} @@ -1045,7 +1045,7 @@ def kernel_een(easting, northing, upward, radius): k_{xxy}(x, y, z) = \frac{ - x + - x }{ (z + \sqrt{x^2 + y^2 + z^2}) \sqrt{x^2 + y^2 + z^2} @@ -1097,7 +1097,7 @@ def kernel_eeu(easting, northing, upward, radius): k_{xxz}(x, y, z) = \frac{ - x + - x }{ (y + \sqrt{x^2 + y^2 + z^2}) \sqrt{x^2 + y^2 + z^2} @@ -1149,7 +1149,7 @@ def kernel_enn(easting, northing, upward, radius): k_{xyy}(x, y, z) = \frac{ - y + - y }{ (z + \sqrt{x^2 + y^2 + z^2}) \sqrt{x^2 + y^2 + z^2} @@ -1201,7 +1201,7 @@ def kernel_nnu(easting, northing, upward, radius): k_{yyz}(x, y, z) = \frac{ - y + - y }{ (x + \sqrt{x^2 + y^2 + z^2}) \sqrt{x^2 + y^2 + z^2} @@ -1253,7 +1253,7 @@ def kernel_euu(easting, northing, upward, radius): k_{xzz}(x, y, z) = \frac{ - z + - z }{ (y + \sqrt{x^2 + y^2 + z^2}) \sqrt{x^2 + y^2 + z^2} @@ -1305,7 +1305,7 @@ def kernel_nuu(easting, northing, upward, radius): k_{yzz}(x, y, z) = \frac{ - z + - z }{ (x + \sqrt{x^2 + y^2 + z^2}) \sqrt{x^2 + y^2 + z^2} @@ -1355,10 +1355,10 @@ def kernel_enu(easting, northing, upward, radius): .. math:: - k_{xyz}(x, y, z) = \frac{1}{\sqrt{x^2 + y^2 + z^2}} + k_{xyz}(x, y, z) = \frac{-1}{\sqrt{x^2 + y^2 + z^2}} """ - return 1 / radius + return -1 / radius @jit(nopython=True) @@ -1372,7 +1372,7 @@ def _kernel_iii(x_i, x_j, x_k, radius): k_{iii}(x_i, x_j, x_k) = \frac{ - x_j x_k (2 x_i^2 + x_j^2 + x_k^2) + - 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} @@ -1409,8 +1409,10 @@ def _kernel_iii(x_i, x_j, x_k, 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) + 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 @@ -1428,7 +1430,7 @@ def _kernel_iij(x_i, x_j, x_k, radius): k_{iij}(x_i, x_j, x_k) = \frac{ - x_i + - x_i }{ (x_k + \sqrt{x_i^2 + x_j^2 + x_k^2}) \sqrt{x_i^2 + x_j^2 + x_k^2} @@ -1470,4 +1472,6 @@ def _kernel_iij(x_i, x_j, x_k, radius): >>> k_euu = _kernel_iij(upward, easting, northing, radius) >>> k_nuu = _kernel_iij(upward, northing, upward, radius) """ - return x_i / ((x_k + radius) * radius) + if x_i == 0 and x_j == 0: + return 0.0 + return -x_i / ((x_k + radius) * radius) From 30591a49182209f1023303a783e1b93ae512e25b Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 26 Sep 2024 16:46:43 -0700 Subject: [PATCH 12/34] Make grad funcs to return nan on singular points Singular points for the grad components should be the same as the ones for the magnetic field components: interior points, edges and vertices. --- choclo/prism/_magnetic.py | 72 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 69 insertions(+), 3 deletions(-) diff --git a/choclo/prism/_magnetic.py b/choclo/prism/_magnetic.py index 86b1162..6e29568 100644 --- a/choclo/prism/_magnetic.py +++ b/choclo/prism/_magnetic.py @@ -1073,7 +1073,29 @@ def magnetic_ee( :func:`choclo.prism.magnetic_nn` :func:`choclo.prism.magnetic_uu` """ - # TODO: Check if observation point falls in a singular point + # 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 # Initialize magnetic gradiometry component b_ee = 0.0 # Iterate over the vertices of the prism @@ -1217,7 +1239,29 @@ def magnetic_nn( :func:`choclo.prism.magnetic_ee` :func:`choclo.prism.magnetic_uu` """ - # TODO: Check if observation point falls in a singular point + # 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 # Initialize magnetic gradiometry component b_nn = 0.0 # Iterate over the vertices of the prism @@ -1361,7 +1405,29 @@ def magnetic_uu( :func:`choclo.prism.magnetic_ee` :func:`choclo.prism.magnetic_nn` """ - # TODO: Check if observation point falls in a singular point + # 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 # Initialize magnetic gradiometry component b_uu = 0.0 # Iterate over the vertices of the prism From 0949c9831ab35555025c093d90acddbfa8fe1fba Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 26 Sep 2024 16:47:36 -0700 Subject: [PATCH 13/34] Fix atol value for finite difference test --- choclo/tests/test_prism_magnetic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index 27235bc..cd69248 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -821,4 +821,4 @@ def test(self, sample_3d_grid, sample_prism, sample_magnetization): forward_func=magnetic_e, delta=self.delta, ) - npt.assert_allclose(b_ee, b_ee_finite_diff, atol=1e-14) + npt.assert_allclose(b_ee, b_ee_finite_diff, atol=1e-11) From c98c4322c7474c70a30e327e39437440a54690d4 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 26 Sep 2024 17:13:35 -0700 Subject: [PATCH 14/34] Extend finite difference tests --- choclo/tests/test_prism_magnetic.py | 36 +++++++++++++++++++++++------ 1 file changed, 29 insertions(+), 7 deletions(-) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index cd69248..653594c 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -810,15 +810,37 @@ class TestMagGradiometryFiniteDifferences: """ delta = 1e-6 # displacement used in the finite difference calculations - - def test(self, sample_3d_grid, sample_prism, sample_magnetization): - b_ee = evaluate(magnetic_ee, sample_3d_grid, sample_prism, sample_magnetization) - b_ee_finite_diff = finite_differences( + rtol, atol = 5e-4, 1e-13 # tolerances used in the comparisons + + @pytest.mark.parametrize("direction_1", ["e", "n", "u"]) + @pytest.mark.parametrize("direction_2", ["e", "n", "u"]) + def test_derivatives_of_magnetic_components( + self, + sample_3d_grid, + sample_prism, + sample_magnetization, + direction_1, + direction_2, + ): + # Get magnetic forward function + mag_func = globals()[f"magnetic_{direction_1}"] + # Get magnetic gradiometry forward function + grad_component = "".join(sorted(f"{direction_1}{direction_2}")) + mag_grad_func = globals()[f"magnetic_{grad_component}"] + # 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="e", - forward_func=magnetic_e, + direction=direction_2, + forward_func=mag_func, delta=self.delta, ) - npt.assert_allclose(b_ee, b_ee_finite_diff, atol=1e-11) + # Compare results + npt.assert_allclose( + mag_grad, mag_grad_finite_diff, rtol=self.rtol, atol=self.atol + ) From a1343e7b769c6b3b1920b26a1cd97893c8444a91 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 26 Sep 2024 17:18:01 -0700 Subject: [PATCH 15/34] Add docstring to test method --- choclo/tests/test_prism_magnetic.py | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index 653594c..a388f75 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -812,20 +812,23 @@ class TestMagGradiometryFiniteDifferences: delta = 1e-6 # displacement used in the finite difference calculations rtol, atol = 5e-4, 1e-13 # tolerances used in the comparisons - @pytest.mark.parametrize("direction_1", ["e", "n", "u"]) - @pytest.mark.parametrize("direction_2", ["e", "n", "u"]) + @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, - direction_1, - direction_2, + self, sample_3d_grid, sample_prism, sample_magnetization, i, j ): + """ + Compare the magnetic gradiometry functions with a finite difference + approximation of it 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 magnetic forward function - mag_func = globals()[f"magnetic_{direction_1}"] + mag_func = globals()[f"magnetic_{i}"] # Get magnetic gradiometry forward function - grad_component = "".join(sorted(f"{direction_1}{direction_2}")) + grad_component = "".join(sorted(f"{i}{j}")) mag_grad_func = globals()[f"magnetic_{grad_component}"] # Evaluate the mag grad function on the sample grid mag_grad = evaluate( @@ -836,7 +839,7 @@ def test_derivatives_of_magnetic_components( sample_3d_grid, sample_prism, sample_magnetization, - direction=direction_2, + direction=j, forward_func=mag_func, delta=self.delta, ) From 22adcfe226c38a716dcf9b4c5fa32c9dbf19adf4 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 26 Sep 2024 17:18:13 -0700 Subject: [PATCH 16/34] Add new functions to the API index --- doc/api/index.rst | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/doc/api/index.rst b/doc/api/index.rst index c7ae9bb..6107e2b 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 From 3f36314ee220d9445daf56659cee766141963684 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 27 Sep 2024 09:51:42 -0700 Subject: [PATCH 17/34] Specify singular points in the docstrings --- choclo/prism/_magnetic.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/choclo/prism/_magnetic.py b/choclo/prism/_magnetic.py index 6e29568..dccce08 100644 --- a/choclo/prism/_magnetic.py +++ b/choclo/prism/_magnetic.py @@ -1026,8 +1026,8 @@ def magnetic_ee( 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 (TODO: specify which are singular points). + Return ``numpy.nan`` if the observation point falls in + a singular point: prism vertices, prism edges or interior points. Notes ----- @@ -1192,8 +1192,8 @@ def magnetic_nn( 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 (TODO: specify which are singular points). + Return ``numpy.nan`` if the observation point falls in + a singular point: prism vertices, prism edges or interior points. Notes ----- @@ -1358,8 +1358,8 @@ def magnetic_uu( 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 (TODO: specify which are singular points). + Return ``numpy.nan`` if the observation point falls in + a singular point: prism vertices, prism edges or interior points. Notes ----- From 3110b573e46e93311057d626e5443aef64cca364 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 27 Sep 2024 10:28:33 -0700 Subject: [PATCH 18/34] Add private function to compute a mag component Reduce repeated code by dumping the evaluation of the kernels on its own private function. --- choclo/prism/_magnetic.py | 400 +++++++++++++++++++------------------- 1 file changed, 204 insertions(+), 196 deletions(-) diff --git a/choclo/prism/_magnetic.py b/choclo/prism/_magnetic.py index dccce08..fd77773 100644 --- a/choclo/prism/_magnetic.py +++ b/choclo/prism/_magnetic.py @@ -514,43 +514,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( @@ -715,43 +696,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( @@ -916,43 +878,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( @@ -1096,43 +1039,24 @@ def magnetic_ee( prism_top, ): return np.nan - # Initialize magnetic gradiometry component - b_ee = 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_eee = kernel_eee(shift_east, shift_north, shift_upward, radius) - k_een = kernel_een(shift_east, shift_north, shift_upward, radius) - k_eeu = kernel_eeu(shift_east, shift_north, shift_upward, radius) - # Compute b_ee using the dot product between the kernel tensor - # and the magnetization vector of the prism - b_ee += (-1) ** (i + j + k) * ( - magnetization_east * k_eee - + magnetization_north * k_een - + magnetization_up * k_eeu - ) + # 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 @@ -1262,43 +1186,24 @@ def magnetic_nn( prism_top, ): return np.nan - # Initialize magnetic gradiometry component - b_nn = 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_enn = kernel_enn(shift_east, shift_north, shift_upward, radius) - k_nnn = kernel_nnn(shift_east, shift_north, shift_upward, radius) - k_nnu = kernel_nnu(shift_east, shift_north, shift_upward, radius) - # Compute b_nn using the dot product between the kernel tensor - # and the magnetization vector of the prism - b_nn += (-1) ** (i + j + k) * ( - magnetization_east * k_enn - + magnetization_north * k_nnn - + magnetization_up * k_nnu - ) + # 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 @@ -1428,8 +1333,111 @@ def magnetic_uu( prism_top, ): return np.nan - # Initialize magnetic gradiometry component - b_uu = 0.0 + # 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 _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, prism_south, prism_north, prism_bottom, prism_top : floats + The 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 @@ -1455,14 +1463,14 @@ def magnetic_uu( # Compute the radius radius = np.sqrt(shift_east_sq + shift_north_sq + shift_upward_sq) # Compute kernel tensor components for the current vertex - k_euu = kernel_euu(shift_east, shift_north, shift_upward, radius) - k_nuu = kernel_nuu(shift_east, shift_north, shift_upward, radius) - k_uuu = kernel_uuu(shift_east, shift_north, shift_upward, radius) - # Compute b_uu using the dot product between the kernel tensor + 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 - b_uu += (-1) ** (i + j + k) * ( - magnetization_east * k_euu - + magnetization_north * k_nuu - + magnetization_up * k_uuu + result += (-1) ** (i + j + k) * ( + magnetization_east * k_e + + magnetization_north * k_n + + magnetization_up * k_u ) - return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_uu + return result From 0fade53ff66933e89044740c6c708c515dbed067 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 27 Sep 2024 11:50:29 -0700 Subject: [PATCH 19/34] Add non-diagonal magnetic component derivatives --- choclo/prism/__init__.py | 3 + choclo/prism/_magnetic.py | 460 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 463 insertions(+) diff --git a/choclo/prism/__init__.py b/choclo/prism/__init__.py index 8f24af6..2a62868 100644 --- a/choclo/prism/__init__.py +++ b/choclo/prism/__init__.py @@ -49,4 +49,7 @@ magnetic_ee, magnetic_nn, magnetic_uu, + magnetic_en, + magnetic_eu, + magnetic_nu, ) diff --git a/choclo/prism/_magnetic.py b/choclo/prism/_magnetic.py index fd77773..e974422 100644 --- a/choclo/prism/_magnetic.py +++ b/choclo/prism/_magnetic.py @@ -27,6 +27,7 @@ kernel_nnu, kernel_euu, kernel_nuu, + kernel_enu, ) from ._utils import ( is_interior_point, @@ -1015,6 +1016,9 @@ def magnetic_ee( :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( @@ -1162,6 +1166,9 @@ def magnetic_nn( :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( @@ -1309,6 +1316,9 @@ def magnetic_uu( :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( @@ -1354,6 +1364,456 @@ def magnetic_uu( 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 (TODO: specify which are singular 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}{\partial i} + \frac{\partial}{\partial j} + \frac{\partial}{\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 (TODO: specify which are singular 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}{\partial i} + \frac{\partial}{\partial j} + \frac{\partial}{\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 (TODO: specify which are singular 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}{\partial i} + \frac{\partial}{\partial j} + \frac{\partial}{\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, From b1776d7220f22bd30c0fb57e056dc10b8a07ee67 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 27 Sep 2024 11:50:50 -0700 Subject: [PATCH 20/34] Update tests --- choclo/tests/test_prism_magnetic.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index a388f75..059875e 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -19,6 +19,9 @@ magnetic_ee, magnetic_nn, magnetic_uu, + magnetic_en, + magnetic_eu, + magnetic_nu, ) @@ -810,7 +813,7 @@ class TestMagGradiometryFiniteDifferences: """ delta = 1e-6 # displacement used in the finite difference calculations - rtol, atol = 5e-4, 1e-13 # tolerances used in the comparisons + 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"]) From 4f215d6f7312a68813d17d63cc3c7f2e02ee371d Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 27 Sep 2024 13:50:37 -0700 Subject: [PATCH 21/34] Minor change on partial derivatives in docstrings --- choclo/prism/_magnetic.py | 36 ++++++++++-------------------------- 1 file changed, 10 insertions(+), 26 deletions(-) diff --git a/choclo/prism/_magnetic.py b/choclo/prism/_magnetic.py index e974422..250d636 100644 --- a/choclo/prism/_magnetic.py +++ b/choclo/prism/_magnetic.py @@ -215,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 @@ -445,8 +444,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 @@ -627,8 +625,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 @@ -809,8 +806,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 @@ -991,9 +987,7 @@ def magnetic_ee( .. math:: u_{ijk} = - \frac{\partial}{\partial i} - \frac{\partial}{\partial j} - \frac{\partial}{\partial k} + \frac{\partial^3}{\partial i \partial j \partial k} \int\limits_R \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} dv @@ -1141,9 +1135,7 @@ def magnetic_nn( .. math:: u_{ijk} = - \frac{\partial}{\partial i} - \frac{\partial}{\partial j} - \frac{\partial}{\partial k} + \frac{\partial^3}{\partial i \partial j \partial k} \int\limits_R \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} dv @@ -1291,9 +1283,7 @@ def magnetic_uu( .. math:: u_{ijk} = - \frac{\partial}{\partial i} - \frac{\partial}{\partial j} - \frac{\partial}{\partial k} + \frac{\partial^3}{\partial i \partial j \partial k} \int\limits_R \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} dv @@ -1441,9 +1431,7 @@ def magnetic_en( .. math:: u_{ijk} = - \frac{\partial}{\partial i} - \frac{\partial}{\partial j} - \frac{\partial}{\partial k} + \frac{\partial^3}{\partial i \partial j \partial k} \int\limits_R \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} dv @@ -1591,9 +1579,7 @@ def magnetic_eu( .. math:: u_{ijk} = - \frac{\partial}{\partial i} - \frac{\partial}{\partial j} - \frac{\partial}{\partial k} + \frac{\partial^3}{\partial i \partial j \partial k} \int\limits_R \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} dv @@ -1741,9 +1727,7 @@ def magnetic_nu( .. math:: u_{ijk} = - \frac{\partial}{\partial i} - \frac{\partial}{\partial j} - \frac{\partial}{\partial k} + \frac{\partial^3}{\partial i \partial j \partial k} \int\limits_R \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} dv From 27c922ae3b8f4d7c900fb2a313a599a1c815003d Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 27 Sep 2024 14:08:53 -0700 Subject: [PATCH 22/34] Predefine components dict to avoid using globals --- choclo/tests/test_prism_magnetic.py | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index 059875e..f17a28c 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -814,6 +814,22 @@ class TestMagGradiometryFiniteDifferences: delta = 1e-6 # displacement used in the finite difference calculations rtol, atol = 5e-4, 5e-12 # tolerances used in the comparisons + 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(self, component: str): + """Return desired magnetic forward function.""" + component = "".join(sorted(component)) + return self.components[component] @pytest.mark.parametrize("i", ["e", "n", "u"]) @pytest.mark.parametrize("j", ["e", "n", "u"]) @@ -828,11 +844,9 @@ def test_derivatives_of_magnetic_components( against a finite difference approximation using the ``magnetic_{i}`` forward function, where ``i`` and ``j`` can be ``e``, ``n`` or ``u``. """ - # Get magnetic forward function - mag_func = globals()[f"magnetic_{i}"] - # Get magnetic gradiometry forward function - grad_component = "".join(sorted(f"{i}{j}")) - mag_grad_func = globals()[f"magnetic_{grad_component}"] + # Get forward functions + mag_func = self.get_forward_function(i) + mag_grad_func = self.get_forward_function(i + j) # Evaluate the mag grad function on the sample grid mag_grad = evaluate( mag_grad_func, sample_3d_grid, sample_prism, sample_magnetization From d6ab667a9d0b34fd1e2ae9eab194c0dd73943c5f Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 27 Sep 2024 14:11:21 -0700 Subject: [PATCH 23/34] Fix format --- choclo/prism/__init__.py | 28 ++++++++++++++-------------- choclo/prism/_magnetic.py | 26 +++++++++++++++----------- choclo/tests/test_prism_magnetic.py | 10 +++++----- 3 files changed, 34 insertions(+), 30 deletions(-) diff --git a/choclo/prism/__init__.py b/choclo/prism/__init__.py index 2a62868..4288e4d 100644 --- a/choclo/prism/__init__.py +++ b/choclo/prism/__init__.py @@ -22,34 +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_eee, - kernel_nnn, kernel_uuu, - kernel_een, - kernel_eeu, - kernel_enn, - kernel_nnu, - kernel_euu, - kernel_nuu, - kernel_enu, ) from ._magnetic import ( magnetic_e, - magnetic_field, - magnetic_n, - magnetic_u, magnetic_ee, - magnetic_nn, - magnetic_uu, magnetic_en, magnetic_eu, + magnetic_field, + magnetic_n, + magnetic_nn, magnetic_nu, + magnetic_u, + magnetic_uu, ) diff --git a/choclo/prism/_magnetic.py b/choclo/prism/_magnetic.py index 250d636..3b263fa 100644 --- a/choclo/prism/_magnetic.py +++ b/choclo/prism/_magnetic.py @@ -13,21 +13,21 @@ from ..constants import VACUUM_MAGNETIC_PERMEABILITY from ._kernels import ( kernel_ee, - kernel_en, - kernel_eu, - kernel_nn, - kernel_nu, - kernel_uu, kernel_eee, - kernel_nnn, - kernel_uuu, kernel_een, kernel_eeu, + kernel_en, kernel_enn, - kernel_nnu, + kernel_enu, + kernel_eu, kernel_euu, + kernel_nn, + kernel_nnn, + kernel_nnu, + kernel_nu, kernel_nuu, - kernel_enu, + kernel_uu, + kernel_uuu, ) from ._utils import ( is_interior_point, @@ -1827,8 +1827,12 @@ def _calculate_component( ---------- easting, northing, upward : floats Coordinates of the observation point. Must be in meters. - prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : floats - The boundaries of the prism. 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}`. diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index f17a28c..b2b27df 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -13,15 +13,15 @@ from ..prism import ( magnetic_e, - magnetic_field, - magnetic_n, - magnetic_u, magnetic_ee, - magnetic_nn, - magnetic_uu, magnetic_en, magnetic_eu, + magnetic_field, + magnetic_n, + magnetic_nn, magnetic_nu, + magnetic_u, + magnetic_uu, ) From 64a835abf88050654d2c8145a7b09dae7b7fa1bc Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Sat, 28 Sep 2024 10:21:23 -0700 Subject: [PATCH 24/34] Continue working on extending tests --- choclo/tests/test_prism_magnetic.py | 191 +++++++++++++++++++++++++--- 1 file changed, 173 insertions(+), 18 deletions(-) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index b2b27df..cafd1a1 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -7,6 +7,7 @@ """ Test magnetic forward modelling functions for rectangular prisms """ +import itertools import numpy as np import numpy.testing as npt import pytest @@ -24,6 +25,24 @@ 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(): @@ -814,22 +833,6 @@ class TestMagGradiometryFiniteDifferences: delta = 1e-6 # displacement used in the finite difference calculations rtol, atol = 5e-4, 5e-12 # tolerances used in the comparisons - 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(self, component: str): - """Return desired magnetic forward function.""" - component = "".join(sorted(component)) - return self.components[component] @pytest.mark.parametrize("i", ["e", "n", "u"]) @pytest.mark.parametrize("j", ["e", "n", "u"]) @@ -845,8 +848,8 @@ def test_derivatives_of_magnetic_components( forward function, where ``i`` and ``j`` can be ``e``, ``n`` or ``u``. """ # Get forward functions - mag_func = self.get_forward_function(i) - mag_grad_func = self.get_forward_function(i + j) + mag_func = get_forward_function(i) + mag_grad_func = get_forward_function(i + j) # Evaluate the mag grad function on the sample grid mag_grad = evaluate( mag_grad_func, sample_3d_grid, sample_prism, sample_magnetization @@ -864,3 +867,155 @@ def test_derivatives_of_magnetic_components( 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 are singular 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 + + def get_coords_on_line(self, prism, direction): + west, east, south, north, bottom, top = prism + distance = 13.4 # define distance from the nearest vertex for obs points + points = [] + if direction == "e": + for northing, upward in itertools.product((south, north), (bottom, top)): + points.append([west - distance, northing, upward]) + points.append([east + distance, northing, upward]) + elif direction == "n": + for easting, upward in itertools.product((west, east), (bottom, top)): + points.append([easting, south - distance, upward]) + points.append([easting, north + distance, upward]) + elif direction == "u": + for easting, northing in itertools.product((west, east), (south, north)): + points.append([easting, northing, bottom - distance]) + points.append([easting, northing, top + distance]) + else: + raise ValueError(f"Invalid direction '{direction}'.") + points = np.array(points).T + easting, northing, upward = points[:] + return easting, northing, upward + + @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.skip + @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. + + 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) + # 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=j, + forward_func=mag_func, + delta=self.delta, + ) + # Compare results + npt.assert_allclose( + mag_grad, mag_grad_finite_diff, rtol=self.rtol, atol=self.atol + ) + + def test_magnetic_ee(self, sample_prism, sample_magnetization): + """ + Compare the magnetic gradiometry functions with a finite difference + approximation of it 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``. + """ + coordinates = self.get_coords_on_line(sample_prism, direction="n") + # Evaluate the mag grad function on the sample grid + mag_grad = evaluate( + magnetic_ee, 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="e", + forward_func=magnetic_e, + delta=self.delta, + ) + # Compare results + npt.assert_allclose( + mag_grad, mag_grad_finite_diff, rtol=self.rtol, atol=self.atol + ) + + def test_magnetic_nn(self, sample_prism, sample_magnetization): + """ + Compare the magnetic gradiometry functions with a finite difference + approximation of it 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``. + """ + coordinates = self.get_coords_on_line(sample_prism, direction="u") + print(coordinates) + # Evaluate the mag grad function on the sample grid + mag_grad = evaluate( + magnetic_nn, 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="n", + forward_func=magnetic_n, + delta=self.delta, + ) + # Compare results + npt.assert_allclose( + mag_grad, mag_grad_finite_diff, rtol=self.rtol, atol=self.atol + ) From 6120d94cae711cdd2dcb6c6be65a377a7f94510d Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Wed, 2 Oct 2024 16:05:33 -0700 Subject: [PATCH 25/34] Run make format --- choclo/tests/test_prism_magnetic.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index f45c18b..55909df 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -8,6 +8,7 @@ Test magnetic forward modelling functions for rectangular prisms """ import itertools + import numpy as np import numpy.testing as npt import pytest @@ -16,15 +17,15 @@ kernel_en, kernel_eu, kernel_nu, - magnetic_field, - magnetic_n, - magnetic_u, magnetic_e, magnetic_ee, magnetic_en, magnetic_eu, + magnetic_field, + magnetic_n, magnetic_nn, magnetic_nu, + magnetic_u, magnetic_uu, ) @@ -1142,4 +1143,4 @@ def test_kernel_nu(self, prism, shift): diff_inline = kernel_inline[0] - kernel_inline[1] diff_shifted = kernel_shifted[0] - kernel_shifted[1] # These two differences should be close enough (not equal!) - npt.assert_allclose(diff_inline, diff_shifted, rtol=1e-7) \ No newline at end of file + npt.assert_allclose(diff_inline, diff_shifted, rtol=1e-7) From d3deb46855af519a9780c46b82623db9fd41b71b Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Wed, 2 Oct 2024 17:17:10 -0700 Subject: [PATCH 26/34] Add tests of mag grad components on faces --- choclo/tests/test_prism_magnetic.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index 55909df..9d14e91 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -828,6 +828,34 @@ 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:]) + class TestMagGradiometryFiniteDifferences: """ From 81620821c21b9c36afa84c7efed532e1586ded80 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Wed, 2 Oct 2024 17:17:41 -0700 Subject: [PATCH 27/34] Remove print statements in tests --- choclo/tests/test_prism_magnetic.py | 1 - 1 file changed, 1 deletion(-) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index 9d14e91..96b1af7 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -1033,7 +1033,6 @@ def test_magnetic_nn(self, sample_prism, sample_magnetization): forward function, where ``i`` and ``j`` can be ``e``, ``n`` or ``u``. """ coordinates = self.get_coords_on_line(sample_prism, direction="u") - print(coordinates) # Evaluate the mag grad function on the sample grid mag_grad = evaluate( magnetic_nn, coordinates, sample_prism, sample_magnetization From a6b70acfe4b97f504f6ad992bb6ca0c0ef330f15 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 3 Oct 2024 13:19:38 -0700 Subject: [PATCH 28/34] Fix minus sign and use atol Fix minus sign applied to a list and set an atol so we can compare values close to zero. --- choclo/tests/test_prism_magnetic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index 96b1af7..22f7782 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -854,7 +854,7 @@ def test_gradiometry_symmetry_on_faces(self, sample_prism, direction, forward_fu forward_func(e, n, u, *sample_prism, *magnetization) for (e, n, u) in zip(easting, northing, upward) ) - npt.assert_allclose(results[0], -results[1:]) + npt.assert_allclose(-results[0], results[1:], atol=1e-23) class TestMagGradiometryFiniteDifferences: From d1fa434ee63d5d1f87935bdb2bbdcdbd3d045c8a Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 3 Oct 2024 16:10:30 -0700 Subject: [PATCH 29/34] Test mag grad components on singular points Test if evaluation of the mag grad components on singular points return nans. --- choclo/tests/test_prism_magnetic.py | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index 22f7782..f71c248 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -676,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): """ @@ -752,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 @@ -767,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): """ @@ -785,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 From 40ba96f0871424a955a0b50fcfaac0ecaeb32ab4 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 4 Oct 2024 12:00:39 -0700 Subject: [PATCH 30/34] Extend lists of See Also in mag functions --- choclo/prism/_magnetic.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/choclo/prism/_magnetic.py b/choclo/prism/_magnetic.py index 3b263fa..f3f8693 100644 --- a/choclo/prism/_magnetic.py +++ b/choclo/prism/_magnetic.py @@ -242,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( @@ -489,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( @@ -670,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( @@ -851,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( From 0dc074f81f35d02c200c7aace5fd66e5f598d617 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 4 Oct 2024 12:03:41 -0700 Subject: [PATCH 31/34] Specify singular points in mag grad components --- choclo/prism/_magnetic.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/choclo/prism/_magnetic.py b/choclo/prism/_magnetic.py index f3f8693..d8ff39c 100644 --- a/choclo/prism/_magnetic.py +++ b/choclo/prism/_magnetic.py @@ -1435,7 +1435,7 @@ def magnetic_en( 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 (TODO: specify which are singular points). + point: prism vertices, prism edges or interior points. Notes ----- @@ -1583,7 +1583,7 @@ def magnetic_eu( 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 (TODO: specify which are singular points). + point: prism vertices, prism edges or interior points. Notes ----- @@ -1731,7 +1731,7 @@ def magnetic_nu( 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 (TODO: specify which are singular points). + point: prism vertices, prism edges or interior points. Notes ----- From d10cf705930e1dcfb9da17926353b766e53968e6 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 4 Oct 2024 12:34:42 -0700 Subject: [PATCH 32/34] Cleanup tests on the limits of mag grad components --- choclo/tests/test_prism_magnetic.py | 93 +++-------------------------- 1 file changed, 8 insertions(+), 85 deletions(-) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index f71c248..44a2c95 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -880,7 +880,7 @@ def test_derivatives_of_magnetic_components( ): """ Compare the magnetic gradiometry functions with a finite difference - approximation of it computed from the magnetic field components. + 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}`` @@ -912,37 +912,15 @@ class TestMagGradiometryLimits: """ Test if limits are well evaluated in magnetic gradiometry kernels. - Most of the third-order kernels are singular 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 + 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 - def get_coords_on_line(self, prism, direction): - west, east, south, north, bottom, top = prism - distance = 13.4 # define distance from the nearest vertex for obs points - points = [] - if direction == "e": - for northing, upward in itertools.product((south, north), (bottom, top)): - points.append([west - distance, northing, upward]) - points.append([east + distance, northing, upward]) - elif direction == "n": - for easting, upward in itertools.product((west, east), (bottom, top)): - points.append([easting, south - distance, upward]) - points.append([easting, north + distance, upward]) - elif direction == "u": - for easting, northing in itertools.product((west, east), (south, north)): - points.append([easting, northing, bottom - distance]) - points.append([easting, northing, top + distance]) - else: - raise ValueError(f"Invalid direction '{direction}'.") - points = np.array(points).T - easting, northing, upward = points[:] - return easting, northing, upward - @pytest.fixture def coordinates(self, sample_prism): """ @@ -967,7 +945,6 @@ def coordinates(self, sample_prism): easting, northing, upward = points[:] return easting, northing, upward - @pytest.mark.skip @pytest.mark.parametrize("i", ["e", "n", "u"]) @pytest.mark.parametrize("j", ["e", "n", "u"]) def test_derivatives_of_magnetic_components_on_edges( @@ -975,7 +952,9 @@ def test_derivatives_of_magnetic_components_on_edges( ): """ Compare the magnetic gradiometry functions with a finite difference - approximation of it computed from the magnetic field components. + 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}`` @@ -983,7 +962,7 @@ def test_derivatives_of_magnetic_components_on_edges( """ # Get forward functions mag_func = get_forward_function(i) - mag_grad_func = get_forward_function(i + j) + mag_grad_func = get_forward_function(f"{i}{j}") # Evaluate the mag grad function on the sample grid mag_grad = evaluate( mag_grad_func, coordinates, sample_prism, sample_magnetization @@ -1002,62 +981,6 @@ def test_derivatives_of_magnetic_components_on_edges( mag_grad, mag_grad_finite_diff, rtol=self.rtol, atol=self.atol ) - def test_magnetic_ee(self, sample_prism, sample_magnetization): - """ - Compare the magnetic gradiometry functions with a finite difference - approximation of it 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``. - """ - coordinates = self.get_coords_on_line(sample_prism, direction="n") - # Evaluate the mag grad function on the sample grid - mag_grad = evaluate( - magnetic_ee, 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="e", - forward_func=magnetic_e, - delta=self.delta, - ) - # Compare results - npt.assert_allclose( - mag_grad, mag_grad_finite_diff, rtol=self.rtol, atol=self.atol - ) - - def test_magnetic_nn(self, sample_prism, sample_magnetization): - """ - Compare the magnetic gradiometry functions with a finite difference - approximation of it 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``. - """ - coordinates = self.get_coords_on_line(sample_prism, direction="u") - # Evaluate the mag grad function on the sample grid - mag_grad = evaluate( - magnetic_nn, 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="n", - forward_func=magnetic_n, - delta=self.delta, - ) - # Compare results - npt.assert_allclose( - mag_grad, mag_grad_finite_diff, rtol=self.rtol, atol=self.atol - ) - class TestBugfixKernelEvaluation: r""" From aaf6f655164ab5958da0d5c984ec31dbb564de04 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 4 Oct 2024 12:37:42 -0700 Subject: [PATCH 33/34] Verbose variable name for direction of finite diff --- choclo/tests/test_prism_magnetic.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index 44a2c95..0b46c3d 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -889,6 +889,7 @@ def test_derivatives_of_magnetic_components( # 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 @@ -898,7 +899,7 @@ def test_derivatives_of_magnetic_components( sample_3d_grid, sample_prism, sample_magnetization, - direction=j, + direction=direction, forward_func=mag_func, delta=self.delta, ) @@ -963,6 +964,7 @@ def test_derivatives_of_magnetic_components_on_edges( # 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 @@ -972,7 +974,7 @@ def test_derivatives_of_magnetic_components_on_edges( coordinates, sample_prism, sample_magnetization, - direction=j, + direction=direction, forward_func=mag_func, delta=self.delta, ) From c397ffce4c4abb216c898fb3c67835cf09884a75 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 7 Oct 2024 13:43:03 -0700 Subject: [PATCH 34/34] Add doctests to new kernels This way I quickly verified that the mathematical expressions in the docstrings match the implementations. --- choclo/prism/_kernels.py | 60 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/choclo/prism/_kernels.py b/choclo/prism/_kernels.py index 168ba0d..b3f8f84 100644 --- a/choclo/prism/_kernels.py +++ b/choclo/prism/_kernels.py @@ -949,6 +949,12 @@ def kernel_eee(easting, northing, upward, radius): \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) @@ -1001,6 +1007,12 @@ def kernel_nnn(easting, northing, upward, radius): \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) @@ -1053,6 +1065,12 @@ def kernel_uuu(easting, northing, upward, radius): \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) @@ -1105,6 +1123,12 @@ def kernel_een(easting, northing, upward, radius): \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) @@ -1157,6 +1181,12 @@ def kernel_eeu(easting, northing, upward, radius): \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) @@ -1209,6 +1239,12 @@ def kernel_enn(easting, northing, upward, radius): \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) @@ -1261,6 +1297,12 @@ def kernel_nnu(easting, northing, upward, radius): \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) @@ -1313,6 +1355,12 @@ def kernel_euu(easting, northing, upward, radius): \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) @@ -1365,6 +1413,12 @@ def kernel_nuu(easting, northing, upward, radius): \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) @@ -1411,6 +1465,12 @@ def kernel_enu(easting, northing, upward, radius): 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