diff --git a/doc/api/index.rst b/doc/api/index.rst index 481d6fa4c..65328119f 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -86,16 +86,12 @@ Magnetic fields: ellipsoid_magnetic prism_magnetic -Ellipsoids: +Geometric bodies: .. autosummary:: :toctree: generated/ - create_ellipsoid - OblateEllipsoid - ProlateEllipsoid - Sphere - TriaxialEllipsoid + Ellipsoid Layers and meshes: @@ -174,4 +170,3 @@ Type hints :toctree: generated/ typing.Coordinates - typing.Ellipsoid diff --git a/harmonica/__init__.py b/harmonica/__init__.py index 72b10625a..3acc064d8 100644 --- a/harmonica/__init__.py +++ b/harmonica/__init__.py @@ -16,15 +16,7 @@ from ._equivalent_sources.spherical import EquivalentSourcesSph from ._euler_deconvolution import EulerDeconvolution from ._forward.dipole import dipole_magnetic -from ._forward.ellipsoids import ( - OblateEllipsoid, - ProlateEllipsoid, - Sphere, - TriaxialEllipsoid, - create_ellipsoid, - ellipsoid_gravity, - ellipsoid_magnetic, -) +from ._forward.ellipsoids import Ellipsoid, ellipsoid_gravity, ellipsoid_magnetic from ._forward.point import point_gravity from ._forward.prism_gravity import prism_gravity from ._forward.prism_layer import DatasetAccessorPrismLayer, prism_layer diff --git a/harmonica/_forward/ellipsoids/__init__.py b/harmonica/_forward/ellipsoids/__init__.py index 4f916405b..9b3a6885c 100644 --- a/harmonica/_forward/ellipsoids/__init__.py +++ b/harmonica/_forward/ellipsoids/__init__.py @@ -8,12 +8,6 @@ Forward modelling of ellipsoids. """ -from .ellipsoids import ( - OblateEllipsoid, - ProlateEllipsoid, - Sphere, - TriaxialEllipsoid, - create_ellipsoid, -) +from .ellipsoids import Ellipsoid from .gravity import ellipsoid_gravity from .magnetic import ellipsoid_magnetic diff --git a/harmonica/_forward/ellipsoids/ellipsoids.py b/harmonica/_forward/ellipsoids/ellipsoids.py index 4b0b4ec43..ea0ef79da 100644 --- a/harmonica/_forward/ellipsoids/ellipsoids.py +++ b/harmonica/_forward/ellipsoids/ellipsoids.py @@ -22,31 +22,14 @@ pyvista = None -def create_ellipsoid( - a: float, - b: float, - c: float, - *, - center: tuple[float, float, float], - yaw: float = 0.0, - pitch: float = 0.0, - roll: float = 0.0, - **kwargs, -) -> "BaseEllipsoid": +class Ellipsoid: """ - Create an ellipsoid given its three semiaxes lenghts. - - This function returns an ellipsoid object of the correct type based on the values of - the three semiaxes lenghts. Semiaxes lengths can be passed in any order. + Ellipsoidal body with arbitrary rotation. Parameters ---------- - a, b, c: float - Semiaxes lenghts in meters. They can be passed in any particular order. They - must all be positive values. - center : tuple of float - Coordinates of the center of the ellipsoid in the following order: `easting`, - `northing`, `upward`. + a, b, c : float + Semi-axis lengths of the ellipsoid in meters. yaw : float, optional Rotation angle about the upward axis, in degrees. pitch : float, optional @@ -56,107 +39,91 @@ def create_ellipsoid( roll : float, optional Rotation angle about the easting axis (after yaw and pitch rotation), in degrees. - **kwargs : dict - Extra arguments passed to the ellipsoid classes. They can be physical properties - like ``density``, ``susceptibility``, and ``remanent_mag``. - - Returns - ------- - ellipsoid : harmonica.typing.Ellipsoid - Instance of one of the ellipsoid classes: :class:`harmonica.OblateEllipsoid`, - :class:`harmonica.ProlateEllipsoid`, :class:`harmonica.Sphere`, or - :class:`harmonica.TriaxialEllipsoid`. - - Examples - -------- - Define a prolate, oblate, triaxial or spherical ellipsoid by passing its semiaxes in - any order: - - >>> ellipsoid = create_ellipsoid(1, 1, 2, center=(1., 2., 3.)) - >>> print(type(ellipsoid).__name__) - ProlateEllipsoid - >>> ellipsoid.a - 2 - >>> ellipsoid.b - 1 - >>> ellipsoid.c - 1 - - >>> ellipsoid = create_ellipsoid(1, 2, 2, center=(1., 2., 3.)) - >>> print(type(ellipsoid).__name__) - OblateEllipsoid - >>> ellipsoid.a - 1 - >>> ellipsoid.b - 2 - >>> ellipsoid.c - 2 - - >>> ellipsoid = create_ellipsoid(1, 2, 3, center=(1., 2., 3.)) - >>> print(type(ellipsoid).__name__) - TriaxialEllipsoid - >>> ellipsoid.a - 3 - >>> ellipsoid.b - 2 - >>> ellipsoid.c - 1 - - We can specify rotation angles while creating the ellipsoid: - - >>> ellipsoid = create_ellipsoid(1, 1, 2, yaw=85, pitch=32, center=(1., 2., 3.)) - >>> ellipsoid.yaw - 85 - >>> ellipsoid.pitch - 32 - - The function can also take physical properties: - - >>> density, susceptibility = -200.0, 0.4 - >>> ellipsoid = create_ellipsoid( - ... 1, 1, 2, center=(1., 2., 3.), - ... density=density, - ... susceptibility=susceptibility, - ... ) + center : tuple of float, optional + Coordinates of the center of the ellipsoid in the following order: `easting`, + `northing`, `upward`, in meters. + density : float or None, optional + Density of the ellipsoid in :math:`kg/m^3`. + susceptibility : float, (3, 3) array or None, optional + Magnetic susceptibility of the ellipsoid in SI units. + A single float represents isotropic susceptibility in the body. + A (3, 3) array represents the susceptibility tensor to account for anisotropy + (defined in the local coordinate system of the ellipsoid). + If None, zero susceptibility will be assigned to the ellipsoid. + remanent_mag : (3,) array or None, optional + Remanent magnetization vector of the ellipsoid 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 ellipsoid. + + Notes + ----- + The three semi-axes ``a``, ``b``, and ``c`` are defined parallel to the ``easting``, + ``northing`` and ``upward`` directions, respectively, before applying any rotation. + + Rotations directed by ``yaw`` and ``roll`` are applied using the right-hand rule + across their respective axes. Pitch rotations are carried out in the opposite + direction, so a positive ``pitch`` *lifts* the side of the ellipsoid pointing in the + easting direction. """ - # Sanity checks - for semiaxis, value in zip(("a", "b", "c"), (a, b, c), strict=True): - if value <= 0: - msg = ( - f"Invalid value of '{semiaxis}' equal to '{value}'. " - "It must be positive." - ) - raise ValueError(msg) - if a == b == c: - return Sphere(a, center=center, **kwargs) + def __init__( + self, + a: float, + b: float, + c: float, + yaw: float = 0.0, + pitch: float = 0.0, + roll: float = 0.0, + center: tuple[float, float, float] = (0.0, 0.0, 0.0), + *, + density: float | None = None, + susceptibility: float | npt.NDArray | None = None, + remanent_mag: npt.NDArray | None = None, + ): + # Sanity checks on semiaxes + for value, semiaxis in zip((a, b, c), ("a", "b", "c"), strict=True): + self._check_positive_semiaxis(value, semiaxis) + self._a, self._b, self._c = a, b, c - # Sort semiaxes so a >= b >= c - c, b, a = sorted((a, b, c)) + # Angles and center + self.yaw, self.pitch, self.roll = yaw, pitch, roll + self.center = center - if a == b: - a = c # set `a` as the smallest semiaxis (`c`) - ellipsoid = OblateEllipsoid(a, b, yaw=yaw, pitch=pitch, center=center, **kwargs) - elif b == c: - ellipsoid = ProlateEllipsoid( - a, b, yaw=yaw, pitch=pitch, center=center, **kwargs - ) - else: - ellipsoid = TriaxialEllipsoid( - a, b, c, yaw=yaw, pitch=pitch, roll=roll, center=center, **kwargs - ) + # Physical properties of the ellipsoid + self.density = density + self.susceptibility = susceptibility + self.remanent_mag = remanent_mag - return ellipsoid + @property + def a(self) -> float: + """Length of the first semiaxis in meters.""" + return self._a + @a.setter + def a(self, value: float) -> None: + self._check_positive_semiaxis(value, "a") + self._a = value -class BaseEllipsoid: - """ - Base class for ellipsoids. + @property + def b(self) -> float: + """Length of the second semiaxis in meters.""" + return self._b + + @b.setter + def b(self, value: float): + self._check_positive_semiaxis(value, "b") + self._b = value - .. important:: + @property + def c(self) -> float: + """Length of the third semiaxis in meters.""" + return self._c - This class is not meant to be instantiated. - """ + @c.setter + def c(self, value: float): + self._check_positive_semiaxis(value, "c") + self._c = value @property def density(self) -> float | None: @@ -321,459 +288,3 @@ def to_pyvista(self, **kwargs): ellipsoid.rotate(rotation=self.rotation_matrix, inplace=True) ellipsoid.translate(self.center, inplace=True) return ellipsoid - - -class TriaxialEllipsoid(BaseEllipsoid): - """ - Triaxial ellipsoid with arbitrary orientation. - - Define a triaxial ellipsoid whose semi-axes lengths are ``a > b > c``. - - Parameters - ---------- - a, b, c : float - Semi-axis lengths of the ellipsoid. Must satisfy the condition: ``a > b > c``. - yaw : float - Rotation angle about the upward axis, in degrees. - pitch : float - Rotation angle about the northing axis (after yaw rotation), in degrees. - A positive pitch angle *lifts* the side of the ellipsoid pointing in easting - direction. - roll : float - Rotation angle about the easting axis (after yaw and pitch rotation), in - degrees. - center : tuple of float - Coordinates of the center of the ellipsoid in the following order: `easting`, - `northing`, `upward`. - density : float or None, optional - Density of the ellipsoid in :math:`kg/m^3`. - susceptibility : float, (3, 3) array or None, optional - Magnetic susceptibility of the ellipsoid 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 ellipsoid. - remanent_mag : (3,) array or None, optional - Remanent magnetization vector of the ellipsoid 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 ellipsoid. - - Notes - ----- - The three semi-axes ``a``, ``b``, and ``c`` are defined parallel to the ``easting``, - ``northing`` and ``upward`` directions, respectively, before applying any rotation. - - Rotations directed by ``yaw`` and ``roll`` are applied using the right-hand rule - across their respective axes. Pitch rotations are carried out in the opposite - direction, so a positive ``pitch`` *lifts* the side of the ellipsoid pointing in the - easting direction. - - """ - - def __init__( - self, - a: float, - b: float, - c: float, - yaw: float, - pitch: float, - roll: float, - center: tuple[float, float, float], - *, - density: float | None = None, - susceptibility: float | npt.NDArray | None = None, - remanent_mag: npt.NDArray | None = None, - ): - # Sanity checks on semiaxes - for value, semiaxis in zip((a, b, c), ("a", "b", "c"), strict=True): - self._check_positive_semiaxis(value, semiaxis) - self._check_semiaxes_lenghts(a, b, c) - self._a, self._b, self._c = a, b, c - - # Angles and center - self.yaw, self.pitch, self.roll = yaw, pitch, roll - self.center = center - - # Physical properties of the ellipsoid - 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._check_semiaxes_lenghts(value, self.b, self.c) - self._a = value - - @property - def b(self) -> float: - """Length of the second semiaxis.""" - return self._b - - @b.setter - def b(self, value: float): - self._check_positive_semiaxis(value, "b") - self._check_semiaxes_lenghts(self.a, value, self.c) - self._b = value - - @property - def c(self) -> float: - """Length of the third semiaxis.""" - return self._c - - @c.setter - def c(self, value: float): - self._check_positive_semiaxis(value, "c") - self._check_semiaxes_lenghts(self.a, self.b, value) - self._c = value - - def _check_semiaxes_lenghts(self, a, b, c): - if not (a > b > c): - msg = ( - "Invalid ellipsoid axis lengths for triaxial ellipsoid: " - f"expected a > b > c but got a = {a}, b = {b}, c = {c}" - ) - raise ValueError(msg) - - -class ProlateEllipsoid(BaseEllipsoid): - """ - Prolate ellipsoid with arbitrary orientation. - - Define a prolate ellipsoid whose semi-axes lengths are ``a > b = c``. - - Parameters - ---------- - a, b : float - Semi-axis lengths of the ellipsoid. Must satisfy the condition: ``a > b = c``. - yaw : float - Rotation angle about the upward axis, in degrees. - pitch : float - Rotation angle about the northing axis (after yaw rotation), in degrees. - A positive pitch angle *lifts* the side of the ellipsoid pointing in easting - direction. - center : tuple of float - Coordinates of the center of the ellipsoid in the following order: `easting`, - `northing`, `upward`. - density : float or None, optional - Density of the ellipsoid in :math:`kg/m^3`. - susceptibility : float or None, optional - Magnetic susceptibility of the ellipsoid 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 ellipsoid. - remanent_mag : (3,) array or None, optional - Remanent magnetization vector of the ellipsoid 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 ellipsoid. - - Notes - ----- - The three semi-axes ``a``, ``b``, and ``c`` are defined parallel to the ``easting``, - ``northing`` and ``upward`` directions, respectively, before applying any rotation. - - Rotations directed by ``yaw`` are applied using the right-hand rule across the - upward axis. Pitch rotations are carried out in the opposite direction, so - a positive ``pitch`` *lifts* the side of the ellipsoid pointing in the easting - direction. - - Roll rotations are not enabled in the prolate ellipsoid, since they don't have any - effect due to symmetry. Hence, ``roll`` is always equal to zero for the - :class:`harmonica.ProlateEllipsoid`. - """ - - def __init__( - self, - a: float, - b: float, - yaw: float, - pitch: float, - center: tuple[float, float, float], - *, - density: float | None = None, - susceptibility: float | npt.NDArray | None = None, - remanent_mag: npt.NDArray | None = None, - ): - # Sanity checks on semiaxes - for value, semiaxis in zip((a, b), ("a", "b"), strict=True): - self._check_positive_semiaxis(value, semiaxis) - self._check_semiaxes_lenghts(a, b) - self._a, self._b = a, b - - # Angles and center - self.yaw, self.pitch = yaw, pitch - self.center = center - - # Physical properties of the ellipsoid - 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): - self._check_positive_semiaxis(value, "a") - self._check_semiaxes_lenghts(value, self.b) - self._a = value - - @property - def b(self) -> float: - """Length of the second semiaxis.""" - return self._b - - @b.setter - def b(self, value: float): - self._check_positive_semiaxis(value, "b") - self._check_semiaxes_lenghts(self.a, value) - self._b = value - - @property - def c(self): - """Length of the third semiaxis, equal to ``b`` by definition.""" - return self.b - - @property - def roll(self): - """Roll angle, equal to zero.""" - return 0.0 - - def _check_semiaxes_lenghts(self, a, b): - if not (a > b): - msg = ( - "Invalid ellipsoid axis lengths for prolate ellipsoid: " - f"expected a > b (= c ) but got a = {a}, b = {b}" - ) - raise ValueError(msg) - - -class OblateEllipsoid(BaseEllipsoid): - """ - Oblate ellipsoid with arbitrary orientation. - - Define a prolate ellipsoid whose semi-axes lengths are ``a < b = c``. - - Parameters - ---------- - a, b : float - Semi-axis lengths of the ellipsoid. Must satisfy the condition: ``a < b = c``. - yaw : float - Rotation angle about the upward axis, in degrees. - pitch : float - Rotation angle about the northing axis (after yaw rotation), in degrees. - A positive pitch angle *lifts* the side of the ellipsoid pointing in easting - direction. - center : tuple of float - Coordinates of the center of the ellipsoid in the following order: `easting`, - `northing`, `upward`. - density : float or None, optional - Density of the ellipsoid in :math:`kg/m^3`. - susceptibility : float or None, optional - Magnetic susceptibility of the ellipsoid 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 ellipsoid. - remanent_mag : (3,) array or None, optional - Remanent magnetization vector of the ellipsoid 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 ellipsoid. - - Notes - ----- - The three semi-axes ``a``, ``b``, and ``c`` are defined parallel to the ``easting``, - ``northing`` and ``upward`` directions, respectively, before applying any rotation. - - Rotations directed by ``yaw`` are applied using the right-hand rule across the - upward axis. Pitch rotations are carried out in the opposite direction, so - a positive ``pitch`` *lifts* the side of the ellipsoid pointing in the easting - direction. - - Roll rotations are not enabled in the prolate ellipsoid, since they don't have any - effect due to symmetry. Hence, ``roll`` is always equal to zero for the - :class:`harmonica.OblateEllipsoid`. - """ - - def __init__( - self, - a: float, - b: float, - yaw: float, - pitch: float, - center: tuple[float, float, float], - *, - density: float | None = None, - susceptibility: float | npt.NDArray | None = None, - remanent_mag: npt.NDArray | None = None, - ): - # Sanity checks on semiaxes - for value, semiaxis in zip((a, b), ("a", "b"), strict=True): - self._check_positive_semiaxis(value, semiaxis) - self._check_semiaxes_lenghts(a, b) - self._a, self._b = a, b - - # Angles and center - self.yaw, self.pitch = yaw, pitch - self.center = center - - # Physical properties of the ellipsoid - 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): - self._check_positive_semiaxis(value, "a") - self._check_semiaxes_lenghts(value, self.b) - self._a = value - - @property - def b(self) -> float: - """Length of the second semiaxis.""" - return self._b - - @b.setter - def b(self, value: float): - self._check_positive_semiaxis(value, "b") - self._check_semiaxes_lenghts(self.a, value) - self._b = value - - @property - def c(self): - """Length of the third semiaxis, equal to ``b`` by definition.""" - return self.b - - @property - def roll(self): - """Roll angle, equal to zero.""" - return 0.0 - - def _check_semiaxes_lenghts(self, a, b): - if not (a < b): - msg = ( - "Invalid ellipsoid axis lengths for oblate ellipsoid: " - 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) diff --git a/harmonica/_forward/ellipsoids/gravity.py b/harmonica/_forward/ellipsoids/gravity.py index c1ecfeeb4..b3baaa986 100644 --- a/harmonica/_forward/ellipsoids/gravity.py +++ b/harmonica/_forward/ellipsoids/gravity.py @@ -16,10 +16,13 @@ from choclo.constants import GRAVITATIONAL_CONST from ...errors import NoPhysicalPropertyWarning -from ...typing import Coordinates, Ellipsoid +from ...typing import Coordinates +from .ellipsoids import Ellipsoid from .utils import ( calculate_lambda, + check_semiaxes_sorted, get_elliptical_integrals, + get_semiaxes_rotation_matrix, is_almost_a_sphere, is_internal, ) @@ -53,10 +56,9 @@ def ellipsoid_gravity( List of arrays containing the ``easting``, ``northing`` and ``upward`` coordinates of the computation points defined on a Cartesian coordinate system. All coordinates should be in meters. - ellipsoid : ellipsoid or list of ellipsoids - Ellipsoidal body represented by an instance of - :class:`harmonica.TriaxialEllipsoid`, :class:`harmonica.ProlateEllipsoid`, - or :class:`harmonica.OblateEllipsoid`, or a list of them. + ellipsoid : harmonica.Ellipsoid or list of harmonica.Ellipsoid + Ellipsoidal body represented by an instance of :class:`harmonica.Ellipsoid` or + a list of them. Returns ------- @@ -97,19 +99,25 @@ def ellipsoid_gravity( origin_e, origin_n, origin_u = ellipsoid.center coords_shifted = (easting - origin_e, northing - origin_n, upward - origin_u) - # Create rotation matrix - r = ellipsoid.rotation_matrix + # Sort the semiaxes (a >= b >= c) + a, b, c = sorted((ellipsoid.a, ellipsoid.b, ellipsoid.c), reverse=True) + + # Get rotation matrix to produce sorted semiaxes in local coordinate system + semiaxes_rotation_matrix = get_semiaxes_rotation_matrix(ellipsoid) + + # Combine the two rotation matrices + rotation = semiaxes_rotation_matrix.T @ ellipsoid.rotation_matrix.T # Rotate observation points - x, y, z = r.T @ np.vstack(coords_shifted) + x, y, z = rotation @ np.vstack(coords_shifted) # Calculate gravity components on local coordinate system gravity_ellipsoid = _compute_gravity_ellipsoid( - x, y, z, ellipsoid.a, ellipsoid.b, ellipsoid.c, ellipsoid.density + x, y, z, a, b, c, ellipsoid.density ) # project onto upward unit vector, axis U - ge_i, gn_i, gu_i = r @ np.vstack(gravity_ellipsoid) + ge_i, gn_i, gu_i = rotation.T @ np.vstack(gravity_ellipsoid) # sum contributions from each ellipsoid ge += ge_i @@ -140,20 +148,21 @@ def _compute_gravity_ellipsoid( Parameters ---------- - x, y, z : arrays + x, y, z : (n,) array Observation coordinates in the local ellipsoid reference frame. - a, b, c : floats - Semiaxis lengths of the ellipsoid. Must conform to the constraints of - the chosen ellipsoid type. + a, b, c : float + Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. density : float Density of the ellipsoidal body in kg/m³. Returns ------- - gx, gy, gz : arrays + gx, gy, gz : (n,) array Gravity acceleration components in the local coordinate system for the ellipsoid. Accelerations are given in SI units (m/s^2). """ + check_semiaxes_sorted(a, b, c) + # Mask internal points internal = is_internal(x, y, z, a, b, c) diff --git a/harmonica/_forward/ellipsoids/magnetic.py b/harmonica/_forward/ellipsoids/magnetic.py index 1075f755d..9159ff6f0 100644 --- a/harmonica/_forward/ellipsoids/magnetic.py +++ b/harmonica/_forward/ellipsoids/magnetic.py @@ -20,12 +20,17 @@ from ..._utils import magnetic_angles_to_vec from ...errors import NoPhysicalPropertyWarning -from ...typing import Coordinates, Ellipsoid +from ...typing import Coordinates +from .ellipsoids import Ellipsoid from .utils import ( calculate_lambda, + check_semiaxes_sorted, get_derivatives_of_elliptical_integrals, get_elliptical_integrals, + get_semiaxes_rotation_matrix, is_almost_a_sphere, + is_almost_oblate, + is_almost_prolate, is_internal, ) @@ -51,10 +56,9 @@ def ellipsoid_magnetic( List of arrays containing the ``easting``, ``northing`` and ``upward`` coordinates of the computation points defined on a Cartesian coordinate system. All coordinates should be in meters. - ellipsoid : ellipsoid or list of ellipsoids - Ellipsoidal body represented by an instance of - :class:`harmonica.TriaxialEllipsoid`, :class:`harmonica.ProlateEllipsoid`, or - :class:`harmonica.OblateEllipsoid`, or a list of them. + ellipsoid : harmonica.Ellipsoid or list of harmonica.Ellipsoid + Ellipsoidal body represented by an instance of :class:`harmonica.Ellipsoid` or + a list of them. external_field : tuple The uniform magnetic field (B) as and array with values of (magnitude, inclination, declination). The magnitude should be in nT, @@ -132,7 +136,8 @@ def _single_ellipsoid_magnetic( ellipsoid : harmonica.typing.Ellipsoid Ellipsoid object for which the magnetic field will be computed. susceptibility : float, (3, 3) array or None - Susceptibility scalar or tensor of the ellipsoid. + Susceptibility scalar or tensor of the ellipsoid + (defined in the local coordinate system of the ellipsoid). remanent_mag : (3,) array or None Remanent magnetization vector of the ellipsoid in SI units. The components of the vector should be in the easting-northing-upward coordinate @@ -153,28 +158,34 @@ def _single_ellipsoid_magnetic( easting, northing, upward = coordinates coords_shifted = (easting - origin_e, northing - origin_n, upward - origin_u) + # Sort the semiaxes (a >= b >= c) + a, b, c = sorted((ellipsoid.a, ellipsoid.b, ellipsoid.c), reverse=True) + + # Get rotation matrix to produce sorted semiaxes in local coordinate system + semiaxes_rotation_matrix = get_semiaxes_rotation_matrix(ellipsoid) + + # Combine the two rotation matrices + rotation = semiaxes_rotation_matrix.T @ ellipsoid.rotation_matrix.T + # Rotate observation points - r_matrix = ellipsoid.rotation_matrix - x, y, z = r_matrix.T @ np.vstack(coords_shifted) + x, y, z = rotation @ np.vstack(coords_shifted) # Build internal demagnetization tensor - n_tensor_internal = get_demagnetization_tensor_internal( - ellipsoid.a, ellipsoid.b, ellipsoid.c - ) + n_tensor_internal = get_demagnetization_tensor_internal(a, b, c) # Cast susceptibility and remanent magnetization into matrix and array, respectively susceptibility = cast_susceptibility(susceptibility) remanent_mag = cast_remanent_magnetization(remanent_mag) # Rotate the external field and the remanent magnetization - h0_field_rotated = r_matrix.T @ h0_field - remnant_mag_rotated = r_matrix.T @ remanent_mag + h0_field_rotated = rotation @ h0_field + remnant_mag_rotated = rotation @ remanent_mag # Get magnetization of the ellipsoid magnetization = get_magnetisation( - ellipsoid.a, - ellipsoid.b, - ellipsoid.c, + a, + b, + c, susceptibility, h0_field_rotated, remnant_mag_rotated, @@ -187,19 +198,19 @@ def _single_ellipsoid_magnetic( b_field = np.zeros((easting.size, 3), dtype=np.float64) # Mask internal observation points - internal = is_internal(x, y, z, ellipsoid.a, ellipsoid.b, ellipsoid.c) + internal = is_internal(x, y, z, a, b, c) # Compute b_field on internal points b_field[internal, :] = mu_0 * (-n_tensor_internal @ magnetization + magnetization) # Compute b_field on external points n_tensors = get_demagnetization_tensor_external( - x[~internal], y[~internal], z[~internal], ellipsoid.a, ellipsoid.b, ellipsoid.c + x[~internal], y[~internal], z[~internal], a, b, c ) b_field[~internal, :] = mu_0 * (-n_tensors @ magnetization) # Rotate the b fields - be, bn, bu = r_matrix @ b_field.T + be, bn, bu = rotation.T @ b_field.T return be, bn, bu @@ -228,10 +239,10 @@ def get_magnetisation( Parameters ---------- - a, b, c : floats - Semi-axes lengths of the ellipsoid. + a, b, c : float + Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. susceptibility : (3, 3) array - Susceptibility tensor. + Susceptibility tensor defined in the local coordinate system of the ellipsoid. h0_field : (3,) array The rotated background field (in local coordinates). remnant_mag : (3,) array @@ -268,6 +279,8 @@ def get_magnetisation( \mathbf{H}(\mathbf{r}) = \mathbf{H}_0 - \mathbf{N}(\mathbf{r}) \mathbf{M}. """ + check_semiaxes_sorted(a, b, c) + if n_tensor is None: n_tensor = get_demagnetization_tensor_internal(a, b, c) eye = np.identity(3) @@ -288,7 +301,8 @@ def cast_susceptibility(susceptibility: float | npt.NDArray | None) -> npt.NDArr susceptibility : float, (3, 3) array or None Magnetic susceptibility value. A single float represents isotropic susceptibility in the body. - A (3, 3) array represents the susceptibility tensor to account for anisotropy. + A (3, 3) array represents the susceptibility tensor to account for anisotropy + (defined in the local coordinate system of the ellipsoid). If None, a susceptibility of zero is assumed. Returns @@ -339,8 +353,8 @@ def get_demagnetization_tensor_internal(a: float, b: float, c: float): Parameters ---------- - a, b, c : floats - Semi-axes lengths of the given ellipsoid. + a, b, c : float + Semi-axes lengths of the given ellipsoid sorted such as ``a >= b >= c``. Returns ------- @@ -360,14 +374,16 @@ def get_demagnetization_tensor_internal(a: float, b: float, c: float): where :math:`\mathbf{M}` is the magnetization vector of the ellipsoid. """ + check_semiaxes_sorted(a, b, c) + if is_almost_a_sphere(a, b, c): n_diagonal = 1 / 3 * np.ones(3) + elif is_almost_prolate(a, b, c): + n_diagonal = _demag_tensor_prolate_internal(a, b) + elif is_almost_oblate(a, b, c): + n_diagonal = _demag_tensor_oblate_internal(b, c) elif a > b > c: n_diagonal = _demag_tensor_triaxial_internal(a, b, c) - elif a > b and b == c: - n_diagonal = _demag_tensor_prolate_internal(a, b) - elif a < b and b == c: - n_diagonal = _demag_tensor_oblate_internal(a, b) else: msg = "Could not determine ellipsoid type for values given." raise ValueError(msg) @@ -382,14 +398,18 @@ def _demag_tensor_triaxial_internal(a: float, b: float, c: float): Parameters ---------- - a, b, c : floats - Semi-axes lengths of the triaxial ellipsoid (a ≥ b ≥ c). + a, b, c : float + Semi-axes lengths of the given ellipsoid sorted such as ``a > b > c``. Returns ------- - nxx, nyy, nzz : floats + nxx, nyy, nzz : float individual diagonal components of the x, y, z matrix. """ + if not a > b > c: + msg = f"Invalid semiaxes length (not a > b > c): a={a}, b={b}, c={c}." + raise ValueError(msg) + # Cache values of E(theta, k) and F(theta, k) so we compute them only once phi = np.arccos(c / a) k = (a**2 - b**2) / (a**2 - c**2) @@ -416,44 +436,44 @@ def _demag_tensor_prolate_internal(a: float, b: float): Parameters ---------- - a, b: floats - Semi-axes lengths of the prolate ellipsoid (a > b = c). + a, b : float + Semi-axes lengths of the prolate ellipsoid (``a > b = c``). Returns ------- - nxx, nyy, nzz : floats - individual diagonal components of the x, y, z matrix. + nxx, nyy, nzz : float + Diagonal components of the internal demagnetization tensor. """ - m = a / b - if not m > 1: - msg = f"Invalid aspect ratio for prolate ellipsoid: a={a}, b={b}, a/b={m}" + if not a > b: + msg = f"Invalid semiaxes for prolate ellipsoid (not a > b): a={a}, b={b}." raise ValueError(msg) + m = a / b sqrt = np.sqrt(m**2 - 1) nxx = (1 / (m**2 - 1)) * (((m / sqrt) * np.log(m + sqrt)) - 1) nyy = nzz = 0.5 * (1 - nxx) return nxx, nyy, nzz -def _demag_tensor_oblate_internal(a: float, b: float): +def _demag_tensor_oblate_internal(b: float, c: float): """ Calculate internal demagnetization factors for oblate case. Parameters ---------- - a, b: floats - Semi-axes lengths of the oblate ellipsoid (a < b = c). + b, c : float + Semi-axes lengths of the oblate ellipsoid (``a == b > c``). Returns ------- - nxx, nyy, nzz : floats - individual diagonal components of the x, y, z matrix. + nxx, nyy, nzz : float + Diagonal components of the internal demagnetization tensor. """ - m = a / b - if not 0 < m < 1: - msg = f"Invalid aspect ratio for oblate ellipsoid: a={a}, b={b}, a/b={m}" + if not b > c: + msg = f"Invalid semiaxes for oblate ellipsoid (not b > c): b={b}, c={c}." raise ValueError(msg) - nxx = 1 / (1 - m**2) * (1 - (m / np.sqrt(1 - m**2)) * np.arccos(m)) - nyy = nzz = 0.5 * (1 - nxx) + m = c / b + nzz = 1 / (1 - m**2) * (1 - (m / np.sqrt(1 - m**2)) * np.arccos(m)) + nxx = nyy = 0.5 * (1 - nzz) return nxx, nyy, nzz @@ -474,7 +494,7 @@ def get_demagnetization_tensor_external( ---------- x, y, z : (n,) array Coordinates of the observation points in the local coordinate system. - a, b, c : floats + a, b, c : float Semi-axes lengths of the given ellipsoid. Returns @@ -551,17 +571,17 @@ def get_demagnetization_tensor_external( deriv_ellip_integrals = get_derivatives_of_elliptical_integrals( a, b, c, lambda_ ) - derivs_lmbda = _spatial_deriv_lambda(x, y, z, a, b, c, lambda_) + derivs_lambda = _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] + derivs_lambda[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] + derivs_lambda[i] * deriv_ellip_integrals[j] * coords[j] ) return demag_tensors @@ -581,16 +601,16 @@ def _spatial_deriv_lambda( Parameters ---------- - x, y, z : floats or (n,) arrays + x, y, z : float or (n,) array Coordinates of the observation points in the local coordinate system. - a, b, c : floats + a, b, c : float Semi-axes lengths of the given ellipsoid. lambda_ : float or (n,) arrays The given lambda value for each point we are considering with this matrix. Returns ------- - derivatives : tuple of floats or tuple of (n,) arrays + derivatives : tuple of float or tuple of (n,) array The spatial derivatives of lambda for each observation point. """ diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index 9d29b3d56..2aa02509e 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -8,21 +8,104 @@ import numpy.typing as npt from scipy.special import ellipeinc, ellipkinc +from ..utils import get_rotation_matrix + # Relative tolerance for two ellipsoid semiaxes to be considered almost equal. # E.g.: two semiaxes a and b are considered almost equal if: # | a - b | < max(a, b) * SEMIAXES_RTOL SEMIAXES_RTOL = 1e-5 +def check_semiaxes_sorted(a, b, c): + """ + Check if ellipsoid's semiaxes are sorted such as ``a >= b >= c``. + + Parameters + ---------- + a, b, c : float + Semiaxes lenghts. + + Raises + ------ + ValueError : if ``not a >= b >= c``. + """ + if not (a >= b >= c): + msg = ( + f"Invalid semiaxes not properly sorted (a >= b >= c): a={a}, b={b}, c={c}." + ) + raise ValueError(msg) + + +def get_semiaxes_rotation_matrix(ellipsoid): + """ + Build extra rotation matrix to align semiaxes in decreasing order. + + Build a 90 degrees rotations matrix that goes from a local coordinate system where: + + - ``x`` points in the direction of ``a``, + - ``y`` points in the direction of ``b``, and + - ``z`` points in the direction of ``c``, + + where ``a >= b >= c``, to a *primed* local coordinate system where: + + - ``x'`` points in the direction of ``ellipsoid.a``, + - ``y'`` points in the direction of ``ellipsoid.b``, and + - ``z'`` points in the direction of ``ellipsoid.c``, + + and ``ellipsoid.a``, ``ellipsoid.b`` and ``ellipsoid.c`` have no particular order. + The ``a``, ``b``, ``c`` are defined as: + + .. code::python + + a, b, c = sorted((ellipsoid.a, ellipsoid.b, ellipsoid.c), reverse=True) + + Parameters + ---------- + ellipsoid : Ellipsoid + Ellipsoid object. + + Returns + ------- + rotation_matrix : (3, 3) np.ndarray + Rotation matrix. + + Notes + ----- + This matrix is not necessarily a permutation matrix, since it can contain -1. + But it acts as a permutation matrix that ensures that ``x``, ``y``, ``z`` form + a right-handed system. + """ + a, b, c = ellipsoid.a, ellipsoid.b, ellipsoid.c + if a >= b >= c: + return np.eye(3, dtype=int) + + if b >= a >= c: + yaw, pitch, roll = 90, 0, 0 + elif c >= b >= a: + yaw, pitch, roll = 0, 90, 0 + elif a >= c >= b: + yaw, pitch, roll = 0, 0, 90 + elif b >= c >= a: + yaw, pitch, roll = 90, 0, 90 + elif c >= a >= b: + yaw, pitch, roll = 90, 90, 0 + else: + msg = f"Invalid semiaxes: a={a}, b={b}, c={c}." + raise ValueError(msg) + + matrix = get_rotation_matrix(yaw, pitch, roll).astype(int) + return matrix + + def is_internal(x, y, z, a, b, c): """ Check if a given point(s) is internal or external to the ellipsoid. Parameters ---------- - x, y, z : (n,) arrays or floats + x, y, z : (n,) array or float Coordinates of the observation point(s) in the local coordinate system. - a, b, c : floats + a, b, c : float Ellipsoid's semiaxes lengths. Returns @@ -36,26 +119,36 @@ def is_almost_a_sphere(a: float, b: float, c: float) -> bool: """ Check if a given ellipsoid approximates a sphere. - Returns True if ellipsoid's semiaxes lenghts are close enough to each other to be + Returns True if ellipsoid's semiaxes lengths are close enough to each other to be approximated by a sphere. + .. important:: + + The semiaxes should be already sorted such as ``a >= b >= c``. + Parameters ---------- - a, b, c: float - Ellipsoid's semiaxes lenghts. + a, b, c : float + Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. Returns ------- bool """ + check_semiaxes_sorted(a, b, c) + # Exactly a sphere if a == b == c: return True - # Prolate or oblate that is almost a sphere + # Prolate that is almost a sphere if b == c and np.abs(a - b) < SEMIAXES_RTOL * max(a, b): return True + # Oblate that is almost a sphere + if a == b and np.abs(b - c) < SEMIAXES_RTOL * max(b, c): + return True + # Triaxial that is almost a sphere if ( # noqa: SIM103 a != b != c @@ -67,6 +160,72 @@ def is_almost_a_sphere(a: float, b: float, c: float) -> bool: return False +def is_almost_prolate(a: float, b: float, c: float) -> bool: + """ + Check if a given ellipsoid approximates a prolate ellipsoid. + + Returns True if ellipsoid's semimiddle and semiminor lengths are close enough to + each other to be approximated by a prolate ellipsoid. + + .. important:: + + The semiaxes should be already sorted such as ``a >= b >= c``. + + Parameters + ---------- + a, b, c : float + Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. + + Returns + ------- + bool + """ + check_semiaxes_sorted(a, b, c) + + # Exactly a prolate + if a > b == c: + return True + + # Triaxial that is almost a prolate + if np.abs(b - c) < SEMIAXES_RTOL * max(b, c): # noqa: SIM103 + return True + + return False + + +def is_almost_oblate(a: float, b: float, c: float) -> bool: + """ + Check if a given ellipsoid approximates an oblate ellipsoid. + + Returns True if ellipsoid's semimajor and semimiddle lengths are close enough to + each other to be approximated by a oblate ellipsoid. + + .. important:: + + The semiaxes should be already sorted such as ``a >= b >= c``. + + Parameters + ---------- + a, b, c : float + Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. + + Returns + ------- + bool + """ + check_semiaxes_sorted(a, b, c) + + # Exactly an oblate + if a == b > c: + return True + + # Triaxial that is almost an oblate + if np.abs(a - b) < SEMIAXES_RTOL * max(a, b): # noqa: SIM103 + return True + + return False + + def calculate_lambda(x, y, z, a, b, c): """ Get the lambda ellipsoidal coordinate for a given ellipsoid and observation points. @@ -77,19 +236,10 @@ def calculate_lambda(x, y, z, a, b, c): Parameters ---------- - x : float or array - X-coordinate(s) of the observation point(s) in the local coordinate - system. - y : float or array - Y-coordinate(s) of the observation point(s). - z : float or array - Z-coordinate(s) of the observation point(s). - a : float - Semi-major axis of the ellipsoid along the x-direction. - b : float - Semi-major axis of the ellipsoid along the y-direction. - c : float - Semi-major axis of the ellipsoid along the z-direction. + x, y, z : float or array + Coordinates of the observation point(s) in the local coordinate system. + a, b, c : float + Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. Returns ------- @@ -97,13 +247,21 @@ def calculate_lambda(x, y, z, a, b, c): The computed value(s) of the lambda parameter. """ - # Solve lambda for prolate and oblate ellipsoids + check_semiaxes_sorted(a, b, c) + + # Solve lambda for prolate (a > b == c ) if b == c: p0 = a**2 * b**2 - b**2 * x**2 - a**2 * (y**2 + z**2) p1 = a**2 + b**2 - x**2 - y**2 - z**2 lambda_ = 0.5 * (np.sqrt(p1**2 - 4 * p0) - p1) - # Solve lambda for triaxial ellipsoids + # Solve lambda for oblate ( a == b > c ) + elif a == b: + p0 = a**2 * c**2 - c**2 * (x**2 + y**2) - a**2 * z**2 + p1 = a**2 + c**2 - x**2 - y**2 - z**2 + lambda_ = 0.5 * (np.sqrt(p1**2 - 4 * p0) - p1) + + # Solve lambda for triaxial (a > b > c) else: p0 = ( a**2 * b**2 * c**2 @@ -127,7 +285,7 @@ def calculate_lambda(x, y, z, a, b, c): # Clip the cos_theta to [-1, 1]. Due to floating point errors its value # could be slightly above 1 or slightly below -1. if isinstance(cos_theta, np.ndarray): - # Run inplace to avoid allocating a new array. + # Run in-place to avoid allocating a new array. np.clip(cos_theta, -1.0, 1.0, out=cos_theta) else: cos_theta = np.clip(cos_theta, -1.0, 1.0) @@ -152,16 +310,21 @@ def get_elliptical_integrals( These integrals are also called :math:`g_1`, :math:`g_2`, and :math:`g_3` in Takahashi et al. (2018). + .. important:: + + The elliptical integrals should not be used with spheres or ellipsoids that + approximate a sphere. + Parameters ---------- - a, b, c : floats - Semi-axes lengths of the given ellipsoid. + a, b, c : float + Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. lambda_ : float or (n,) array The given lambda value for the point we are considering. Returns ------- - A, B, C : floats or tuple of (n,) arrays + A, B, C : float or tuple of (n,) array The elliptical integrals evaluated for the given ellipsoid and observation point. @@ -198,34 +361,43 @@ def get_elliptical_integrals( Expressions of the elliptic integrals vary for each type of ellipsoid (triaxial, oblate and prolate). """ - if a > b > c: - g1, g2, g3 = _get_elliptical_integrals_triaxial(a, b, c, lambda_) - elif a > b and b == c: + check_semiaxes_sorted(a, b, c) + + if is_almost_a_sphere(a, b, c): + msg = ( + "Invalid semiaxes that create (almost) a spherical ellipsoid: " + f"a={a}, b={b}, c={c}." + ) + raise ValueError(msg) + + if is_almost_prolate(a, b, c): g1, g2, g3 = _get_elliptical_integrals_prolate(a, b, lambda_) - elif a < b and b == c: - g1, g2, g3 = _get_elliptical_integrals_oblate(a, b, lambda_) + elif is_almost_oblate(a, b, c): + g1, g2, g3 = _get_elliptical_integrals_oblate(b, c, lambda_) + elif a > b > c: + g1, g2, g3 = _get_elliptical_integrals_triaxial(a, b, c, lambda_) else: - msg = f"Invalid semiaxis lenghts: a={a}, b={b}, c={c}." + msg = f"Invalid semiaxes lengths: a={a}, b={b}, c={c}." raise ValueError(msg) return g1, g2, g3 -def _get_elliptical_integrals_triaxial(a, b, c, lmbda): +def _get_elliptical_integrals_triaxial(a, b, c, lambda_): r""" Compute elliptical integrals for a triaxial ellipsoid. Parameters ---------- - a, b, c : floats - Semi-axes lengths of the given ellipsoid. - lmbda : float + a, b, c : float + Semi-axes lengths of the ellipsoid sorted such as ``a > b > c``. + lambda_ : float or (n,) array The given lambda value for the point we are considering. Returns ------- - floats + A, B, C : float or tuple of (n,) array The elliptical integrals evaluated for the given ellipsoid and observation - point. + point(s). Notes ----- @@ -282,8 +454,12 @@ def _get_elliptical_integrals_triaxial(a, b, c, lmbda): from Takahashi et al. (2018). """ + if not a > b > c: + msg = f"Invalid semiaxes length (not a > b > c): a={a}, b={b}, c={c}." + raise ValueError(msg) + # Compute phi and kappa - int_arcsin = np.sqrt((a**2 - c**2) / (a**2 + lmbda)) + int_arcsin = np.sqrt((a**2 - c**2) / (a**2 + lambda_)) phi = np.arcsin(int_arcsin) k = (a**2 - b**2) / (a**2 - c**2) @@ -298,7 +474,7 @@ def _get_elliptical_integrals_triaxial(a, b, c, lmbda): g2_multiplier = (2 * np.sqrt(a**2 - c**2)) / ((a**2 - b**2) * (b**2 - c**2)) g2_elliptics = ellipe - ((b**2 - c**2) / (a**2 - c**2)) * ellipk g2_last_term = ((a**2 - b**2) / np.sqrt(a**2 - c**2)) * np.sqrt( - (c**2 + lmbda) / ((a**2 + lmbda) * (b**2 + lmbda)) + (c**2 + lambda_) / ((a**2 + lambda_) * (b**2 + lambda_)) ) g2 = g2_multiplier * (g2_elliptics - g2_last_term) @@ -307,29 +483,29 @@ def _get_elliptical_integrals_triaxial(a, b, c, lmbda): # (the minus sign is missing in Takahashi (2018)). g3_term_1 = -(2 / ((b**2 - c**2) * np.sqrt(a**2 - c**2))) * ellipe g3_term_2 = (2 / (b**2 - c**2)) * np.sqrt( - (b**2 + lmbda) / ((a**2 + lmbda) * (c**2 + lmbda)) + (b**2 + lambda_) / ((a**2 + lambda_) * (c**2 + lambda_)) ) g3 = g3_term_1 + g3_term_2 return g1, g2, g3 -def _get_elliptical_integrals_prolate(a, b, lmbda): +def _get_elliptical_integrals_prolate(a, b, lambda_): r""" Compute elliptical integrals for a prolate ellipsoid. Parameters ---------- - a, b : floats - Semi-axes lengths of the given ellipsoid. - lmbda : float + a, b : float + Semi-axes lengths of the given ellipsoid, where ``a > b``. + lambda_ : float or (n,) array The given lambda value for the point we are considering. Returns ------- - floats + A, B, C : float or tuple of (n,) array The elliptical integrals evaluated for the given ellipsoid and observation - point. + point(s). Notes ----- @@ -380,34 +556,38 @@ def _get_elliptical_integrals_prolate(a, b, lmbda): C(\lambda) = B(\lambda) """ + if not a > b: + msg = f"Invalid semiaxes length (not a > b): a={a}, b={b}." + raise ValueError(msg) + # Cache some reused variables e2 = a**2 - b**2 sqrt_e = np.sqrt(e2) - sqrt_l1 = np.sqrt(a**2 + lmbda) - sqrt_l2 = np.sqrt(b**2 + lmbda) + sqrt_l1 = np.sqrt(a**2 + lambda_) + sqrt_l2 = np.sqrt(b**2 + lambda_) log = np.log((sqrt_e + sqrt_l1) / sqrt_l2) g1 = (2 / (sqrt_e**3)) * (log - sqrt_e / sqrt_l1) - g2 = (1 / (sqrt_e**3)) * ((sqrt_e * sqrt_l1) / (b**2 + lmbda) - log) + g2 = (1 / (sqrt_e**3)) * ((sqrt_e * sqrt_l1) / (b**2 + lambda_) - log) return g1, g2, g2 -def _get_elliptical_integrals_oblate(a, b, lmbda): +def _get_elliptical_integrals_oblate(b, c, lambda_): r""" Compute elliptical integrals for a oblate ellipsoid. Parameters ---------- - a, b : floats - Semi-axes lengths of the given ellipsoid. - lmbda : float + b, c : float + Semi-axes lengths of the given ellipsoid, where ``b > c``. + lambda_ : float or (n,) array The given lambda value for the point we are considering. Returns ------- - floats + A, B, C : float or tuple of (n,) array The elliptical integrals evaluated for the given ellipsoid and observation - point. + point(s). Notes ----- @@ -418,46 +598,56 @@ def _get_elliptical_integrals_oblate(a, b, lmbda): .. math:: A(\lambda) = - \frac{ 2 }{ (b^2 - a^2)^{\frac{3}{2}} } + \frac{ 1 }{ (b^2 - c^2)^{\frac{3}{2}} } \left\{ - \sqrt{ \frac{b^2 - a^2}{a^2 + \lambda} } + \arctan \left[ \sqrt{ \frac{b^2 - c^2}{c^2 + \lambda} } \right] - - \arctan \left[ \sqrt{ \frac{b^2 - a^2}{a^2 + \lambda} } \right] + \frac{ + \sqrt{ (b^2 - c^2) (c^2 + \lambda) } + }{ + b^2 + \lambda + } \right\} .. math:: - B(\lambda) = - \frac{ 1 }{ (b^2 - a^2)^{\frac{3}{2}} } + C(\lambda) = + \frac{ 2 }{ (b^2 - c^2)^{\frac{3}{2}} } \left\{ - \arctan \left[ \sqrt{ \frac{b^2 - a^2}{a^2 + \lambda} } \right] + \sqrt{ \frac{b^2 - c^2}{c^2 + \lambda} } - - \frac{ - \sqrt{ (b^2 - a^2) (a^2 + \lambda) } - }{ - b^2 + \lambda - } + \arctan \left[ \sqrt{ \frac{b^2 - c^2}{c^2 + \lambda} } \right] \right\} + and .. math:: - C(\lambda) = B(\lambda) + B(\lambda) = A(\lambda) + + .. important:: + + These equations are modified versions of the ones in Takahashi (2018), adapted to + any oblate ellipsoid defined as: ``a = b > c``. """ - arctan = np.arctan(np.sqrt((b**2 - a**2) / (a**2 + lmbda))) + if not b > c: + msg = f"Invalid semiaxes length (not b > c): b={b}, c={c}." + raise ValueError(msg) + + arctan = np.arctan(np.sqrt((b**2 - c**2) / (c**2 + lambda_))) g1 = ( - 2 - / ((b**2 - a**2) ** (3 / 2)) - * ((np.sqrt((b**2 - a**2) / (a**2 + lmbda))) - arctan) - ) - g2 = ( 1 - / ((b**2 - a**2) ** (3 / 2)) - * (arctan - (np.sqrt((b**2 - a**2) * (a**2 + lmbda))) / (b**2 + lmbda)) + / ((b**2 - c**2) ** (3 / 2)) + * (arctan - (np.sqrt((b**2 - c**2) * (c**2 + lambda_))) / (b**2 + lambda_)) ) - return g1, g2, g2 + g3 = ( + 2 + / ((b**2 - c**2) ** (3 / 2)) + * ((np.sqrt((b**2 - c**2) / (c**2 + lambda_))) - arctan) + ) + return g1, g1, g3 def get_derivatives_of_elliptical_integrals( @@ -471,16 +661,17 @@ def get_derivatives_of_elliptical_integrals( Parameters ---------- - a, b, c : floats - Semi-axes lengths of the given ellipsoid. - lambda_ : float + a, b, c : float + Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. + lambda_ : float or (n,) array The given lambda value for the point we are considering. Returns ------- - hx, hy, hz : tuple of floats + hx, hy, hz : tuple of float or tuple of (n,) array The h values for the given observation point. """ + check_semiaxes_sorted(a, b, c) r = np.sqrt((a**2 + lambda_) * (b**2 + lambda_) * (c**2 + lambda_)) hx, hy, hz = tuple(-1 / (e**2 + lambda_) / r for e in (a, b, c)) return hx, hy, hz diff --git a/harmonica/tests/ellipsoids/test_ellipsoids.py b/harmonica/tests/ellipsoids/test_ellipsoids.py index 0d629bf24..b069bda5e 100644 --- a/harmonica/tests/ellipsoids/test_ellipsoids.py +++ b/harmonica/tests/ellipsoids/test_ellipsoids.py @@ -8,20 +8,13 @@ Test ellipsoid classes. """ -import itertools import re from unittest.mock import patch import numpy as np import pytest -from harmonica import ( - OblateEllipsoid, - ProlateEllipsoid, - Sphere, - TriaxialEllipsoid, - create_ellipsoid, -) +from harmonica import Ellipsoid try: import pyvista @@ -29,161 +22,8 @@ pyvista = None -class TestProlateEllipsoid: - """Test the ProlateEllipsoid class.""" - - @pytest.mark.parametrize(("a", "b"), [(35.0, 50.0), (50.0, 50.0)]) - def test_invalid_semiaxes(self, a, b): - """Test error if not a > b.""" - msg = re.escape("Invalid ellipsoid axis lengths for prolate ellipsoid") - with pytest.raises(ValueError, match=msg): - ProlateEllipsoid(a, b, yaw=0, pitch=0, center=(0, 0, 0)) - - @pytest.mark.parametrize("semiaxis", ["a", "b"]) - def test_non_positive_semiaxis(self, semiaxis): - """Test error after non-positive semiaxis.""" - match semiaxis: - case "a": - a, b = -1.0, 40.0 - case "b": - a, b = 50.0, -1.0 - case _: - raise ValueError() - msg = re.escape(f"Invalid value of '{semiaxis}' equal to '{-1.0}'") - with pytest.raises(ValueError, match=msg): - ProlateEllipsoid(a, b, yaw=0, pitch=0, center=(0, 0, 0)) - - @pytest.mark.parametrize("semiaxis", ["a", "b"]) - def test_non_positive_semiaxis_setter(self, semiaxis): - """Test error after non-positive semiaxis when using the setter.""" - a, b = 50.0, 35.0 - ellipsoid = ProlateEllipsoid(a, b, yaw=0, pitch=0, center=(0, 0, 0)) - msg = re.escape(f"Invalid value of '{semiaxis}' equal to '{-1.0}'") - with pytest.raises(ValueError, match=msg): - setattr(ellipsoid, semiaxis, -1.0) - - def test_semiaxes_setter(self): - """Test setters for semiaxes.""" - a, b = 50.0, 35.0 - ellipsoid = ProlateEllipsoid(a, b, yaw=0, pitch=0, center=(0, 0, 0)) - # Test setter of a - new_a = a + 1 - ellipsoid.a = new_a - assert ellipsoid.a == new_a - # Test setter of b - new_b = b + 1 - ellipsoid.b = new_b - assert ellipsoid.b == new_b - - def test_invalid_semiaxes_setter(self): - """Test error if not a > b when using the setter.""" - a, b = 50.0, 35.0 - ellipsoid = ProlateEllipsoid(a, b, yaw=0, pitch=0, center=(0, 0, 0)) - msg = re.escape("Invalid ellipsoid axis lengths for prolate ellipsoid") - with pytest.raises(ValueError, match=msg): - ellipsoid.a = 20.0 - with pytest.raises(ValueError, match=msg): - ellipsoid.b = 70.0 - - def test_value_of_c(self): - """Test if c is always equal to b.""" - a, b = 50.0, 35.0 - ellipsoid = ProlateEllipsoid(a, b, yaw=0, pitch=0, center=(0, 0, 0)) - assert ellipsoid.b == ellipsoid.c - # Update the value of b and check again - ellipsoid.b = 45.0 - assert ellipsoid.b == ellipsoid.c - - def test_roll_equal_to_zero(self): - """Test if roll is always equal to zero.""" - a, b = 50.0, 35.0 - ellipsoid = ProlateEllipsoid(a, b, yaw=0, pitch=0, center=(0, 0, 0)) - assert ellipsoid.roll == 0.0 - - -class TestOblateEllipsoid: - """Test the OblateEllipsoid class.""" - - @pytest.mark.parametrize(("a", "b"), [(50.0, 35.0), (50.0, 50.0)]) - def test_invalid_semiaxes(self, a, b): - """Test error if not a < b.""" - msg = re.escape("Invalid ellipsoid axis lengths for oblate ellipsoid") - with pytest.raises(ValueError, match=msg): - OblateEllipsoid(a, b, yaw=0, pitch=0, center=(0, 0, 0)) - - @pytest.mark.parametrize("semiaxis", ["a", "b"]) - def test_non_positive_semiaxis(self, semiaxis): - """Test error after non-positive semiaxis.""" - match semiaxis: - case "a": - a, b = -1.0, 40.0 - case "b": - a, b = 50.0, -1.0 - case _: - raise ValueError() - msg = re.escape(f"Invalid value of '{semiaxis}' equal to '{-1.0}'") - with pytest.raises(ValueError, match=msg): - OblateEllipsoid(a, b, yaw=0, pitch=0, center=(0, 0, 0)) - - @pytest.mark.parametrize("semiaxis", ["a", "b"]) - def test_non_positive_semiaxis_setter(self, semiaxis): - """Test error after non-positive semiaxis when using the setter.""" - a, b = 35.0, 50.0 - ellipsoid = OblateEllipsoid(a, b, yaw=0, pitch=0, center=(0, 0, 0)) - msg = re.escape(f"Invalid value of '{semiaxis}' equal to '{-1.0}'") - with pytest.raises(ValueError, match=msg): - setattr(ellipsoid, semiaxis, -1.0) - - def test_semiaxes_setter(self): - """Test setters for semiaxes.""" - a, b = 35.0, 50.0 - ellipsoid = OblateEllipsoid(a, b, yaw=0, pitch=0, center=(0, 0, 0)) - # Test setter of a - new_a = a + 1 - ellipsoid.a = new_a - assert ellipsoid.a == new_a - # Test setter of b - new_b = b + 1 - ellipsoid.b = new_b - assert ellipsoid.b == new_b - - def test_invalid_semiaxes_setter(self): - """Test error if not a < b when using the setter.""" - a, b = 35.0, 50.0 - ellipsoid = OblateEllipsoid(a, b, yaw=0, pitch=0, center=(0, 0, 0)) - msg = re.escape("Invalid ellipsoid axis lengths for oblate ellipsoid") - with pytest.raises(ValueError, match=msg): - ellipsoid.a = 70.0 - with pytest.raises(ValueError, match=msg): - ellipsoid.b = 20.0 - - def test_value_of_c(self): - """Test if c is always equal to b.""" - a, b = 35.0, 50.0 - ellipsoid = OblateEllipsoid(a, b, yaw=0, pitch=0, center=(0, 0, 0)) - assert ellipsoid.b == ellipsoid.c - # Update the value of b and check again - ellipsoid.b = 45.0 - assert ellipsoid.b == ellipsoid.c - - def test_roll_equal_to_zero(self): - """Test if roll is always equal to zero.""" - a, b = 35.0, 50.0 - ellipsoid = OblateEllipsoid(a, b, yaw=0, pitch=0, center=(0, 0, 0)) - assert ellipsoid.roll == 0.0 - - -class TestTriaxialEllipsoid: - """Test the TriaxialEllipsoid class.""" - - @pytest.mark.parametrize( - ("a", "b", "c"), [(50.0, 35.0, 45.0), (50.0, 50.0, 50.0), (60.0, 50.0, 50.0)] - ) - def test_invalid_semiaxes(self, a, b, c): - """Test error if not a > b > c.""" - msg = re.escape("Invalid ellipsoid axis lengths for triaxial ellipsoid") - with pytest.raises(ValueError, match=msg): - TriaxialEllipsoid(a, b, c, yaw=0, pitch=0, roll=0, center=(0, 0, 0)) +class TestEllipsoid: + """Test the Ellipsoid class.""" @pytest.mark.parametrize("semiaxis", ["a", "b", "c"]) def test_non_positive_semiaxis(self, semiaxis): @@ -199,13 +39,13 @@ def test_non_positive_semiaxis(self, semiaxis): raise ValueError() msg = re.escape(f"Invalid value of '{semiaxis}' equal to '{-1.0}'") with pytest.raises(ValueError, match=msg): - TriaxialEllipsoid(a, b, c, yaw=0, pitch=0, roll=0, center=(0, 0, 0)) + Ellipsoid(a, b, c) @pytest.mark.parametrize("semiaxis", ["a", "b", "c"]) def test_non_positive_semiaxis_setter(self, semiaxis): """Test error after non-positive semiaxis when using the setter.""" a, b, c = 50.0, 40.0, 35.0 - ellipsoid = TriaxialEllipsoid(a, b, c, yaw=0, pitch=0, roll=0, center=(0, 0, 0)) + ellipsoid = Ellipsoid(a, b, c) msg = re.escape(f"Invalid value of '{semiaxis}' equal to '{-1.0}'") with pytest.raises(ValueError, match=msg): setattr(ellipsoid, semiaxis, -1.0) @@ -213,7 +53,7 @@ def test_non_positive_semiaxis_setter(self, semiaxis): def test_semiaxes_setter(self): """Test setters for semiaxes.""" a, b, c = 50.0, 40.0, 35.0 - ellipsoid = TriaxialEllipsoid(a, b, c, yaw=0, pitch=0, roll=0, center=(0, 0, 0)) + ellipsoid = Ellipsoid(a, b, c) # Test setter of a new_a = a + 1 ellipsoid.a = new_a @@ -227,149 +67,56 @@ def test_semiaxes_setter(self): ellipsoid.c = new_c assert ellipsoid.c == new_c - def test_invalid_semiaxes_setter(self): - """Test error if not a > b > c when using the setter.""" - a, b, c = 50.0, 40.0, 30.0 - ellipsoid = TriaxialEllipsoid(a, b, c, yaw=0, pitch=0, roll=0, center=(0, 0, 0)) - msg = re.escape("Invalid ellipsoid axis lengths for triaxial ellipsoid") - with pytest.raises(ValueError, match=msg): - ellipsoid.a = 30.0 - with pytest.raises(ValueError, match=msg): - ellipsoid.b = 20.0 - with pytest.raises(ValueError, match=msg): - ellipsoid.c = 70.0 - - -class TestSphere: - """Test the Sphere class.""" - - def test_non_positive_semiaxis(self): - """Test error after non-positive a semiaxis.""" - a = -1.0 - msg = re.escape(f"Invalid value of 'a' equal to '{a}'") - with pytest.raises(ValueError, match=msg): - Sphere(a=a, center=(0, 0, 0)) - - def test_non_positive_semiaxis_setter(self): - """Test error after non-positive semiaxis when using the setter.""" - sphere = Sphere(a=10.0, center=(0, 0, 0)) - new_value = -1.0 - msg = re.escape(f"Invalid value of 'a' equal to '{new_value}'") - with pytest.raises(ValueError, match=msg): - sphere.a = new_value - - def test_semiaxes_setter(self): - """Test setters for semiaxes.""" - a = 10.0 - sphere = Sphere(a, center=(0, 0, 0)) - new_a = a + 1 - sphere.a = new_a - assert sphere.a == new_a - - @pytest.mark.parametrize("semiaxis", ["b", "c"]) - def test_value_of_b_and_c(self, semiaxis): - """Test if b and c are always equal to a.""" - a = 10.0 - sphere = Sphere(a, center=(0, 0, 0)) - assert sphere.a == getattr(sphere, semiaxis) - # Update the value of the semiaxis and check again - sphere.a = 45.0 - assert sphere.a == getattr(sphere, semiaxis) - - def test_radius(self): - """Test that the radius is always equal to a.""" - a = 10.0 - sphere = Sphere(a, center=(0, 0, 0)) - assert sphere.a == sphere.radius - # Update the value of the semiaxis and check again - sphere.a = 45.0 - assert sphere.a == sphere.radius - - def test_angles_equal_to_zero(self): - """Test if all angles are always equal to zero.""" - sphere = Sphere(10.0, center=(0, 0, 0)) - assert sphere.yaw == 0.0 - assert sphere.pitch == 0.0 - assert sphere.roll == 0.0 - - def test_rotation_matrix(self): + def test_rotation_matrix_of_sphere(self): """Make sure the rotation matrix of a sphere is always the identity.""" - sphere = Sphere(10.0, center=(0, 0, 0)) + a = 10.0 + sphere = Ellipsoid(a, a, a) np.testing.assert_array_equal( sphere.rotation_matrix, np.eye(3, dtype=np.float64) ) -@pytest.mark.parametrize( - "ellipsoid_class", [OblateEllipsoid, ProlateEllipsoid, TriaxialEllipsoid, Sphere] -) class TestPhysicalProperties: @pytest.fixture - def ellipsoid_args(self, ellipsoid_class): - if ellipsoid_class is OblateEllipsoid: - args = { - "a": 20.0, - "b": 50.0, - "pitch": 0.0, - "yaw": 0.0, - "center": (0, 0, 0), - } - elif ellipsoid_class is ProlateEllipsoid: - args = { - "a": 50.0, - "b": 20.0, - "pitch": 0.0, - "yaw": 0.0, - "center": (0, 0, 0), - } - elif ellipsoid_class is TriaxialEllipsoid: - args = { - "a": 50.0, - "b": 20.0, - "c": 10.0, - "pitch": 0.0, - "yaw": 0.0, - "roll": 0.0, - "center": (0, 0, 0), - } - elif ellipsoid_class is Sphere: - args = {"a": 50.0, "center": (0, 0, 0)} - else: - raise TypeError() - return args + def semiaxes(self): + """Ellipsoid semiaxes.""" + return 50.0, 35.0, 25.0 - def test_density(self, ellipsoid_class, ellipsoid_args): + def test_density(self, semiaxes): """ Test assigning density to the ellipsoid. """ + a, b, c = semiaxes density = 3.0 - ellipsoid = ellipsoid_class(**ellipsoid_args, density=density) + ellipsoid = Ellipsoid(a, b, c, density=density) assert ellipsoid.density == density # Check overwriting it new_density = -4.0 ellipsoid.density = new_density assert ellipsoid.density == new_density # Check density as None - ellipsoid = ellipsoid_class(**ellipsoid_args, density=None) + ellipsoid = Ellipsoid(a, b, c, density=None) assert ellipsoid.density is None - def test_invalid_density(self, ellipsoid_class, ellipsoid_args): + def test_invalid_density(self, semiaxes): """ Test errors after invalid density. """ + a, b, c = semiaxes density = [1.0, 3.0] msg = re.escape("Invalid 'density' of type 'list'") with pytest.raises(TypeError, match=msg): - ellipsoid_class(**ellipsoid_args, density=density) + Ellipsoid(a, b, c, density=density) @pytest.mark.parametrize("susceptibility", [0.1, "tensor", None]) - def test_susceptibility(self, ellipsoid_class, ellipsoid_args, susceptibility): + def test_susceptibility(self, semiaxes, susceptibility): """ Test assigning susceptibility to the ellipsoid. """ + a, b, c = semiaxes if susceptibility == "tensor": susceptibility = np.random.default_rng(seed=42).uniform(size=(3, 3)) - ellipsoid = ellipsoid_class(**ellipsoid_args, susceptibility=susceptibility) + ellipsoid = Ellipsoid(a, b, c, susceptibility=susceptibility) # Check if it was correctly assigned if susceptibility is None: @@ -384,29 +131,32 @@ def test_susceptibility(self, ellipsoid_class, ellipsoid_args, susceptibility): ellipsoid.susceptibility = new_sus assert ellipsoid.susceptibility == new_sus - def test_invalid_susceptibility_type(self, ellipsoid_class, ellipsoid_args): + def test_invalid_susceptibility_type(self, semiaxes): """ Test errors after invalid susceptibility type. """ + a, b, c = semiaxes susceptibility = [1.0, 3.0] msg = re.escape("Invalid 'susceptibility' of type 'list'") with pytest.raises(TypeError, match=msg): - ellipsoid_class(**ellipsoid_args, susceptibility=susceptibility) + Ellipsoid(a, b, c, susceptibility=susceptibility) - def test_invalid_susceptibility_shape(self, ellipsoid_class, ellipsoid_args): + def test_invalid_susceptibility_shape(self, semiaxes): """ Test errors after invalid susceptibility shape. """ + a, b, c = semiaxes susceptibility = np.array([[1, 2, 3], [4, 5, 6]]) msg = re.escape("Invalid 'susceptibility' as an array with shape '(2, 3)'") with pytest.raises(ValueError, match=msg): - ellipsoid_class(**ellipsoid_args, susceptibility=susceptibility) + Ellipsoid(a, b, c, susceptibility=susceptibility) @pytest.mark.parametrize("remanent_mag_type", ["array", "list", None]) - def test_remanent_mag(self, ellipsoid_class, ellipsoid_args, remanent_mag_type): + def test_remanent_mag(self, semiaxes, remanent_mag_type): """ Test assigning susceptibility to the ellipsoid. """ + a, b, c = semiaxes if remanent_mag_type == "array": remanent_mag = np.array([1.0, 2.0, 3.0]) elif remanent_mag_type == "list": @@ -414,7 +164,7 @@ def test_remanent_mag(self, ellipsoid_class, ellipsoid_args, remanent_mag_type): else: remanent_mag = None - ellipsoid = ellipsoid_class(**ellipsoid_args, remanent_mag=remanent_mag) + ellipsoid = Ellipsoid(a, b, c, remanent_mag=remanent_mag) # Check if it was correctly assigned if remanent_mag is None: @@ -427,152 +177,40 @@ def test_remanent_mag(self, ellipsoid_class, ellipsoid_args, remanent_mag_type): ellipsoid.remanent_mag = new_remanent_mag np.testing.assert_almost_equal(ellipsoid.remanent_mag, new_remanent_mag) - def test_invalid_remanent_mag_type(self, ellipsoid_class, ellipsoid_args): + def test_invalid_remanent_mag_type(self, semiaxes): """ Test errors after invalid remanent_mag type. """ class Dummy: ... + a, b, c = semiaxes remanent_mag = Dummy() msg = re.escape("Invalid 'remanent_mag' of type 'Dummy'") with pytest.raises(TypeError, match=msg): - ellipsoid_class(**ellipsoid_args, remanent_mag=remanent_mag) + Ellipsoid(a, b, c, remanent_mag=remanent_mag) - def test_invalid_remanent_mag_shape(self, ellipsoid_class, ellipsoid_args): + def test_invalid_remanent_mag_shape(self, semiaxes): """ Test errors after invalid remanent_mag shape. """ + a, b, c = semiaxes remanent_mag = np.array([1.0, 2.0]) msg = re.escape("Invalid 'remanent_mag' with shape '(2,)'") with pytest.raises(ValueError, match=msg): - ellipsoid_class(**ellipsoid_args, remanent_mag=remanent_mag) - - -class TestCreateEllipsoid: - """ - Test the ``create_ellipsoid`` function. - """ - - def test_sphere(self): - """Test ``create_ellipsoid`` when expecting a ``Sphere``.""" - a = 50 - center = (1.0, 2.0, 3.0) - sphere = create_ellipsoid(a, a, a, center=center) - assert isinstance(sphere, Sphere) - assert sphere.a == a - np.testing.assert_almost_equal(sphere.center, center) - - @pytest.mark.parametrize( - ("a", "b", "c"), [(2.0, 1.0, 1.0), (1.0, 2.0, 1.0), (1.0, 1.0, 2.0)] - ) - def test_prolate(self, a, b, c): - """Test ``create_ellipsoid`` when expecting a ``ProlateEllipsoid``.""" - center = (1.0, 2.0, 3.0) - yaw, pitch = 30.0, -17.0 - ellipsoid = create_ellipsoid(a, b, c, center=center, yaw=yaw, pitch=pitch) - assert isinstance(ellipsoid, ProlateEllipsoid) - assert ellipsoid.a == 2.0 - assert ellipsoid.b == 1.0 - assert ellipsoid.c == 1.0 - assert ellipsoid.yaw == yaw - assert ellipsoid.pitch == pitch - np.testing.assert_almost_equal(ellipsoid.center, center) - - @pytest.mark.parametrize( - ("a", "b", "c"), [(2.0, 3.0, 3.0), (3.0, 2.0, 3.0), (3.0, 3.0, 2.0)] - ) - def test_oblate(self, a, b, c): - """Test ``create_ellipsoid`` when expecting a ``OblateEllipsoid``.""" - center = (1.0, 2.0, 3.0) - yaw, pitch = 30.0, -17.0 - ellipsoid = create_ellipsoid(a, b, c, center=center, yaw=yaw, pitch=pitch) - assert isinstance(ellipsoid, OblateEllipsoid) - assert ellipsoid.a == 2.0 - assert ellipsoid.b == 3.0 - assert ellipsoid.c == 3.0 - assert ellipsoid.yaw == yaw - assert ellipsoid.pitch == pitch - np.testing.assert_almost_equal(ellipsoid.center, center) - - @pytest.mark.parametrize(("a", "b", "c"), itertools.permutations((1.0, 2.0, 3.0))) - def test_triaxial(self, a, b, c): - """Test ``create_ellipsoid`` when expecting a ``TriaxialEllipsoid``.""" - center = (1.0, 2.0, 3.0) - yaw, pitch, roll = 30.0, -17.0, 45.0 - ellipsoid = create_ellipsoid( - a, b, c, center=center, yaw=yaw, pitch=pitch, roll=roll - ) - assert isinstance(ellipsoid, TriaxialEllipsoid) - assert ellipsoid.a == 3.0 - assert ellipsoid.b == 2.0 - assert ellipsoid.c == 1.0 - assert ellipsoid.yaw == yaw - assert ellipsoid.pitch == pitch - assert ellipsoid.roll == roll - np.testing.assert_almost_equal(ellipsoid.center, center) - - @pytest.mark.parametrize( - ("a", "b", "c"), - [(-1.0, 2.0, 3.0), (1.0, -2.0, 3.0), (1.0, 2.0, -3.0), (0.0, 1.0, 2.0)], - ids=["a", "b", "c", "a-zero"], - ) - def test_invalid_semiaxes(self, a, b, c): - """Test error when semiaxes are non-positive.""" - msg = ( - r"Invalid value of '[abc]' equal to '-?[0-9]+\.[0-9]+'\. " - r"It must be positive\." - ) - with pytest.raises(ValueError, match=msg): - create_ellipsoid(a, b, c, center=(0, 1, 2)) - - @pytest.mark.parametrize( - ("a", "b", "c"), - [(1.0, 1.0, 1.0), (1.0, 1.0, 2.0), (2.0, 2.0, 1.0), (1.0, 2.0, 3.0)], - ids=["sphere", "prolate", "oblate", "triaxial"], - ) - @pytest.mark.parametrize( - "physical_property", - ["density", "susceptibility", "remanent_mag"], - ) - def test_physical_properties(self, a, b, c, physical_property): - """Test if physical properties are present in the objects.""" - match physical_property: - case "density": - value = 200.0 - case "susceptibility": - value = 0.1 - case "remanent_mag": - value = np.array([1.0, 2.0, 3.0]) - case _: - raise ValueError() - kwargs = {physical_property: value} - center = (1.0, 2.0, 3.0) - ellipsoid = create_ellipsoid(a, b, c, center=center, **kwargs) - np.testing.assert_allclose(getattr(ellipsoid, physical_property), value) + Ellipsoid(a, b, c, remanent_mag=remanent_mag) @pytest.mark.skipif(pyvista is None, reason="requires pyvista") class TestToPyvista: """Test exporting ellipsoids to PyVista objects.""" - @pytest.fixture(params=["oblate", "prolate", "triaxial", "sphere"]) - def ellipsoid(self, request): + @pytest.fixture + def ellipsoid(self): a, b, c = 3.0, 2.0, 1.0 yaw, pitch, roll = 73.0, 14.0, -35.0 center = (43.0, -72.0, 105) - match request.param: - case "oblate": - ellipsoid = OblateEllipsoid(b, a, yaw, pitch, center=center) - case "prolate": - ellipsoid = ProlateEllipsoid(a, b, yaw, pitch, center=center) - case "triaxial": - ellipsoid = TriaxialEllipsoid(a, b, c, yaw, pitch, roll, center=center) - case "sphere": - ellipsoid = Sphere(a, center=center) - case _: - raise ValueError() - return ellipsoid + return Ellipsoid(a, b, c, yaw=yaw, pitch=pitch, roll=roll, center=center) @patch("harmonica._forward.ellipsoids.ellipsoids.pyvista", None) def test_pyvista_missing_error(self, ellipsoid): @@ -604,64 +242,21 @@ class TestString: yaw, pitch, roll = 73, 14, -35 center = (43.0, -72.0, 105) - def build_ellipsoid(self, ellipsoid_type): - match ellipsoid_type: - case "oblate": - ellipsoid = OblateEllipsoid( - self.b, self.a, self.yaw, self.pitch, center=self.center - ) - case "prolate": - ellipsoid = ProlateEllipsoid( - self.a, self.b, self.yaw, self.pitch, center=self.center - ) - case "triaxial": - ellipsoid = TriaxialEllipsoid( - self.a, - self.b, - self.c, - self.yaw, - self.pitch, - self.roll, - center=self.center, - ) - case "sphere": - ellipsoid = Sphere(self.a, center=self.center) - case _: - raise ValueError() - return ellipsoid - - def test_prolate(self): - ellipsoid = self.build_ellipsoid("prolate") - expected = ( - "ProlateEllipsoid:" - "\n • a: 3.0 m" - "\n • b: 2.0 m" - "\n • c: 2.0 m" - "\n • center: (43.0, -72.0, 105.0) m" - "\n • yaw: 73.0" - "\n • pitch: 14.0" - "\n • roll: 0.0" - ) - assert expected == str(ellipsoid) - - def test_oblate(self): - ellipsoid = self.build_ellipsoid("oblate") - expected = ( - "OblateEllipsoid:" - "\n • a: 2.0 m" - "\n • b: 3.0 m" - "\n • c: 3.0 m" - "\n • center: (43.0, -72.0, 105.0) m" - "\n • yaw: 73.0" - "\n • pitch: 14.0" - "\n • roll: 0.0" + @pytest.fixture + def ellipsoid(self): + return Ellipsoid( + self.a, + self.b, + self.c, + yaw=self.yaw, + pitch=self.pitch, + roll=self.roll, + center=self.center, ) - assert expected == str(ellipsoid) - def test_triaxial(self): - ellipsoid = self.build_ellipsoid("triaxial") + def test_triaxial(self, ellipsoid): expected = ( - "TriaxialEllipsoid:" + "Ellipsoid:" "\n • a: 3.0 m" "\n • b: 2.0 m" "\n • c: 1.0 m" @@ -672,25 +267,7 @@ def test_triaxial(self): ) assert expected == str(ellipsoid) - def test_sphere(self): - ellipsoid = self.build_ellipsoid("sphere") - expected = ( - "Sphere:" - "\n • a: 3.0 m" - "\n • b: 3.0 m" - "\n • c: 3.0 m" - "\n • center: (43.0, -72.0, 105.0) m" - "\n • yaw: 0.0" - "\n • pitch: 0.0" - "\n • roll: 0.0" - ) - assert expected == str(ellipsoid) - - @pytest.mark.parametrize( - "ellipsoid_type", ["prolate", "oblate", "sphere", "triaxial"] - ) - def test_density(self, ellipsoid_type): - ellipsoid = self.build_ellipsoid(ellipsoid_type) + def test_density(self, ellipsoid): ellipsoid.density = -400 (density_line,) = [ line for line in str(ellipsoid).split("\n") if "density" in line @@ -698,11 +275,7 @@ def test_density(self, ellipsoid_type): expected_line = " • density: -400.0 kg/m³" assert density_line == expected_line - @pytest.mark.parametrize( - "ellipsoid_type", ["prolate", "oblate", "sphere", "triaxial"] - ) - def test_susceptibility(self, ellipsoid_type): - ellipsoid = self.build_ellipsoid(ellipsoid_type) + def test_susceptibility(self, ellipsoid): ellipsoid.susceptibility = 0.3 (sus_line,) = [ line for line in str(ellipsoid).splitlines() if "susceptibility" in line @@ -710,11 +283,7 @@ def test_susceptibility(self, ellipsoid_type): expected_line = " • susceptibility: 0.3" assert sus_line == expected_line - @pytest.mark.parametrize( - "ellipsoid_type", ["prolate", "oblate", "sphere", "triaxial"] - ) - def test_remanent_mag(self, ellipsoid_type): - ellipsoid = self.build_ellipsoid(ellipsoid_type) + def test_remanent_mag(self, ellipsoid): ellipsoid.remanent_mag = (12, -43, 59) (rem_line,) = [ line for line in str(ellipsoid).splitlines() if "remanent_mag" in line @@ -722,11 +291,7 @@ def test_remanent_mag(self, ellipsoid_type): expected_line = " • remanent_mag: (12.0, -43.0, 59.0) A/m" assert rem_line == expected_line - @pytest.mark.parametrize( - "ellipsoid_type", ["prolate", "oblate", "sphere", "triaxial"] - ) - def test_susceptibility_tensor(self, ellipsoid_type): - ellipsoid = self.build_ellipsoid(ellipsoid_type) + def test_susceptibility_tensor(self, ellipsoid): sus = np.array([[1, 0, 0], [0, 2, 0], [0, 0, 3]]) ellipsoid.susceptibility = sus # Grab the susceptibility tensor lines diff --git a/harmonica/tests/ellipsoids/test_gravity.py b/harmonica/tests/ellipsoids/test_gravity.py index cb626dc84..68a55d6e0 100644 --- a/harmonica/tests/ellipsoids/test_gravity.py +++ b/harmonica/tests/ellipsoids/test_gravity.py @@ -22,12 +22,8 @@ import verde as vd from harmonica import ellipsoid_gravity, point_gravity -from harmonica._forward.ellipsoids import ( - OblateEllipsoid, - ProlateEllipsoid, - Sphere, - TriaxialEllipsoid, -) +from harmonica._forward.ellipsoids import Ellipsoid +from harmonica._forward.ellipsoids.utils import SEMIAXES_RTOL from harmonica.errors import NoPhysicalPropertyWarning @@ -46,49 +42,21 @@ def build_ellipsoid(ellipsoid_type, *, center=(0, 0, 0), density=None): match ellipsoid_type: case "triaxial": a, b, c = 3.2, 2.1, 1.3 - ellipsoid = TriaxialEllipsoid( - a=a, b=b, c=c, yaw=0, pitch=0, roll=0, center=center, density=density - ) case "prolate": a, b = 3.2, 2.1 - ellipsoid = ProlateEllipsoid( - a=a, b=b, yaw=0, pitch=0, center=center, density=density - ) + c = b case "oblate": a, b = 2.2, 3.1 - ellipsoid = OblateEllipsoid( - a=a, b=b, yaw=0, pitch=0, center=center, density=density - ) + c = b case "sphere": a = 3.2 - ellipsoid = Sphere(a=a, center=center, density=density) + b, c = a, a case _: msg = f"Invalid ellipsoid type: {ellipsoid_type}" raise ValueError(msg) - return ellipsoid - - -def test_degenerate_ellipsoid_cases(): - """ - Test cases where the ellipsoid axes lengths are close to the boundary of - accepted values. - - """ - # ellipsoids take (a, b, #c, yaw, pitch, #roll, center) - a, b, c = 5, 4.99999999, 4.99999998 - yaw, pitch, roll = 0, 0, 0 - center = (0, 0, 0) - density = 2000 - tri = TriaxialEllipsoid(a, b, c, yaw, pitch, roll, center, density=density) - pro = ProlateEllipsoid(a, b, yaw, pitch, center, density=density) - obl = OblateEllipsoid(b, a, yaw, pitch, center, density=density) - coordinates = vd.grid_coordinates( - region=(-20, 20, -20, 20), spacing=0.5, extra_coords=5 - ) - _, _, _ = ellipsoid_gravity(coordinates, tri) - _, _, _ = ellipsoid_gravity(coordinates, pro) - _, _, _ = ellipsoid_gravity(coordinates, obl) + ellipsoid = Ellipsoid(a, b, c, center=center, density=density) + return ellipsoid def test_opposite_planes(): @@ -99,12 +67,8 @@ def test_opposite_planes(): """ a, b, c = (4, 3, 2) # triaxial ellipsoid - yaw, pitch, roll = 90, 0, 0 - center = (0, 0, 0) density = 2000 - triaxial_example = TriaxialEllipsoid( - a, b, c, yaw, pitch, roll, center, density=density - ) + triaxial_example = Ellipsoid(a, b, c, density=density) # define observation points (2D grid) at surface height (z axis, # 'Upward') = 5 @@ -130,9 +94,7 @@ def test_int_ext_boundary(): # compare a set value apart a, b, c = (5, 4, 3) - ellipsoid = TriaxialEllipsoid( - a, b, c, yaw=0, pitch=0, roll=0, center=(0, 0, 0), density=2000 - ) + ellipsoid = Ellipsoid(a, b, c, density=2000) e = np.array([[4.9999999, 5.00000001]]) n = np.array([[0.0, 0.0]]) @@ -330,7 +292,7 @@ def test_sphere_vs_point_source(self): radius = 50.0 center = (10, -29, 105) density = 200.0 - sphere = Sphere(radius, center=center, density=density) + sphere = Ellipsoid(radius, radius, radius, center=center, density=density) # Build a 3d grid of observation points centered in (0, 0, 0) n = 51 @@ -374,55 +336,43 @@ class TestEllipsoidVsSphere: center = (28, 19, -50) density = 200.0 - # Difference between ellipsoid's semiaxes. - # It should be small compared to the sphere radius, so the ellipsoid approximates - # a sphere. - delta = 1e-4 + # Ratio between ellipsoid's semiaxes, defined as ratio = | a - b | / max(a, b). + # Make it small enough so ellipsoids approximate a sphere, but not too small that + # might trigger the fixes for numerical uncertainties. + ratio = 1e-4 + assert ratio > SEMIAXES_RTOL @pytest.fixture def sphere(self): """Sphere used to compare gravity fields.""" - return Sphere(self.radius, self.center, density=self.density) + return Ellipsoid( + self.radius, + self.radius, + self.radius, + center=self.center, + density=self.density, + ) @pytest.fixture(params=["oblate", "prolate", "triaxial"]) def ellipsoid(self, request): """ Ellipsoid that approximates a sphere. """ - yaw, pitch, roll = 0, 0, 0 a = self.radius match request.param: case "oblate": - ellipsoid = OblateEllipsoid( - a=a, - b=a + self.delta, - yaw=yaw, - pitch=pitch, - center=self.center, - density=self.density, - ) + a = b = self.radius + c = (1 - self.ratio) * b case "prolate": - ellipsoid = ProlateEllipsoid( - a=a, - b=a - self.delta, - yaw=yaw, - pitch=pitch, - center=self.center, - density=self.density, - ) + a = self.radius + b = c = (1 - self.ratio) * a case "triaxial": - ellipsoid = TriaxialEllipsoid( - a=a, - b=a - self.delta, - c=a - 2 * self.delta, - yaw=yaw, - pitch=pitch, - roll=roll, - center=self.center, - density=self.density, - ) + a = self.radius + b = (1 - self.ratio) * a + c = (1 - 2 * self.ratio) * a case _: raise ValueError() + ellipsoid = Ellipsoid(a, b, c, center=self.center, density=self.density) return ellipsoid @pytest.fixture @@ -448,9 +398,7 @@ def test_ellipsoid_vs_sphere(self, coordinates, ellipsoid, sphere): g_ellipsoid = ellipsoid_gravity(coordinates, ellipsoid) # Compare the two fields - maxabs = vd.maxabs(*g_sphere, *g_ellipsoid) - atol = maxabs * 1e-5 - np.testing.assert_allclose(g_sphere, g_ellipsoid, atol=atol) + np.testing.assert_allclose(g_sphere, g_ellipsoid, rtol=7e-4) class TestSymmetryOnRotations: @@ -467,8 +415,7 @@ def flip_ellipsoid(self, ellipsoid): """ ellipsoid.yaw += 180 ellipsoid.pitch *= -1 - if isinstance(ellipsoid, TriaxialEllipsoid): - ellipsoid.roll *= -1 + ellipsoid.roll *= -1 return ellipsoid @pytest.fixture(params=["oblate", "prolate", "triaxial"]) @@ -482,41 +429,24 @@ def ellipsoid(self, request): density = 238 match ellipsoid_type: case "oblate": - ellipsoid = OblateEllipsoid( - a=semiminor, - b=semimajor, - yaw=yaw, - pitch=pitch, - center=center, - density=density, - ) + a, b = semiminor, semimajor + c = b case "prolate": - ellipsoid = ProlateEllipsoid( - a=semimajor, - b=semiminor, - yaw=yaw, - pitch=pitch, - center=center, - density=density, - ) + a, b = semimajor, semiminor + c = b case "triaxial": - ellipsoid = TriaxialEllipsoid( - a=semimajor, - b=semimiddle, - c=semiminor, - yaw=yaw, - pitch=pitch, - roll=roll, - center=center, - density=density, - ) + a, b, c = semimajor, semimiddle, semiminor + c = b case _: raise ValueError() + ellipsoid = Ellipsoid( + a, b, c, yaw=yaw, pitch=pitch, roll=roll, center=center, density=density + ) return ellipsoid def test_symmetry_when_flipping(self, ellipsoid): """ - Test symmetry of magnetic field when flipping the ellipsoid. + Test symmetry of gravity field when flipping the ellipsoid. Rotate the ellipsoid so the geometry is preserved. The gravity field generated by the ellipsoid should be the same as before the rotation. @@ -529,7 +459,7 @@ def test_symmetry_when_flipping(self, ellipsoid): # Generate a flipped ellipsoid ellipsoid_flipped = self.flip_ellipsoid(copy(ellipsoid)) - # Compute magnetic fields + # Compute gravity fields g_field, g_field_flipped = tuple( ellipsoid_gravity(coordinates, ell) for ell in (ellipsoid, ellipsoid_flipped) @@ -558,23 +488,25 @@ def coordinates(self): def ellipsoids(self): """Sample ellipsoids.""" ellipsoids = [ - OblateEllipsoid( + Ellipsoid( a=20, b=60, + c=60, yaw=30.2, pitch=-23, center=(-10.0, 20.0, -10.0), density=200.0, ), - ProlateEllipsoid( + Ellipsoid( a=40, b=15, + c=15, yaw=170.2, pitch=71, center=(15.0, 0.0, -40.0), density=-400.0, ), - TriaxialEllipsoid( + Ellipsoid( a=60, b=18, c=15, @@ -607,52 +539,15 @@ def test_multiple_ellipsoids(self, coordinates, ellipsoids): np.testing.assert_allclose(gz, gz_expected) -@pytest.mark.parametrize( - "ellipsoid_class", [OblateEllipsoid, ProlateEllipsoid, TriaxialEllipsoid] -) class TestNoneDensity: """Test warning when ellipsoid has no density.""" - @pytest.fixture - def ellipsoid_args(self, ellipsoid_class): - if ellipsoid_class is OblateEllipsoid: - args = { - "a": 20.0, - "b": 50.0, - "pitch": 0.0, - "yaw": 0.0, - "center": (0, 0, 0), - } - elif ellipsoid_class is ProlateEllipsoid: - args = { - "a": 50.0, - "b": 20.0, - "pitch": 0.0, - "yaw": 0.0, - "center": (0, 0, 0), - } - elif ellipsoid_class is TriaxialEllipsoid: - args = { - "a": 50.0, - "b": 20.0, - "c": 10.0, - "pitch": 0.0, - "yaw": 0.0, - "roll": 0.0, - "center": (0, 0, 0), - } - elif ellipsoid_class is Sphere: - args = {"a": 50.0, "center": (0, 0, 0)} - else: - raise TypeError() - return args - - def test_warning(self, ellipsoid_class, ellipsoid_args): + def test_warning(self): """ Test warning about ellipsoid with no density being skipped. """ coordinates = (0.0, 0.0, 0.0) - ellipsoid = ellipsoid_class(**ellipsoid_args) + ellipsoid = build_ellipsoid("triaxial") msg = re.escape( f"Ellipsoid {ellipsoid} doesn't have a density value. It will be skipped." @@ -701,7 +596,8 @@ def get_coordinates(self, azimuth=45, polar=45): @pytest.fixture def sphere(self): - return Sphere(self.radius, center=self.center, density=self.density) + a = self.radius + return Ellipsoid(a, a, a, center=self.center, density=self.density) def test_gravity_prolate(self, sphere): """ @@ -714,9 +610,7 @@ def test_gravity_prolate(self, sphere): a = self.radius b = (1 - ratio) * a - ellipsoid = ProlateEllipsoid( - a, b, yaw=0, pitch=0, center=self.center, density=self.density - ) + ellipsoid = Ellipsoid(a, b, b, center=self.center, density=self.density) ge_sphere, gn_sphere, gu_sphere = ellipsoid_gravity(coordinates, sphere) ge_ell, gn_ell, gu_ell = ellipsoid_gravity(coordinates, ellipsoid) @@ -736,9 +630,7 @@ def test_gravity_oblate(self, sphere): a = self.radius b = (1 + ratio) * a - ellipsoid = OblateEllipsoid( - a, b, yaw=0, pitch=0, center=self.center, density=self.density - ) + ellipsoid = Ellipsoid(a, a, b, center=self.center, density=self.density) ge_sphere, gn_sphere, gu_sphere = ellipsoid_gravity(coordinates, sphere) ge_ell, gn_ell, gu_ell = ellipsoid_gravity(coordinates, ellipsoid) @@ -759,9 +651,7 @@ def test_gravity_triaxial(self, sphere): a = self.radius b = (1 - ratio) * a c = (1 - 2 * ratio) * a - ellipsoid = TriaxialEllipsoid( - a, b, c, yaw=0, pitch=0, roll=0, center=self.center, density=self.density - ) + ellipsoid = Ellipsoid(a, b, c, center=self.center, density=self.density) ge_sphere, gn_sphere, gu_sphere = ellipsoid_gravity(coordinates, sphere) ge_ell, gn_ell, gu_ell = ellipsoid_gravity(coordinates, ellipsoid) @@ -769,3 +659,256 @@ def test_gravity_triaxial(self, sphere): np.testing.assert_allclose(ge_ell, ge_sphere, rtol=rtol, atol=atol) np.testing.assert_allclose(gn_ell, gn_sphere, rtol=rtol, atol=atol) np.testing.assert_allclose(gu_ell, gu_sphere, rtol=rtol, atol=atol) + + +class TestTriaxialOnLimits: + """ + Test a triaxial ellipsoid vs oblate and prolate ellipsoids. + + Test the gravity fields of a triaxial ellipsoid that approximates an oblate and + a prolate one against their gravity fields. + """ + + semimajor = 50.0 + atol_ratio = 1e-5 + rtol = 1e-5 + + def get_coordinates(self, azimuth=45, polar=45): + """ + Generate coordinates of observation points along a certain direction. + """ + r = np.linspace(self.semimajor, 5 * self.semimajor, 501) + azimuth, polar = np.rad2deg(azimuth), np.rad2deg(polar) + easting = r * np.cos(azimuth) * np.cos(polar) + northing = r * np.sin(azimuth) * np.cos(polar) + upward = r * np.sin(polar) + return (easting, northing, upward) + + def test_triaxial_vs_prolate(self): + """Compare triaxial with prolate ellipsoid.""" + coordinates = self.get_coordinates() + a, b = self.semimajor, 20.0 + c = b - 1e-4 + density = 200.0 + triaxial = Ellipsoid(a, b, c, density=density) + prolate = Ellipsoid(a, b, b, density=density) + + g_triaxial, g_prolate = tuple( + ellipsoid_gravity(coordinates, ell) for ell in (triaxial, prolate) + ) + + for gi_triaxial, gi_prolate in zip(g_triaxial, g_prolate, strict=True): + atol = self.atol_ratio * vd.maxabs(gi_prolate) + np.testing.assert_allclose( + gi_triaxial, gi_prolate, atol=atol, rtol=self.rtol + ) + + def test_triaxial_vs_oblate(self): + """Compare triaxial with oblate ellipsoid.""" + coordinates = self.get_coordinates() + a = self.semimajor + b = a - 1e-4 + c = 20.0 + density = 200.0 + triaxial = Ellipsoid(a, b, c, density=density) + oblate = Ellipsoid(a, a, c, density=density) + + g_triaxial, g_oblate = tuple( + ellipsoid_gravity(coordinates, ell) for ell in (triaxial, oblate) + ) + + for gi_triaxial, gi_oblate in zip(g_triaxial, g_oblate, strict=True): + atol = self.atol_ratio * vd.maxabs(gi_oblate) + np.testing.assert_allclose( + gi_triaxial, gi_oblate, atol=atol, rtol=self.rtol + ) + + +class TestSemiaxesArbitraryOrder: + """ + Test gravity fields when defining the same ellipsoids with different orders of the + semiaxes plus needed rotation angles. + """ + + atol_ratio = 0.0 + rtol = 1e-7 + + def get_coordinates(self, semimajor, azimuth=45, polar=45): + """ + Generate coordinates of observation points along a certain direction. + """ + r = np.linspace(0.5 * semimajor, 5.5 * semimajor, 501) + azimuth, polar = np.rad2deg(azimuth), np.rad2deg(polar) + easting = r * np.cos(azimuth) * np.cos(polar) + northing = r * np.sin(azimuth) * np.cos(polar) + upward = r * np.sin(polar) + return (easting, northing, upward) + + @pytest.mark.parametrize( + ("semiaxes", "yaw", "pitch", "roll"), + [ + ((30, 20, 10), 0, 0, 0), + ((20, 30, 10), 90, 0, 0), + ((10, 20, 30), 0, 90, 0), + ((30, 10, 20), 0, 0, 90), + ((20, 10, 30), 90, 90, 0), + ((10, 30, 20), 90, 0, 90), + ], + ) + def test_triaxial(self, semiaxes, yaw, pitch, roll): + a, b, c = semiaxes + semiaxes_sorted = sorted(semiaxes, reverse=True) + density = 200.0 + ellipsoid = Ellipsoid(a, b, c, density=density) + ellipsoid_rotated = Ellipsoid( + *semiaxes_sorted, yaw=yaw, pitch=pitch, roll=roll, density=density + ) + + coordinates = self.get_coordinates(semiaxes_sorted[0]) + g_field = ellipsoid_gravity(coordinates, ellipsoid) + g_rotated = ellipsoid_gravity(coordinates, ellipsoid_rotated) + + np.testing.assert_allclose(g_field, g_rotated, rtol=self.rtol) + + @pytest.mark.parametrize( + ("semiaxes", "yaw", "pitch", "roll"), + [ + ((30, 10, 10), 0, 0, 0), + ((10, 30, 10), 90, 0, 0), + ((10, 10, 30), 0, 90, 0), + ], + ) + def test_prolate(self, semiaxes, yaw, pitch, roll): + a, b, c = semiaxes + semiaxes_sorted = sorted(semiaxes, reverse=True) + density = 200.0 + ellipsoid = Ellipsoid(a, b, c, density=density) + ellipsoid_rotated = Ellipsoid( + *semiaxes_sorted, yaw=yaw, pitch=pitch, roll=roll, density=density + ) + + coordinates = self.get_coordinates(semiaxes_sorted[0]) + g_field = ellipsoid_gravity(coordinates, ellipsoid) + g_rotated = ellipsoid_gravity(coordinates, ellipsoid_rotated) + + np.testing.assert_allclose(g_field, g_rotated, rtol=self.rtol) + + @pytest.mark.parametrize( + ("semiaxes", "yaw", "pitch", "roll"), + [ + ((20, 20, 10), 0, 0, 0), + ((20, 10, 20), 0, 0, 90), + ((10, 20, 20), 0, 90, 0), + ], + ) + def test_oblate(self, semiaxes, yaw, pitch, roll): + a, b, c = semiaxes + semiaxes_sorted = sorted(semiaxes, reverse=True) + density = 200.0 + ellipsoid = Ellipsoid(a, b, c, density=density) + ellipsoid_rotated = Ellipsoid( + *semiaxes_sorted, yaw=yaw, pitch=pitch, roll=roll, density=density + ) + + coordinates = self.get_coordinates(semiaxes_sorted[0]) + g_field = ellipsoid_gravity(coordinates, ellipsoid) + g_rotated = ellipsoid_gravity(coordinates, ellipsoid_rotated) + + np.testing.assert_allclose(g_field, g_rotated, rtol=self.rtol) + + +class TestNumericalInstabilitiesTriaxial: + """ + Test fix for numerical instabilities when triaxial approximates a prolate or oblate. + + When two of the three semiaxes of a triaxial ellipsoid are almost equal to each + other (in the order of machine precision), analytic solutions might fall under + singularities, triggering numerical instabilities. + + Given two semiaxes a and b, define a ratio as: + + .. code:: + + ratio = | a - b | / max(a, b) + + These test functions will build the ellipsoids using very small values of this + ratio. + """ + + # Properties of the sphere + semimajor, semiminor = 50.0, 30.0 + center = (0, 0, 0) + density = 200.0 + + def get_coordinates(self, azimuth=45, polar=45): + """ + Generate coordinates of observation points along a certain direction. + """ + r = np.linspace(0.5 * self.semimajor, 5.5 * self.semimajor, 501) + azimuth, polar = np.rad2deg(azimuth), np.rad2deg(polar) + easting = r * np.cos(azimuth) * np.cos(polar) + northing = r * np.sin(azimuth) * np.cos(polar) + upward = r * np.sin(polar) + return (easting, northing, upward) + + @pytest.fixture + def prolate(self): + return Ellipsoid( + self.semimajor, + self.semiminor, + self.semiminor, + center=self.center, + density=self.density, + ) + + @pytest.fixture + def oblate(self): + return Ellipsoid( + self.semimajor, + self.semimajor, + self.semiminor, + center=self.center, + density=self.density, + ) + + def test_triaxial_vs_prolate(self, prolate): + """ + Test gravity field of triaxial ellipsoid that is almost a prolate. + """ + coordinates = self.get_coordinates() + + # Semiaxes ratio. Sufficiently small to trigger the numerical instabilities + ratio = 1e-14 + + a = self.semimajor + b = self.semiminor + c = (1 - ratio) * b + ellipsoid = Ellipsoid(a, b, c, center=self.center, density=self.density) + ge_prolate, gn_prolate, gz_prolate = ellipsoid_gravity(coordinates, prolate) + ge_ell, gn_ell, gz_ell = ellipsoid_gravity(coordinates, ellipsoid) + + rtol, atol = 1e-7, 1e-8 + np.testing.assert_allclose(ge_ell, ge_prolate, rtol=rtol, atol=atol) + np.testing.assert_allclose(gn_ell, gn_prolate, rtol=rtol, atol=atol) + np.testing.assert_allclose(gz_ell, gz_prolate, rtol=rtol, atol=atol) + + def test_triaxial_vs_oblate(self, oblate): + """ + Test gravity field of triaxial ellipsoid that is almost a oblate. + """ + coordinates = self.get_coordinates() + + # Semiaxes ratio. Sufficiently small to trigger the numerical instabilities + ratio = 1e-14 + + a = self.semimajor + b = (1 - ratio) * a + c = self.semiminor + ellipsoid = Ellipsoid(a, b, c, center=self.center, density=self.density) + ge_oblate, gn_oblate, gz_oblate = ellipsoid_gravity(coordinates, oblate) + ge_ell, gn_ell, gz_ell = ellipsoid_gravity(coordinates, ellipsoid) + + rtol, atol = 1e-7, 1e-8 + np.testing.assert_allclose(ge_ell, ge_oblate, rtol=rtol, atol=atol) + np.testing.assert_allclose(gn_ell, gn_oblate, rtol=rtol, atol=atol) + np.testing.assert_allclose(gz_ell, gz_oblate, rtol=rtol, atol=atol) diff --git a/harmonica/tests/ellipsoids/test_magnetic.py b/harmonica/tests/ellipsoids/test_magnetic.py index 8fab168c5..c94e8e4b5 100644 --- a/harmonica/tests/ellipsoids/test_magnetic.py +++ b/harmonica/tests/ellipsoids/test_magnetic.py @@ -24,16 +24,12 @@ import harmonica as hm from harmonica import ellipsoid_magnetic -from harmonica._forward.ellipsoids import ( - OblateEllipsoid, - ProlateEllipsoid, - Sphere, - TriaxialEllipsoid, -) +from harmonica._forward.ellipsoids import Ellipsoid from harmonica._forward.ellipsoids.magnetic import ( get_demagnetization_tensor_internal, get_magnetisation, ) +from harmonica._forward.ellipsoids.utils import SEMIAXES_RTOL from harmonica._forward.utils import get_rotation_matrix from harmonica.errors import NoPhysicalPropertyWarning @@ -47,22 +43,12 @@ def test_euler_returns(): def test_magnetic_symmetry(): """ - Check the symmetry of magentic calculations at surfaces above and below + Check the symmetry of magnetic calculations at surfaces above and below the body. """ - a, b, c = (4, 3, 2) # triaxial ellipsoid - yaw, pitch, roll = (0, 0, 0) - external_field = (10_000, 0, 0) susceptibility = 0.1 - triaxial_example = TriaxialEllipsoid( - a, b, c, yaw, pitch, roll, (0, 0, 0), susceptibility=susceptibility - ) - triaxial_example2 = TriaxialEllipsoid( - a, b, c, yaw, pitch, roll, (0, 0, 0), susceptibility=susceptibility - ) + ellipsoid = Ellipsoid(4, 3, 2, susceptibility=susceptibility) - # define observation points (2D grid) at surface height (z axis, - # 'Upward') = 5 coordinates = vd.grid_coordinates( region=(-20, 20, -20, 20), spacing=0.5, extra_coords=5 ) @@ -70,8 +56,9 @@ def test_magnetic_symmetry(): region=(-20, 20, -20, 20), spacing=0.5, extra_coords=-5 ) - be1, bn1, bu1 = ellipsoid_magnetic(coordinates, triaxial_example, external_field) - be2, bn2, bu2 = ellipsoid_magnetic(coordinates2, triaxial_example2, external_field) + external_field = (10_000, 0, 0) + be1, bn1, bu1 = ellipsoid_magnetic(coordinates, ellipsoid, external_field) + be2, bn2, bu2 = ellipsoid_magnetic(coordinates2, ellipsoid, external_field) np.testing.assert_allclose(np.abs(be1), np.flip(np.abs(be2))) np.testing.assert_allclose(np.abs(bn1), np.flip(np.abs(bn2))) @@ -83,25 +70,18 @@ def test_flipped_h0(): Check that reversing the magentising field produces the same (reversed) field. """ - - a, b = (2, 4) # triaxial ellipsoid - yaw, pitch = 0, 0 - - external_field1 = np.array((55_000, 0.0, 90.0)) - external_field2 = -external_field1 + a, c = (2, 4) susceptibility = 0.1 - oblate_example = OblateEllipsoid( - a, b, yaw, pitch, (0, 0, 0), susceptibility=susceptibility - ) + oblate = Ellipsoid(a, a, c, susceptibility=susceptibility) - # define observation points (2D grid) at surface height (z axis, - # 'Upward') = 5 coordinates = vd.grid_coordinates( region=(-20, 20, -20, 20), spacing=0.5, extra_coords=5 ) - be1, bn1, bu1 = ellipsoid_magnetic(coordinates, oblate_example, external_field1) - be2, bn2, bu2 = ellipsoid_magnetic(coordinates, oblate_example, external_field2) + external_field1 = np.array((55_000, 0.0, 90.0)) + external_field2 = -external_field1 + be1, bn1, bu1 = ellipsoid_magnetic(coordinates, oblate, external_field1) + be2, bn2, bu2 = ellipsoid_magnetic(coordinates, oblate, external_field2) np.testing.assert_allclose(np.abs(be1), np.abs(be2)) np.testing.assert_allclose(np.abs(bn1), np.abs(bn2)) @@ -112,12 +92,8 @@ def test_zero_susceptibility(): """ Test for the case of 0 susceptibility == inducing field. """ - - a, b = 1, 2 susceptibility = 0 - ellipsoid = OblateEllipsoid( - a, b, yaw=0, pitch=0, center=(0, 0, 0), susceptibility=susceptibility - ) + ellipsoid = Ellipsoid(2, 2, 1, susceptibility=susceptibility) coordinates = vd.grid_coordinates( region=(-10, 10, -10, 10), spacing=1.0, extra_coords=5 ) @@ -134,18 +110,13 @@ def test_zero_field(): """ Test that zero field produces zero anomalies. """ - - a, b = 1, 2 - external_field = np.array([0, 0, 0]) susceptibility = 0.01 - - ellipsoid = OblateEllipsoid( - a, b, yaw=0, pitch=0, center=(0, 0, 0), susceptibility=susceptibility - ) + ellipsoid = Ellipsoid(2, 1, 1, susceptibility=susceptibility) coordinates = vd.grid_coordinates( region=(-10, 10, -10, 10), spacing=1.0, extra_coords=5 ) + external_field = (0, 0, 0) be, bn, bu = ellipsoid_magnetic(coordinates, ellipsoid, external_field) np.testing.assert_allclose(be[0], 0) @@ -159,13 +130,11 @@ def test_mag_ext_int_boundary(): consistent. """ - a, b = 50, 60 + a, c = 50, 60 external_field = (55_000, 0.0, 90.0) susceptibility = 0.01 - ellipsoid = OblateEllipsoid( - a, b, yaw=0, pitch=0, center=(0, 0, 0), susceptibility=susceptibility - ) + ellipsoid = Ellipsoid(a, a, c, susceptibility=susceptibility) e = np.array([49.99, 50.00]) n = np.array([0.0, 0.0]) @@ -187,10 +156,10 @@ def test_mag_flipped_ellipsoid(): external_field = (10_000, 0, 0) susceptibility = 0.01 - triaxial_example = TriaxialEllipsoid( + triaxial_example = Ellipsoid( a, b, c, yaw=0, pitch=0, roll=0, center=(0, 0, 0), susceptibility=susceptibility ) - triaxial_example2 = TriaxialEllipsoid( + triaxial_example2 = Ellipsoid( a, b, c, @@ -245,68 +214,27 @@ def check_rotation_equivalence(base_ellipsoid, rotated_ellipsoids): np.testing.assert_allclose(np.abs(bu), np.abs(base_bu), rtol=1e-4) # triaxial cases - base_tri = TriaxialEllipsoid( - a, b, c, yaw=0, pitch=0, roll=0, center=(0, 0, 0), susceptibility=susceptibility - ) + base_tri = Ellipsoid(a, b, c, susceptibility=susceptibility) tri_rotated = [ - TriaxialEllipsoid( - a, - b, - c, - yaw=360, - pitch=0, - roll=0, - center=(0, 0, 0), - susceptibility=susceptibility, - ), - TriaxialEllipsoid( - a, - b, - c, - yaw=0, - pitch=180, - roll=0, - center=(0, 0, 0), - susceptibility=susceptibility, - ), - TriaxialEllipsoid( - a, - b, - c, - yaw=0, - pitch=360, - roll=360, - center=(0, 0, 0), - susceptibility=susceptibility, - ), + Ellipsoid(a, b, c, yaw=360, susceptibility=susceptibility), + Ellipsoid(a, b, c, pitch=180, susceptibility=susceptibility), + Ellipsoid(a, b, c, pitch=360, roll=360, susceptibility=susceptibility), ] check_rotation_equivalence(base_tri, tri_rotated) # prolate cases - base_pro = ProlateEllipsoid( - a, b, yaw=0, pitch=0, center=(0, 0, 0), susceptibility=susceptibility - ) + base_pro = Ellipsoid(a, b, b, susceptibility=susceptibility) pro_rotated = [ - ProlateEllipsoid( - a, b, yaw=360, pitch=0, center=(0, 0, 0), susceptibility=susceptibility - ), - ProlateEllipsoid( - a, b, yaw=0, pitch=180, center=(0, 0, 0), susceptibility=susceptibility - ), + Ellipsoid(a, b, b, yaw=360, susceptibility=susceptibility), + Ellipsoid(a, b, b, pitch=180, susceptibility=susceptibility), ] check_rotation_equivalence(base_pro, pro_rotated) # oblate cases - base_obl = OblateEllipsoid( - b, a, yaw=0, pitch=0, center=(0, 0, 0), susceptibility=susceptibility - ) + base_obl = Ellipsoid(a, a, c, susceptibility=susceptibility) obl_rotated = [ - OblateEllipsoid( - b, a, yaw=360, pitch=0, center=(0, 0, 0), susceptibility=susceptibility - ), - OblateEllipsoid( - b, a, yaw=0, pitch=180, center=(0, 0, 0), susceptibility=susceptibility - ), + Ellipsoid(a, a, c, yaw=360, susceptibility=susceptibility), + Ellipsoid(a, a, c, pitch=180, susceptibility=susceptibility), ] check_rotation_equivalence(base_obl, obl_rotated) @@ -321,11 +249,11 @@ def ellipsoid_semiaxes(self, request): ellipsoid_type = request.param match ellipsoid_type: case "oblate": - a, b = 50.0, 60.0 - c = b + a = b = 60.0 + c = 50.0 case "prolate": - a, b = 60.0, 50.0 - c = b + a = 60.0 + b = c = 50.0 case "triaxial": a, b, c = 70.0, 60.0, 50.0 case _: @@ -361,7 +289,7 @@ def test_demagnetization(self, ellipsoid_semiaxes): @pytest.mark.parametrize( ("a", "b", "c"), [ - (50.0, 60.0, 60.0), # oblate + (60.0, 60.0, 50.0), # oblate (60.0, 50.0, 50.0), # prolate (70.0, 60.0, 50.0), # triaxial (70.0, 70.0, 70.0), # sphere @@ -411,6 +339,12 @@ class TestMagnetizationVersusSphere: Test if ellipsoid's magnetization approximates the one of the sphere. """ + # Ratio between ellipsoid's semiaxes, defined as ratio = | a - b | / max(a, b). + # Make it small enough so ellipsoids approximate a sphere, but not too small that + # might trigger the fixes for numerical uncertainties. + ratio = 1e-4 + assert ratio > SEMIAXES_RTOL + @pytest.fixture def radius(self): """ @@ -423,16 +357,18 @@ def ellipsoid_semiaxes(self, radius, request): """ Ellipsoid's semiaxes that approximate a sphere. """ - a = radius ellipsoid_type = request.param match ellipsoid_type: case "oblate": - b = c = a + 1e-2 + a = b = radius + c = (1 - self.ratio) * b case "prolate": - b = c = a - 1e-2 + a = radius + b = c = (1 - self.ratio) * a case "triaxial": - b = a - 1e-3 - c = a - 1e-2 + a = radius + b = (1 - self.ratio) * a + c = (1 - 2 * self.ratio) * a case _: raise ValueError() return a, b, c @@ -472,15 +408,15 @@ class TestMagneticFieldVersusSphere: Test if magnetic field of ellipsoid approximates the one of the sphere. """ - # Sphere radius, center, and susceptibility. + # Sphere radius and susceptibility. radius = 50.0 - center = (0, 0, 0) susceptibility = 0.5 - # Difference between ellipsoid's semiaxes. - # It should be small compared to the sphere radius, so the ellipsoid approximates - # a sphere. - delta = 0.001 + # Ratio between ellipsoid's semiaxes, defined as ratio = | a - b | / max(a, b). + # Make it small enough so ellipsoids approximate a sphere, but not too small that + # might trigger the fixes for numerical uncertainties. + ratio = 1e-4 + assert ratio > SEMIAXES_RTOL # Define external field external_field = (55_123.0, 32.0, -28.9) @@ -505,42 +441,22 @@ def get_ellipsoid(self, ellipsoid_type: str): Returns ------- - OblateEllipsoid, ProlateEllipsoid or TriaxialEllipsoid + Ellipsoid """ - yaw, pitch, roll = 0, 0, 0 - a = self.radius match ellipsoid_type: case "oblate": - ellipsoid = OblateEllipsoid( - a=a, - b=a + self.delta, - yaw=yaw, - pitch=pitch, - center=self.center, - susceptibility=self.susceptibility, - ) + a = b = self.radius + c = (1 - self.ratio) * b case "prolate": - ellipsoid = ProlateEllipsoid( - a=a, - b=a - self.delta, - yaw=yaw, - pitch=pitch, - center=self.center, - susceptibility=self.susceptibility, - ) + a = self.radius + b = c = (1 - self.ratio) * a case "triaxial": - ellipsoid = TriaxialEllipsoid( - a=a, - b=a - self.delta, - c=a - 2 * self.delta, - yaw=yaw, - pitch=pitch, - roll=roll, - center=self.center, - susceptibility=self.susceptibility, - ) + a = self.radius + b = (1 - self.ratio) * a + c = (1 - 2 * self.ratio) * a case _: raise ValueError() + ellipsoid = Ellipsoid(a, b, c, susceptibility=self.susceptibility) return ellipsoid @pytest.mark.parametrize("ellipsoid_type", ["oblate", "prolate", "triaxial"]) @@ -548,8 +464,8 @@ def test_magnetic_field_vs_sphere(self, coordinates, ellipsoid_type): """ Test magnetic field of ellipsoids against the one for a sphere. """ - sphere = Sphere( - self.radius, center=self.center, susceptibility=self.susceptibility + sphere = Ellipsoid( + self.radius, self.radius, self.radius, susceptibility=self.susceptibility ) b_sphere = ellipsoid_magnetic(coordinates, sphere, self.external_field) @@ -559,7 +475,7 @@ def test_magnetic_field_vs_sphere(self, coordinates, ellipsoid_type): rtol = 5e-4 for bi_sphere, bi_ellipsoid in zip(b_sphere, b_ellipsoid, strict=True): maxabs = vd.maxabs(bi_sphere, bi_ellipsoid) - atol = maxabs * 1e-4 + atol = maxabs * 5e-4 np.testing.assert_allclose(bi_sphere, bi_ellipsoid, atol=atol, rtol=rtol) @@ -573,10 +489,11 @@ class TestMagneticFieldVersusDipole: center = (0, 0, 0) susceptibility = 0.001 # use small sus to reduce demag effects - # Difference between ellipsoid's semiaxes. - # It should be small compared to the sphere radius, so the ellipsoid approximates - # a sphere. - delta = 0.001 + # Ratio between ellipsoid's semiaxes, defined as ratio = | a - b | / max(a, b). + # Make it small enough so ellipsoids approximate a sphere, but not too small that + # might trigger the fixes for numerical uncertainties. + ratio = 1e-4 + assert ratio > SEMIAXES_RTOL # Define external field external_field = (55_123.0, 32.0, -28.9) @@ -601,46 +518,24 @@ def get_ellipsoid(self, ellipsoid_type: str): Returns ------- - OblateEllipsoid, ProlateEllipsoid, Sphere or TriaxialEllipsoid + Ellipsoid """ - yaw, pitch, roll = 0, 0, 0 - a = self.radius match ellipsoid_type: case "oblate": - ellipsoid = OblateEllipsoid( - a=a, - b=a + self.delta, - yaw=yaw, - pitch=pitch, - center=self.center, - susceptibility=self.susceptibility, - ) + a = b = self.radius + c = (1 - self.ratio) * b case "prolate": - ellipsoid = ProlateEllipsoid( - a=a, - b=a - self.delta, - yaw=yaw, - pitch=pitch, - center=self.center, - susceptibility=self.susceptibility, - ) + a = self.radius + b = c = (1 - self.ratio) * a case "triaxial": - ellipsoid = TriaxialEllipsoid( - a=a, - b=a - self.delta, - c=a - 2 * self.delta, - yaw=yaw, - pitch=pitch, - roll=roll, - center=self.center, - susceptibility=self.susceptibility, - ) + a = self.radius + b = (1 - self.ratio) * a + c = (1 - 2 * self.ratio) * a case "sphere": - ellipsoid = Sphere( - a=a, center=self.center, susceptibility=self.susceptibility - ) + a = b = c = self.radius case _: raise ValueError() + ellipsoid = Ellipsoid(a, b, c, susceptibility=self.susceptibility) return ellipsoid def get_dipole_moment(self, ellipsoid): @@ -674,8 +569,8 @@ def test_magnetic_field_vs_dipole(self, coordinates, ellipsoid_type): rtol = 5e-4 for bi_dipole, bi_ellipsoid in zip(b_dipole, b_ellipsoid, strict=True): maxabs = vd.maxabs(bi_dipole, bi_ellipsoid) - atol = maxabs * 1e-4 - np.testing.assert_allclose(bi_dipole, bi_ellipsoid, atol=atol, rtol=rtol) + atol = maxabs * 5e-4 + np.testing.assert_allclose(bi_dipole, bi_ellipsoid, rtol=rtol, atol=atol) class TestSymmetryOnRotations: @@ -692,8 +587,7 @@ def flip_ellipsoid(self, ellipsoid): """ ellipsoid.yaw += 180 ellipsoid.pitch *= -1 - if isinstance(ellipsoid, TriaxialEllipsoid): - ellipsoid.roll *= -1 + ellipsoid.roll *= -1 return ellipsoid @pytest.fixture(params=["oblate", "prolate", "triaxial"]) @@ -706,25 +600,17 @@ def ellipsoid(self, request): yaw, pitch, roll = 62.3, 48.2, 14.9 match ellipsoid_type: case "oblate": - ellipsoid = OblateEllipsoid( - a=semiminor, b=semimajor, yaw=yaw, pitch=pitch, center=center - ) + a = b = semimajor + c = semiminor case "prolate": - ellipsoid = ProlateEllipsoid( - a=semimajor, b=semiminor, yaw=yaw, pitch=pitch, center=center - ) + a = semimajor + b = c = semiminor case "triaxial": - ellipsoid = TriaxialEllipsoid( - a=semimajor, - b=semimiddle, - c=semiminor, - yaw=yaw, - pitch=pitch, - roll=roll, - center=center, - ) + a, b, c = semimajor, semimiddle, semiminor case _: raise ValueError() + + ellipsoid = Ellipsoid(a, b, c, yaw=yaw, pitch=pitch, roll=roll, center=center) return ellipsoid @pytest.mark.parametrize("magnetization_type", ["induced", "remanent", "both"]) @@ -786,21 +672,23 @@ def coordinates(self): def ellipsoids(self): """Sample ellipsoids.""" ellipsoids = [ - OblateEllipsoid( - a=20, + Ellipsoid( + a=60, b=60, + c=20, yaw=30.2, pitch=-23, center=(-10.0, 20.0, -10.0), ), - ProlateEllipsoid( + Ellipsoid( a=40, b=15, + c=15, yaw=170.2, pitch=71, center=(15.0, 0.0, -40.0), ), - TriaxialEllipsoid( + Ellipsoid( a=60, b=18, c=15, @@ -896,52 +784,16 @@ def test_multiple_ellipsoids_remanence(self, coordinates, ellipsoids): np.testing.assert_allclose(bz, bz_expected) -@pytest.mark.parametrize( - "ellipsoid_class", [OblateEllipsoid, ProlateEllipsoid, TriaxialEllipsoid, Sphere] -) class TestNoMagnetic: """Test warning when ellipsoid has no susceptibility nor remanent magnetization.""" - @pytest.fixture - def ellipsoid_args(self, ellipsoid_class): - if ellipsoid_class is OblateEllipsoid: - args = { - "a": 20.0, - "b": 50.0, - "pitch": 0.0, - "yaw": 0.0, - "center": (0, 0, 0), - } - elif ellipsoid_class is ProlateEllipsoid: - args = { - "a": 50.0, - "b": 20.0, - "pitch": 0.0, - "yaw": 0.0, - "center": (0, 0, 0), - } - elif ellipsoid_class is TriaxialEllipsoid: - args = { - "a": 50.0, - "b": 20.0, - "c": 10.0, - "pitch": 0.0, - "yaw": 0.0, - "roll": 0.0, - "center": (0, 0, 0), - } - elif ellipsoid_class is Sphere: - args = {"a": 50.0, "center": (0, 0, 0)} - else: - raise TypeError() - return args - - def test_warning(self, ellipsoid_class, ellipsoid_args): + def test_warning(self): """ Test warning about ellipsoid with no susceptibility nor remanence being skipped. """ coordinates = (0.0, 0.0, 0.0) - ellipsoid = ellipsoid_class(**ellipsoid_args) + a, b, c = 50.0, 35.0, 25.0 + ellipsoid = Ellipsoid(a, b, c) external_field = (55_000.0, 13, 71) msg = re.escape( @@ -951,6 +803,312 @@ def test_warning(self, ellipsoid_class, ellipsoid_args): with pytest.warns(NoPhysicalPropertyWarning, match=msg): bx, by, bz = ellipsoid_magnetic(coordinates, ellipsoid, external_field) - # Check the gravity acceleration components are zero + # Check the magnetic field components are zero for b_component in (bx, by, bz): assert b_component == 0.0 + + +class TestTriaxialOnLimits: + """ + Test a triaxial ellipsoid vs oblate and prolate ellipsoids. + + Test the magnetic fields of a triaxial ellipsoid that approximates an oblate and + a prolate one against their magnetic fields. + """ + + semimajor = 50.0 + external_field = (55_000, 71, 18) + atol_ratio = 1e-6 + rtol = 1e-5 + + def get_coordinates(self, azimuth=45, polar=45): + """ + Generate coordinates of observation points along a certain direction. + """ + r = np.linspace(self.semimajor, 5 * self.semimajor, 501) + azimuth, polar = np.rad2deg(azimuth), np.rad2deg(polar) + easting = r * np.cos(azimuth) * np.cos(polar) + northing = r * np.sin(azimuth) * np.cos(polar) + upward = r * np.sin(polar) + return (easting, northing, upward) + + def test_triaxial_vs_prolate(self): + """Compare triaxial with prolate ellipsoid.""" + coordinates = self.get_coordinates() + a, b = self.semimajor, 20.0 + c = b - 1e-4 + susceptibility = 0.2 + triaxial = Ellipsoid(a, b, c, susceptibility=susceptibility) + prolate = Ellipsoid(a, b, b, susceptibility=susceptibility) + + b_triaxial, b_prolate = tuple( + ellipsoid_magnetic(coordinates, ell, external_field=self.external_field) + for ell in (triaxial, prolate) + ) + + for bi_triaxial, bi_prolate in zip(b_triaxial, b_prolate, strict=True): + atol = self.atol_ratio * vd.maxabs(bi_prolate) + np.testing.assert_allclose( + bi_triaxial, bi_prolate, atol=atol, rtol=self.rtol + ) + + def test_triaxial_vs_oblate(self): + """Compare triaxial with oblate ellipsoid.""" + coordinates = self.get_coordinates() + a = self.semimajor + b = a - 1e-4 + c = 20.0 + susceptibility = 0.2 + triaxial = Ellipsoid(a, b, c, susceptibility=susceptibility) + oblate = Ellipsoid(a, a, c, susceptibility=susceptibility) + + b_triaxial, b_oblate = tuple( + ellipsoid_magnetic(coordinates, ell, external_field=self.external_field) + for ell in (triaxial, oblate) + ) + + for bi_triaxial, bi_oblate in zip(b_triaxial, b_oblate, strict=True): + atol = self.atol_ratio * vd.maxabs(bi_oblate) + np.testing.assert_allclose( + bi_triaxial, bi_oblate, atol=atol, rtol=self.rtol + ) + + +class TestSemiaxesArbitraryOrder: + """ + Test magnetic fields when defining the same ellipsoids with different orders of the + semiaxes plus needed rotation angles. + """ + + atol_ratio = 0.0 + rtol = 1e-7 + external_field = (55_000, 17, -21) + + def get_coordinates(self, semimajor, azimuth=45, polar=45): + """ + Generate coordinates of observation points along a certain direction. + """ + r = np.linspace(0.5 * semimajor, 5.5 * semimajor, 501) + azimuth, polar = np.rad2deg(azimuth), np.rad2deg(polar) + easting = r * np.cos(azimuth) * np.cos(polar) + northing = r * np.sin(azimuth) * np.cos(polar) + upward = r * np.sin(polar) + return (easting, northing, upward) + + @pytest.mark.parametrize( + ("semiaxes", "yaw", "pitch", "roll"), + [ + ((30, 20, 10), 0, 0, 0), + ((20, 30, 10), 90, 0, 0), + ((10, 20, 30), 0, 90, 0), + ((30, 10, 20), 0, 0, 90), + ((20, 10, 30), 90, 90, 0), + ((10, 30, 20), 90, 0, 90), + ], + ) + def test_triaxial(self, semiaxes, yaw, pitch, roll): + a, b, c = semiaxes + semiaxes_sorted = sorted(semiaxes, reverse=True) + susceptibility = 0.2 + rem_mag = np.array([1.0, 2.0, -3.0]) + ellipsoid = Ellipsoid( + a, b, c, susceptibility=susceptibility, remanent_mag=rem_mag + ) + ellipsoid_rotated = Ellipsoid( + *semiaxes_sorted, + yaw=yaw, + pitch=pitch, + roll=roll, + susceptibility=susceptibility, + remanent_mag=rem_mag, + ) + + coordinates = self.get_coordinates(semiaxes_sorted[0]) + g_field = ellipsoid_magnetic( + coordinates, ellipsoid, external_field=self.external_field + ) + g_rotated = ellipsoid_magnetic( + coordinates, ellipsoid_rotated, external_field=self.external_field + ) + + np.testing.assert_allclose(g_field, g_rotated, rtol=self.rtol) + + @pytest.mark.parametrize( + ("semiaxes", "yaw", "pitch", "roll"), + [ + ((30, 10, 10), 0, 0, 0), + ((10, 30, 10), 90, 0, 0), + ((10, 10, 30), 0, 90, 0), + ], + ) + def test_prolate(self, semiaxes, yaw, pitch, roll): + a, b, c = semiaxes + semiaxes_sorted = sorted(semiaxes, reverse=True) + susceptibility = 0.2 + rem_mag = np.array([1.0, 2.0, -3.0]) + ellipsoid = Ellipsoid( + a, b, c, susceptibility=susceptibility, remanent_mag=rem_mag + ) + ellipsoid_rotated = Ellipsoid( + *semiaxes_sorted, + yaw=yaw, + pitch=pitch, + roll=roll, + susceptibility=susceptibility, + remanent_mag=rem_mag, + ) + + coordinates = self.get_coordinates(semiaxes_sorted[0]) + g_field = ellipsoid_magnetic( + coordinates, ellipsoid, external_field=self.external_field + ) + g_rotated = ellipsoid_magnetic( + coordinates, ellipsoid_rotated, external_field=self.external_field + ) + + np.testing.assert_allclose(g_field, g_rotated, rtol=self.rtol) + + @pytest.mark.parametrize( + ("semiaxes", "yaw", "pitch", "roll"), + [ + ((20, 20, 10), 0, 0, 0), + ((20, 10, 20), 0, 0, 90), + ((10, 20, 20), 0, 90, 0), + ], + ) + def test_oblate(self, semiaxes, yaw, pitch, roll): + a, b, c = semiaxes + semiaxes_sorted = sorted(semiaxes, reverse=True) + susceptibility = 0.2 + rem_mag = np.array([1.0, 2.0, -3.0]) + ellipsoid = Ellipsoid( + a, b, c, susceptibility=susceptibility, remanent_mag=rem_mag + ) + ellipsoid_rotated = Ellipsoid( + *semiaxes_sorted, + yaw=yaw, + pitch=pitch, + roll=roll, + susceptibility=susceptibility, + remanent_mag=rem_mag, + ) + + coordinates = self.get_coordinates(semiaxes_sorted[0]) + g_field = ellipsoid_magnetic( + coordinates, ellipsoid, external_field=self.external_field + ) + g_rotated = ellipsoid_magnetic( + coordinates, ellipsoid_rotated, external_field=self.external_field + ) + + np.testing.assert_allclose(g_field, g_rotated, rtol=self.rtol) + + +class TestNumericalInstabilitiesTriaxial: + """ + Test fix for numerical instabilities when triaxial approximates a prolate or oblate. + + When two of the three semiaxes of a triaxial ellipsoid are almost equal to each + other (in the order of machine precision), analytic solutions might fall under + singularities, triggering numerical instabilities. + + Given two semiaxes a and b, define a ratio as: + + .. code:: + + ratio = | a - b | / max(a, b) + + These test functions will build the ellipsoids using very small values of this + ratio. + """ + + # Properties of the sphere + semimajor, semiminor = 50.0, 30.0 + center = (0, 0, 0) + susceptibility = 0.1 + external_field = (55_000, 17, 21) + + def get_coordinates(self, azimuth=45, polar=45): + """ + Generate coordinates of observation points along a certain direction. + """ + r = np.linspace(1e-3 * self.semimajor, 5.5 * self.semimajor, 501) + azimuth, polar = np.rad2deg(azimuth), np.rad2deg(polar) + easting = r * np.cos(azimuth) * np.cos(polar) + northing = r * np.sin(azimuth) * np.cos(polar) + upward = r * np.sin(polar) + return (easting, northing, upward) + + @pytest.fixture + def prolate(self): + return Ellipsoid( + self.semimajor, + self.semiminor, + self.semiminor, + center=self.center, + susceptibility=self.susceptibility, + ) + + @pytest.fixture + def oblate(self): + return Ellipsoid( + self.semimajor, + self.semimajor, + self.semiminor, + center=self.center, + susceptibility=self.susceptibility, + ) + + def test_triaxial_vs_prolate(self, prolate): + """ + Test magnetic field of triaxial ellipsoid that is almost a prolate. + """ + coordinates = self.get_coordinates() + + # Semiaxes ratio. Sufficiently small to trigger the numerical instabilities + ratio = 1e-14 + + a = self.semimajor + b = self.semiminor + c = (1 - ratio) * b + ellipsoid = Ellipsoid( + a, b, c, center=self.center, susceptibility=self.susceptibility + ) + be_prolate, bn_prolate, bu_prolate = ellipsoid_magnetic( + coordinates, prolate, self.external_field + ) + be_ell, bn_ell, bu_ell = ellipsoid_magnetic( + coordinates, ellipsoid, self.external_field + ) + + rtol, atol = 1e-7, 1e-8 + np.testing.assert_allclose(be_ell, be_prolate, rtol=rtol, atol=atol) + np.testing.assert_allclose(bn_ell, bn_prolate, rtol=rtol, atol=atol) + np.testing.assert_allclose(bu_ell, bu_prolate, rtol=rtol, atol=atol) + + def test_triaxial_vs_oblate(self, oblate): + """ + Test magnetic field of triaxial ellipsoid that is almost a oblate. + """ + coordinates = self.get_coordinates() + + # Semiaxes ratio. Sufficiently small to trigger the numerical instabilities + ratio = 1e-14 + + a = self.semimajor + b = (1 - ratio) * a + c = self.semiminor + ellipsoid = Ellipsoid( + a, b, c, center=self.center, susceptibility=self.susceptibility + ) + be_oblate, bn_oblate, bu_oblate = ellipsoid_magnetic( + coordinates, oblate, self.external_field + ) + be_ell, bn_ell, bu_ell = ellipsoid_magnetic( + coordinates, ellipsoid, self.external_field + ) + + rtol, atol = 1e-7, 1e-8 + np.testing.assert_allclose(be_ell, be_oblate, rtol=rtol, atol=atol) + np.testing.assert_allclose(bn_ell, bn_oblate, rtol=rtol, atol=atol) + np.testing.assert_allclose(bu_ell, bu_oblate, rtol=rtol, atol=atol) diff --git a/harmonica/tests/ellipsoids/test_utils.py b/harmonica/tests/ellipsoids/test_utils.py index 914427083..3a97957cb 100644 --- a/harmonica/tests/ellipsoids/test_utils.py +++ b/harmonica/tests/ellipsoids/test_utils.py @@ -4,15 +4,23 @@ # # This code is part of the Fatiando a Terra project (https://www.fatiando.org) # +import itertools +import re + import numpy as np import pytest -from harmonica._forward.ellipsoids.utils import calculate_lambda +from harmonica import Ellipsoid +from harmonica._forward.ellipsoids.utils import ( + calculate_lambda, + get_semiaxes_rotation_matrix, + is_almost_a_sphere, +) @pytest.mark.parametrize( ("a", "b", "c"), - [(6.0, 5.0, 5.0), (4.0, 5.0, 5.0), (6.0, 5.0, 4.0)], + [(6.0, 5.0, 5.0), (5.0, 5.0, 4.0), (6.0, 5.0, 4.0)], ids=["prolate", "oblate", "triaxial"], ) def test_lambda(a, b, c): @@ -39,9 +47,20 @@ def test_lambda(a, b, c): ) +def test_lambda_unsorted(): + """ + Test error if semiaxes are not sorted. + """ + a, b, c = 1.0, 2.0, 3.0 + x, y, z = np.meshgrid(*[np.linspace(-10, 10, 41) for _ in range(3)]) + msg = re.escape("Invalid semiaxes not properly sorted") + with pytest.raises(ValueError, match=msg): + calculate_lambda(x, y, z, a, b, c) + + @pytest.mark.parametrize( ("a", "b", "c"), - [(3.0, 2.0, 2.0), (1.0, 2.0, 2.0), (3.0, 2.0, 1.0)], + [(3.0, 2.0, 2.0), (2.0, 2.0, 1.0), (3.0, 2.0, 1.0)], ids=["prolate", "oblate", "triaxial"], ) def test_zero_cases_for_lambda(a, b, c): @@ -50,25 +69,25 @@ def test_zero_cases_for_lambda(a, b, c): may throw zero errors. """ x, y, z = 0, 0, 5 - lmbda = calculate_lambda(x, y, z, a, b, c) - lmbda_eqn = z**2 - c**2 - np.testing.assert_allclose(lmbda, lmbda_eqn) + lambda_ = calculate_lambda(x, y, z, a, b, c) + lambda_eqn = z**2 - c**2 + np.testing.assert_allclose(lambda_, lambda_eqn) x, y, z = 0, 5, 0 - lmbda = calculate_lambda(x, y, z, a, b, c) - lmbda_eqn = y**2 - b**2 - np.testing.assert_allclose(lmbda, lmbda_eqn) + lambda_ = calculate_lambda(x, y, z, a, b, c) + lambda_eqn = y**2 - b**2 + np.testing.assert_allclose(lambda_, lambda_eqn) x, y, z = 5, 0, 0 - lmbda = calculate_lambda(x, y, z, a, b, c) - lmbda_eqn = x**2 - a**2 - np.testing.assert_allclose(lmbda, lmbda_eqn) + lambda_ = calculate_lambda(x, y, z, a, b, c) + lambda_eqn = x**2 - a**2 + np.testing.assert_allclose(lambda_, lambda_eqn) @pytest.mark.parametrize("zero_coord", ["x", "y", "z"]) @pytest.mark.parametrize( ("a", "b", "c"), - [(3.4, 2.2, 2.2), (1.1, 2.8, 2.8), (3.4, 2.2, 1.1)], + [(3.4, 2.2, 2.2), (2.8, 2.8, 1.1), (3.4, 2.2, 1.1)], ids=["prolate", "oblate", "triaxial"], ) def test_second_order_equations(a, b, c, zero_coord): @@ -94,3 +113,70 @@ def test_second_order_equations(a, b, c, zero_coord): p0 = b**2 * c**2 - y**2 * c**2 - z**2 * b**2 expected_lambda = 0.5 * (np.sqrt(p1**2 - 4 * p0) - p1) np.testing.assert_allclose(expected_lambda, lambda_) + + +class TestSemiaxesRotationMatrix: + """ + Test the ``get_semiaxes_rotation_matrix`` function. + """ + + @pytest.mark.parametrize(("a", "b", "c"), itertools.permutations((3.0, 2.0, 1.0))) + def test_sorting_of_semiaxes(self, a, b, c): + """ + Test if the transposed matrix correctly sorts the semiaxes. + + When we apply the transposed matrix to unsorted semiaxes, we should recover the + semiaxes in the right order. Although some of them might have a minus sign due + to the rotation, so we'll compare with absolute values. + + For example: + + .. code::python + + matrix.T @ np.array([1.0, 2.0, 3.0]) + + Should return: + + .. code::python + + np.array([3.0, 2.0, -1.0]) + """ + ellipsoid = Ellipsoid(a, b, c) + matrix = get_semiaxes_rotation_matrix(ellipsoid) + semiaxes = np.array([a, b, c]) + expected = sorted((a, b, c), reverse=True) + np.testing.assert_allclose(np.abs(matrix.T @ semiaxes), expected) + + def test_example(self): + a, b, c = 1, 2, 3 + ellipsoid = Ellipsoid(a, b, c) + matrix = get_semiaxes_rotation_matrix(ellipsoid) + semiaxes = np.array([a, b, c]) + expected = [c, b, -a] + np.testing.assert_allclose(matrix.T @ semiaxes, expected) + + +@pytest.mark.parametrize( + ("a", "b", "c", "expected"), + [ + (1, 1, 1, True), # exact sphere + (1.00001, 1, 1, True), # prolate as sphere + (1.00001, 1.00001, 1, True), # oblate as sphere + (3.00002, 3.00001, 3, True), # triaxial as sphere + (2, 1, 1, False), # non-spherical prolate + (2, 2, 1, False), # non-spherical oblate + (3, 2, 1, False), # non-spherical triaxial + ], + ids=[ + "exact sphere", + "prolate as sphere", + "oblate as sphere", + "triaxial as sphere", + "non-spherical prolate", + "non-spherical oblate", + "non-spherical triaxial", + ], +) +def test_is_almost_a_sphere(a, b, c, expected): + """Test the ``is_almost_a_sphere`` function.""" + assert is_almost_a_sphere(a, b, c) == expected diff --git a/harmonica/typing.py b/harmonica/typing.py index f554e29da..a5c81dbcc 100644 --- a/harmonica/typing.py +++ b/harmonica/typing.py @@ -13,68 +13,8 @@ import numpy as np import numpy.typing as npt -from typing_extensions import Protocol Coordinates: TypeAlias = Sequence[npt.NDArray[np.float64]] | Sequence[float] """ Type alias to represent 3D coordinates. """ - - -class Ellipsoid(Protocol): - """Protocol to define an Ellipsoid object.""" - - @property - def a(self) -> float: - """Length of the first semiaxis.""" - raise NotImplementedError - - @property - def b(self) -> float: - """Length of the second semiaxis.""" - raise NotImplementedError - - @property - def c(self) -> float: - """Length of the third semiaxis.""" - raise NotImplementedError - - @property - def yaw(self) -> float: - """Yaw angle in degrees.""" - raise NotImplementedError - - @property - def pitch(self) -> float: - """Pitch angle in degrees.""" - raise NotImplementedError - - @property - def roll(self) -> float: - """Roll angle in degrees.""" - raise NotImplementedError - - @property - def center(self) -> tuple[float, float, float]: - """Coordinates of ellipsoid's center.""" - raise NotImplementedError - - @property - def density(self) -> float | None: - """Density of the ellipsoid in :math:`kg/m^3`.""" - raise NotImplementedError - - @property - def susceptibility(self) -> float | npt.NDArray | None: - """Magnetic susceptibility of the ellipsoid in SI units.""" - raise NotImplementedError - - @property - def remanent_mag(self) -> npt.NDArray | None: - """Remanent magnetization of the ellipsoid in A/m.""" - raise NotImplementedError - - @property - def rotation_matrix(self) -> npt.NDArray: - """Rotation matrix for the ellipsoid.""" - raise NotImplementedError