From 8d50e5af630b4604646deb965254dedc532a5dd0 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 May 2024 14:42:13 -0700 Subject: [PATCH 01/19] Fix style errors in docstrings (#85) Fix style for `density` argument in prism functions. Capitalize See Also sections in prism functions. Remove double empty lines in docstrings. --- choclo/prism/_gravity.py | 20 ++++++++++---------- choclo/prism/_kernels.py | 1 - choclo/prism/_magnetic.py | 9 ++++----- 3 files changed, 14 insertions(+), 16 deletions(-) diff --git a/choclo/prism/_gravity.py b/choclo/prism/_gravity.py index 241fac48..d03840ed 100644 --- a/choclo/prism/_gravity.py +++ b/choclo/prism/_gravity.py @@ -72,7 +72,7 @@ def gravity_pot( The bottom boundary of the prism. Must be in meters. prism_top : float The top boundary of the prism. Must be in meters. - density: float + density : float Density of the rectangular prism in kilograms per cubic meter. Returns @@ -218,7 +218,7 @@ def gravity_e( The bottom boundary of the prism. Must be in meters. prism_top : float The top boundary of the prism. Must be in meters. - density: float + density : float Density of the rectangular prism in kilograms per cubic meter. Returns @@ -365,7 +365,7 @@ def gravity_n( The bottom boundary of the prism. Must be in meters. prism_top : float The top boundary of the prism. Must be in meters. - density: float + density : float Density of the rectangular prism in kilograms per cubic meter. Returns @@ -512,7 +512,7 @@ def gravity_u( The bottom boundary of the prism. Must be in meters. prism_top : float The top boundary of the prism. Must be in meters. - density: float + density : float Density of the rectangular prism in kilograms per cubic meter. Returns @@ -665,7 +665,7 @@ def gravity_ee( The bottom boundary of the prism. Must be in meters. prism_top : float The top boundary of the prism. Must be in meters. - density: float + density : float Density of the rectangular prism in kilograms per cubic meter. Returns @@ -836,7 +836,7 @@ def gravity_nn( The bottom boundary of the prism. Must be in meters. prism_top : float The top boundary of the prism. Must be in meters. - density: float + density : float Density of the rectangular prism in kilograms per cubic meter. Returns @@ -1007,7 +1007,7 @@ def gravity_uu( The bottom boundary of the prism. Must be in meters. prism_top : float The top boundary of the prism. Must be in meters. - density: float + density : float Density of the rectangular prism in kilograms per cubic meter. Returns @@ -1178,7 +1178,7 @@ def gravity_en( The bottom boundary of the prism. Must be in meters. prism_top : float The top boundary of the prism. Must be in meters. - density: float + density : float Density of the rectangular prism in kilograms per cubic meter. Returns @@ -1325,7 +1325,7 @@ def gravity_eu( The bottom boundary of the prism. Must be in meters. prism_top : float The top boundary of the prism. Must be in meters. - density: float + density : float Density of the rectangular prism in kilograms per cubic meter. Returns @@ -1473,7 +1473,7 @@ def gravity_nu( The bottom boundary of the prism. Must be in meters. prism_top : float The top boundary of the prism. Must be in meters. - density: float + density : float Density of the rectangular prism in kilograms per cubic meter. Returns diff --git a/choclo/prism/_kernels.py b/choclo/prism/_kernels.py index 94032f1b..9646c72e 100644 --- a/choclo/prism/_kernels.py +++ b/choclo/prism/_kernels.py @@ -366,7 +366,6 @@ def kernel_u(easting, northing, upward, radius): by [Nagy2000]_ in order to compute the numerical kernel for the **upward** component instead for the downward one. - References ---------- - [Nagy2000]_ diff --git a/choclo/prism/_magnetic.py b/choclo/prism/_magnetic.py index 7a946db9..7431b622 100644 --- a/choclo/prism/_magnetic.py +++ b/choclo/prism/_magnetic.py @@ -205,7 +205,6 @@ def magnetic_field( \mathbf{B}(\mathbf{p}) = \frac{\mu_0}{4\pi} \mathbf{U} \cdot \mathbf{M} - References ---------- - [Blakely1995]_ @@ -214,7 +213,7 @@ def magnetic_field( - [Nagy2002]_ - [Fukushima2020]_ - See also + See Also -------- :func:`choclo.prism.magnetic_e` :func:`choclo.prism.magnetic_n` @@ -462,7 +461,7 @@ def magnetic_e( - [Nagy2002]_ - [Fukushima2020]_ - See also + See Also -------- :func:`choclo.prism.magnetic_field` :func:`choclo.prism.magnetic_n` @@ -663,7 +662,7 @@ def magnetic_n( - [Nagy2002]_ - [Fukushima2020]_ - See also + See Also -------- :func:`choclo.prism.magnetic_field` :func:`choclo.prism.magnetic_e` @@ -864,7 +863,7 @@ def magnetic_u( - [Nagy2002]_ - [Fukushima2020]_ - See also + See Also -------- :func:`choclo.prism.magnetic_field` :func:`choclo.prism.magnetic_e` From 403808bf7af4a7f9e4a93fbab93da2110d28eebe Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 30 May 2024 08:31:00 -0700 Subject: [PATCH 02/19] Improve docstrings of `magnetic_field` functions (#87) Improve the description of the outputs of the `magnetic_field` functions for dipoles and prisms. List each element of the tuple that is returned as separate returns. --- choclo/dipole/_forward.py | 14 +++++++++----- choclo/prism/_magnetic.py | 22 +++++++++++++++------- 2 files changed, 24 insertions(+), 12 deletions(-) diff --git a/choclo/dipole/_forward.py b/choclo/dipole/_forward.py index 90527257..c30f9857 100644 --- a/choclo/dipole/_forward.py +++ b/choclo/dipole/_forward.py @@ -61,11 +61,15 @@ def magnetic_field( Returns ------- - b : array - Array containing the three components of the magnetic field generated - by the dipole on the observation point in :math:`\text{T}`. - The components are returned in the following order: ``b_e``, ``b_n``, - ``b_u``. + b_e : float + Easting component of the magnetic field generated by the dipole + on the observation point in :math:`\text{T}`. + b_n : float + Northing component of the magnetic field generated by the dipole + on the observation point in :math:`\text{T}`. + b_u : float + Upward component of the magnetic field generated by the dipole + on the observation point in :math:`\text{T}`. Notes ----- diff --git a/choclo/prism/_magnetic.py b/choclo/prism/_magnetic.py index 7431b622..3d358d7a 100644 --- a/choclo/prism/_magnetic.py +++ b/choclo/prism/_magnetic.py @@ -81,13 +81,21 @@ def magnetic_field( Returns ------- - b_e, b_n, b_u : floats - The three components of the magnetic field generated by the prism on - the observation point in :math:`\text{T}`. - The components are returned in the following order: ``b_e``, ``b_n``, - ``b_u``. - Return a tuple full of ``numpy.nan`` if the observation point falls in - a singular point: prism vertices, prism edges or interior points. + b_e : float + Easting component of the magnetic field generated by the prism + on the observation point in :math:`\text{T}`. + It will be ``numpy.nan`` if the observation point falls in a singular + point: prism vertices, prism edges or interior points. + b_n : float + Northing component of the magnetic field generated by the prism + on the observation point in :math:`\text{T}`. + It will be ``numpy.nan`` if the observation point falls in a singular + point: prism vertices, prism edges or interior points. + b_u : float + Upward component of the magnetic field generated by the prism + on the observation point in :math:`\text{T}`. + It will be ``numpy.nan`` if the observation point falls in a singular + point: prism vertices, prism edges or interior points. Notes ----- From c3f2188d30d6d471e2c3609e0ea4227bf32d671e Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Tue, 18 Jun 2024 11:00:42 -0300 Subject: [PATCH 03/19] Bump pypa/gh-action-pypi-publish from 1.8.14 to 1.9.0 (#89) --- .github/workflows/pypi.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/pypi.yml b/.github/workflows/pypi.yml index 3e375e46..36c2d286 100644 --- a/.github/workflows/pypi.yml +++ b/.github/workflows/pypi.yml @@ -102,7 +102,7 @@ jobs: - name: Publish to Test PyPI # Only publish to TestPyPI when a PR is merged (pushed to main) if: success() && github.event_name == 'push' - uses: pypa/gh-action-pypi-publish@v1.8.14 + uses: pypa/gh-action-pypi-publish@v1.9.0 with: repository_url: https://test.pypi.org/legacy/ # Allow existing releases on test PyPI without errors. @@ -112,4 +112,4 @@ jobs: - name: Publish to PyPI # Only publish to PyPI when a release triggers the build if: success() && github.event_name == 'release' - uses: pypa/gh-action-pypi-publish@v1.8.14 + uses: pypa/gh-action-pypi-publish@v1.9.0 From b23a0281247c43a7d6aa59954f75622ac694d61c Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 21 Jun 2024 09:34:16 -0700 Subject: [PATCH 04/19] Update how output variables are stored in Actions (#90) Replace the deprecated `set-output` in GitHub Actions for concatenating name and values into the `GITHUB_OUTPUT` env variable. --- .github/workflows/docs.yml | 2 +- .github/workflows/test.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 1aa7fabb..2ea13841 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -79,7 +79,7 @@ jobs: - name: Get the pip cache folder id: pip-cache run: | - echo "::set-output name=dir::$(pip cache dir)" + echo "dir=$(pip cache dir)" >> $GITHUB_OUTPUT - name: Setup caching for pip packages uses: actions/cache@v4 diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 06b3bb64..c0ef1d93 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -115,7 +115,7 @@ jobs: - name: Get the pip cache folder id: pip-cache run: | - echo "::set-output name=dir::$(pip cache dir)" + echo "dir=$(pip cache dir)" >> $GITHUB_OUTPUT - name: Setup caching for pip packages uses: actions/cache@v4 From 5f0b36c377ac416156760278d38f771209a11c3f Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Aug 2024 12:27:57 -0300 Subject: [PATCH 05/19] Move push to codecov to its own job in Actions (#88) Remove the push to codecov step from the `test` job into a new job that depends on the test job. Upload the coverage reports as artifacts after testing, and reuse the artifacts in the new job. Upload all coverage reports in a single push to Codecov to minimize the number of hits. --- .github/workflows/test.yml | 35 +++++++++++++++++++++++++++++++---- 1 file changed, 31 insertions(+), 4 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index c0ef1d93..6cbd3440 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -28,7 +28,7 @@ defaults: jobs: ############################################################################# - # Run tests and upload to codecov + # Run tests test: name: ${{ matrix.os }} python=${{ matrix.python }} dependencies=${{ matrix.dependencies }} runs-on: ${{ matrix.os }} @@ -146,11 +146,38 @@ jobs: - name: Convert coverage report to XML for codecov run: coverage xml - - name: Upload coverage to Codecov + - name: Upload coverage as artifact + uses: actions/upload-artifact@v4 + with: + name: coverage_${{ matrix.os }}_${{ matrix.dependencies }} + path: ./coverage.xml + + + ############################################################################# + # Upload coverage report to codecov + codecov-upload: + runs-on: ubuntu-latest + needs: test + + steps: + + - name: Download coverage report artifacts + # Download coverage reports from every runner. + # Maximum coverage is achieved by combining reports from every runner. + # Each coverage file will live in its own folder with the same name as + # the artifact. + uses: actions/download-artifact@v4 + with: + pattern: coverage_* + + - name: List all downloaded artifacts + run: ls -l -R . + + - name: Upload coverage reports to Codecov uses: codecov/codecov-action@v4 with: - file: ./coverage.xml - env_vars: OS,PYTHON,DEPENDENCIES + # Upload all coverage report files + files: ./coverage_*/coverage.xml # Fail the job so we know coverage isn't being updated. Otherwise it # can silently drop and we won't know. fail_ci_if_error: true From f933989865c2826a27f177a95669416ca752c664 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Tue, 3 Sep 2024 09:39:42 -0700 Subject: [PATCH 06/19] Bump pypa/gh-action-pypi-publish from 1.9.0 to 1.10.0 (#92) --- .github/workflows/pypi.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/pypi.yml b/.github/workflows/pypi.yml index 36c2d286..ace6b8b4 100644 --- a/.github/workflows/pypi.yml +++ b/.github/workflows/pypi.yml @@ -102,7 +102,7 @@ jobs: - name: Publish to Test PyPI # Only publish to TestPyPI when a PR is merged (pushed to main) if: success() && github.event_name == 'push' - uses: pypa/gh-action-pypi-publish@v1.9.0 + uses: pypa/gh-action-pypi-publish@v1.10.0 with: repository_url: https://test.pypi.org/legacy/ # Allow existing releases on test PyPI without errors. @@ -112,4 +112,4 @@ jobs: - name: Publish to PyPI # Only publish to PyPI when a release triggers the build if: success() && github.event_name == 'release' - uses: pypa/gh-action-pypi-publish@v1.9.0 + uses: pypa/gh-action-pypi-publish@v1.10.0 From 92319f9d39fe9e2509faaa810ac384f95b0dbf12 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 5 Sep 2024 15:12:45 -0700 Subject: [PATCH 07/19] Replace `build` for `python-build` in environment.yml (#91) The `build` package in `conda-forge` is outdated, `python-build` should be used instead. --- environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index 5274b2a9..ca823e77 100644 --- a/environment.yml +++ b/environment.yml @@ -9,7 +9,7 @@ dependencies: - numpy - numba # Build - - build + - python-build - twine # Test - pytest From eaf80e78329fabe52f45c047950d76d6fadb2b35 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Wed, 11 Sep 2024 14:56:20 -0700 Subject: [PATCH 08/19] Bump pypa/gh-action-pypi-publish from 1.10.0 to 1.10.1 (#93) --- .github/workflows/pypi.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/pypi.yml b/.github/workflows/pypi.yml index ace6b8b4..14597091 100644 --- a/.github/workflows/pypi.yml +++ b/.github/workflows/pypi.yml @@ -102,7 +102,7 @@ jobs: - name: Publish to Test PyPI # Only publish to TestPyPI when a PR is merged (pushed to main) if: success() && github.event_name == 'push' - uses: pypa/gh-action-pypi-publish@v1.10.0 + uses: pypa/gh-action-pypi-publish@v1.10.1 with: repository_url: https://test.pypi.org/legacy/ # Allow existing releases on test PyPI without errors. @@ -112,4 +112,4 @@ jobs: - name: Publish to PyPI # Only publish to PyPI when a release triggers the build if: success() && github.event_name == 'release' - uses: pypa/gh-action-pypi-publish@v1.10.0 + uses: pypa/gh-action-pypi-publish@v1.10.1 From a11b4082d06f7f1850b4046f0d49e6fa67deed6b Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 23 Sep 2024 16:32:02 -0700 Subject: [PATCH 09/19] Bump pypa/gh-action-pypi-publish from 1.10.1 to 1.10.2 (#95) --- .github/workflows/pypi.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/pypi.yml b/.github/workflows/pypi.yml index 14597091..2f1add19 100644 --- a/.github/workflows/pypi.yml +++ b/.github/workflows/pypi.yml @@ -102,7 +102,7 @@ jobs: - name: Publish to Test PyPI # Only publish to TestPyPI when a PR is merged (pushed to main) if: success() && github.event_name == 'push' - uses: pypa/gh-action-pypi-publish@v1.10.1 + uses: pypa/gh-action-pypi-publish@v1.10.2 with: repository_url: https://test.pypi.org/legacy/ # Allow existing releases on test PyPI without errors. @@ -112,4 +112,4 @@ jobs: - name: Publish to PyPI # Only publish to PyPI when a release triggers the build if: success() && github.event_name == 'release' - uses: pypa/gh-action-pypi-publish@v1.10.1 + uses: pypa/gh-action-pypi-publish@v1.10.2 From 0542ab8dc1e504b9f8cc843546377905f4ab80dc Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 24 Sep 2024 10:51:47 -0700 Subject: [PATCH 10/19] Simplify tests for prism magnetic forward funcs (#96) Simplify some of the tests for the magnetic forward modelling functions to reduce complexity, improve readability and set the ground to extend them. --- choclo/tests/test_prism_magnetic.py | 549 +++++++++------------------- 1 file changed, 181 insertions(+), 368 deletions(-) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index 9859e6a1..fcecdddb 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -68,24 +68,98 @@ def get_prism_center(prism): return easting, northing, upward +def evaluate(forward_func, coordinates, prism, magnetization): + """ + Evaluate a forward function on a set of observation points. + + Parameters + ---------- + forward_func : callable + Forward modelling function to evaluate. + coordinates : tuple of (n_data) arrays + Coordinates of the observation points. + prism : (6) array + Boundaries of the prism. + magnetization : (3) array + Magnetization vector of the prism. + + Returns + ------- + array + Array with the result of evaluating the ``forward_func`` on every + observation point. + """ + coordinates = tuple(c.ravel() for c in coordinates) + result = np.array( + list( + forward_func(e, n, u, *prism, *magnetization) + for e, n, u in zip(*coordinates) + ) + ) + return result + + +def finite_differences( + coordinates, prism, magnetization, direction, forward_func, delta=1e-6 +): + """ + Compute spatial derivatives through finite differences. + + Parameters + ---------- + coordinates : tuple of (n_data) arrays + Coordinates of the observation points. + prism : (6) array + Boundaries of the prism. + magnetization : (3) array + Magnetization vector of the prism. + direction : {"e", "n", "u"} + Direction along which take the derivative. + forward_func : callable + Forward modelling function to use to compute the finite difference. + delta : float, optional + Displacement for the finite differences. + + Returns + ------- + (n_data) array + Array with the derivatives of the ``forward_func`` on each observation + point calculated using finite differences. + """ + # Get original coordinates + easting, northing, upward = coordinates + if direction == "e": + shifted_coords = (easting + delta, northing, upward) + elif direction == "n": + shifted_coords = (easting, northing + delta, upward) + elif direction == "u": + shifted_coords = (easting, northing, upward + delta) + else: + ValueError(f"Invalid direction '{direction}'") + # Compute field on original and shifted coordinates + field = evaluate(forward_func, coordinates, prism, magnetization) + field_shifted = evaluate(forward_func, shifted_coords, prism, magnetization) + # Compute spatial derivative + spatial_derivative = (field_shifted - field) / delta + return spatial_derivative + + class TestSymmetryBe: """ Test symmetry of easting component of the magnetic field """ atol = 1e-17 # absolute tolerance for values near zero - - @pytest.mark.parametrize( - "magnetization", - [ - (500, 0, 0), - (-500, 0, 0), - (0, 500, 0), - (0, -500, 0), - (0, 0, 500), - (0, 0, -500), - ], - ) + MAGNETIZATIONS = [ + (500, 0, 0), + (-500, 0, 0), + (0, 500, 0), + (0, -500, 0), + (0, 0, 500), + (0, 0, -500), + ] + + @pytest.mark.parametrize("magnetization", MAGNETIZATIONS) def test_symmetry_across_easting_northing( self, sample_3d_grid, sample_prism, magnetization ): @@ -105,21 +179,11 @@ def test_symmetry_across_easting_northing( upward_bottom = 2 * prism_center_u - upward_top # Compute magnetic_e on every observation point magnetization = np.array(magnetization) - b_e_top = np.array( - list( - magnetic_e(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting.ravel(), northing.ravel(), upward_top.ravel() - ) - ) + b_e_top = evaluate( + magnetic_e, (easting, northing, upward_top), sample_prism, magnetization ) - b_e_bottom = np.array( - list( - magnetic_e(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting.ravel(), northing.ravel(), upward_bottom.ravel() - ) - ) + b_e_bottom = evaluate( + magnetic_e, (easting, northing, upward_bottom), sample_prism, magnetization ) # Check symmetry between top and bottom if magnetization[2] != 0: @@ -127,17 +191,7 @@ def test_symmetry_across_easting_northing( else: npt.assert_allclose(b_e_top, b_e_bottom, atol=self.atol) - @pytest.mark.parametrize( - "magnetization", - [ - (500, 0, 0), - (-500, 0, 0), - (0, 500, 0), - (0, -500, 0), - (0, 0, 500), - (0, 0, -500), - ], - ) + @pytest.mark.parametrize("magnetization", MAGNETIZATIONS) def test_symmetry_across_easting_upward( self, sample_3d_grid, sample_prism, magnetization ): @@ -157,21 +211,11 @@ def test_symmetry_across_easting_upward( northing_south = 2 * prism_center_n - northing_north # Compute magnetic_e on every observation point magnetization = np.array(magnetization) - b_e_north = np.array( - list( - magnetic_e(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting.ravel(), northing_north.ravel(), upward.ravel() - ) - ) + b_e_north = evaluate( + magnetic_e, (easting, northing_north, upward), sample_prism, magnetization ) - b_e_south = np.array( - list( - magnetic_e(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting.ravel(), northing_south.ravel(), upward.ravel() - ) - ) + b_e_south = evaluate( + magnetic_e, (easting, northing_south, upward), sample_prism, magnetization ) # Check symmetry between south and north if magnetization[1] != 0: @@ -179,17 +223,7 @@ def test_symmetry_across_easting_upward( else: npt.assert_allclose(b_e_north, b_e_south, atol=self.atol) - @pytest.mark.parametrize( - "magnetization", - [ - (500, 0, 0), - (-500, 0, 0), - (0, 500, 0), - (0, -500, 0), - (0, 0, 500), - (0, 0, -500), - ], - ) + @pytest.mark.parametrize("magnetization", MAGNETIZATIONS) def test_symmetry_across_northing_upward( self, sample_3d_grid, sample_prism, magnetization ): @@ -209,21 +243,11 @@ def test_symmetry_across_northing_upward( easting_west = 2 * prism_center_e - easting_east # Compute magnetic_e on every observation point magnetization = np.array(magnetization) - b_e_east = np.array( - list( - magnetic_e(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting_east.ravel(), northing.ravel(), upward.ravel() - ) - ) + b_e_east = evaluate( + magnetic_e, (easting_east, northing, upward), sample_prism, magnetization ) - b_e_west = np.array( - list( - magnetic_e(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting_west.ravel(), northing.ravel(), upward.ravel() - ) - ) + b_e_west = evaluate( + magnetic_e, (easting_west, northing, upward), sample_prism, magnetization ) # Check symmetry between west and east if magnetization[0] != 0: @@ -247,17 +271,11 @@ def test_symmetry_when_flipping(self, sample_3d_grid, sample_prism): magnetization_west = np.array([-500.0, 0, 0]) # Compute the magnetic field generated by each moment on the # observation points - b_e_east = np.array( - list( - magnetic_e(e, n, u, *sample_prism, *magnetization_east) - for e, n, u in zip(easting, northing, upward) - ) + b_e_east = evaluate( + magnetic_e, (easting, northing, upward), sample_prism, magnetization_east ) - b_e_west = np.array( - list( - magnetic_e(e, n, u, *sample_prism, *magnetization_west) - for e, n, u in zip(easting, northing, upward) - ) + b_e_west = evaluate( + magnetic_e, (easting, northing, upward), sample_prism, magnetization_west ) # Check if the sign gets inverted npt.assert_allclose(b_e_east, -b_e_west) @@ -269,18 +287,16 @@ class TestSymmetryBn: """ atol = 1e-17 # absolute tolerance for values near zero - - @pytest.mark.parametrize( - "magnetization", - [ - (500, 0, 0), - (-500, 0, 0), - (0, 500, 0), - (0, -500, 0), - (0, 0, 500), - (0, 0, -500), - ], - ) + MAGNETIZATIONS = [ + (500, 0, 0), + (-500, 0, 0), + (0, 500, 0), + (0, -500, 0), + (0, 0, 500), + (0, 0, -500), + ] + + @pytest.mark.parametrize("magnetization", MAGNETIZATIONS) def test_symmetry_across_easting_northing( self, sample_3d_grid, sample_prism, magnetization ): @@ -300,21 +316,11 @@ def test_symmetry_across_easting_northing( upward_bottom = 2 * prism_center_u - upward_top # Compute magnetic_n on every observation point magnetization = np.array(magnetization) - b_n_top = np.array( - list( - magnetic_n(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting.ravel(), northing.ravel(), upward_top.ravel() - ) - ) + b_n_top = evaluate( + magnetic_n, (easting, northing, upward_top), sample_prism, magnetization ) - b_n_bottom = np.array( - list( - magnetic_n(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting.ravel(), northing.ravel(), upward_bottom.ravel() - ) - ) + b_n_bottom = evaluate( + magnetic_n, (easting, northing, upward_bottom), sample_prism, magnetization ) # Check symmetry between top and bottom if magnetization[2] != 0: @@ -322,17 +328,7 @@ def test_symmetry_across_easting_northing( else: npt.assert_allclose(b_n_top, b_n_bottom, atol=self.atol) - @pytest.mark.parametrize( - "magnetization", - [ - (500, 0, 0), - (-500, 0, 0), - (0, 500, 0), - (0, -500, 0), - (0, 0, 500), - (0, 0, -500), - ], - ) + @pytest.mark.parametrize("magnetization", MAGNETIZATIONS) def test_symmetry_across_easting_upward( self, sample_3d_grid, sample_prism, magnetization ): @@ -352,21 +348,11 @@ def test_symmetry_across_easting_upward( northing_south = 2 * prism_center_n - northing_north # Compute magnetic_n on every observation point magnetization = np.array(magnetization) - b_n_north = np.array( - list( - magnetic_n(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting.ravel(), northing_north.ravel(), upward.ravel() - ) - ) + b_n_north = evaluate( + magnetic_n, (easting, northing_north, upward), sample_prism, magnetization ) - b_n_south = np.array( - list( - magnetic_n(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting.ravel(), northing_south.ravel(), upward.ravel() - ) - ) + b_n_south = evaluate( + magnetic_n, (easting, northing_south, upward), sample_prism, magnetization ) # Check symmetry between south and north if magnetization[1] != 0: @@ -374,17 +360,7 @@ def test_symmetry_across_easting_upward( else: npt.assert_allclose(b_n_north, -b_n_south, atol=self.atol) - @pytest.mark.parametrize( - "magnetization", - [ - (500, 0, 0), - (-500, 0, 0), - (0, 500, 0), - (0, -500, 0), - (0, 0, 500), - (0, 0, -500), - ], - ) + @pytest.mark.parametrize("magnetization", MAGNETIZATIONS) def test_symmetry_across_northing_upward( self, sample_3d_grid, sample_prism, magnetization ): @@ -404,21 +380,11 @@ def test_symmetry_across_northing_upward( easting_west = 2 * prism_center_e - easting_east # Compute magnetic_n on every observation point magnetization = np.array(magnetization) - b_n_east = np.array( - list( - magnetic_n(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting_east.ravel(), northing.ravel(), upward.ravel() - ) - ) + b_n_east = evaluate( + magnetic_n, (easting_east, northing, upward), sample_prism, magnetization ) - b_n_west = np.array( - list( - magnetic_n(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting_west.ravel(), northing.ravel(), upward.ravel() - ) - ) + b_n_west = evaluate( + magnetic_n, (easting_west, northing, upward), sample_prism, magnetization ) # Check symmetry between west and east if magnetization[0] != 0: @@ -442,17 +408,11 @@ def test_symmetry_when_flipping(self, sample_3d_grid, sample_prism): magnetization_south = np.array([0, -500.0, 0]) # Compute the magnetic field generated by each moment on the # observation points - b_n_north = np.array( - list( - magnetic_n(e, n, u, *sample_prism, *magnetization_north) - for e, n, u in zip(easting, northing, upward) - ) + b_n_north = evaluate( + magnetic_n, (easting, northing, upward), sample_prism, magnetization_north ) - b_n_south = np.array( - list( - magnetic_n(e, n, u, *sample_prism, *magnetization_south) - for e, n, u in zip(easting, northing, upward) - ) + b_n_south = evaluate( + magnetic_n, (easting, northing, upward), sample_prism, magnetization_south ) # Check if the sign gets inverted npt.assert_allclose(b_n_north, -b_n_south) @@ -464,18 +424,16 @@ class TestSymmetryBu: """ atol = 1e-17 # absolute tolerance for values near zero - - @pytest.mark.parametrize( - "magnetization", - [ - (500, 0, 0), - (-500, 0, 0), - (0, 500, 0), - (0, -500, 0), - (0, 0, 500), - (0, 0, -500), - ], - ) + MAGNETIZATIONS = [ + (500, 0, 0), + (-500, 0, 0), + (0, 500, 0), + (0, -500, 0), + (0, 0, 500), + (0, 0, -500), + ] + + @pytest.mark.parametrize("magnetization", MAGNETIZATIONS) def test_symmetry_across_easting_northing( self, sample_3d_grid, sample_prism, magnetization ): @@ -495,21 +453,11 @@ def test_symmetry_across_easting_northing( upward_bottom = 2 * prism_center_u - upward_top # Compute magnetic_u on every observation point magnetization = np.array(magnetization) - b_u_top = np.array( - list( - magnetic_u(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting.ravel(), northing.ravel(), upward_top.ravel() - ) - ) + b_u_top = evaluate( + magnetic_u, (easting, northing, upward_top), sample_prism, magnetization ) - b_u_bottom = np.array( - list( - magnetic_u(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting.ravel(), northing.ravel(), upward_bottom.ravel() - ) - ) + b_u_bottom = evaluate( + magnetic_u, (easting, northing, upward_bottom), sample_prism, magnetization ) # Check symmetry between top and bottom if magnetization[2] != 0: @@ -517,17 +465,7 @@ def test_symmetry_across_easting_northing( else: npt.assert_allclose(b_u_top, -b_u_bottom, atol=self.atol) - @pytest.mark.parametrize( - "magnetization", - [ - (500, 0, 0), - (-500, 0, 0), - (0, 500, 0), - (0, -500, 0), - (0, 0, 500), - (0, 0, -500), - ], - ) + @pytest.mark.parametrize("magnetization", MAGNETIZATIONS) def test_symmetry_across_easting_upward( self, sample_3d_grid, sample_prism, magnetization ): @@ -547,21 +485,11 @@ def test_symmetry_across_easting_upward( northing_south = 2 * prism_center_n - northing_north # Compute magnetic_u on every observation point magnetization = np.array(magnetization) - b_u_north = np.array( - list( - magnetic_u(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting.ravel(), northing_north.ravel(), upward.ravel() - ) - ) + b_u_north = evaluate( + magnetic_u, (easting, northing_north, upward), sample_prism, magnetization ) - b_u_south = np.array( - list( - magnetic_u(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting.ravel(), northing_south.ravel(), upward.ravel() - ) - ) + b_u_south = evaluate( + magnetic_u, (easting, northing_south, upward), sample_prism, magnetization ) # Check symmetry between south and north if magnetization[1] != 0: @@ -569,17 +497,7 @@ def test_symmetry_across_easting_upward( else: npt.assert_allclose(b_u_north, b_u_south, atol=self.atol) - @pytest.mark.parametrize( - "magnetization", - [ - (500, 0, 0), - (-500, 0, 0), - (0, 500, 0), - (0, -500, 0), - (0, 0, 500), - (0, 0, -500), - ], - ) + @pytest.mark.parametrize("magnetization", MAGNETIZATIONS) def test_symmetry_across_northing_upward( self, sample_3d_grid, sample_prism, magnetization ): @@ -599,21 +517,11 @@ def test_symmetry_across_northing_upward( easting_west = 2 * prism_center_e - easting_east # Compute magnetic_u on every observation point magnetization = np.array(magnetization) - b_u_east = np.array( - list( - magnetic_u(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting_east.ravel(), northing.ravel(), upward.ravel() - ) - ) + b_u_east = evaluate( + magnetic_u, (easting_east, northing, upward), sample_prism, magnetization ) - b_u_west = np.array( - list( - magnetic_u(e, n, u, *sample_prism, *magnetization) - for e, n, u in zip( - easting_west.ravel(), northing.ravel(), upward.ravel() - ) - ) + b_u_west = evaluate( + magnetic_u, (easting_west, northing, upward), sample_prism, magnetization ) # Check symmetry between east and west if magnetization[0] != 0: @@ -637,17 +545,11 @@ def test_symmetry_when_flipping(self, sample_3d_grid, sample_prism): magnetization_down = np.array([0, 0, -500.0]) # Compute the magnetic field generated by each moment on the # observation points - b_u_up = np.array( - list( - magnetic_u(e, n, u, *sample_prism, *magnetization_up) - for e, n, u in zip(easting, northing, upward) - ) + b_u_up = evaluate( + magnetic_u, (easting, northing, upward), sample_prism, magnetization_up ) - b_u_down = np.array( - list( - magnetic_u(e, n, u, *sample_prism, *magnetization_down) - for e, n, u in zip(easting, northing, upward) - ) + b_u_down = evaluate( + magnetic_u, (easting, northing, upward), sample_prism, magnetization_down ) # Check if the sign gets inverted npt.assert_allclose(b_u_up, -b_u_down) @@ -667,31 +569,17 @@ def test_magnetic_field(self, sample_3d_grid, sample_prism, sample_magnetization Test magnetic_field against each one of the other functions """ # Compute all components of B using magnetic_field - b = np.array( - list( - magnetic_field(e, n, u, *sample_prism, *sample_magnetization) - for e, n, u in zip(*sample_3d_grid) - ) - ) + b = evaluate(magnetic_field, sample_3d_grid, sample_prism, sample_magnetization) b_e, b_n, b_u = tuple(b[:, i] for i in range(3)) # Computed the individual fields - b_e_expected = np.array( - list( - magnetic_e(e, n, u, *sample_prism, *sample_magnetization) - for e, n, u in zip(*sample_3d_grid) - ) + b_e_expected = evaluate( + magnetic_e, sample_3d_grid, sample_prism, sample_magnetization ) - b_n_expected = np.array( - list( - magnetic_n(e, n, u, *sample_prism, *sample_magnetization) - for e, n, u in zip(*sample_3d_grid) - ) + b_n_expected = evaluate( + magnetic_n, sample_3d_grid, sample_prism, sample_magnetization ) - b_u_expected = np.array( - list( - magnetic_u(e, n, u, *sample_prism, *sample_magnetization) - for e, n, u in zip(*sample_3d_grid) - ) + b_u_expected = evaluate( + magnetic_u, sample_3d_grid, sample_prism, sample_magnetization ) npt.assert_allclose(b_e, b_e_expected) npt.assert_allclose(b_n, b_n_expected) @@ -700,100 +588,25 @@ def test_magnetic_field(self, sample_3d_grid, sample_prism, sample_magnetization class TestDivergenceOfB: """ - Test if the divergence of the magnetic field is equal to zero - - Compute the derivatives of B through finite differences + Test if the divergence of the magnetic field is equal to zero. """ - # Displacement used in the finite differences - delta = 1e-6 - - def get_b_ee_finite_differences(self, coordinates, prism, magnetization): - """ - Compute b_ee using finite differences - """ - # Get original coordinates - easting, northing, upward = coordinates - # Shift coordinates using delta - easting_shifted = easting + self.delta - # Compute b_e on original and shifted coordinates - b_e = np.array( - list( - magnetic_e(e, n, u, *prism, *magnetization) - for e, n, u in zip(easting, northing, upward) - ) - ) - b_e_shifted = np.array( - list( - magnetic_e(e, n, u, *prism, *magnetization) - for e, n, u in zip(easting_shifted, northing, upward) - ) - ) - # Compute b_ee - b_ee = (b_e_shifted - b_e) / self.delta - return b_ee - - def get_b_nn_finite_differences(self, coordinates, prism, magnetization): - """ - Compute b_nn using finite differences - """ - # Get original coordinates - easting, northing, upward = coordinates - # Shift coordinates using delta - northing_shifted = northing + self.delta - # Compute b_e on original and shifted coordinates - b_n = np.array( - list( - magnetic_n(e, n, u, *prism, *magnetization) - for e, n, u in zip(easting, northing, upward) - ) - ) - b_n_shifted = np.array( - list( - magnetic_n(e, n, u, *prism, *magnetization) - for e, n, u in zip(easting, northing_shifted, upward) - ) - ) - # Compute b_nn - b_nn = (b_n_shifted - b_n) / self.delta - return b_nn - - def get_b_uu_finite_differences(self, coordinates, prism, magnetization): + def test_divergence_of_b_finite_differences( + self, sample_3d_grid, sample_prism, sample_magnetization + ): """ - Compute b_uu using finite differences + Test div of B computing its derivatives through finite differences. """ - # Get original coordinates - easting, northing, upward = coordinates - # Shift coordinates using delta - upward_shifted = upward + self.delta - # Compute b_e on original and shifted coordinates - b_u = np.array( - list( - magnetic_u(e, n, u, *prism, *magnetization) - for e, n, u in zip(easting, northing, upward) - ) - ) - b_u_shifted = np.array( - list( - magnetic_u(e, n, u, *prism, *magnetization) - for e, n, u in zip(easting, northing, upward_shifted) - ) - ) - # Compute b_uu - b_uu = (b_u_shifted - b_u) / self.delta - return b_uu - - def test_divergence_of_b(self, sample_3d_grid, sample_prism, sample_magnetization): - # Compute b_ee, b_nn and b_uu using finite differences - b_ee = self.get_b_ee_finite_differences( - sample_3d_grid, sample_prism, sample_magnetization - ) - b_nn = self.get_b_nn_finite_differences( - sample_3d_grid, sample_prism, sample_magnetization - ) - b_uu = self.get_b_uu_finite_differences( - sample_3d_grid, sample_prism, sample_magnetization + delta = 1e-6 + kwargs = dict( + coordinates=sample_3d_grid, + prism=sample_prism, + magnetization=sample_magnetization, + delta=delta, ) + b_ee = finite_differences(direction="e", forward_func=magnetic_e, **kwargs) + b_nn = finite_differences(direction="n", forward_func=magnetic_n, **kwargs) + b_uu = finite_differences(direction="u", forward_func=magnetic_u, **kwargs) # Check if the divergence of B is zero npt.assert_allclose(-b_uu, b_ee + b_nn, atol=1e-11) From 7bb37dc76c5e6e7dd0e07f884dfce87d20146a94 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Wed, 2 Oct 2024 16:01:26 -0700 Subject: [PATCH 11/19] Fix bug on `safe_log` and solves discontinuous magnetic fields (#100) Fix how `_safe_log` checks if `r == |x|`. Add the other two shifted coordinates (`y` and `z`) as arguments for `_safe_log` so we can use them to assess that the observation point is inline with a pair of vertices or not. Make the check by evaluating if `y == 0 and z == 0` instead of `-x == r` (with `x < 0`), since `-x` and `r` could be the same value (up to machine precision) while `y` and `z` are non-zero (specially when $|x| \gg |y|$ and/or $|x| \gg |z|$). Update calls to `_safe_log`, update its docstrings, and update docstrings of kernels that mentioned the safe log function. Add tests that catches the bug, and that pass with the bugfix. --- choclo/prism/_kernels.py | 176 ++++++++++++++++++---------- choclo/tests/test_prism_magnetic.py | 131 ++++++++++++++++++++- 2 files changed, 245 insertions(+), 62 deletions(-) diff --git a/choclo/prism/_kernels.py b/choclo/prism/_kernels.py index 9646c72e..4b1eea65 100644 --- a/choclo/prism/_kernels.py +++ b/choclo/prism/_kernels.py @@ -55,9 +55,9 @@ def kernel_pot(easting, northing, upward, radius): .. math:: k_V(x, y, z) &= - x y \, \operatorname{safe-ln} (z + r) - + y z \, \operatorname{safe-ln} (x + r) - + z x \, \operatorname{safe-ln} (y + r) \\ + x y \, \operatorname{safe\_ln} (z, r) + + y z \, \operatorname{safe\_ln} (x, r) + + z x \, \operatorname{safe\_ln} (y, r) \\ & - \frac{x^2}{2} \operatorname{safe-arctan} \left( yz, xr \right) - \frac{y^2}{2} \operatorname{safe-arctan} \left( zx, yr \right) - \frac{z^2}{2} \operatorname{safe-arctan} \left( xy, zr \right) @@ -66,10 +66,12 @@ def kernel_pot(easting, northing, upward, radius): .. math:: - \operatorname{safe-ln}(x) = + \operatorname{safe\_ln}(x, r) = \begin{cases} - 0 & |x| < 10^{-10} \\ - \ln (x) + 0 & r = 0 \\ + \ln(x + r) & x \ge 0 \\ + \ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\ + -\ln(-2 x) & x < 0, r = |x| \end{cases} and @@ -91,9 +93,9 @@ def kernel_pot(easting, northing, upward, radius): - [Fukushima2020]_ """ kernel = ( - easting * northing * _safe_log(upward, radius) - + northing * upward * _safe_log(easting, radius) - + easting * upward * _safe_log(northing, radius) + easting * northing * _safe_log(upward, easting, northing, radius) + + northing * upward * _safe_log(easting, northing, upward, radius) + + easting * upward * _safe_log(northing, upward, easting, radius) - 0.5 * easting**2 * _safe_atan2(upward * northing, easting * radius) - 0.5 * northing**2 * _safe_atan2(upward * easting, northing * radius) - 0.5 * upward**2 * _safe_atan2(easting * northing, upward * radius) @@ -147,8 +149,8 @@ def kernel_e(easting, northing, upward, radius): k_x(x, y, z) = -\left[ - y \, \operatorname{safe-ln} (z + r) - + z \, \operatorname{safe-ln} (y + r) + y \, \operatorname{safe\_ln} (z, r) + + z \, \operatorname{safe\_ln} (y, r) - x \, \operatorname{safe-arctan} \left( yz, xr \right) \right] @@ -156,10 +158,12 @@ def kernel_e(easting, northing, upward, radius): .. math:: - \operatorname{safe-ln}(x) = + \operatorname{safe\_ln}(x, r) = \begin{cases} - 0 & |x| < 10^{-10} \\ - \ln (x) + 0 & r = 0 \\ + \ln(x + r) & x \ge 0 \\ + \ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\ + -\ln(-2 x) & x < 0, r = |x| \end{cases} and @@ -187,8 +191,8 @@ def kernel_e(easting, northing, upward, radius): - [Fukushima2020]_ """ kernel = -( - northing * _safe_log(upward, radius) - + upward * _safe_log(northing, radius) + northing * _safe_log(upward, easting, northing, radius) + + upward * _safe_log(northing, upward, easting, radius) - easting * _safe_atan2(northing * upward, easting * radius) ) return kernel @@ -240,8 +244,8 @@ def kernel_n(easting, northing, upward, radius): k_y(x, y, z) = -\left[ - z \, \operatorname{safe-ln} (x + r) - + x \, \operatorname{safe-ln} (z + r) + z \, \operatorname{safe\_ln} (x, r) + + x \, \operatorname{safe\_ln} (z, r) - y \, \operatorname{safe-arctan} \left( zx, yr \right) \right] @@ -249,10 +253,12 @@ def kernel_n(easting, northing, upward, radius): .. math:: - \operatorname{safe-ln}(x) = + \operatorname{safe\_ln}(x, r) = \begin{cases} - 0 & |x| < 10^{-10} \\ - \ln (x) + 0 & r = 0 \\ + \ln(x + r) & x \ge 0 \\ + \ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\ + -\ln(-2 x) & x < 0, r = |x| \end{cases} and @@ -280,8 +286,8 @@ def kernel_n(easting, northing, upward, radius): - [Fukushima2020]_ """ kernel = -( - upward * _safe_log(easting, radius) - + easting * _safe_log(upward, radius) + upward * _safe_log(easting, northing, upward, radius) + + easting * _safe_log(upward, easting, northing, radius) - northing * _safe_atan2(easting * upward, northing * radius) ) return kernel @@ -333,8 +339,8 @@ def kernel_u(easting, northing, upward, radius): k_z(x, y, z) = - \left[ - x \, \operatorname{safe-ln} (y + r) - + y \, \operatorname{safe-ln} (x + r) + x \, \operatorname{safe\_ln} (y, r) + + y \, \operatorname{safe\_ln} (x, r) - z \, \operatorname{safe-arctan} \left( xy, zr \right) \right] @@ -342,10 +348,12 @@ def kernel_u(easting, northing, upward, radius): .. math:: - \operatorname{safe-ln}(x) = + \operatorname{safe\_ln}(x, r) = \begin{cases} - 0 & |x| < 10^{-10} \\ - \ln (x) + 0 & r = 0 \\ + \ln(x + r) & x \ge 0 \\ + \ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\ + -\ln(-2 x) & x < 0, r = |x| \end{cases} and @@ -375,8 +383,8 @@ def kernel_u(easting, northing, upward, radius): # The minus sign is to return the kernel for the upward component instead # of the downward one. kernel = -( - easting * _safe_log(northing, radius) - + northing * _safe_log(easting, radius) + easting * _safe_log(northing, upward, easting, radius) + + northing * _safe_log(easting, northing, upward, radius) - upward * _safe_atan2(easting * northing, upward * radius) ) return kernel @@ -619,16 +627,18 @@ def kernel_en(easting, northing, upward, radius): .. math:: - k_{xy}(x, y, z) = \operatorname{safe-ln} \left( z + r \right) + k_{xy}(x, y, z) = \operatorname{safe\_ln} \left(z, r \right) where .. math:: - \operatorname{safe-ln}(x) = + \operatorname{safe\_ln}(x, r) = \begin{cases} - 0 & |x| < 10^{-10} \\ - \ln (x) + 0 & r = 0 \\ + \ln(x + r) & x \ge 0 \\ + \ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\ + -\ln(-2 x) & x < 0, r = |x| \end{cases} References @@ -637,7 +647,7 @@ def kernel_en(easting, northing, upward, radius): - [Nagy2002]_ - [Fukushima2020]_ """ - return _safe_log(upward, radius) + return _safe_log(upward, easting, northing, radius) @jit(nopython=True) @@ -682,16 +692,18 @@ def kernel_eu(easting, northing, upward, radius): .. math:: - k_{xz}(x, y, z) = \operatorname{safe-ln} \left( y + r \right) + k_{xz}(x, y, z) = \operatorname{safe\_ln} \left(y, r \right) where .. math:: - \operatorname{safe-ln}(x) = + \operatorname{safe\_ln}(x, r) = \begin{cases} - 0 & |x| < 10^{-10} \\ - \ln (x) + 0 & r = 0 \\ + \ln(x + r) & x \ge 0 \\ + \ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\ + -\ln(-2 x) & x < 0, r = |x| \end{cases} References @@ -700,7 +712,7 @@ def kernel_eu(easting, northing, upward, radius): - [Nagy2002]_ - [Fukushima2020]_ """ - return _safe_log(northing, radius) + return _safe_log(northing, easting, upward, radius) @jit(nopython=True) @@ -745,16 +757,18 @@ def kernel_nu(easting, northing, upward, radius): .. math:: - k_{yz}(x, y, z) = \operatorname{safe-ln} \left( x + r \right) + k_{yz}(x, y, z) = \operatorname{safe\_ln} \left(x, r \right) where .. math:: - \operatorname{safe-ln}(x) = + \operatorname{safe\_ln}(x, r) = \begin{cases} - 0 & |x| < 10^{-10} \\ - \ln (x) + 0 & r = 0 \\ + \ln(x + r) & x \ge 0 \\ + \ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\ + -\ln(-2 x) & x < 0, r = |x| \end{cases} References @@ -763,7 +777,7 @@ def kernel_nu(easting, northing, upward, radius): - [Nagy2002]_ - [Fukushima2020]_ """ - return _safe_log(easting, radius) + return _safe_log(easting, northing, upward, radius) @jit(nopython=True) @@ -806,42 +820,82 @@ def _safe_atan2(y, x): @jit(nopython=True) -def _safe_log(x, r): +def _safe_log(x, y, z, r): r""" - Safe log function to use in the prism kernels + Safe log function to use in the prism kernels. - Evaluates the :math:`\ln{x + r}` where :math:`x` is one of the shifted + Evaluates the :math:`\ln(x + r)` where :math:`x` is one of the shifted coordinate of the prism vertex and :math:`r` is the Euclidean distance (always non-negative) from the vertex to the observation point. + Parameters + ---------- + x : float + Shifted coordinate of the prism vertex that will be used for evaluating + the log. + y, z : floats + The other two shifted coordinates of the prism vertex. They are used to + determine if ``abs(x) == r`` with more accuracy, and to evaluate the + log function to avoid floating point errors when :math:`x < 0` and + :math:`r \ne |x|`. + r : float + Euclidean distance from the vertex to the observation point. We require + this argument because this quantity is usually precomputed before + evaluating the kernel, thus saving computation time. + + Returns + ------- + float + Evaluation of the :math:`\ln{x + r}` to be used on kernels. + + Notes + ----- + The :math:`\text{safe\_ln}(x, r)` function is defined as: + .. math:: - \text{safe_ln}(x, r) = + \text{safe\_ln}(x, r) = \begin{cases} - 0 & x = 0, r = 0 \\ + 0 & r = 0 \\ \ln(x + r) & x \ge 0 \\ - \ln((r^2 - x^2) / (r - x)) & x < 0, r \ne -x \\ - -\ln(-2 x) & x < 0, r = -x + \ln((y^2 + z^2) / (r - x)) & x < 0, r \ne |x| \\ + -\ln(-2 x) & x < 0, r = |x| \end{cases} - This function returns 0 when the observation point is located on the vertex - of the prism (:math:`r=0`); and two modified versions in case that - :math:`x` is negative: if :math:`x = -r` then the :math:`\ln{x + r}` can be - replaced by :math:`-\ln{|x| + r} = -\ln(-2x)`, and for any other negative - value of :math:`x` it returns :math:`\ln((r^2 - x^2) / (r - x))` which - helps by reducing the floating point errors ([Fukushima2020_]). - This modified version was inspired by [Nagy2000] and [Fukushima2020_]: + The function evaluates :math:`\ln(x + r)` when :math:`x \ge 0`. Otherwise, + the returned value takes into account the limit cases: + + * If :math:`r = 0` (observation point is on one of the vertices of the + prism), it returns zero. This accounts for the limits of the + zeroth-order and first-order kernels when evaluated on the vertices of + the prism. The second-order kernels are not defined on the vertices. + * If :math:`x < 0, r \ne |x|`, then it evaluates + :math:`\ln((y^2 + z^2) / (r - x))` to reduce floating point errors + ([Fukushima2020]_). + * If :math:`x < 0, r = |x|`, then it returns one of the terms of the limit + of the second-order kernels (when the observation point is inline with + two nodes): :math:`-\ln(-2x)`. + + When checking if :math:`r = |x|`, we will evaluate if ``y == 0.0 and z == + 0.0``, to avoid issues with floating point errors in cases in which + :math:`|x| \gg |y|` and/or :math:`|x| \gg |z|`. In such cases, the ``r`` + value could be exactly the same as ``x`` (up to machine precision) and + not reflect that ``y`` and/or ``z`` are not-null. + + This modified version of the log function was inspired by [Nagy2000]_ and + [Fukushima2020]_. References ---------- + - [Nagy2000]_ - [Fukushima2020]_ """ if r == 0: return 0 if x < 0: - if r == -x: + if y == 0.0 and z == 0.0: result = -np.log(-2 * x) else: - result = np.log((r**2 - x**2) / (r - x)) + result = np.log((y**2 + z**2) / (r - x)) return result return np.log(x + r) diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index fcecdddb..224cb9c1 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 ( + kernel_en, + kernel_eu, + kernel_nu, + magnetic_e, + magnetic_field, + magnetic_n, + magnetic_u, +) @pytest.fixture(name="sample_prism") @@ -778,3 +786,124 @@ 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 TestBugfixKernelEvaluation: + r""" + Tests that ensure that the bug related to wrong evaluation of the + non-diagonal second-order kernels was fixed. + + This bug is due to wrongly evaluating the safe_log function on two of the + vertices of the prism. + """ + + @pytest.fixture + def prism(self): + return np.array([-10.0, 10.0, -12.0, 12.0, -20.0, -5.0]) + + def evaluate_kernel(self, coordinates, vertex, kernel): + """ + Evaluate a given kernel. + + Parameters + ---------- + coordinates : tuple of floats + Coordinates of the observation point. + vertex : tuple of floats + Coordinates of the prism vertex. + kernel : callable + Kernel function to evaluate. + + Returns + ------- + float + """ + shifted_east, shifted_north, shifted_upward = tuple( + vertex[i] - coordinates[i] for i in range(3) + ) + radius = np.sqrt(shifted_east**2 + shifted_north**2 + shifted_upward**2) + return kernel(shifted_east, shifted_north, shifted_upward, radius) + + @pytest.mark.parametrize("shift", ("easting", "northing")) + def test_kernel_en(self, prism, shift): + """ + Test bugfix on kernel_en. + """ + # Get the coordinates of the prism and collect only two of the vertices + _, east, south, _, bottom, top = prism[:] + vertices = [[east, south, top], [east, south, bottom]] + # Build observation points: + # one above the vertices, one slightly shifted + delta = 1e-6 + if shift == "easting": + coords = [(east, south, top + 50.0), (east + delta, south, top + 50.0)] + else: + coords = [(east, south, top + 50.0), (east, south + delta, top + 50.0)] + # Evaluate kernel on the two nodes + kernel_above = [ + self.evaluate_kernel(coords[0], vertex, kernel_en) for vertex in vertices + ] + kernel_shifted = [ + self.evaluate_kernel(coords[1], vertex, kernel_en) for vertex in vertices + ] + # Compute the differences between the evaluations on both nodes + diff_above = kernel_above[0] - kernel_above[1] + diff_shifted = kernel_shifted[0] - kernel_shifted[1] + # These two differences should be close enough (not equal!) + npt.assert_allclose(diff_above, diff_shifted, rtol=1e-7) + + @pytest.mark.parametrize("shift", ("easting", "upward")) + def test_kernel_eu(self, prism, shift): + """ + Test bugfix on kernel_eu. + """ + # Get the coordinates of the prism and collect only two of the vertices + _, east, south, north, _, top = prism[:] + vertices = [[east, north, top], [east, south, top]] + # Build observation points: + # one to the north of the vertices, one slightly shifted. + delta = 1e-6 + if shift == "easting": + coords = [(east, north + 50.0, top), (east + delta, north + 50.0, top)] + else: + coords = [(east, north + 50.0, top), (east, north + 50.0, top + delta)] + # Evaluate kernel on the two nodes + kernel_inline = [ + self.evaluate_kernel(coords[0], vertex, kernel_eu) for vertex in vertices + ] + kernel_shifted = [ + self.evaluate_kernel(coords[1], vertex, kernel_eu) for vertex in vertices + ] + # Compute the differences between the evaluations on both nodes + 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) + + @pytest.mark.parametrize("shift", ("northing", "upward")) + def test_kernel_nu(self, prism, shift): + """ + Test bugfix on kernel_nu. + """ + # Get the coordinates of the prism and collect only two of the vertices + west, east, _, north, _, top = prism[:] + vertices = [[east, north, top], [west, north, top]] + # Build observation points: + # one to the east of the vertices, one slightly shifted. + delta = 1e-6 + if shift == "northing": + coords = [(east + 50.0, north, top), (east + 50.0, north + delta, top)] + else: + coords = [(east + 50.0, north, top), (east + 50.0, north, top + delta)] + # Evaluate kernel on the two nodes + kernel_inline = [ + self.evaluate_kernel(coords[0], vertex, kernel_nu) for vertex in vertices + ] + kernel_shifted = [ + self.evaluate_kernel(coords[1], vertex, kernel_nu) for vertex in vertices + ] + # Compute the differences between the evaluations on both nodes + 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) From 74aa55c0a63e31550c6067881eef8140a6dd33fa Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 7 Oct 2024 10:59:04 -0700 Subject: [PATCH 12/19] Add some more tests for _safe_log (#101) Add some additional tests for the `_safe_log` function. --- choclo/tests/test_prism_safe_log.py | 92 +++++++++++++++++++++++++++++ 1 file changed, 92 insertions(+) create mode 100644 choclo/tests/test_prism_safe_log.py diff --git a/choclo/tests/test_prism_safe_log.py b/choclo/tests/test_prism_safe_log.py new file mode 100644 index 00000000..c146700b --- /dev/null +++ b/choclo/tests/test_prism_safe_log.py @@ -0,0 +1,92 @@ +# Copyright (c) 2022 The Choclo Developers. +# Distributed under the terms of the BSD 3-Clause License. +# SPDX-License-Identifier: BSD-3-Clause +# +# This code is part of the Fatiando a Terra project (https://www.fatiando.org) +# +""" +Test the _safe_log private function for prisms. +""" + +import numpy as np +import numpy.testing as npt +import pytest + +from ..prism._kernels import _safe_log + + +class TestSafeLogFixedX: + r""" + Test some properties of ``_safe_log`` when ``x`` is fixed. + + We define: + + .. math:: + + f(x, y, z) = \ln(x + \sqrt{x^2 + y^2 + z^2}) + + When fixing ``x``, this function has some continuity, symmetry and monotony + properties. + + For simplicity, let's define: + + .. math:: + + g(z) = f(x=x_0, y=0, z) = \ln(x_0 + \sqrt{x_0^2 + z^2}) + + with :math:`x_0 = \text{constant}`. + """ + + def eval_safe_log(self, z, x_0=0.0): + r = np.sqrt(x_0**2 + z**2) + return _safe_log(x_0, 0.0, z, r) + + @pytest.mark.parametrize( + "x_0", (-10.0, 0.0, 10.0), ids=("x_0=-10.0", "x_0=0.0", "x_0=10.0") + ) + def test_monotony_z_negative(self, x_0): + """ + Test monotony of the safe_log function with negative z values. + + Function :math:`g(x, z)` with :math:`x` fixed monotonically decreases + on :math:`z < 0`. Without fixes to reduce the floating point errors, + the safe_log function won't be monotone. + """ + z = np.linspace(-1e-3, 0, 1001) + z = z[z < 0] + results = np.array([self.eval_safe_log(z_i, x_0) for z_i in z]) + assert (results[1:] < results[:-1]).all() + + @pytest.mark.parametrize( + "x_0", (-10.0, 0.0, 10.0), ids=("x_0=-10.0", "x_0=0.0", "x_0=10.0") + ) + def test_monotony_z_positive(self, x_0): + """ + Test monotony of the safe_log function with positive z values. + + Function :math:`g(x, z)` with :math:`x` fixed monotonically increases + on :math:`z > 0`. Without fixes to reduce the floating point errors, + the safe_log function won't be monotone. + """ + z = np.linspace(0, 1e-3, 1001) + z = z[z > 0] + results = np.array([self.eval_safe_log(z_i, x_0) for z_i in z]) + assert (results[1:] > results[:-1]).all() + + @pytest.mark.parametrize( + "x_0", (-10.0, 0.0, 10.0), ids=("x_0=-10.0", "x_0=0.0", "x_0=10.0") + ) + def test_symmetry(self, x_0): + """ + Test symmetry of the function with respect to zero. + + Function :math:`g(x, z)` with :math:`x` fixed is symmetric with respect + to :math:`z = 0`. + """ + vmax = 1e-3 + z = np.linspace(-vmax, vmax, 1001) + z = z[z != 0] + results = np.array([self.eval_safe_log(z_i, x_0) for z_i in z]) + results_z_negative = results[z < 0] + results_z_positive = results[z > 0] + npt.assert_allclose(results_z_negative, results_z_positive[::-1]) From 36b36e16cb346545f01f8a95aae8df01a5d64fee Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 7 Oct 2024 14:02:10 -0700 Subject: [PATCH 13/19] Add magnetic gradiometry forwards for prisms (#97) Add forward modelling functions for the magnetic gradiometry components due to rectangular prisms: spatial derivatives of the magnetic field components. Include third-order kernel functions that are needed to compute the derivatives of the magnetic field components. Create private function to evaluate the kernels on the vertices of the prism provided its magnetization so we can reduce the amount of repeated code in the magnetic forward modelling functions. Add tests for the new functions. --- choclo/prism/__init__.py | 23 +- choclo/prism/_kernels.py | 690 +++++++++++++++ choclo/prism/_magnetic.py | 1234 ++++++++++++++++++++++++--- choclo/tests/test_prism_magnetic.py | 214 ++++- doc/api/index.rst | 16 + 5 files changed, 2047 insertions(+), 130 deletions(-) diff --git a/choclo/prism/__init__.py b/choclo/prism/__init__.py index 04f37b8a..4288e4d4 100644 --- a/choclo/prism/__init__.py +++ b/choclo/prism/__init__.py @@ -22,13 +22,34 @@ from ._kernels import ( kernel_e, kernel_ee, + kernel_eee, + kernel_een, + kernel_eeu, kernel_en, + kernel_enn, + kernel_enu, kernel_eu, + kernel_euu, kernel_n, kernel_nn, + kernel_nnn, + kernel_nnu, kernel_nu, + kernel_nuu, kernel_pot, kernel_u, kernel_uu, + kernel_uuu, +) +from ._magnetic import ( + magnetic_e, + magnetic_ee, + magnetic_en, + magnetic_eu, + magnetic_field, + magnetic_n, + magnetic_nn, + magnetic_nu, + magnetic_u, + magnetic_uu, ) -from ._magnetic import magnetic_e, magnetic_field, magnetic_n, magnetic_u diff --git a/choclo/prism/_kernels.py b/choclo/prism/_kernels.py index 4b1eea65..b3f8f845 100644 --- a/choclo/prism/_kernels.py +++ b/choclo/prism/_kernels.py @@ -899,3 +899,693 @@ def _safe_log(x, y, z, r): result = np.log((y**2 + z**2) / (r - x)) return result return np.log(x + r) + + +@jit(nopython=True) +def kernel_eee(easting, northing, upward, radius): + r""" + Easting-easting-easting kernel due to a prism. + + Evaluates the integration kernel for the easting-easting-easting component + of the grad tensor generated by a prism [Nagy2000]_ on a single vertex of + the prism. The coordinates that must be passed are shifted coordinates: the + coordinates of the vertex from a Cartesian coordinate system whose origin + is located in the observation point. + + Parameters + ---------- + easting : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate of the vertex of the prism. Must be in + meters. + radius : float + Square root of the sum of the squares of the ``easting``, ``northing`` + and ``upward`` shifted coordinates. + + Returns + ------- + kernel : float + Value of the kernel function for the easting-easting-easting component + of the grad tensor due to a rectangular prism evaluated on a single + vertex. + + Notes + ----- + Computes the following numerical kernel on the passed *shifted + coordinates*: + + .. math:: + + k_{xxx}(x, y, z) = + \frac{ + - y z (2 x^2 + y^2 + z^2) + }{ + (x^2 + y^2)(x^2 + z^2) + \sqrt{x^2 + y^2 + z^2} + } + + Examples + -------- + >>> x, y, z = 3.1, 5.2, -3.0 + >>> r = np.sqrt(x**2 + y**2 + z**2) + >>> float(kernel_eee(x, y, z, r)) # doctest: +NUMBER + 0.18706595 + """ + return _kernel_iii(easting, northing, upward, radius) + + +@jit(nopython=True) +def kernel_nnn(easting, northing, upward, radius): + r""" + Northing-northing-northing kernel due to a prism. + + Evaluates the integration kernel for the northing-northing-northing + component of the grad tensor generated by a prism [Nagy2000]_ on a single + vertex of the prism. The coordinates that must be passed are shifted + coordinates: the coordinates of the vertex from a Cartesian coordinate + system whose origin is located in the observation point. + + Parameters + ---------- + easting : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate of the vertex of the prism. Must be in + meters. + radius : float + Square root of the sum of the squares of the ``easting``, ``northing`` + and ``upward`` shifted coordinates. + + Returns + ------- + kernel : float + Value of the kernel function for the northing-northing-northing + component of the grad tensor due to a rectangular prism evaluated on + a single vertex. + + Notes + ----- + Computes the following numerical kernel on the passed *shifted + coordinates*: + + .. math:: + + k_{yyy}(x, y, z) = + \frac{ + - x z (x^2 + 2 y^2 + z^2) + }{ + (x^2 + y^2)(y^2 + z^2) + \sqrt{x^2 + y^2 + z^2} + } + + Examples + -------- + >>> x, y, z = 3.1, 5.2, -3.0 + >>> r = np.sqrt(x**2 + y**2 + z**2) + >>> float(kernel_nnn(x, y, z, r)) # doctest: +NUMBER + 0.07574927 + """ + return _kernel_iii(northing, upward, easting, radius) + + +@jit(nopython=True) +def kernel_uuu(easting, northing, upward, radius): + r""" + Upward-upward-upward kernel due to a prism. + + Evaluates the integration kernel for the upward-upward-upward + component of the grad tensor generated by a prism [Nagy2000]_ on a single + vertex of the prism. The coordinates that must be passed are shifted + coordinates: the coordinates of the vertex from a Cartesian coordinate + system whose origin is located in the observation point. + + Parameters + ---------- + easting : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate of the vertex of the prism. Must be in + meters. + radius : float + Square root of the sum of the squares of the ``easting``, ``northing`` + and ``upward`` shifted coordinates. + + Returns + ------- + kernel : float + Value of the kernel function for the upward-upward-upward component of + the grad tensor due to a rectangular prism evaluated on a single + vertex. + + Notes + ----- + Computes the following numerical kernel on the passed *shifted + coordinates*: + + .. math:: + + k_{zzz}(x, y, z) = + \frac{ + - x y (x^2 + y^2 + 2 z^2) + }{ + (x^2 + z^2)(y^2 + z^2) + \sqrt{x^2 + y^2 + z^2} + } + + Examples + -------- + >>> x, y, z = 3.1, 5.2, -3.0 + >>> r = np.sqrt(x**2 + y**2 + z**2) + >>> float(kernel_uuu(x, y, z, r)) # doctest: +NUMBER + -0.19440331 + """ + return _kernel_iii(upward, northing, easting, radius) + + +@jit(nopython=True) +def kernel_een(easting, northing, upward, radius): + r""" + Easting-easting-northing kernel due to a prism. + + Evaluates the integration kernel for the easting-easting-northing + component of the grad tensor generated by a prism [Nagy2000]_ on a single + vertex of the prism. The coordinates that must be passed are shifted + coordinates: the coordinates of the vertex from a Cartesian coordinate + system whose origin is located in the observation point. + + Parameters + ---------- + easting : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate of the vertex of the prism. Must be in + meters. + radius : float + Square root of the sum of the squares of the ``easting``, ``northing`` + and ``upward`` shifted coordinates. + + Returns + ------- + kernel : float + Value of the kernel function for the easting-easting-northing component + of the grad tensor due to a rectangular prism evaluated on a single + vertex. + + Notes + ----- + Computes the following numerical kernel on the passed *shifted + coordinates*: + + .. math:: + + k_{xxy}(x, y, z) = + \frac{ + - x + }{ + (z + \sqrt{x^2 + y^2 + z^2}) + \sqrt{x^2 + y^2 + z^2} + } + + Examples + -------- + >>> x, y, z = 3.1, 5.2, -3.0 + >>> r = np.sqrt(x**2 + y**2 + z**2) + >>> float(kernel_een(x, y, z, r)) # doctest: +NUMBER + -0.12214070 + """ + return _kernel_iij(easting, northing, upward, radius) + + +@jit(nopython=True) +def kernel_eeu(easting, northing, upward, radius): + r""" + Easting-easting-upward kernel due to a prism. + + Evaluates the integration kernel for the easting-easting-upward + component of the grad tensor generated by a prism [Nagy2000]_ on a single + vertex of the prism. The coordinates that must be passed are shifted + coordinates: the coordinates of the vertex from a Cartesian coordinate + system whose origin is located in the observation point. + + Parameters + ---------- + easting : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate of the vertex of the prism. Must be in + meters. + radius : float + Square root of the sum of the squares of the ``easting``, ``northing`` + and ``upward`` shifted coordinates. + + Returns + ------- + kernel : float + Value of the kernel function for the easting-easting-upward component + of the grad tensor due to a rectangular prism evaluated on a single + vertex. + + Notes + ----- + Computes the following numerical kernel on the passed *shifted + coordinates*: + + .. math:: + + k_{xxz}(x, y, z) = + \frac{ + - x + }{ + (y + \sqrt{x^2 + y^2 + z^2}) + \sqrt{x^2 + y^2 + z^2} + } + + Examples + -------- + >>> x, y, z = 3.1, 5.2, -3.0 + >>> r = np.sqrt(x**2 + y**2 + z**2) + >>> float(kernel_eeu(x, y, z, r)) # doctest: +NUMBER + -0.03837408 + """ + return _kernel_iij(easting, upward, northing, radius) + + +@jit(nopython=True) +def kernel_enn(easting, northing, upward, radius): + r""" + Easting-northing-northing kernel due to a prism. + + Evaluates the integration kernel for the easting-northing-northing + component of the grad tensor generated by a prism [Nagy2000]_ on a single + vertex of the prism. The coordinates that must be passed are shifted + coordinates: the coordinates of the vertex from a Cartesian coordinate + system whose origin is located in the observation point. + + Parameters + ---------- + easting : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate of the vertex of the prism. Must be in + meters. + radius : float + Square root of the sum of the squares of the ``easting``, ``northing`` + and ``upward`` shifted coordinates. + + Returns + ------- + kernel : float + Value of the kernel function for the easting-northing-northing + component of the grad tensor due to a rectangular prism evaluated on + a single vertex. + + Notes + ----- + Computes the following numerical kernel on the passed *shifted + coordinates*: + + .. math:: + + k_{xyy}(x, y, z) = + \frac{ + - y + }{ + (z + \sqrt{x^2 + y^2 + z^2}) + \sqrt{x^2 + y^2 + z^2} + } + + Examples + -------- + >>> x, y, z = 3.1, 5.2, -3.0 + >>> r = np.sqrt(x**2 + y**2 + z**2) + >>> float(kernel_enn(x, y, z, r)) # doctest: +NUMBER + -0.20488118 + """ + return _kernel_iij(northing, easting, upward, radius) + + +@jit(nopython=True) +def kernel_nnu(easting, northing, upward, radius): + r""" + Northing-northing-upward kernel due to a prism. + + Evaluates the integration kernel for the northing-northing-upward + component of the grad tensor generated by a prism [Nagy2000]_ on a single + vertex of the prism. The coordinates that must be passed are shifted + coordinates: the coordinates of the vertex from a Cartesian coordinate + system whose origin is located in the observation point. + + Parameters + ---------- + easting : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate of the vertex of the prism. Must be in + meters. + radius : float + Square root of the sum of the squares of the ``easting``, ``northing`` + and ``upward`` shifted coordinates. + + Returns + ------- + kernel : float + Value of the kernel function for the northing-northing-upward + component of the grad tensor due to a rectangular prism evaluated on + a single vertex. + + Notes + ----- + Computes the following numerical kernel on the passed *shifted + coordinates*: + + .. math:: + + k_{yyz}(x, y, z) = + \frac{ + - y + }{ + (x + \sqrt{x^2 + y^2 + z^2}) + \sqrt{x^2 + y^2 + z^2} + } + + Examples + -------- + >>> x, y, z = 3.1, 5.2, -3.0 + >>> r = np.sqrt(x**2 + y**2 + z**2) + >>> float(kernel_nnu(x, y, z, r)) # doctest: +NUMBER + -0.07808384 + """ + return _kernel_iij(northing, upward, easting, radius) + + +@jit(nopython=True) +def kernel_euu(easting, northing, upward, radius): + r""" + Easting-upward-upward kernel due to a prism. + + Evaluates the integration kernel for the easting-upward-upward + component of the grad tensor generated by a prism [Nagy2000]_ on a single + vertex of the prism. The coordinates that must be passed are shifted + coordinates: the coordinates of the vertex from a Cartesian coordinate + system whose origin is located in the observation point. + + Parameters + ---------- + easting : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate of the vertex of the prism. Must be in + meters. + radius : float + Square root of the sum of the squares of the ``easting``, ``northing`` + and ``upward`` shifted coordinates. + + Returns + ------- + kernel : float + Value of the kernel function for the easting-upward-upward + component of the grad tensor due to a rectangular prism evaluated on + a single vertex. + + Notes + ----- + Computes the following numerical kernel on the passed *shifted + coordinates*: + + .. math:: + + k_{xzz}(x, y, z) = + \frac{ + - z + }{ + (y + \sqrt{x^2 + y^2 + z^2}) + \sqrt{x^2 + y^2 + z^2} + } + + Examples + -------- + >>> x, y, z = 3.1, 5.2, -3.0 + >>> r = np.sqrt(x**2 + y**2 + z**2) + >>> float(kernel_euu(x, y, z, r)) # doctest: +NUMBER + 0.03713621 + """ + return _kernel_iij(upward, easting, northing, radius) + + +@jit(nopython=True) +def kernel_nuu(easting, northing, upward, radius): + r""" + Northing-upward-upward kernel due to a prism. + + Evaluates the integration kernel for the northing-upward-upward + component of the grad tensor generated by a prism [Nagy2000]_ on a single + vertex of the prism. The coordinates that must be passed are shifted + coordinates: the coordinates of the vertex from a Cartesian coordinate + system whose origin is located in the observation point. + + Parameters + ---------- + easting : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate of the vertex of the prism. Must be in + meters. + radius : float + Square root of the sum of the squares of the ``easting``, ``northing`` + and ``upward`` shifted coordinates. + + Returns + ------- + kernel : float + Value of the kernel function for the northing-upward-upward + component of the grad tensor due to a rectangular prism evaluated on + a single vertex. + + Notes + ----- + Computes the following numerical kernel on the passed *shifted + coordinates*: + + .. math:: + + k_{yzz}(x, y, z) = + \frac{ + - z + }{ + (x + \sqrt{x^2 + y^2 + z^2}) + \sqrt{x^2 + y^2 + z^2} + } + + Examples + -------- + >>> x, y, z = 3.1, 5.2, -3.0 + >>> r = np.sqrt(x**2 + y**2 + z**2) + >>> float(kernel_nuu(x, y, z, r)) # doctest: +NUMBER + 0.04504837 + """ + return _kernel_iij(upward, northing, easting, radius) + + +@jit(nopython=True) +def kernel_enu(easting, northing, upward, radius): + r""" + Easting-northing-upward kernel due to a prism. + + Evaluates the integration kernel for the easting-northing-upward + component of the grad tensor generated by a prism [Nagy2000]_ on a single + vertex of the prism. The coordinates that must be passed are shifted + coordinates: the coordinates of the vertex from a Cartesian coordinate + system whose origin is located in the observation point. + + Parameters + ---------- + easting : float + Shifted easting coordinate of the vertex of the prism. Must be in + meters. + northing : float + Shifted northing coordinate of the vertex of the prism. Must be in + meters. + upward : float + Shifted upward coordinate of the vertex of the prism. Must be in + meters. + radius : float + Square root of the sum of the squares of the ``easting``, ``northing`` + and ``upward`` shifted coordinates. + + Returns + ------- + kernel : float + Value of the kernel function for the easting-northing-upward component + of the grad tensor due to a rectangular prism evaluated on a single + vertex. + + Notes + ----- + Computes the following numerical kernel on the passed *shifted + coordinates*: + + .. math:: + + k_{xyz}(x, y, z) = \frac{-1}{\sqrt{x^2 + y^2 + z^2}} + + Examples + -------- + >>> x, y, z = 3.1, 5.2, -3.0 + >>> r = np.sqrt(x**2 + y**2 + z**2) + >>> float(kernel_enu(x, y, z, r)) # doctest: +NUMBER + -0.1480061 + """ + return -1 / radius + + +@jit(nopython=True) +def _kernel_iii(x_i, x_j, x_k, radius): + r""" + Diagonal 3rd order kernels of a prism. + + Evaluates the following integration kernel: + + .. math:: + + k_{iii}(x_i, x_j, x_k) = + \frac{ + - x_j x_k (2 x_i^2 + x_j^2 + x_k^2) + }{ + (x_i^2 + x_j^2)(x_i^2 + x_k^2) + \sqrt{x_i^2 + x_j^2 + x_k^2} + } + + This function allows us to evaluate the third order kernels + :math:`k_{xxx}(x, y, z)`, :math:`k_{yyy}(x, y, z)`, and :math:`k_{zzz}(x, + y, z)` by cycling through the coordinates of the nodes of the prism. + + Parameters + ---------- + x_i, x_j, x_k : floats + Shifted coordinates of the vertex of the prism. Must be in meters. + radius : float + Square root of the sum of the squares of the ``x_i``, ``x_j``, and + ``x_k`` shifted coordinates. + + Returns + ------- + kernel : float + Value of the kernel function. + + Examples + -------- + Given the shifted coordinates of a prism vertex ``easting``, ``northing`` + and ``upward``, we can use this function to compute the diagonal third + order kernels by cycling the order of the shifted coordinates of the + vertex: + + >>> import numpy as np + >>> easting, northing, upward = 1.7, 2.8, 3.9 + >>> radius = np.sqrt(easting**2 + northing**2 + upward**2) + >>> k_eee = _kernel_iii(easting, northing, upward, radius) + >>> k_nnn = _kernel_iii(northing, upward, easting, radius) + >>> k_uuu = _kernel_iii(upward, easting, northing, radius) + """ + if (x_i == 0 and x_j == 0) or (x_i == 0 and x_k == 0): + return 0.0 + x_i_sq, x_j_sq, x_k_sq = x_i**2, x_j**2, x_k**2 + numerator = -x_j * x_k * (2 * x_i_sq + x_j_sq + x_k_sq) + denominator = (x_i_sq + x_j_sq) * (x_i_sq + x_k_sq) * radius + return numerator / denominator + + +@jit(nopython=True) +def _kernel_iij(x_i, x_j, x_k, radius): + r""" + Non-diagonal 3rd order kernels of a prism. + + Evaluates the following integration kernel: + + .. math:: + + easting / ((upward + radius) * radius) + + k_{iij}(x_i, x_j, x_k) = + \frac{ + - x_i + }{ + (x_k + \sqrt{x_i^2 + x_j^2 + x_k^2}) + \sqrt{x_i^2 + x_j^2 + x_k^2} + } + + This function allows us to evaluate the non-diagonal third order kernels + :math:`k_{xxy}(x, y, z)`, :math:`k_{xxz}(x, y, z)`, + :math:`k_{xyy}(x, y, z)`, :math:`k_{yyz}(x, y, z)`, + :math:`k_{xzz}(x, y, z)`, and :math:`k_{yzz}(x, y, z)` by cycling + through the coordinates of the nodes of the prism. + + Parameters + ---------- + x_i, x_j, x_k : floats + Shifted coordinates of the vertex of the prism. Must be in meters. + radius : float + Square root of the sum of the squares of the ``x_i``, ``x_j``, and + ``x_k`` shifted coordinates. + + Returns + ------- + kernel : float + Value of the kernel function. + + Examples + -------- + Given the shifted coordinates of a prism vertex ``easting``, ``northing`` + and ``upward``, we can use this function to compute the non-diagonal third + order kernels by changing the order of the shifted coordinates of the + vertex: + + >>> import numpy as np + >>> easting, northing, upward = 1.7, 2.8, 3.9 + >>> radius = np.sqrt(easting**2 + northing**2 + upward**2) + >>> k_een = _kernel_iij(easting, northing, upward, radius) + >>> k_eeu = _kernel_iij(easting, upward, northing, radius) + >>> k_enn = _kernel_iij(northing, easting, upward, radius) + >>> k_nnu = _kernel_iij(northing, upward, easting, radius) + >>> k_euu = _kernel_iij(upward, easting, northing, radius) + >>> k_nuu = _kernel_iij(upward, northing, upward, radius) + """ + if x_i == 0 and x_j == 0: + return 0.0 + return -x_i / ((x_k + radius) * radius) diff --git a/choclo/prism/_magnetic.py b/choclo/prism/_magnetic.py index 3d358d7a..d8ff39c9 100644 --- a/choclo/prism/_magnetic.py +++ b/choclo/prism/_magnetic.py @@ -11,7 +11,24 @@ from numba import jit from ..constants import VACUUM_MAGNETIC_PERMEABILITY -from ._kernels import kernel_ee, kernel_en, kernel_eu, kernel_nn, kernel_nu, kernel_uu +from ._kernels import ( + kernel_ee, + kernel_eee, + kernel_een, + kernel_eeu, + kernel_en, + kernel_enn, + kernel_enu, + kernel_eu, + kernel_euu, + kernel_nn, + kernel_nnn, + kernel_nnu, + kernel_nu, + kernel_nuu, + kernel_uu, + kernel_uuu, +) from ._utils import ( is_interior_point, is_point_on_east_face, @@ -198,8 +215,7 @@ def magnetic_field( .. math:: u_{ij} = - \frac{\partial}{\partial i} - \frac{\partial}{\partial j} + \frac{\partial^2}{\partial i \partial j} \int\limits_R \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} dv @@ -226,6 +242,12 @@ def magnetic_field( :func:`choclo.prism.magnetic_e` :func:`choclo.prism.magnetic_n` :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` """ # Check if observation point falls in a singular point if is_point_on_edge( @@ -428,8 +450,7 @@ def magnetic_e( .. math:: u_{ij} = - \frac{\partial}{\partial i} - \frac{\partial}{\partial j} + \frac{\partial^2}{\partial i \partial j} \int\limits_R \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} dv @@ -474,6 +495,12 @@ def magnetic_e( :func:`choclo.prism.magnetic_field` :func:`choclo.prism.magnetic_n` :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` """ # Check if observation point falls in a singular point if is_point_on_edge( @@ -498,43 +525,24 @@ def magnetic_e( prism_top, ): return np.nan - # Initialize magnetic field vector component - b_e = 0.0 - # Iterate over the vertices of the prism - for i in range(2): - # Compute shifted easting coordinate - if i == 0: - shift_east = prism_east - easting - else: - shift_east = prism_west - easting - shift_east_sq = shift_east**2 - for j in range(2): - # Compute shifted northing coordinate - if j == 0: - shift_north = prism_north - northing - else: - shift_north = prism_south - northing - shift_north_sq = shift_north**2 - for k in range(2): - # Compute shifted upward coordinate - if k == 0: - shift_upward = prism_top - upward - else: - shift_upward = prism_bottom - upward - shift_upward_sq = shift_upward**2 - # Compute the radius - radius = np.sqrt(shift_east_sq + shift_north_sq + shift_upward_sq) - # Compute kernel tensor components for the current vertex - k_ee = kernel_ee(shift_east, shift_north, shift_upward, radius) - k_en = kernel_en(shift_east, shift_north, shift_upward, radius) - k_eu = kernel_eu(shift_east, shift_north, shift_upward, radius) - # Compute b_e using the dot product between the kernel tensor - # and the magnetization vector of the prism - b_e += (-1) ** (i + j + k) * ( - magnetization_east * k_ee - + magnetization_north * k_en - + magnetization_up * k_eu - ) + # Compute magnetic field vector component + b_e = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_ee, + kernel_en, + kernel_eu, + ) # Add 4 pi to Be if computing on the eastmost face, to correctly evaluate # the limit approaching from outside (approaching from the east) if is_point_on_east_face( @@ -629,8 +637,7 @@ def magnetic_n( .. math:: u_{ij} = - \frac{\partial}{\partial i} - \frac{\partial}{\partial j} + \frac{\partial^2}{\partial i \partial j} \int\limits_R \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} dv @@ -675,6 +682,12 @@ def magnetic_n( :func:`choclo.prism.magnetic_field` :func:`choclo.prism.magnetic_e` :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` """ # Check if observation point falls in a singular point if is_point_on_edge( @@ -699,43 +712,24 @@ def magnetic_n( prism_top, ): return np.nan - # Initialize magnetic field vector component - b_n = 0.0 - # Iterate over the vertices of the prism - for i in range(2): - # Compute shifted easting coordinate - if i == 0: - shift_east = prism_east - easting - else: - shift_east = prism_west - easting - shift_east_sq = shift_east**2 - for j in range(2): - # Compute shifted northing coordinate - if j == 0: - shift_north = prism_north - northing - else: - shift_north = prism_south - northing - shift_north_sq = shift_north**2 - for k in range(2): - # Compute shifted upward coordinate - if k == 0: - shift_upward = prism_top - upward - else: - shift_upward = prism_bottom - upward - shift_upward_sq = shift_upward**2 - # Compute the radius - radius = np.sqrt(shift_east_sq + shift_north_sq + shift_upward_sq) - # Compute kernel tensor components for the current vertex - k_nn = kernel_nn(shift_east, shift_north, shift_upward, radius) - k_en = kernel_en(shift_east, shift_north, shift_upward, radius) - k_nu = kernel_nu(shift_east, shift_north, shift_upward, radius) - # Compute b_n using the dot product between the kernel tensor - # and the magnetization vector of the prism - b_n += (-1) ** (i + j + k) * ( - magnetization_east * k_en - + magnetization_north * k_nn - + magnetization_up * k_nu - ) + # Compute magnetic field vector component + b_n = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_en, + kernel_nn, + kernel_nu, + ) # Add 4 pi to Bn if computing on the northmost face, to correctly evaluate # the limit approaching from outside (approaching from the north) if is_point_on_north_face( @@ -830,8 +824,7 @@ def magnetic_u( .. math:: u_{ij} = - \frac{\partial}{\partial i} - \frac{\partial}{\partial j} + \frac{\partial^2}{\partial i \partial j} \int\limits_R \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} dv @@ -876,6 +869,12 @@ def magnetic_u( :func:`choclo.prism.magnetic_field` :func:`choclo.prism.magnetic_e` :func:`choclo.prism.magnetic_n` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` """ # Check if observation point falls in a singular point if is_point_on_edge( @@ -900,43 +899,24 @@ def magnetic_u( prism_top, ): return np.nan - # Initialize magnetic field vector component - b_u = 0.0 - # Iterate over the vertices of the prism - for i in range(2): - # Compute shifted easting coordinate - if i == 0: - shift_east = prism_east - easting - else: - shift_east = prism_west - easting - shift_east_sq = shift_east**2 - for j in range(2): - # Compute shifted northing coordinate - if j == 0: - shift_north = prism_north - northing - else: - shift_north = prism_south - northing - shift_north_sq = shift_north**2 - for k in range(2): - # Compute shifted upward coordinate - if k == 0: - shift_upward = prism_top - upward - else: - shift_upward = prism_bottom - upward - shift_upward_sq = shift_upward**2 - # Compute the radius - radius = np.sqrt(shift_east_sq + shift_north_sq + shift_upward_sq) - # Compute kernel tensor components for the current vertex - k_uu = kernel_uu(shift_east, shift_north, shift_upward, radius) - k_eu = kernel_eu(shift_east, shift_north, shift_upward, radius) - k_nu = kernel_nu(shift_east, shift_north, shift_upward, radius) - # Compute b_n using the dot product between the kernel tensor - # and the magnetization vector of the prism - b_u += (-1) ** (i + j + k) * ( - magnetization_east * k_eu - + magnetization_north * k_nu - + magnetization_up * k_uu - ) + # Compute magnetic field vector component + b_u = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_eu, + kernel_nu, + kernel_uu, + ) # Add 4 pi to Bu if computing on the north face, to correctly evaluate the # limit approaching from outside (approaching from the top) if is_point_on_top_face( @@ -952,3 +932,1017 @@ def magnetic_u( ): b_u += 4 * np.pi * magnetization_up return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_u + + +@jit(nopython=True) +def magnetic_ee( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, +): + r""" + Easting derivative of the easting component of the magnetic field. + + Returns the easting derivative of the easting component of the magnetic + field due to a single rectangular prism on a single computation point. + + Parameters + ---------- + easting : float + Easting coordinate of the observation point. Must be in meters. + northing : float + Northing coordinate of the observation point. Must be in meters. + upward : float + Upward coordinate of the observation point. Must be in meters. + prism_west : float + The West boundary of the prism. Must be in meters. + prism_east : float + The East boundary of the prism. Must be in meters. + prism_south : float + The South boundary of the prism. Must be in meters. + prism_north : float + The North boundary of the prism. Must be in meters. + prism_bottom : float + The bottom boundary of the prism. Must be in meters. + prism_top : float + The top boundary of the prism. Must be in meters. + magnetization_east : float + The East component of the magnetization vector of the prism. Must be in + :math:`A m^{-1}`. + magnetization_north : float + The North component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + magnetization_up : float + The upward component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + + Returns + ------- + b_ee : float + Easting derivative of the easting component of the magnetic field + generated by the prism on the observation point in :math:`\text{T}`. + Return ``numpy.nan`` if the observation point falls in + a singular point: prism vertices, prism edges or interior points. + + Notes + ----- + Computes the easting derivative of the easting component of the magnetic + field :math:`\mathbf{B}(\mathbf{p})` generated by a rectangular prism + :math:`R` with a magnetization vector :math:`M` on the observation point + :math:`\mathbf{p}` as follows: + + .. math:: + + B_{xx}(\mathbf{p}) = + \frac{\mu_0}{4\pi} + \left( M_x u_{xxx} + M_y u_{xxy} + M_z u_{xxz} \right) + + Where :math:`u_{ijk}` are: + + .. math:: + + u_{ijk} = + \frac{\partial^3}{\partial i \partial j \partial k} + \int\limits_R + \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} + dv + + with :math:`i,j,k \in \{x, y, z\}`. + + References + ---------- + - [Blakely1995]_ + - [Oliveira2015]_ + - [Nagy2000]_ + - [Nagy2002]_ + - [Fukushima2020]_ + + See Also + -------- + :func:`choclo.prism.magnetic_field` + :func:`choclo.prism.magnetic_e` + :func:`choclo.prism.magnetic_n` + :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` + """ + # Check if observation point falls in a singular point + if is_point_on_edge( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ) or is_interior_point( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ): + return np.nan + # Compute magnetic gradiometry component + b_ee = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_eee, + kernel_een, + kernel_eeu, + ) + return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_ee + + +@jit(nopython=True) +def magnetic_nn( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, +): + r""" + Northing derivative of the northing component of the magnetic field. + + Returns the northing derivative of the northing component of the magnetic + field due to a single rectangular prism on a single computation point. + + Parameters + ---------- + easting : float + Easting coordinate of the observation point. Must be in meters. + northing : float + Northing coordinate of the observation point. Must be in meters. + upward : float + Upward coordinate of the observation point. Must be in meters. + prism_west : float + The West boundary of the prism. Must be in meters. + prism_east : float + The East boundary of the prism. Must be in meters. + prism_south : float + The South boundary of the prism. Must be in meters. + prism_north : float + The North boundary of the prism. Must be in meters. + prism_bottom : float + The bottom boundary of the prism. Must be in meters. + prism_top : float + The top boundary of the prism. Must be in meters. + magnetization_east : float + The East component of the magnetization vector of the prism. Must be in + :math:`A m^{-1}`. + magnetization_north : float + The North component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + magnetization_up : float + The upward component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + + Returns + ------- + b_nn : float + Northing derivative of the northing component of the magnetic field + generated by the prism on the observation point in :math:`\text{T}`. + Return ``numpy.nan`` if the observation point falls in + a singular point: prism vertices, prism edges or interior points. + + Notes + ----- + Computes the northing derivative of the northing component of the magnetic + field :math:`\mathbf{B}(\mathbf{p})` generated by a rectangular prism + :math:`R` with a magnetization vector :math:`M` on the observation point + :math:`\mathbf{p}` as follows: + + .. math:: + + B_{yy}(\mathbf{p}) = + \frac{\mu_0}{4\pi} + \left( M_x u_{xyy} + M_y u_{yyy} + M_z u_{yyz} \right) + + Where :math:`u_{ijk}` are: + + .. math:: + + u_{ijk} = + \frac{\partial^3}{\partial i \partial j \partial k} + \int\limits_R + \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} + dv + + with :math:`i,j,k \in \{x, y, z\}`. + + References + ---------- + - [Blakely1995]_ + - [Oliveira2015]_ + - [Nagy2000]_ + - [Nagy2002]_ + - [Fukushima2020]_ + + See Also + -------- + :func:`choclo.prism.magnetic_field` + :func:`choclo.prism.magnetic_e` + :func:`choclo.prism.magnetic_n` + :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` + """ + # Check if observation point falls in a singular point + if is_point_on_edge( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ) or is_interior_point( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ): + return np.nan + # Compute magnetic gradiometry component + b_nn = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_enn, + kernel_nnn, + kernel_nnu, + ) + return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_nn + + +@jit(nopython=True) +def magnetic_uu( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, +): + r""" + Upward derivative of the upward component of the magnetic field. + + Returns the upward derivative of the upward component of the magnetic + field due to a single rectangular prism on a single computation point. + + Parameters + ---------- + easting : float + Easting coordinate of the observation point. Must be in meters. + northing : float + Northing coordinate of the observation point. Must be in meters. + upward : float + Upward coordinate of the observation point. Must be in meters. + prism_west : float + The West boundary of the prism. Must be in meters. + prism_east : float + The East boundary of the prism. Must be in meters. + prism_south : float + The South boundary of the prism. Must be in meters. + prism_north : float + The North boundary of the prism. Must be in meters. + prism_bottom : float + The bottom boundary of the prism. Must be in meters. + prism_top : float + The top boundary of the prism. Must be in meters. + magnetization_east : float + The East component of the magnetization vector of the prism. Must be in + :math:`A m^{-1}`. + magnetization_north : float + The North component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + magnetization_up : float + The upward component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + + Returns + ------- + b_uu : float + Upward derivative of the upward component of the magnetic field + generated by the prism on the observation point in :math:`\text{T}`. + Return ``numpy.nan`` if the observation point falls in + a singular point: prism vertices, prism edges or interior points. + + Notes + ----- + Computes the upward derivative of the upward component of the magnetic + field :math:`\mathbf{B}(\mathbf{p})` generated by a rectangular prism + :math:`R` with a magnetization vector :math:`M` on the observation point + :math:`\mathbf{p}` as follows: + + .. math:: + + B_{zz}(\mathbf{p}) = + \frac{\mu_0}{4\pi} + \left( M_x u_{xzz} + M_y u_{yzz} + M_z u_{zzz} \right) + + Where :math:`u_{ijk}` are: + + .. math:: + + u_{ijk} = + \frac{\partial^3}{\partial i \partial j \partial k} + \int\limits_R + \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} + dv + + with :math:`i,j,k \in \{x, y, z\}`. + + References + ---------- + - [Blakely1995]_ + - [Oliveira2015]_ + - [Nagy2000]_ + - [Nagy2002]_ + - [Fukushima2020]_ + + See Also + -------- + :func:`choclo.prism.magnetic_field` + :func:`choclo.prism.magnetic_e` + :func:`choclo.prism.magnetic_n` + :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` + """ + # Check if observation point falls in a singular point + if is_point_on_edge( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ) or is_interior_point( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ): + return np.nan + # Compute magnetic gradiometry component + b_uu = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_euu, + kernel_nuu, + kernel_uuu, + ) + return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_uu + + +@jit(nopython=True) +def magnetic_en( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, +): + r""" + Northing derivative of the easting component of the magnetic field. + + Returns the northing derivative of the easting component of the magnetic + field due to a single rectangular prism on a single computation point. + + Parameters + ---------- + easting : float + Easting coordinate of the observation point. Must be in meters. + northing : float + Northing coordinate of the observation point. Must be in meters. + upward : float + Upward coordinate of the observation point. Must be in meters. + prism_west : float + The West boundary of the prism. Must be in meters. + prism_east : float + The East boundary of the prism. Must be in meters. + prism_south : float + The South boundary of the prism. Must be in meters. + prism_north : float + The North boundary of the prism. Must be in meters. + prism_bottom : float + The bottom boundary of the prism. Must be in meters. + prism_top : float + The top boundary of the prism. Must be in meters. + magnetization_east : float + The East component of the magnetization vector of the prism. Must be in + :math:`A m^{-1}`. + magnetization_north : float + The North component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + magnetization_up : float + The upward component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + + Returns + ------- + b_en : float + Northing derivative of the easting component of the magnetic field + generated by the prism on the observation point in :math:`\text{T}`. + Return ``numpy.nan`` if the observation point falls in a singular + point: prism vertices, prism edges or interior points. + + Notes + ----- + Computes the northing derivative of the easting component of the magnetic + field :math:`\mathbf{B}(\mathbf{p})` generated by a rectangular prism + :math:`R` with a magnetization vector :math:`M` on the observation point + :math:`\mathbf{p}` as follows: + + .. math:: + + B_{xy}(\mathbf{p}) = + \frac{\mu_0}{4\pi} + \left( M_x u_{xxy} + M_y u_{xyy} + M_z u_{xyz} \right) + + Where :math:`u_{ijk}` are: + + .. math:: + + u_{ijk} = + \frac{\partial^3}{\partial i \partial j \partial k} + \int\limits_R + \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} + dv + + with :math:`i,j,k \in \{x, y, z\}`. + + References + ---------- + - [Blakely1995]_ + - [Oliveira2015]_ + - [Nagy2000]_ + - [Nagy2002]_ + - [Fukushima2020]_ + + See Also + -------- + :func:`choclo.prism.magnetic_field` + :func:`choclo.prism.magnetic_e` + :func:`choclo.prism.magnetic_n` + :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_eu` + :func:`choclo.prism.magnetic_nu` + """ + # Check if observation point falls in a singular point + if is_point_on_edge( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ) or is_interior_point( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ): + return np.nan + # Compute magnetic gradiometry component + b_en = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_een, + kernel_enn, + kernel_enu, + ) + return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_en + + +@jit(nopython=True) +def magnetic_eu( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, +): + r""" + Upward derivative of the easting component of the magnetic field. + + Returns the upward derivative of the easting component of the magnetic + field due to a single rectangular prism on a single computation point. + + Parameters + ---------- + easting : float + Easting coordinate of the observation point. Must be in meters. + northing : float + Northing coordinate of the observation point. Must be in meters. + upward : float + Upward coordinate of the observation point. Must be in meters. + prism_west : float + The West boundary of the prism. Must be in meters. + prism_east : float + The East boundary of the prism. Must be in meters. + prism_south : float + The South boundary of the prism. Must be in meters. + prism_north : float + The North boundary of the prism. Must be in meters. + prism_bottom : float + The bottom boundary of the prism. Must be in meters. + prism_top : float + The top boundary of the prism. Must be in meters. + magnetization_east : float + The East component of the magnetization vector of the prism. Must be in + :math:`A m^{-1}`. + magnetization_north : float + The North component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + magnetization_up : float + The upward component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + + Returns + ------- + b_eu : float + Upward derivative of the easting component of the magnetic field + generated by the prism on the observation point in :math:`\text{T}`. + Return ``numpy.nan`` if the observation point falls in a singular + point: prism vertices, prism edges or interior points. + + Notes + ----- + Computes the northing derivative of the easting component of the magnetic + field :math:`\mathbf{B}(\mathbf{p})` generated by a rectangular prism + :math:`R` with a magnetization vector :math:`M` on the observation point + :math:`\mathbf{p}` as follows: + + .. math:: + + B_{xz}(\mathbf{p}) = + \frac{\mu_0}{4\pi} + \left( M_x u_{xxz} + M_y u_{xyz} + M_z u_{xzz} \right) + + Where :math:`u_{ijk}` are: + + .. math:: + + u_{ijk} = + \frac{\partial^3}{\partial i \partial j \partial k} + \int\limits_R + \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} + dv + + with :math:`i,j,k \in \{x, y, z\}`. + + References + ---------- + - [Blakely1995]_ + - [Oliveira2015]_ + - [Nagy2000]_ + - [Nagy2002]_ + - [Fukushima2020]_ + + See Also + -------- + :func:`choclo.prism.magnetic_field` + :func:`choclo.prism.magnetic_e` + :func:`choclo.prism.magnetic_n` + :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_nu` + """ + # Check if observation point falls in a singular point + if is_point_on_edge( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ) or is_interior_point( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ): + return np.nan + # Compute magnetic gradiometry component + b_eu = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_eeu, + kernel_enu, + kernel_euu, + ) + return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_eu + + +@jit(nopython=True) +def magnetic_nu( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, +): + r""" + Upward derivative of the northing component of the magnetic field. + + Returns the upward derivative of the northing component of the magnetic + field due to a single rectangular prism on a single computation point. + + Parameters + ---------- + easting : float + Easting coordinate of the observation point. Must be in meters. + northing : float + Northing coordinate of the observation point. Must be in meters. + upward : float + Upward coordinate of the observation point. Must be in meters. + prism_west : float + The West boundary of the prism. Must be in meters. + prism_east : float + The East boundary of the prism. Must be in meters. + prism_south : float + The South boundary of the prism. Must be in meters. + prism_north : float + The North boundary of the prism. Must be in meters. + prism_bottom : float + The bottom boundary of the prism. Must be in meters. + prism_top : float + The top boundary of the prism. Must be in meters. + magnetization_east : float + The East component of the magnetization vector of the prism. Must be in + :math:`A m^{-1}`. + magnetization_north : float + The North component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + magnetization_up : float + The upward component of the magnetization vector of the prism. Must be + in :math:`A m^{-1}`. + + Returns + ------- + b_nu : float + Upward derivative of the northing component of the magnetic field + generated by the prism on the observation point in :math:`\text{T}`. + Return ``numpy.nan`` if the observation point falls in a singular + point: prism vertices, prism edges or interior points. + + Notes + ----- + Computes the northing derivative of the easting component of the magnetic + field :math:`\mathbf{B}(\mathbf{p})` generated by a rectangular prism + :math:`R` with a magnetization vector :math:`M` on the observation point + :math:`\mathbf{p}` as follows: + + .. math:: + + B_{yz}(\mathbf{p}) = + \frac{\mu_0}{4\pi} + \left( M_x u_{xyz} + M_y u_{yyz} + M_z u_{yzz} \right) + + Where :math:`u_{ijk}` are: + + .. math:: + + u_{ijk} = + \frac{\partial^3}{\partial i \partial j \partial k} + \int\limits_R + \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} + dv + + with :math:`i,j,k \in \{x, y, z\}`. + + References + ---------- + - [Blakely1995]_ + - [Oliveira2015]_ + - [Nagy2000]_ + - [Nagy2002]_ + - [Fukushima2020]_ + + See Also + -------- + :func:`choclo.prism.magnetic_field` + :func:`choclo.prism.magnetic_e` + :func:`choclo.prism.magnetic_n` + :func:`choclo.prism.magnetic_u` + :func:`choclo.prism.magnetic_ee` + :func:`choclo.prism.magnetic_nn` + :func:`choclo.prism.magnetic_uu` + :func:`choclo.prism.magnetic_en` + :func:`choclo.prism.magnetic_eu` + """ + # Check if observation point falls in a singular point + if is_point_on_edge( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ) or is_interior_point( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + ): + return np.nan + # Compute magnetic gradiometry component + b_nu = _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_enu, + kernel_nnu, + kernel_nuu, + ) + return VACUUM_MAGNETIC_PERMEABILITY / 4 / np.pi * b_nu + + +@jit(nopython=True) +def _calculate_component( + easting, + northing, + upward, + prism_west, + prism_east, + prism_south, + prism_north, + prism_bottom, + prism_top, + magnetization_east, + magnetization_north, + magnetization_up, + kernel_i, + kernel_j, + kernel_k, +): + r""" + Calculate field component for a single prism and observation point. + + Evaluate the provided kernels on the shifted coordinates of prism vertices + to compute the magnetic field component given a magnetization vector of the + prism. + + Parameters + ---------- + easting, northing, upward : floats + Coordinates of the observation point. Must be in meters. + prism_west, prism_east : floats + The easting boundaries of the prism. Must be in meters. + prism_south, prism_north : floats + The northing boundaries of the prism. Must be in meters. + prism_bottom, prism_top : floats + The upward boundaries of the prism. Must be in meters. + magnetization_east, magnetization_north, magnetization_up : floats + The components of the magnetization vector of the prism. Must be in + :math:`A m^{-1}`. + kernel_i, kernel_j, kernel_k : callables + Kernel functions to be evaluated on each vertex of the prism. + + Returns + ------- + float + + Notes + ----- + Given the kernels :math:`k_i(\hat{x}, \hat{y}, \hat{z})`, + :math:`k_j(\hat{x}, \hat{y}, \hat{z})`, and :math:`k_k(\hat{x}, \hat{y}, + \hat{z})`; a prism provided by its boundaries :math:`x_1`, :math:`x_2`, + :math:`y_1`, :math:`y_2`, :math:`z_1`, and :math:`z_2`; a magnetization + vector :math:`\mathbf{M}=(M_x, M_y, M_z)`; and an observation point + :math:`\mathbf{p}=(x, y, z)`, this function returns: + + .. math:: + + M_x u_x(x, y, z) + M_y u_y(x, y, z) + M_z u_z(x, y, z), + + where + + .. math:: + + u_x(x, y, z) = + \Bigg\lvert\Bigg\lvert\Bigg\lvert + k_i(\hat{x}, \hat{y}, \hat{z}) + \Bigg\rvert_{x_1 - x}^{x_2 - x} + \Bigg\rvert_{y_1 - y}^{y_2 - y} + \Bigg\rvert_{z_1 - z}^{z_2 - z} + + .. math:: + + u_y(x, y, z) = + \Bigg\lvert\Bigg\lvert\Bigg\lvert + k_j(\hat{x}, \hat{y}, \hat{z}) + \Bigg\rvert_{x_1 - x}^{x_2 - x} + \Bigg\rvert_{y_1 - y}^{y_2 - y} + \Bigg\rvert_{z_1 - z}^{z_2 - z} + + .. math:: + + u_z(x, y, z) = + \Bigg\lvert\Bigg\lvert\Bigg\lvert + k_k(\hat{x}, \hat{y}, \hat{z}) + \Bigg\rvert_{x_1 - x}^{x_2 - x} + \Bigg\rvert_{y_1 - y}^{y_2 - y} + \Bigg\rvert_{z_1 - z}^{z_2 - z} + """ + result = 0.0 + # Iterate over the vertices of the prism + for i in range(2): + # Compute shifted easting coordinate + if i == 0: + shift_east = prism_east - easting + else: + shift_east = prism_west - easting + shift_east_sq = shift_east**2 + for j in range(2): + # Compute shifted northing coordinate + if j == 0: + shift_north = prism_north - northing + else: + shift_north = prism_south - northing + shift_north_sq = shift_north**2 + for k in range(2): + # Compute shifted upward coordinate + if k == 0: + shift_upward = prism_top - upward + else: + shift_upward = prism_bottom - upward + shift_upward_sq = shift_upward**2 + # Compute the radius + radius = np.sqrt(shift_east_sq + shift_north_sq + shift_upward_sq) + # Compute kernel tensor components for the current vertex + k_e = kernel_i(shift_east, shift_north, shift_upward, radius) + k_n = kernel_j(shift_east, shift_north, shift_upward, radius) + k_u = kernel_k(shift_east, shift_north, shift_upward, radius) + # Compute b_en using the dot product between the kernel tensor + # and the magnetization vector of the prism + result += (-1) ** (i + j + k) * ( + magnetization_east * k_e + + magnetization_north * k_n + + magnetization_up * k_u + ) + return result diff --git a/choclo/tests/test_prism_magnetic.py b/choclo/tests/test_prism_magnetic.py index 224cb9c1..0b46c3d6 100644 --- a/choclo/tests/test_prism_magnetic.py +++ b/choclo/tests/test_prism_magnetic.py @@ -7,6 +7,8 @@ """ Test magnetic forward modelling functions for rectangular prisms """ +import itertools + import numpy as np import numpy.testing as npt import pytest @@ -16,11 +18,35 @@ kernel_eu, kernel_nu, magnetic_e, + magnetic_ee, + magnetic_en, + magnetic_eu, magnetic_field, magnetic_n, + magnetic_nn, + magnetic_nu, magnetic_u, + magnetic_uu, ) +COMPONENTS = { + "e": magnetic_e, + "n": magnetic_n, + "u": magnetic_u, + "ee": magnetic_ee, + "nn": magnetic_nn, + "uu": magnetic_uu, + "en": magnetic_en, + "eu": magnetic_eu, + "nu": magnetic_nu, +} + + +def get_forward_function(component: str): + """Return desired magnetic forward function.""" + component = "".join(sorted(component)) + return COMPONENTS[component] + @pytest.fixture(name="sample_prism") def fixture_sample_prism(): @@ -599,6 +625,21 @@ class TestDivergenceOfB: Test if the divergence of the magnetic field is equal to zero. """ + def test_divergence_of_b(self, sample_3d_grid, sample_prism, sample_magnetization): + """ + Test div of B through magnetic gradiometry analytical functions. + """ + kwargs = dict( + coordinates=sample_3d_grid, + prism=sample_prism, + magnetization=sample_magnetization, + ) + b_ee = evaluate(magnetic_ee, **kwargs) + b_nn = evaluate(magnetic_nn, **kwargs) + b_uu = evaluate(magnetic_uu, **kwargs) + # Check if the divergence of B is zero + npt.assert_allclose(-b_uu, b_ee + b_nn, atol=1e-11) + def test_divergence_of_b_finite_differences( self, sample_3d_grid, sample_prism, sample_magnetization ): @@ -635,6 +676,19 @@ class TestMagneticFieldSingularities: we approach from outside of the prism. """ + COMPONENTS = ( + magnetic_field, + magnetic_e, + magnetic_n, + magnetic_u, + magnetic_ee, + magnetic_en, + magnetic_eu, + magnetic_nn, + magnetic_nu, + magnetic_uu, + ) + @pytest.fixture() def sample_prism(self): """ @@ -711,9 +765,7 @@ def get_interior_points(self, prism): coordinates = tuple(c.ravel() for c in np.meshgrid(easting, northing, upward)) return coordinates - @pytest.mark.parametrize( - "forward_func", (magnetic_field, magnetic_e, magnetic_n, magnetic_u) - ) + @pytest.mark.parametrize("forward_func", COMPONENTS) def test_on_vertices(self, sample_prism, forward_func): """ Test if magnetic field components on vertices are equal to NaN @@ -726,9 +778,7 @@ def test_on_vertices(self, sample_prism, forward_func): ) assert np.isnan(results).all() - @pytest.mark.parametrize( - "forward_func", (magnetic_field, magnetic_e, magnetic_n, magnetic_u) - ) + @pytest.mark.parametrize("forward_func", COMPONENTS) @pytest.mark.parametrize("direction", ("easting", "northing", "upward")) def test_on_edges(self, sample_prism, direction, forward_func): """ @@ -744,9 +794,7 @@ def test_on_edges(self, sample_prism, direction, forward_func): ) assert np.isnan(results).all() - @pytest.mark.parametrize( - "forward_func", (magnetic_field, magnetic_e, magnetic_n, magnetic_u) - ) + @pytest.mark.parametrize("forward_func", COMPONENTS) def test_on_interior_points(self, sample_prism, forward_func): """ Test if magnetic field components are NaN on internal points @@ -787,6 +835,154 @@ def test_symmetry_on_faces(self, sample_prism, direction, forward_func): ) npt.assert_allclose(results, results[0]) + @pytest.mark.parametrize( + "forward_func", + ( + magnetic_ee, + magnetic_nn, + magnetic_uu, + magnetic_en, + magnetic_eu, + magnetic_nu, + ), + ) + @pytest.mark.parametrize("direction", ("easting", "northing", "upward")) + def test_gradiometry_symmetry_on_faces(self, sample_prism, direction, forward_func): + """ + Tests symmetry of magnetic gradiometry components on the center of + faces normal to the component direction + + For example, check if ``magnetic_ee`` has the opposite value on points + in the top face and points of the bottom face. + """ + easting, northing, upward = self.get_faces_centers(sample_prism, direction) + magnetization = np.array([1.0, 1.0, 1.0]) + results = list( + forward_func(e, n, u, *sample_prism, *magnetization) + for (e, n, u) in zip(easting, northing, upward) + ) + npt.assert_allclose(-results[0], results[1:], atol=1e-23) + + +class TestMagGradiometryFiniteDifferences: + """ + Test the magnetic gradiometry components by comparing them to numerical + derivatives computed through finite differences of the magnetic components. + """ + + delta = 1e-6 # displacement used in the finite difference calculations + rtol, atol = 5e-4, 5e-12 # tolerances used in the comparisons + + @pytest.mark.parametrize("i", ["e", "n", "u"]) + @pytest.mark.parametrize("j", ["e", "n", "u"]) + def test_derivatives_of_magnetic_components( + self, sample_3d_grid, sample_prism, sample_magnetization, i, j + ): + """ + Compare the magnetic gradiometry functions with a finite difference + approximation of them computed from the magnetic field components. + + With this function we can compare the result of ``magnetic_{i}{j}`` + against a finite difference approximation using the ``magnetic_{i}`` + forward function, where ``i`` and ``j`` can be ``e``, ``n`` or ``u``. + """ + # Get forward functions + mag_func = get_forward_function(i) + mag_grad_func = get_forward_function(i + j) + direction = j # direction along which to apply the finite differences + # Evaluate the mag grad function on the sample grid + mag_grad = evaluate( + mag_grad_func, sample_3d_grid, sample_prism, sample_magnetization + ) + # Compute the mag grad function using finite differences + mag_grad_finite_diff = finite_differences( + sample_3d_grid, + sample_prism, + sample_magnetization, + direction=direction, + forward_func=mag_func, + delta=self.delta, + ) + # Compare results + npt.assert_allclose( + mag_grad, mag_grad_finite_diff, rtol=self.rtol, atol=self.atol + ) + + +class TestMagGradiometryLimits: + """ + Test if limits are well evaluated in magnetic gradiometry kernels. + + Most of the third-order kernels have singular points along the lines that + matches with the edges of the prism, but the magnetic gradiometry functions + are defined (and finite) in those regions. These tests ensure that these + limits are well resolved. + """ + + delta = 1e-6 # displacement used in the finite difference calculations + rtol, atol = 5e-4, 5e-12 # tolerances used in the comparisons + + @pytest.fixture + def coordinates(self, sample_prism): + """ + Return a set of coordinates located along the lines of the prism edges. + """ + west, east, south, north, bottom, top = sample_prism + distance = 13.4 # define distance from the nearest vertex for obs points + points = [] + # Define points on lines parallel to upward edges + for easting, northing in itertools.product((west, east), (south, north)): + points.append([easting, northing, bottom - distance]) + points.append([easting, northing, top + distance]) + # Define points on lines parallel to easting edges + for northing, upward in itertools.product((south, north), (bottom, top)): + points.append([west - distance, northing, upward]) + points.append([east + distance, northing, upward]) + # Define points on lines parallel to northing edges + for easting, upward in itertools.product((west, east), (bottom, top)): + points.append([easting, south - distance, upward]) + points.append([easting, north + distance, upward]) + points = np.array(points).T + easting, northing, upward = points[:] + return easting, northing, upward + + @pytest.mark.parametrize("i", ["e", "n", "u"]) + @pytest.mark.parametrize("j", ["e", "n", "u"]) + def test_derivatives_of_magnetic_components_on_edges( + self, coordinates, sample_prism, sample_magnetization, i, j + ): + """ + Compare the magnetic gradiometry functions with a finite difference + approximation of it computed from the magnetic field components on + observation points that are located along the same lines of the prism + edges. + + With this function we can compare the result of ``magnetic_{i}{j}`` + against a finite difference approximation using the ``magnetic_{i}`` + forward function, where ``i`` and ``j`` can be ``e``, ``n`` or ``u``. + """ + # Get forward functions + mag_func = get_forward_function(i) + mag_grad_func = get_forward_function(f"{i}{j}") + direction = j # direction along which to apply the finite differences + # Evaluate the mag grad function on the sample grid + mag_grad = evaluate( + mag_grad_func, coordinates, sample_prism, sample_magnetization + ) + # Compute the mag grad function using finite differences + mag_grad_finite_diff = finite_differences( + coordinates, + sample_prism, + sample_magnetization, + direction=direction, + forward_func=mag_func, + delta=self.delta, + ) + # Compare results + npt.assert_allclose( + mag_grad, mag_grad_finite_diff, rtol=self.rtol, atol=self.atol + ) + class TestBugfixKernelEvaluation: r""" diff --git a/doc/api/index.rst b/doc/api/index.rst index c7ae9bbf..6107e2bd 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -95,6 +95,12 @@ Magnetic prism.magnetic_e prism.magnetic_n prism.magnetic_u + prism.magnetic_ee + prism.magnetic_nn + prism.magnetic_uu + prism.magnetic_en + prism.magnetic_eu + prism.magnetic_nu Kernels ^^^^^^^ @@ -112,6 +118,16 @@ Kernels prism.kernel_en prism.kernel_eu prism.kernel_nu + prism.kernel_eee + prism.kernel_nnn + prism.kernel_uuu + prism.kernel_een + prism.kernel_eeu + prism.kernel_enn + prism.kernel_nnu + prism.kernel_euu + prism.kernel_nuu + prism.kernel_enu Utilities From 6aee0a42fa41a2e286cba5246699d5a84d77273a Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Tue, 8 Oct 2024 08:54:40 -0700 Subject: [PATCH 14/19] Bump pypa/gh-action-pypi-publish from 1.10.2 to 1.10.3 (#102) --- .github/workflows/pypi.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/pypi.yml b/.github/workflows/pypi.yml index 2f1add19..3ed3ce78 100644 --- a/.github/workflows/pypi.yml +++ b/.github/workflows/pypi.yml @@ -102,7 +102,7 @@ jobs: - name: Publish to Test PyPI # Only publish to TestPyPI when a PR is merged (pushed to main) if: success() && github.event_name == 'push' - uses: pypa/gh-action-pypi-publish@v1.10.2 + uses: pypa/gh-action-pypi-publish@v1.10.3 with: repository_url: https://test.pypi.org/legacy/ # Allow existing releases on test PyPI without errors. @@ -112,4 +112,4 @@ jobs: - name: Publish to PyPI # Only publish to PyPI when a release triggers the build if: success() && github.event_name == 'release' - uses: pypa/gh-action-pypi-publish@v1.10.2 + uses: pypa/gh-action-pypi-publish@v1.10.3 From 6714a2289e59ddc98a4b77313291273338b68c44 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 8 Oct 2024 12:46:39 -0700 Subject: [PATCH 15/19] Add changelog for v0.3.0 (#104) Add changelog for Choclo v0.3.0. Include link to new docs in `versions.rst`. --- doc/changes.rst | 36 ++++++++++++++++++++++++++++++++++++ doc/versions.rst | 1 + 2 files changed, 37 insertions(+) diff --git a/doc/changes.rst b/doc/changes.rst index c567fccc..736debef 100644 --- a/doc/changes.rst +++ b/doc/changes.rst @@ -3,6 +3,42 @@ Changelog ========= +Version 0.3.0 +------------- + +Released on: 2024/10/08 + +doi: https://doi.org/10.5281/zenodo.13905447 + +Bug fixes: + +- Fix bug on ``safe_log`` and solves discontinuous magnetic fields (`#100 `__) + +New features: + +- Add forward modelling functions for the magnetic gradiometry components of prisms (`#97 `__) + +Maintenance: + +- Run tests with oldest dependencies on x86 macos (`#83 `__) +- Replace ``_version_generated.py`` for ``_version.py`` in Makefile (`#82 `__) +- Update how output variables are stored in Actions (`#90 `__) +- Move push to codecov to its own job in Actions (`#88 `__) +- Replace ``build`` for ``python-build`` in ``environment.yml`` (`#91 `__) +- Simplify tests for prism magnetic forward funcs (`#96 `__) +- Add some more tests for ``_safe_log`` (`#101 `__) + +Documentation: + +- Replace Sphinx napoleon for numpydoc (`#84 `__) +- Fix style errors in docstrings (`#85 `__) +- Improve docstrings of ``magnetic_field`` functions (`#87 `__) + +This release contains contributions from: + +- Santiago Soler + + Version 0.2.0 ------------- diff --git a/doc/versions.rst b/doc/versions.rst index 9d94cccb..a9fbe402 100644 --- a/doc/versions.rst +++ b/doc/versions.rst @@ -7,6 +7,7 @@ Use the links below to access documentation for specific versions * `Latest release `__ * `Development `__ (reflects the current development branch on GitHub) +* `v0.3.0 `__ * `v0.2.0 `__ * `v0.1.0 `__ * `v0.0.1 `__ From 1392845c087888444516e8d0f2d9d02632ae9a19 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Tue, 12 Nov 2024 10:25:30 -0800 Subject: [PATCH 16/19] Bump pypa/gh-action-pypi-publish from 1.10.3 to 1.12.2 (#107) --- .github/workflows/pypi.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/pypi.yml b/.github/workflows/pypi.yml index 3ed3ce78..ff04a9e3 100644 --- a/.github/workflows/pypi.yml +++ b/.github/workflows/pypi.yml @@ -102,7 +102,7 @@ jobs: - name: Publish to Test PyPI # Only publish to TestPyPI when a PR is merged (pushed to main) if: success() && github.event_name == 'push' - uses: pypa/gh-action-pypi-publish@v1.10.3 + uses: pypa/gh-action-pypi-publish@v1.12.2 with: repository_url: https://test.pypi.org/legacy/ # Allow existing releases on test PyPI without errors. @@ -112,4 +112,4 @@ jobs: - name: Publish to PyPI # Only publish to PyPI when a release triggers the build if: success() && github.event_name == 'release' - uses: pypa/gh-action-pypi-publish@v1.10.3 + uses: pypa/gh-action-pypi-publish@v1.12.2 From f405e25a2ec684a985a8a5162bcbe30716df5ea9 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 14 Nov 2024 13:04:51 -0800 Subject: [PATCH 17/19] Group similar parameters in docstrings (#105) Group similar parameters together in the docstrings of forward and kernel functions for points, dipoles and prisms. This makes the docstrings of the functions shorter and avoid too many repetitions. Configure flake8 for maximum line length in docstrings to 88. --- .flake8 | 2 +- choclo/dipole/_forward.py | 128 ++++++------------- choclo/point/_forward.py | 170 ++++++++----------------- choclo/point/_kernels.py | 160 ++++++----------------- choclo/prism/_gravity.py | 253 ++++++++----------------------------- choclo/prism/_kernels.py | 246 ++++++++++-------------------------- choclo/prism/_magnetic.py | 258 +++++++++----------------------------- choclo/prism/_utils.py | 184 ++++++--------------------- choclo/utils.py | 39 +++--- 9 files changed, 355 insertions(+), 1085 deletions(-) diff --git a/.flake8 b/.flake8 index 7540f326..b7fbefdd 100644 --- a/.flake8 +++ b/.flake8 @@ -2,7 +2,7 @@ [flake8] max-line-length = 88 -max-doc-length = 79 +max-doc-length = 88 ignore = # Too many leading '#' for block comment E266, diff --git a/choclo/dipole/_forward.py b/choclo/dipole/_forward.py index c30f9857..526a2b39 100644 --- a/choclo/dipole/_forward.py +++ b/choclo/dipole/_forward.py @@ -37,39 +37,20 @@ def magnetic_field( Parameters ---------- - easting_p : float - Easting coordinate of the observation point in meters. - northing_p : float - Northing coordinate of the observation point in meters. - upward_p : float - Upward coordinate of the observation point in meters. - easting_q : float - Easting coordinate of the dipole in meters. - northing_q : float - Northing coordinate of the dipole in meters. - upward_q : float - Upward coordinate of the dipole in meters. - magnetic_moment_east : float - The East component of the magnetic moment vector of the dipole. Must be - in :math:`A m^2`. - magnetic_moment_north : float - The North component of the magnetic moment vector of the dipole. Must - be in :math:`A m^2`. - magnetic_moment_up : float - The upward component of the magnetic moment vector of the dipole. Must - be in :math:`A m^2`. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of the observation point in + meters. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of the dipole in meters. + magnetic_moment_east, magnetic_moment_north, magnetic_moment_up : float + The east, north and upward component of the magnetic moment vector of + the dipole. Must be in :math:`A m^2`. Returns ------- - b_e : float - Easting component of the magnetic field generated by the dipole - on the observation point in :math:`\text{T}`. - b_n : float - Northing component of the magnetic field generated by the dipole - on the observation point in :math:`\text{T}`. - b_u : float - Upward component of the magnetic field generated by the dipole - on the observation point in :math:`\text{T}`. + b_e, b_n, b_u : float + Easting, northing and upward components of the magnetic field generated + by the dipole on the observation point in :math:`\text{T}`. Notes ----- @@ -141,27 +122,14 @@ def magnetic_e( Parameters ---------- - easting_p : float - Easting coordinate of the observation point in meters. - northing_p : float - Northing coordinate of the observation point in meters. - upward_p : float - Upward coordinate of the observation point in meters. - easting_q : float - Easting coordinate of the dipole in meters. - northing_q : float - Northing coordinate of the dipole in meters. - upward_q : float - Upward coordinate of the dipole in meters. - magnetic_moment_east : float - The East component of the magnetic moment vector of the dipole. Must be - in :math:`A m^2`. - magnetic_moment_north : float - The North component of the magnetic moment vector of the dipole. Must - be in :math:`A m^2`. - magnetic_moment_up : float - The upward component of the magnetic moment vector of the dipole. Must - be in :math:`A m^2`. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of the observation point in + meters. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of the dipole in meters. + magnetic_moment_east, magnetic_moment_north, magnetic_moment_up : float + The east, north and upward component of the magnetic moment vector of + the dipole. Must be in :math:`A m^2`. Returns ------- @@ -232,27 +200,14 @@ def magnetic_n( Parameters ---------- - easting_p : float - Easting coordinate of the observation point in meters. - northing_p : float - Northing coordinate of the observation point in meters. - upward_p : float - Upward coordinate of the observation point in meters. - easting_q : float - Easting coordinate of the dipole in meters. - northing_q : float - Northing coordinate of the dipole in meters. - upward_q : float - Upward coordinate of the dipole in meters. - magnetic_moment_east : float - The East component of the magnetic moment vector of the dipole. Must be - in :math:`A m^2`. - magnetic_moment_north : float - The North component of the magnetic moment vector of the dipole. Must - be in :math:`A m^2`. - magnetic_moment_up : float - The upward component of the magnetic moment vector of the dipole. Must - be in :math:`A m^2`. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of the observation point in + meters. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of the dipole in meters. + magnetic_moment_east, magnetic_moment_north, magnetic_moment_up : float + The east, north and upward component of the magnetic moment vector of + the dipole. Must be in :math:`A m^2`. Returns ------- @@ -323,27 +278,14 @@ def magnetic_u( Parameters ---------- - easting_p : float - Easting coordinate of the observation point in meters. - northing_p : float - Northing coordinate of the observation point in meters. - upward_p : float - Upward coordinate of the observation point in meters. - easting_q : float - Easting coordinate of the dipole in meters. - northing_q : float - Northing coordinate of the dipole in meters. - upward_q : float - Upward coordinate of the dipole in meters. - magnetic_moment_east : float - The East component of the magnetic moment vector of the dipole. Must be - in :math:`A m^2`. - magnetic_moment_north : float - The North component of the magnetic moment vector of the dipole. Must - be in :math:`A m^2`. - magnetic_moment_up : float - The upward component of the magnetic moment vector of the dipole. Must - be in :math:`A m^2`. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of the observation point in + meters. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of the dipole in meters. + magnetic_moment_east, magnetic_moment_north, magnetic_moment_up : float + The east, north and upward component of the magnetic moment vector of + the dipole. Must be in :math:`A m^2`. Returns ------- diff --git a/choclo/point/_forward.py b/choclo/point/_forward.py index 4332e8f2..c722db26 100644 --- a/choclo/point/_forward.py +++ b/choclo/point/_forward.py @@ -35,18 +35,11 @@ def gravity_pot(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q Parameters ---------- - easting_p : float - Easting coordinate of the observation point in meters. - northing_p : float - Northing coordinate of the observation point in meters. - upward_p : float - Upward coordinate of the observation point in meters. - easting_q : float - Easting coordinate of the point source in meters. - northing_q : float - Northing coordinate of the point source in meters. - upward_q : float - Upward coordinate of the point source in meters. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of the observation point in + meters. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of the point source in meters. mass : float Mass of the point source in kilograms. @@ -90,18 +83,11 @@ def gravity_e(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, Parameters ---------- - easting_p : float - Easting coordinate of the observation point in meters. - northing_p : float - Northing coordinate of the observation point in meters. - upward_p : float - Upward coordinate of the observation point in meters. - easting_q : float - Easting coordinate of the point source in meters. - northing_q : float - Northing coordinate of the point source in meters. - upward_q : float - Upward coordinate of the point source in meters. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of the observation point in + meters. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of the point source in meters. mass : float Mass of the point source in kilograms. @@ -150,18 +136,11 @@ def gravity_n(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, Parameters ---------- - easting_p : float - Easting coordinate of the observation point in meters. - northing_p : float - Northing coordinate of the observation point in meters. - upward_p : float - Upward coordinate of the observation point in meters. - easting_q : float - Easting coordinate of the point source in meters. - northing_q : float - Northing coordinate of the point source in meters. - upward_q : float - Upward coordinate of the point source in meters. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of the observation point in + meters. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of the point source in meters. mass : float Mass of the point source in kilograms. @@ -210,18 +189,11 @@ def gravity_u(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, Parameters ---------- - easting_p : float - Easting coordinate of the observation point in meters. - northing_p : float - Northing coordinate of the observation point in meters. - upward_p : float - Upward coordinate of the observation point in meters. - easting_q : float - Easting coordinate of the point source in meters. - northing_q : float - Northing coordinate of the point source in meters. - upward_q : float - Upward coordinate of the point source in meters. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of the observation point in + meters. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of the point source in meters. mass : float Mass of the point source in kilograms. @@ -270,18 +242,11 @@ def gravity_ee(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, Parameters ---------- - easting_p : float - Easting coordinate of the observation point in meters. - northing_p : float - Northing coordinate of the observation point in meters. - upward_p : float - Upward coordinate of the observation point in meters. - easting_q : float - Easting coordinate of the point source in meters. - northing_q : float - Northing coordinate of the point source in meters. - upward_q : float - Upward coordinate of the point source in meters. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of the observation point in + meters. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of the point source in meters. mass : float Mass of the point source in kilograms. @@ -337,18 +302,11 @@ def gravity_nn(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, Parameters ---------- - easting_p : float - Easting coordinate of the observation point in meters. - northing_p : float - Northing coordinate of the observation point in meters. - upward_p : float - Upward coordinate of the observation point in meters. - easting_q : float - Easting coordinate of the point source in meters. - northing_q : float - Northing coordinate of the point source in meters. - upward_q : float - Upward coordinate of the point source in meters. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of the observation point in + meters. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of the point source in meters. mass : float Mass of the point source in kilograms. @@ -404,18 +362,11 @@ def gravity_uu(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, Parameters ---------- - easting_p : float - Easting coordinate of the observation point in meters. - northing_p : float - Northing coordinate of the observation point in meters. - upward_p : float - Upward coordinate of the observation point in meters. - easting_q : float - Easting coordinate of the point source in meters. - northing_q : float - Northing coordinate of the point source in meters. - upward_q : float - Upward coordinate of the point source in meters. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of the observation point in + meters. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of the point source in meters. mass : float Mass of the point source in kilograms. @@ -471,18 +422,11 @@ def gravity_en(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, Parameters ---------- - easting_p : float - Easting coordinate of the observation point in meters. - northing_p : float - Northing coordinate of the observation point in meters. - upward_p : float - Upward coordinate of the observation point in meters. - easting_q : float - Easting coordinate of the point source in meters. - northing_q : float - Northing coordinate of the point source in meters. - upward_q : float - Upward coordinate of the point source in meters. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of the observation point in + meters. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of the point source in meters. mass : float Mass of the point source in kilograms. @@ -533,18 +477,11 @@ def gravity_eu(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, Parameters ---------- - easting_p : float - Easting coordinate of the observation point in meters. - northing_p : float - Northing coordinate of the observation point in meters. - upward_p : float - Upward coordinate of the observation point in meters. - easting_q : float - Easting coordinate of the point source in meters. - northing_q : float - Northing coordinate of the point source in meters. - upward_q : float - Upward coordinate of the point source in meters. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of the observation point in + meters. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of the point source in meters. mass : float Mass of the point source in kilograms. @@ -595,18 +532,11 @@ def gravity_nu(easting_p, northing_p, upward_p, easting_q, northing_q, upward_q, Parameters ---------- - easting_p : float - Easting coordinate of the observation point in meters. - northing_p : float - Northing coordinate of the observation point in meters. - upward_p : float - Upward coordinate of the observation point in meters. - easting_q : float - Easting coordinate of the point source in meters. - northing_q : float - Northing coordinate of the point source in meters. - upward_q : float - Upward coordinate of the point source in meters. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of the observation point in + meters. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of the point source in meters. mass : float Mass of the point source in kilograms. diff --git a/choclo/point/_kernels.py b/choclo/point/_kernels.py index dc061d5a..faa815b4 100644 --- a/choclo/point/_kernels.py +++ b/choclo/point/_kernels.py @@ -24,18 +24,10 @@ def kernel_pot( Parameters ---------- - easting_p : float - Easting coordinate of point :math:`\mathbf{p}`. - northing_p : float - Northing coordinate of point :math:`\mathbf{p}`. - upward_p : float - Upward coordinate of point :math:`\mathbf{p}`. - easting_q : float - Easting coordinate of point :math:`\mathbf{q}`. - northing_q : float - Northing coordinate of point :math:`\mathbf{q}`. - upward_q : float - Upward coordinate of point :math:`\mathbf{q}`. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of point :math:`\mathbf{p}`. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of point :math:`\mathbf{q}`. distance : float Euclidean distance between points :math:`\mathbf{p}` and :math:`\mathbf{q}`. @@ -80,18 +72,10 @@ def kernel_e( Parameters ---------- - easting_p : float - Easting coordinate of point :math:`\mathbf{p}`. - northing_p : float - Northing coordinate of point :math:`\mathbf{p}`. - upward_p : float - Upward coordinate of point :math:`\mathbf{p}`. - easting_q : float - Easting coordinate of point :math:`\mathbf{q}`. - northing_q : float - Northing coordinate of point :math:`\mathbf{q}`. - upward_q : float - Upward coordinate of point :math:`\mathbf{q}`. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of point :math:`\mathbf{p}`. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of point :math:`\mathbf{q}`. distance : float Euclidean distance between points :math:`\mathbf{p}` and :math:`\mathbf{q}`. @@ -141,18 +125,10 @@ def kernel_n( Parameters ---------- - easting_p : float - Easting coordinate of point :math:`\mathbf{p}`. - northing_p : float - Northing coordinate of point :math:`\mathbf{p}`. - upward_p : float - Upward coordinate of point :math:`\mathbf{p}`. - easting_q : float - Easting coordinate of point :math:`\mathbf{q}`. - northing_q : float - Northing coordinate of point :math:`\mathbf{q}`. - upward_q : float - Upward coordinate of point :math:`\mathbf{q}`. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of point :math:`\mathbf{p}`. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of point :math:`\mathbf{q}`. distance : float Euclidean distance between points :math:`\mathbf{p}` and :math:`\mathbf{q}`. @@ -202,18 +178,10 @@ def kernel_u( Parameters ---------- - easting_p : float - Easting coordinate of point :math:`\mathbf{p}`. - northing_p : float - Northing coordinate of point :math:`\mathbf{p}`. - upward_p : float - Upward coordinate of point :math:`\mathbf{p}`. - easting_q : float - Easting coordinate of point :math:`\mathbf{q}`. - northing_q : float - Northing coordinate of point :math:`\mathbf{q}`. - upward_q : float - Upward coordinate of point :math:`\mathbf{q}`. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of point :math:`\mathbf{p}`. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of point :math:`\mathbf{q}`. distance : float Euclidean distance between points :math:`\mathbf{p}` and :math:`\mathbf{q}`. @@ -263,18 +231,10 @@ def kernel_ee( Parameters ---------- - easting_p : float - Easting coordinate of point :math:`\mathbf{p}`. - northing_p : float - Northing coordinate of point :math:`\mathbf{p}`. - upward_p : float - Upward coordinate of point :math:`\mathbf{p}`. - easting_q : float - Easting coordinate of point :math:`\mathbf{q}`. - northing_q : float - Northing coordinate of point :math:`\mathbf{q}`. - upward_q : float - Upward coordinate of point :math:`\mathbf{q}`. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of point :math:`\mathbf{p}`. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of point :math:`\mathbf{q}`. distance : float Euclidean distance between points :math:`\mathbf{p}` and :math:`\mathbf{q}`. @@ -329,18 +289,10 @@ def kernel_nn( Parameters ---------- - easting_p : float - Easting coordinate of point :math:`\mathbf{p}`. - northing_p : float - Northing coordinate of point :math:`\mathbf{p}`. - upward_p : float - Upward coordinate of point :math:`\mathbf{p}`. - easting_q : float - Easting coordinate of point :math:`\mathbf{q}`. - northing_q : float - Northing coordinate of point :math:`\mathbf{q}`. - upward_q : float - Upward coordinate of point :math:`\mathbf{q}`. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of point :math:`\mathbf{p}`. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of point :math:`\mathbf{q}`. distance : float Euclidean distance between points :math:`\mathbf{p}` and :math:`\mathbf{q}`. @@ -395,18 +347,10 @@ def kernel_uu( Parameters ---------- - easting_p : float - Easting coordinate of point :math:`\mathbf{p}`. - northing_p : float - Northing coordinate of point :math:`\mathbf{p}`. - upward_p : float - Upward coordinate of point :math:`\mathbf{p}`. - easting_q : float - Easting coordinate of point :math:`\mathbf{q}`. - northing_q : float - Northing coordinate of point :math:`\mathbf{q}`. - upward_q : float - Upward coordinate of point :math:`\mathbf{q}`. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of point :math:`\mathbf{p}`. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of point :math:`\mathbf{q}`. distance : float Euclidean distance between points :math:`\mathbf{p}` and :math:`\mathbf{q}`. @@ -461,18 +405,10 @@ def kernel_en( Parameters ---------- - easting_p : float - Easting coordinate of point :math:`\mathbf{p}`. - northing_p : float - Northing coordinate of point :math:`\mathbf{p}`. - upward_p : float - Upward coordinate of point :math:`\mathbf{p}`. - easting_q : float - Easting coordinate of point :math:`\mathbf{q}`. - northing_q : float - Northing coordinate of point :math:`\mathbf{q}`. - upward_q : float - Upward coordinate of point :math:`\mathbf{q}`. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of point :math:`\mathbf{p}`. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of point :math:`\mathbf{q}`. distance : float Euclidean distance between points :math:`\mathbf{p}` and :math:`\mathbf{q}`. @@ -522,18 +458,10 @@ def kernel_eu( Parameters ---------- - easting_p : float - Easting coordinate of point :math:`\mathbf{p}`. - northing_p : float - Northing coordinate of point :math:`\mathbf{p}`. - upward_p : float - Upward coordinate of point :math:`\mathbf{p}`. - easting_q : float - Easting coordinate of point :math:`\mathbf{q}`. - northing_q : float - Northing coordinate of point :math:`\mathbf{q}`. - upward_q : float - Upward coordinate of point :math:`\mathbf{q}`. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of point :math:`\mathbf{p}`. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of point :math:`\mathbf{q}`. distance : float Euclidean distance between points :math:`\mathbf{p}` and :math:`\mathbf{q}`. @@ -583,18 +511,10 @@ def kernel_nu( Parameters ---------- - easting_p : float - Easting coordinate of point :math:`\mathbf{p}`. - northing_p : float - Northing coordinate of point :math:`\mathbf{p}`. - upward_p : float - Upward coordinate of point :math:`\mathbf{p}`. - easting_q : float - Easting coordinate of point :math:`\mathbf{q}`. - northing_q : float - Northing coordinate of point :math:`\mathbf{q}`. - upward_q : float - Upward coordinate of point :math:`\mathbf{q}`. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of point :math:`\mathbf{p}`. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of point :math:`\mathbf{q}`. distance : float Euclidean distance between points :math:`\mathbf{p}` and :math:`\mathbf{q}`. diff --git a/choclo/prism/_gravity.py b/choclo/prism/_gravity.py index d03840ed..03239639 100644 --- a/choclo/prism/_gravity.py +++ b/choclo/prism/_gravity.py @@ -54,24 +54,11 @@ def gravity_pot( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. density : float Density of the rectangular prism in kilograms per cubic meter. @@ -200,24 +187,11 @@ def gravity_e( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. density : float Density of the rectangular prism in kilograms per cubic meter. @@ -347,24 +321,11 @@ def gravity_n( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. density : float Density of the rectangular prism in kilograms per cubic meter. @@ -494,24 +455,11 @@ def gravity_u( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. density : float Density of the rectangular prism in kilograms per cubic meter. @@ -647,24 +595,11 @@ def gravity_ee( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. density : float Density of the rectangular prism in kilograms per cubic meter. @@ -818,24 +753,11 @@ def gravity_nn( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. density : float Density of the rectangular prism in kilograms per cubic meter. @@ -989,24 +911,11 @@ def gravity_uu( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. density : float Density of the rectangular prism in kilograms per cubic meter. @@ -1160,24 +1069,11 @@ def gravity_en( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. density : float Density of the rectangular prism in kilograms per cubic meter. @@ -1307,24 +1203,11 @@ def gravity_eu( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. density : float Density of the rectangular prism in kilograms per cubic meter. @@ -1455,24 +1338,11 @@ def gravity_nu( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. density : float Density of the rectangular prism in kilograms per cubic meter. @@ -1600,24 +1470,11 @@ def _evaluate_kernel( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. kernel : callable Kernel function that will be evaluated on each one of the shifted vertices of the prism. diff --git a/choclo/prism/_kernels.py b/choclo/prism/_kernels.py index b3f8f845..2d4126fc 100644 --- a/choclo/prism/_kernels.py +++ b/choclo/prism/_kernels.py @@ -28,15 +28,9 @@ def kernel_pot(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -120,15 +114,9 @@ def kernel_e(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -215,15 +203,9 @@ def kernel_n(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -310,15 +292,9 @@ def kernel_u(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -406,15 +382,9 @@ def kernel_ee(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -471,15 +441,9 @@ def kernel_nn(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -536,15 +500,9 @@ def kernel_uu(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -601,15 +559,9 @@ def kernel_en(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -666,15 +618,9 @@ def kernel_eu(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -731,15 +677,9 @@ def kernel_nu(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -833,7 +773,7 @@ def _safe_log(x, y, z, r): x : float Shifted coordinate of the prism vertex that will be used for evaluating the log. - y, z : floats + y, z : float The other two shifted coordinates of the prism vertex. They are used to determine if ``abs(x) == r`` with more accuracy, and to evaluate the log function to avoid floating point errors when :math:`x < 0` and @@ -914,15 +854,9 @@ def kernel_eee(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -972,15 +906,9 @@ def kernel_nnn(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -1030,15 +958,9 @@ def kernel_uuu(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -1088,15 +1010,9 @@ def kernel_een(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -1146,15 +1062,9 @@ def kernel_eeu(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -1204,15 +1114,9 @@ def kernel_enn(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -1262,15 +1166,9 @@ def kernel_nnu(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -1320,15 +1218,9 @@ def kernel_euu(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -1378,15 +1270,9 @@ def kernel_nuu(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -1436,15 +1322,9 @@ def kernel_enu(easting, northing, upward, radius): 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. + easting, northing, upward : float + Shifted easting, northing and upward coordinates 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. @@ -1498,7 +1378,7 @@ def _kernel_iii(x_i, x_j, x_k, radius): Parameters ---------- - x_i, x_j, x_k : floats + x_i, x_j, x_k : float 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 @@ -1558,7 +1438,7 @@ def _kernel_iij(x_i, x_j, x_k, radius): Parameters ---------- - x_i, x_j, x_k : floats + x_i, x_j, x_k : float 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 diff --git a/choclo/prism/_magnetic.py b/choclo/prism/_magnetic.py index d8ff39c9..598e84a3 100644 --- a/choclo/prism/_magnetic.py +++ b/choclo/prism/_magnetic.py @@ -68,24 +68,11 @@ def magnetic_field( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries 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}`. @@ -98,19 +85,9 @@ def magnetic_field( Returns ------- - b_e : float - Easting component of the magnetic field generated by the prism - on the observation point in :math:`\text{T}`. - It will be ``numpy.nan`` if the observation point falls in a singular - point: prism vertices, prism edges or interior points. - b_n : float - Northing component of the magnetic field generated by the prism - on the observation point in :math:`\text{T}`. - It will be ``numpy.nan`` if the observation point falls in a singular - point: prism vertices, prism edges or interior points. - b_u : float - Upward component of the magnetic field generated by the prism - on the observation point in :math:`\text{T}`. + b_e, b_n, b_u : float + Easting, northing and upward component of the magnetic field generated + by the prism on the observation point in :math:`\text{T}`. It will be ``numpy.nan`` if the observation point falls in a singular point: prism vertices, prism edges or interior points. @@ -396,24 +373,11 @@ def magnetic_e( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries 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}`. @@ -583,24 +547,11 @@ def magnetic_n( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries 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}`. @@ -770,24 +721,11 @@ def magnetic_u( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries 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}`. @@ -957,24 +895,11 @@ def magnetic_ee( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries 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}`. @@ -1105,24 +1030,11 @@ def magnetic_nn( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries 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}`. @@ -1253,24 +1165,11 @@ def magnetic_uu( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries 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}`. @@ -1401,24 +1300,11 @@ def magnetic_en( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries 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}`. @@ -1549,24 +1435,11 @@ def magnetic_eu( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries 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}`. @@ -1697,24 +1570,11 @@ def magnetic_nu( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries 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}`. @@ -1849,15 +1709,11 @@ def _calculate_component( Parameters ---------- - easting, northing, upward : floats + easting, northing, upward : float Coordinates of the observation point. Must be in meters. - prism_west, prism_east : floats - The easting boundaries of the prism. Must be in meters. - prism_south, prism_north : floats - The northing boundaries of the prism. Must be in meters. - prism_bottom, prism_top : floats - The upward boundaries of the prism. Must be in meters. - magnetization_east, magnetization_north, magnetization_up : floats + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. + magnetization_east, magnetization_north, magnetization_up : float The components of the magnetization vector of the prism. Must be in :math:`A m^{-1}`. kernel_i, kernel_j, kernel_k : callables diff --git a/choclo/prism/_utils.py b/choclo/prism/_utils.py index 30998639..2c43b172 100644 --- a/choclo/prism/_utils.py +++ b/choclo/prism/_utils.py @@ -31,24 +31,11 @@ def is_interior_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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. Returns ------- @@ -83,24 +70,11 @@ def is_point_on_edge( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. Returns ------- @@ -167,24 +141,11 @@ def is_point_on_easting_edge( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. Returns ------- @@ -222,24 +183,11 @@ def is_point_on_northing_edge( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. Returns ------- @@ -277,24 +225,11 @@ def is_point_on_upward_edge( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. Returns ------- @@ -332,24 +267,11 @@ def is_point_on_east_face( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. Returns ------- @@ -388,24 +310,11 @@ def is_point_on_north_face( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. Returns ------- @@ -443,24 +352,11 @@ def is_point_on_top_face( 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. + easting, northing, upward : float + Easting, northing and upward coordinates of the observation point. Must + be in meters. + prism_west, prism_east, prism_south, prism_north, prism_bottom, prism_top : float + The boundaries of the prism. Must be in meters. Returns ------- diff --git a/choclo/utils.py b/choclo/utils.py index 1a07325f..76cdaeae 100644 --- a/choclo/utils.py +++ b/choclo/utils.py @@ -24,18 +24,10 @@ def distance_cartesian( Parameters ---------- - easting_p : float - Easting coordinate of point :math:`\mathbf{p}`. - northing_p : float - Northing coordinate of point :math:`\mathbf{p}`. - upward_p : float - Upward coordinate of point :math:`\mathbf{p}`. - easting_q : float - Easting coordinate of point :math:`\mathbf{q}`. - northing_q : float - Northing coordinate of point :math:`\mathbf{q}`. - upward_q : float - Upward coordinate of point :math:`\mathbf{q}`. + easting_p, northing_p, upward_p : float + Easting, northing and upward coordinates of point :math:`\mathbf{p}`. + easting_q, northing_q, upward_q : float + Easting, northing and upward coordinates of point :math:`\mathbf{q}`. Returns ------- @@ -75,23 +67,20 @@ def distance_spherical( Parameters ---------- - longitude_p : float - Longitude coordinate of point :math:`\mathbf{p}` in degrees. - latitude_p : float - Latitude coordinate of point :math:`\mathbf{p}` in degrees. - radius_p : float - Radial coordinate of point :math:`\mathbf{p}` in meters. - longitude_q : float - Longitude coordinate of point :math:`\mathbf{q}` in degrees. - latitude_q : float - Latitude coordinate of point :math:`\mathbf{q}` in degrees. - radius_q : float - Radial coordinate of point :math:`\mathbf{q}` in meters. + longitude_p, latitude_p, radius_p : float + Longitude, latitude and radial coordinates of point :math:`\mathbf{p}`. + Longitude and latitude must be in degrees. Radial coordinate should be + in meters. + longitude_q, latitude_q, radius_q : float + Longitude, latitude and radial coordinates of point :math:`\mathbf{q}`. + Longitude and latitude must be in degrees. Radial coordinate should be + in meters. Returns ------- distance : float - Euclidean distance between ``point_p`` and ``point_q``. + Euclidean distance between point :math:`\mathbf{p}` and point + :math:`\mathbf{q}`. Notes ----- From a81aaf4f93212473721e04fd062f603c3f83e613 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 18 Nov 2024 15:29:41 -0800 Subject: [PATCH 18/19] Bump codecov/codecov-action from 4 to 5 (#109) --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 6cbd3440..ca458c30 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -174,7 +174,7 @@ jobs: run: ls -l -R . - name: Upload coverage reports to Codecov - uses: codecov/codecov-action@v4 + uses: codecov/codecov-action@v5 with: # Upload all coverage report files files: ./coverage_*/coverage.xml From e38581158f6a0679d95115c41d719bc8cc5a3819 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Wed, 20 Nov 2024 15:31:22 -0800 Subject: [PATCH 19/19] Make prism kernels return `np.nan` on vertices (#108) Make the 2nd and 3rd order kernel functions for prisms to return `np.nan` when the observation points matches any vertex of the prisms. Fixes division-by-zero error when evaluating `kernel_enu` when `radius` is zero, and makes the behavior consistent with the other third-order kernels. --- choclo/prism/_kernels.py | 58 ++++++++++++++++++----- choclo/tests/test_prism_kernels.py | 76 ++++++++++++++++++++++++++++++ 2 files changed, 122 insertions(+), 12 deletions(-) create mode 100644 choclo/tests/test_prism_kernels.py diff --git a/choclo/prism/_kernels.py b/choclo/prism/_kernels.py index 2d4126fc..14b95b12 100644 --- a/choclo/prism/_kernels.py +++ b/choclo/prism/_kernels.py @@ -394,6 +394,7 @@ def kernel_ee(easting, northing, upward, radius): kernel : float Value of the kernel function for the easting-easting component of the tensor due to a rectangular prism evaluated on a single vertex. + Return ``np.nan`` if ``radius`` is zero. Notes ----- @@ -422,6 +423,8 @@ def kernel_ee(easting, northing, upward, radius): - [Nagy2002]_ - [Fukushima2020]_ """ + if radius == 0.0: + return np.nan return -_safe_atan2(northing * upward, easting * radius) @@ -453,6 +456,7 @@ def kernel_nn(easting, northing, upward, radius): kernel : float Value of the kernel function for the northing-northing component of the tensor due to a rectangular prism evaluated on a single vertex. + Return ``np.nan`` if ``radius`` is zero. Notes ----- @@ -481,6 +485,8 @@ def kernel_nn(easting, northing, upward, radius): - [Nagy2002]_ - [Fukushima2020]_ """ + if radius == 0.0: + return np.nan return -_safe_atan2(easting * upward, northing * radius) @@ -512,6 +518,7 @@ def kernel_uu(easting, northing, upward, radius): kernel : float Value of the kernel function for the upward-upward component of the tensor due to a rectangular prism evaluated on a single vertex. + Return ``np.nan`` if ``radius`` is zero. Notes ----- @@ -540,6 +547,8 @@ def kernel_uu(easting, northing, upward, radius): - [Nagy2002]_ - [Fukushima2020]_ """ + if radius == 0.0: + return np.nan return -_safe_atan2(easting * northing, upward * radius) @@ -571,6 +580,7 @@ def kernel_en(easting, northing, upward, radius): kernel : float Value of the kernel function for the easting-northing component of the tensor due to a rectangular prism evaluated on a single vertex. + Return ``np.nan`` if ``radius`` is zero. Notes ----- @@ -599,6 +609,8 @@ def kernel_en(easting, northing, upward, radius): - [Nagy2002]_ - [Fukushima2020]_ """ + if radius == 0.0: + return np.nan return _safe_log(upward, easting, northing, radius) @@ -630,6 +642,7 @@ def kernel_eu(easting, northing, upward, radius): kernel : float Value of the kernel function for the easting-upward component of the tensor due to a rectangular prism evaluated on a single vertex. + Return ``np.nan`` if ``radius`` is zero. Notes ----- @@ -658,6 +671,8 @@ def kernel_eu(easting, northing, upward, radius): - [Nagy2002]_ - [Fukushima2020]_ """ + if radius == 0.0: + return np.nan return _safe_log(northing, easting, upward, radius) @@ -689,6 +704,7 @@ def kernel_nu(easting, northing, upward, radius): kernel : float Value of the kernel function for the northing-upward component of the tensor due to a rectangular prism evaluated on a single vertex. + Return ``np.nan`` if ``radius`` is zero. Notes ----- @@ -717,6 +733,8 @@ def kernel_nu(easting, northing, upward, radius): - [Nagy2002]_ - [Fukushima2020]_ """ + if radius == 0.0: + return np.nan return _safe_log(easting, northing, upward, radius) @@ -866,7 +884,8 @@ def kernel_eee(easting, northing, upward, radius): 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. + vertex. Return ``np.nan`` if ``radius`` is zero (i.e. observation point + is on the vertex). Notes ----- @@ -918,7 +937,8 @@ def kernel_nnn(easting, northing, upward, radius): 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. + a single vertex. Return ``np.nan`` if ``radius`` is zero (i.e. + observation point is on the vertex). Notes ----- @@ -970,7 +990,8 @@ def kernel_uuu(easting, northing, upward, radius): 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. + vertex. Return ``np.nan`` if ``radius`` is zero (i.e. observation point + is on the vertex). Notes ----- @@ -1022,7 +1043,8 @@ def kernel_een(easting, northing, upward, radius): 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. + vertex. Return ``np.nan`` if ``radius`` is zero (i.e. observation point + is on the vertex). Notes ----- @@ -1074,7 +1096,8 @@ def kernel_eeu(easting, northing, upward, radius): 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. + vertex. Return ``np.nan`` if ``radius`` is zero (i.e. observation point + is on the vertex). Notes ----- @@ -1126,7 +1149,8 @@ def kernel_enn(easting, northing, upward, radius): 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. + a single vertex. Return ``np.nan`` if ``radius`` is zero (i.e. + observation point is on the vertex). Notes ----- @@ -1178,7 +1202,8 @@ def kernel_nnu(easting, northing, upward, radius): 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. + a single vertex. Return ``np.nan`` if ``radius`` is zero (i.e. + observation point is on the vertex). Notes ----- @@ -1230,7 +1255,8 @@ def kernel_euu(easting, northing, upward, radius): 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. + a single vertex. Return ``np.nan`` if ``radius`` is zero (i.e. + observation point is on the vertex). Notes ----- @@ -1282,7 +1308,8 @@ def kernel_nuu(easting, northing, upward, radius): 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. + a single vertex. Return ``np.nan`` if ``radius`` is zero (i.e. + observation point is on the vertex). Notes ----- @@ -1334,7 +1361,8 @@ def kernel_enu(easting, northing, upward, radius): 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. + vertex. Return ``np.nan`` if ``radius`` is zero (i.e. observation point + is on the vertex). Notes ----- @@ -1352,6 +1380,8 @@ def kernel_enu(easting, northing, upward, radius): >>> float(kernel_enu(x, y, z, r)) # doctest: +NUMBER -0.1480061 """ + if radius == 0.0: + return np.nan return -1 / radius @@ -1387,7 +1417,7 @@ def _kernel_iii(x_i, x_j, x_k, radius): Returns ------- kernel : float - Value of the kernel function. + Value of the kernel function. Return ``np.nan`` if ``radius`` is zero. Examples -------- @@ -1403,6 +1433,8 @@ 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 radius == 0.0: + return np.nan 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 @@ -1447,7 +1479,7 @@ def _kernel_iij(x_i, x_j, x_k, radius): Returns ------- kernel : float - Value of the kernel function. + Value of the kernel function. Return ``np.nan`` if ``radius`` is zero. Examples -------- @@ -1466,6 +1498,8 @@ 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) """ + if radius == 0.0: + return np.nan if x_i == 0 and x_j == 0: return 0.0 return -x_i / ((x_k + radius) * radius) diff --git a/choclo/tests/test_prism_kernels.py b/choclo/tests/test_prism_kernels.py new file mode 100644 index 00000000..b0f1072c --- /dev/null +++ b/choclo/tests/test_prism_kernels.py @@ -0,0 +1,76 @@ +# Copyright (c) 2022 The Choclo Developers. +# Distributed under the terms of the BSD 3-Clause License. +# SPDX-License-Identifier: BSD-3-Clause +# +# This code is part of the Fatiando a Terra project (https://www.fatiando.org) +# +""" +Additional tests to prism kernels. +""" + +import numpy as np +import pytest + +from choclo.prism import ( + kernel_ee, + kernel_eee, + kernel_een, + kernel_eeu, + kernel_en, + kernel_enn, + kernel_enu, + kernel_eu, + kernel_euu, + kernel_nn, + kernel_nnn, + kernel_nnu, + kernel_nu, + kernel_nuu, + kernel_uu, + kernel_uuu, +) + + +class TestThirdOrderKernelsOnVertex: + """ + Test if third-order kernels evaluated on a vertex return np.nan. + + This functionality is meant to be supported for the public usage of this + kernels, although the forward modelling functions within Choclo won't use + it since they filter out singular points before they call kernels. + """ + + SECOND_ORDER_KERNELS = ( + kernel_ee, + kernel_nn, + kernel_uu, + kernel_en, + kernel_eu, + kernel_nu, + ) + THIRD_ORDER_KERNELS = ( + kernel_eee, + kernel_nnn, + kernel_uuu, + kernel_een, + kernel_eeu, + kernel_enn, + kernel_nnu, + kernel_euu, + kernel_nuu, + kernel_enu, + ) + + @pytest.mark.parametrize("kernel", THIRD_ORDER_KERNELS) + def test_third_order_kernels_on_vertex(self, kernel): + easting, northing, upward = 0.0, 0.0, 0.0 + radius = 0.0 + result = kernel(easting, northing, upward, radius) + assert np.isnan(result) + + @pytest.mark.parametrize("kernel", SECOND_ORDER_KERNELS) + def test_second_order_kernels_on_vertex(self, kernel): + easting, northing, upward = 0.0, 0.0, 0.0 + radius = 0.0 + result = kernel(easting, northing, upward, radius) + assert np.isnan(result)