Add magnetic field forward modelling of rectangular prisms

doc/api/index.rst

@@ -81,6 +81,8 @@ Magnetic fields:
 
  dipole_magnetic
  dipole_magnetic_component
+ prism_magnetic
+ prism_magnetic_component
 
 Layers and meshes:
 

doc/user_guide/forward_modelling/prism.rst

@@ -49,7 +49,7 @@ computation point:
 
 
 Gravitational fields
-^^^^^^^^^^^^^^^^^^^^
+--------------------
 
 The :func:`harmonica.prism_gravity` is able to compute the gravitational
 potential (``"potential"``), the acceleration components (``"g_e"``, ``"g_n"``,
@@ -139,8 +139,9 @@ Plot the results:
  plt.colorbar(tmp, ax=ax, label="Eotvos", orientation="horizontal", pad=0.08)
  plt.show()
 
+
 Passing multiple prisms
-^^^^^^^^^^^^^^^^^^^^^^^
+~~~~~~~~~~~~~~~~~~~~~~~
 
 We can compute the gravitational field of a set of prisms by passing a list of
 them, where each prism is defined as mentioned before, and then making
@@ -200,10 +201,123 @@ Lets plot this gravitational field:
  fig.show()
 
 
+Magnetic fields
+---------------
+
+The :func:`harmonica.prism_magnetic` and
+:func:`harmonica.prism_magnetic_component` functions allows us to calculate the
+magnetic fields generated by rectangular prisms on a set of observation points.
+Each rectangular prism is defined in the same way we did for the gravity
+forward modelling, although we now need to define a magnetization vector for
+each one of them.
+
+For example:
+
+.. jupyter-execute::
+
+ import numpy as np
+
+ prisms = [
+ [-5e3, -3e3, -5e3, -2e3, -10e3, -1e3],
+ [3e3, 4e3, 4e3, 5e3, -9e3, -1e3],
+ ]
+
+ magnetization = np.array([
+ [0.5, 0.5, -0.78989],
+ [-0.4, 0.3, 0.2],
+ ])
+
+Each row of the ``magnetization`` 2D vector corresponds to the easting,
+northing and upward components of the magnetization vector of each prism in
+``prisms``, provided in :math:`Am^2`.
+
+With the :func:`harmonica.prism_magnetic` function we can compute the three
+components of the magnetic field the prisms generates on any set of observation
+points.
+
+.. jupyter-execute::
+
+ region = (-10e3, 10e3, -10e3, 10e3)
+ shape = (51, 51)
+ height = 10
+ coordinates = vd.grid_coordinates(region, shape=shape, extra_coords=height)
+
+.. jupyter-execute::
+
+ b_e, b_n, b_u = hm.prism_magnetic(coordinates, prisms, magnetization)
+
+.. jupyter-execute::
+
+ fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 6))
+
+ for ax, mag_component, title in zip(axes, (b_e, b_n, b_u), ("Be", "Bn", "Bu")):
+ maxabs = vd.maxabs(mag_component)
+ tmp = ax.pcolormesh(
+ coordinates[0],
+ coordinates[1],
+ mag_component,
+ vmin=-maxabs,
+ vmax=maxabs,
+ cmap="RdBu_r",
+ )
+ ax.contour(
+ coordinates[0],
+ coordinates[1],
+ mag_component,
+ colors="k",
+ linewidths=0.5,
+ )
+ ax.set_title(title)
+ ax.set_aspect("equal")
+ plt.colorbar(
+ tmp,
+ ax=ax,
+ orientation="horizontal",
+ label="nT",
+ pad=0.08,
+ aspect=42,
+ shrink=0.8,
+ )
+
+ plt.show()
+
+
+Alternatively, we can use :func:`harmonica.prism_magnetic_component` to
+calculate only a single component of the magnetic field.
+
+.. important::
+
+ Using :func:`harmonica.prism_magnetic_component` to compute several magnetic
+ components is less efficient that using :func:`harmonica.prism_magnetic`.
+ Use the former only when a single component is needed.
+
+For example, we can calculate only the upward component of the magnetic field
+generated by these two prisms:
+
+.. jupyter-execute::
+
+ b_u = hm.prism_magnetic_component(
+ coordinates, prisms, magnetization, component="upward"
+ )
+
+.. jupyter-execute::
+
+ maxabs = vd.maxabs(b_u)
+
+ tmp = plt.pcolormesh(
+ coordinates[0], coordinates[1], b_u, vmin=-maxabs, vmax=maxabs, cmap="RdBu_r"
+ )
+ plt.contour(coordinates[0], coordinates[1], b_u, colors="k", linewidths=0.5)
+ plt.title("Bu")
+ plt.gca().set_aspect("equal")
+ plt.colorbar(tmp, label="nT", pad=0.03, aspect=42, shrink=0.8)
+ plt.show()
+
+
 .. _prism_layer:
 
 Prism layer
-^^^^^^^^^^^
+-----------
 
 One of the most common usage of prisms is to model geologic structures.
 Harmonica offers the possibility to define a layer of prisms through the
@@ -241,8 +355,6 @@ sample a trigonometric function for this simple example:
 
 .. jupyter-execute::
 
- import numpy as np
-
  wavelength = 24 * spacing
  surface = np.abs(np.sin(easting * 2 * np.pi / wavelength))
 

harmonica/__init__.py

@@ -11,8 +11,9 @@
 from ._equivalent_sources.spherical import EquivalentSourcesSph
 from ._forward.dipole import dipole_magnetic, dipole_magnetic_component
 from ._forward.point import point_gravity
-from ._forward.prism import prism_gravity
+from ._forward.prism_gravity import prism_gravity
 from ._forward.prism_layer import DatasetAccessorPrismLayer, prism_layer
+from ._forward.prism_magnetic import prism_magnetic, prism_magnetic_component
 from ._forward.tesseroid import tesseroid_gravity
 from ._forward.tesseroid_layer import DatasetAccessorTesseroidLayer, tesseroid_layer
 from ._gravity_corrections import bouguer_correction

harmonica/_forward/_tesseroid_variable_density.py

@@ -165,6 +165,7 @@ def _density_based_discretization(tesseroid, density):
  tesseroids : list
  List containing the boundaries of discretized tesseroids.
  """
+
  # Define normalized density
  def normalized_density(radius):
  return (density(radius) - density_min) / (density_max - density_min)

harmonica/_forward/prism_layer.py

@@ -14,7 +14,7 @@
 import xarray as xr
 
 from ..visualization import prism_to_pyvista
-from .prism import prism_gravity
+from .prism_gravity import prism_gravity
 
 
 def prism_layer(