From 351795274812abc2050522f8b94fd9c62cfb6515 Mon Sep 17 00:00:00 2001 From: Leonardo Uieda Date: Fri, 22 Aug 2025 15:04:27 -0300 Subject: [PATCH] Adapt forward modeling for multiple components If we want to work with the diamond directions, we'll need to start calculating the full magnetic field vector of the dipoles. Adapt the functions to have a `field` argument like in Harmonica. Use xarray directly to make the grid instead of Verde. Move forward modeling functions to a `_forward.py` module. --- src/magali/_forward.py | 122 +++++++++++++++++++++++++++++++++++++++ src/magali/_inversion.py | 32 +++++----- src/magali/_synthetic.py | 72 ----------------------- 3 files changed, 139 insertions(+), 87 deletions(-) create mode 100644 src/magali/_forward.py diff --git a/src/magali/_forward.py b/src/magali/_forward.py new file mode 100644 index 00000000..9ea08f1f --- /dev/null +++ b/src/magali/_forward.py @@ -0,0 +1,122 @@ +# Copyright (c) 2024 The Magali 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) +# +""" +Forward modeling code to calculate magnetic fields of finite sources. +""" + +import harmonica as hm +import verde as vd +import xarray as xr + +from ._units import coordinates_micrometer_to_meter + + +def dipole_magnetic(coordinates, dipole_coordinates, dipole_moments, field): + """ + Compute the magnetic field produced by one or more dipoles. + + Uses :func:`harmonica.dipole_magnetic` for the calculations and handles the + necessary unit versions from micrometers to meters. + + Parameters + ---------- + coordinates : tuple of float + Observation point coordinates in micrometers (μm). + dipole_coordinates : tuple of float + Dipole location coordinates in micrometers (μm). + dipole_moments : tuple of float + Dipole moment components (Am²). + field : str + Magnetic field that will be computed. The available fields are: + + - The full magnetic vector: ``b`` + - x component of the magnetic vector: ``b_x`` + - y component of the magnetic vector: ``b_y`` + - z (upward) component of the magnetic vector: ``b_z`` + + Returns + ------- + magnetic_field : array or tuple of arrays + Computed magnetic field on every observation point in :math:`nT`. + If ``field`` is set to a single component, then a single array with the + computed magnetic field component will be returned. + If ``field`` is set to ``"b"``, then a tuple containing the three + components of the magnetic vector will be returned in the following + order: ``b_x``, ``b_y``, ``b_z``. + """ + coordinates_m = coordinates_micrometer_to_meter(coordinates) + dipole_coordinates_m = coordinates_micrometer_to_meter(dipole_coordinates) + translator = {"b_x": "b_e", "b_y": "b_n", "b_z": "b_u", "b": "b"} + if field not in translator: + message = ( + f"Invalid field argument '{field}'. Must be one of {translator.keys()}." + ) + raise ValueError(message) + return hm.dipole_magnetic( + coordinates_m, dipole_coordinates_m, dipole_moments, field=translator[field] + ) + + +def dipole_magnetic_grid( + region, spacing, sensor_sample_distance, dipole_coordinates, dipole_moments, field +): + """ + Generate a grid of the magnetic field of one or more dipoles. + + Uses :func:`magali.dipole_magnetic` to generate a regular grid of the + magnetic field produced by a dipole model. The grid is returned as an + :mod:`xarray` container type. + + Parameters + ---------- + region : tuple = (x_min, x_max, y_min, y_max) + The bounding box of the grid in micrometers (μm). + spacing : float + Grid spacing in micrometers (μm). + sensor_sample_distance : float + Distance of the sample from the grid in micrometers (μm). + dipole_coordinates : tuple = (x, y, z) + Arrays with the x, y, and z coordinates of the dipoles in micrometers + (μm). The arrays can have any shape as long as they all have the same + shape. + dipole_moments : tuple = (mx, my, mz) + Arrays with the x, y, and z components of the dipole moment of each + dipole in Am². The arrays should have a shape compatible with the + dipole coordinate arrays. + field : str + Magnetic field that will be computed. The available fields are: + + - The full magnetic vector: ``b`` + - x component of the magnetic vector: ``b_x`` + - y component of the magnetic vector: ``b_y`` + - z (upward) component of the magnetic vector: ``b_z`` + + Returns + ------- + grid : :class:`xarray.DataArray` or :class:`xarray.Dataset` + Gridded magnetic field of the dipole model given. If *field* is + a single component, the grid is returned as + a :class:`xarray.DataArray`. If ``field == "b"``, the three components + of the magnetic field will be returned as a :class:`xarray.Dataset`. + """ + coordinates = vd.grid_coordinates( + region=region, # µm + spacing=spacing, # µm + extra_coords=sensor_sample_distance, + ) + data = dipole_magnetic(coordinates, dipole_coordinates, dipole_moments, field) + dims = ("y", "x") + coords = { + "x": coordinates[0][0, :], + "y": coordinates[1][:, 0], + "z": (dims, coordinates[2]), + } + + grid.x.attrs = {"units": "µm"} + grid.y.attrs = {"units": "µm"} + grid.bz.attrs = {"long_name": "vertical magnetic field", "units": "nT"} + return data.bz diff --git a/src/magali/_inversion.py b/src/magali/_inversion.py index 1063e93b..5a944f12 100644 --- a/src/magali/_inversion.py +++ b/src/magali/_inversion.py @@ -57,9 +57,9 @@ def fit(self, coordinates, data): Parameters ---------- coordinates : tuple = (x, y, z) - Arrays with the x, y, and z coordinates of the observations points. - The arrays can have any shape as long as they all have the same - shape. + Arrays with the x, y, and z coordinates of the observations points + in micrometers (μm). The arrays can have any shape as long as they + all have the same shape. data : array Array with the observed Bz component of the magnetic field (in nT) at the locations provided in *coordinates*. Must have the same @@ -87,9 +87,9 @@ def jacobian(self, coordinates): Parameters ---------- coordinates : tuple = (x, y, z) - Arrays with the x, y, and z coordinates of the observations points. - The arrays can have any shape as long as they all have the same - shape. + Arrays with the x, y, and z coordinates of the observations points + in micrometers (μm). The arrays can have any shape as long as they + all have the same shape. Returns ------- @@ -185,9 +185,10 @@ def predict(self, coordinates): Parameters ---------- - coordinates : tuple of array-like - Arrays with the x, y, z coordinates of the observation points, - in µm. The arrays must have the same shape. + coordinates : tuple = (x, y, z) + Arrays with the x, y, and z coordinates of the computation points + in micrometers (μm). The arrays can have any shape as long as they + all have the same shape. Returns ------- @@ -221,9 +222,9 @@ def jacobian(self, coordinates, location, moment, jacobian): Parameters ---------- coordinates : tuple = (x, y, z) - Arrays with the x, y, and z coordinates of the observations points. - The arrays can have any shape as long as they all have the same - shape. + Arrays with the x, y, and z coordinates of the observation points + in micrometers (μm). The arrays can have any shape as long as they + all have the same shape. location : array-like Dipole location (x, y, z), in µm. moment : array-like @@ -286,9 +287,10 @@ def fit(self, coordinates, data): Parameters ---------- - coordinates : tuple of array-like - Arrays with the x, y, z coordinates of the observation points. - The arrays can have any shape as long as they all have the same shape. + coordinates : tuple = (x, y, z) + Arrays with the x, y, and z coordinates of the observation points + in micrometers (μm). The arrays can have any shape as long as they + all have the same shape. data : array-like Observed Bz component of the magnetic field (in nT) at the observation points. Must have the same shape as the coordinate arrays. diff --git a/src/magali/_synthetic.py b/src/magali/_synthetic.py index 1e7d3a97..dd1a9ffd 100644 --- a/src/magali/_synthetic.py +++ b/src/magali/_synthetic.py @@ -10,9 +10,6 @@ import harmonica as hm import numpy as np -import verde as vd - -from ._units import coordinates_micrometer_to_meter def random_directions( @@ -145,72 +142,3 @@ def _rotate_vector(x, y, z, inclination, declination): y_r = np.sin(phi) * (x * np.cos(theta) + z * np.sin(theta)) + y * np.cos(phi) z_r = -x * np.sin(theta) + z * np.cos(theta) return x_r, y_r, z_r - - -def dipole_bz(coordinates, dipole_coordinates, dipole_moments): - """ - Compute the vertical component of the magnetic field (Bz) produced by a magnetic dipole. - - Parameters - ---------- - coordinates : tuple of float - Observation point coordinates in micrometers (μm). - dipole_coordinates : tuple of float - Dipole location coordinates in micrometers (μm). - dipole_moments : tuple of float - Dipole moment components (Am²). - - Returns - ------- - bz : float - The vertical component of the magnetic field (Bz) at the given observation point. - """ - coordinates_m = coordinates_micrometer_to_meter(coordinates) - dipole_coordinates_m = coordinates_micrometer_to_meter(dipole_coordinates) - return hm.dipole_magnetic( - coordinates_m, dipole_coordinates_m, dipole_moments, field="b_u" - ) - - -def dipole_bz_grid( - region, spacing, sensor_sample_distance, dipole_coordinates, dipole_moments -): - """ - Generate a grid of the vertical component of the magnetic field of a dipole. - - Parameters - ---------- - region : tuple of float - The spatial region for the grid in micrometers (μm), defined as - (x_min, x_max, y_min, y_max). - spacing : float - Grid spacing in micrometers (μm). - sensor_sample_distance : float - Distance of the sensor from the grid in micrometers (μm). - dipole_coordinates : tuple of float - Dipole location coordinates in micrometers (μm). - dipole_moments : tuple of float - Dipole moment components (Am²). - - Returns - ------- - data : xarray.DataArray - Gridded dataset containing the vertical component of the magnetic - field (Bz). - The dataset includes: - - "bz" : vertical magnetic field (nT) - - "x" and "y" coordinates with units in micrometers (μm) - """ - coordinates = vd.grid_coordinates( - region=region, # µm - spacing=spacing, # µm - extra_coords=sensor_sample_distance, - ) - bz = dipole_bz(coordinates, dipole_coordinates, dipole_moments) - data = vd.make_xarray_grid( - coordinates, bz, data_names=["bz"], dims=("y", "x"), extra_coords_names="z" - ) - data.x.attrs = {"units": "µm"} - data.y.attrs = {"units": "µm"} - data.bz.attrs = {"long_name": "vertical magnetic field", "units": "nT"} - return data.bz