Add magnetic field forward modelling of rectangular prisms (#369)

Add forward modelling functions that calculate the magnetic fields of
rectangular prisms. One function efficiently computes the three
components of the magnetic field, while the other one calculates just
one component. The functions make use of Choclo's forward modelling
functions, while providing a nice interface, sanity checks and
parallelization through Numba. Add examples and tests.


**Relevant issues/PRs:**

Related to #406

---------

Co-authored-by: Lu Li <[email protected]>

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^{-1}`.
+
+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(