Add a new Sphere class as a special case of ellipsoids (#602)

Add a new `Sphere` class as a special case of ellipsoids with `a == b ==
c`. Implement gravity and magnetic forward modelling for the sphere,
both inside and outside of the body. Extend tests to `Sphere` objects.
Use the forward modelling of the `Sphere` to compare the fields of
ellipsoids against it.

Author: santisoler

doc/api/index.rst

@@ -93,6 +93,7 @@ Ellipsoids:
 
     OblateEllipsoid
     ProlateEllipsoid
+    Sphere
     TriaxialEllipsoid
 
 Layers and meshes:

harmonica/__init__.py

@@ -21,6 +21,7 @@
 from ._forward.ellipsoids import (
     OblateEllipsoid,
     ProlateEllipsoid,
+    Sphere,
     TriaxialEllipsoid,
 )
 from ._forward.point import point_gravity

harmonica/_forward/ellipsoid_gravity.py

@@ -152,6 +152,22 @@ def _compute_gravity_ellipsoid(
     # Mask internal points
     internal = is_internal(x, y, z, a, b, c)
 
+    if a == b == c:
+        # Fallback to sphere equations which are simpler
+        factor = -4 / 3 * np.pi * a**3 * GRAVITATIONAL_CONST * density
+        gx, gy, gz = tuple(np.zeros_like(x) for _ in range(3))
+
+        gx[internal] = factor * x[internal] / a**3
+        gy[internal] = factor * y[internal] / a**3
+        gz[internal] = factor * z[internal] / a**3
+
+        r = np.sqrt(x[~internal] ** 2 + y[~internal] ** 2 + z[~internal] ** 2)
+        gx[~internal] = factor * x[~internal] / r**3
+        gy[~internal] = factor * y[~internal] / r**3
+        gz[~internal] = factor * z[~internal] / r**3
+
+        return gx, gy, gz
+
     # Compute lambda on all observation points:
     # Make it zero on internal points, calculate it for external points.
     lambda_ = np.zeros_like(x, dtype=np.float64)

harmonica/_forward/ellipsoid_magnetics.py

@@ -365,6 +365,8 @@ def get_demagnetization_tensor_internal(a: float, b: float, c: float):
         n_diagonal = _demag_tensor_prolate_internal(a, b)
     elif a < b and b == c:
         n_diagonal = _demag_tensor_oblate_internal(a, b)
+    elif a == b == c:
+        n_diagonal = 1 / 3 * np.ones(3)
     else:
         msg = "Could not determine ellipsoid type for values given."
         raise ValueError(msg)
@@ -528,23 +530,38 @@ def get_demagnetization_tensor_external(
     # Allocate array for all demagnetization tensors (one for each observation point)
     demag_tensors = np.empty((x.size, 3, 3), dtype=np.float64)
 
-    # Calculate lambda and other quantities needed to build the tensors
-    lambda_ = calculate_lambda(x, y, z, a, b, c)
-    ellip_integrals = get_elliptical_integrals(a, b, c, lambda_)
-    deriv_ellip_integrals = get_derivatives_of_elliptical_integrals(a, b, c, lambda_)
-    derivs_lmbda = _spatial_deriv_lambda(x, y, z, a, b, c, lambda_)
-
     coords = (x, y, z)
-    for i, j in itertools.product(range(3), range(3)):
-        if i == j:
-            demag_tensors[:, i, i] = ((a * b * c) / 2) * (
-                derivs_lmbda[i] * deriv_ellip_integrals[i] * coords[i]
-                + ellip_integrals[i]
-            )
-        else:
-            demag_tensors[:, i, j] = ((a * b * c) / 2) * (
-                derivs_lmbda[i] * deriv_ellip_integrals[j] * coords[j]
-            )
+
+    if a == b == c:
+        r = np.sqrt(x**2 + y**2 + z**2)
+        r5, r3 = r**5, r**3
+        factor = -(a**3 / 3)
+
+        for i, j in itertools.product(range(3), range(3)):
+            if i == j:
+                demag_tensors[:, i, i] = factor * (3 * coords[i] ** 2 / r5 - 1 / r3)
+            else:
+                demag_tensors[:, i, j] = factor * (3 * coords[i] * coords[j] / r5)
+
+    else:
+        # Calculate lambda and other quantities needed to build the tensors
+        lambda_ = calculate_lambda(x, y, z, a, b, c)
+        ellip_integrals = get_elliptical_integrals(a, b, c, lambda_)
+        deriv_ellip_integrals = get_derivatives_of_elliptical_integrals(
+            a, b, c, lambda_
+        )
+        derivs_lmbda = _spatial_deriv_lambda(x, y, z, a, b, c, lambda_)
+
+        for i, j in itertools.product(range(3), range(3)):
+            if i == j:
+                demag_tensors[:, i, i] = ((a * b * c) / 2) * (
+                    derivs_lmbda[i] * deriv_ellip_integrals[i] * coords[i]
+                    + ellip_integrals[i]
+                )
+            else:
+                demag_tensors[:, i, j] = ((a * b * c) / 2) * (
+                    derivs_lmbda[i] * deriv_ellip_integrals[j] * coords[j]
+                )
 
     return demag_tensors
 

harmonica/_forward/ellipsoids.py

@@ -26,10 +26,6 @@ class BaseEllipsoid:
         This class is not meant to be instantiated.
     """
 
-    yaw: float
-    pitch: float
-    roll: float
-
     @property
     def density(self) -> float | None:
         """Density of the ellipsoid in :math:`kg/m^3`."""
@@ -482,3 +478,111 @@ def _check_semiaxes_lenghts(self, a, b):
                 f"expected a < b (= c ) but got a = {a}, b = {b}"
             )
             raise ValueError(msg)
+
+
+class Sphere(BaseEllipsoid):
+    """
+    Homogeneous sphere.
+
+    Represent a sphere as a particular type of ellipsoid where ``a == b == c``.
+
+    Parameters
+    ----------
+    a : float
+        Radius or semiaxes lengths in meters.
+    center : tuple of float
+        Coordinates of the center of the sphere in the following order: `easting`,
+        `northing`, `upward`. All must be in meters.
+    density : float or None, optional
+        Density of the sphere in :math:`kg/m^3`.
+    susceptibility : float, (3, 3) array or None, optional
+        Magnetic susceptibility of the sphere in SI units.
+        A single float represents isotropic susceptibility in the body.
+        A (3, 3) array represents the susceptibility tensor to account for anisotropy.
+        If None, zero susceptibility will be assigned to the sphere.
+    remanent_mag : (3,) array or None, optional
+        Remanent magnetization vector of the sphere in A/m units. Its components
+        are defined in the easting-northing-upward coordinate system and should be
+        passed in that order.
+        If None, no remanent magnetization will be assigned to the sphere.
+
+    Notes
+    -----
+    All semiaxes (``a``, ``b`` and ``c``) are equal to each other.
+    All rotation angles (``yaw``, ``pitch`` and ``roll``) are equal to zero for the
+    sphere, since sphere is invariant to rotations.
+    """
+
+    def __init__(
+        self,
+        a: float,
+        center: tuple[float, float, float],
+        *,
+        density: float | None = None,
+        susceptibility: float | npt.NDArray | None = None,
+        remanent_mag: npt.NDArray | None = None,
+    ):
+        self._check_positive_semiaxis(a, "a")
+        self._a = a
+        self.center = center
+
+        # Physical properties of the sphere
+        self.density = density
+        self.susceptibility = susceptibility
+        self.remanent_mag = remanent_mag
+
+    @property
+    def a(self) -> float:
+        """Length of the first semiaxis."""
+        return self._a
+
+    @a.setter
+    def a(self, value: float) -> None:
+        self._check_positive_semiaxis(value, "a")
+        self._a = value
+
+    @property
+    def b(self) -> float:
+        """Length of the second semiaxis."""
+        return self._a
+
+    @property
+    def c(self) -> float:
+        """Length of the third semiaxis."""
+        return self._a
+
+    @property
+    def radius(self) -> float:
+        """Sphere radius."""
+        return self._a
+
+    @property
+    def yaw(self):
+        """Yaw angle, equal to zero."""
+        return 0.0
+
+    @property
+    def pitch(self):
+        """Pitch angle, equal to zero."""
+        return 0.0
+
+    @property
+    def roll(self):
+        """Roll angle, equal to zero."""
+        return 0.0
+
+    @property
+    def rotation_matrix(self) -> npt.NDArray:
+        """
+        Create a rotation matrix for the sphere.
+
+        .. important:
+
+            The rotation matrix for the sphere is always the identity matrix.
+
+        Returns
+        -------
+        rotation_matrix : (3, 3) array
+            Identity matrix.
+        """
+        return np.eye(3, dtype=np.float64)