From da1c361442affa3a26a581096fc401609f4daae1 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 27 Nov 2025 13:56:20 -0800 Subject: [PATCH 01/49] Start drafting a single `Ellipsoid` class Start implementing gravity field with a permutation matrix to sort out the semiaxes. --- harmonica/__init__.py | 11 +- harmonica/_forward/ellipsoids/__init__.py | 11 +- harmonica/_forward/ellipsoids/ellipsoids.py | 1379 +++++++++++-------- harmonica/_forward/ellipsoids/gravity.py | 15 +- harmonica/_forward/ellipsoids/utils.py | 95 ++ 5 files changed, 926 insertions(+), 585 deletions(-) diff --git a/harmonica/__init__.py b/harmonica/__init__.py index 72b10625a..5f9e7a7d1 100644 --- a/harmonica/__init__.py +++ b/harmonica/__init__.py @@ -17,11 +17,12 @@ from ._euler_deconvolution import EulerDeconvolution from ._forward.dipole import dipole_magnetic from ._forward.ellipsoids import ( - OblateEllipsoid, - ProlateEllipsoid, - Sphere, - TriaxialEllipsoid, - create_ellipsoid, + Ellipsoid, + # OblateEllipsoid, + # ProlateEllipsoid, + # Sphere, + # TriaxialEllipsoid, + # create_ellipsoid, ellipsoid_gravity, ellipsoid_magnetic, ) diff --git a/harmonica/_forward/ellipsoids/__init__.py b/harmonica/_forward/ellipsoids/__init__.py index 4f916405b..e8a37c7b1 100644 --- a/harmonica/_forward/ellipsoids/__init__.py +++ b/harmonica/_forward/ellipsoids/__init__.py @@ -9,11 +9,12 @@ """ from .ellipsoids import ( - OblateEllipsoid, - ProlateEllipsoid, - Sphere, - TriaxialEllipsoid, - create_ellipsoid, + # OblateEllipsoid, + # ProlateEllipsoid, + # Sphere, + # TriaxialEllipsoid, + # create_ellipsoid, + 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 9dc18e1f5..7163b9fb1 100644 --- a/harmonica/_forward/ellipsoids/ellipsoids.py +++ b/harmonica/_forward/ellipsoids/ellipsoids.py @@ -22,141 +22,107 @@ 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`. - yaw : float, optional + a, b, c : float + Semi-axis lengths of the ellipsoid in meters. + yaw : float Rotation angle about the upward axis, in degrees. - pitch : float, optional + 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, optional + roll : float 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 + 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. + 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, + 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._a, self._b, self._c = a, b, c + + # Angles and center + self.yaw, self.pitch, self.roll = yaw, pitch, roll + self.center = center - # Sort semiaxes so a >= b >= c - c, b, a = sorted((a, b, c)) + # Physical properties of the ellipsoid + self.density = density + self.susceptibility = susceptibility + self.remanent_mag = remanent_mag - 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 - ) + @property + def a(self) -> float: + """Length of the first semiaxis in meters.""" + return self._a - return ellipsoid + @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 in meters.""" + return self._b -class BaseEllipsoid: - """ - Base class for ellipsoids. + @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: @@ -295,457 +261,730 @@ def to_pyvista(self, **kwargs): 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) +# 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": +# """ +# 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. +# +# 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`. +# yaw : float, optional +# Rotation angle about the upward axis, in degrees. +# pitch : float, optional +# 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, 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, +# ... ) +# """ +# # 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) +# +# # Sort semiaxes so a >= b >= c +# c, b, a = sorted((a, b, c)) +# +# 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 +# ) +# +# return ellipsoid +# +# +# class BaseEllipsoid: +# """ +# Base class for ellipsoids. +# +# .. important:: +# +# This class is not meant to be instantiated. +# """ +# +# @property +# def density(self) -> float | None: +# """Density of the ellipsoid in :math:`kg/m^3`.""" +# return self._density +# +# @density.setter +# def density(self, value) -> None: +# if not isinstance(value, Real) and value is not None: +# msg = ( +# f"Invalid 'density' of type '{type(value).__name__}'. " +# "It must be a float or None." +# ) +# raise TypeError(msg) +# if isinstance(value, Real): +# value = float(value) +# self._density = value +# +# @property +# def susceptibility(self) -> float | npt.NDArray | None: +# """Magnetic susceptibility of the ellipsoid in SI units.""" +# return self._susceptibility +# +# @susceptibility.setter +# def susceptibility(self, value) -> None: +# if isinstance(value, np.ndarray): +# if value.shape != (3, 3): +# msg = ( +# f"Invalid 'susceptibility' as an array with shape '{value.shape}'. " +# "It must be a (3, 3) array, a single float or None." +# ) +# raise ValueError(msg) +# elif not isinstance(value, Real) and value is not None: +# msg = ( +# f"Invalid 'susceptibility' of type '{type(value).__name__}'. " +# "It must be a (3, 3) array, a single float or None." +# ) +# raise TypeError(msg) +# if isinstance(value, Real): +# value = float(value) +# self._susceptibility = value +# +# @property +# def remanent_mag(self) -> npt.NDArray | None: +# """Remanent magnetization of the ellipsoid in A/m.""" +# return self._remanent_mag +# +# @remanent_mag.setter +# def remanent_mag(self, value) -> None: +# if isinstance(value, Sequence): +# value = np.asarray(value) +# if isinstance(value, np.ndarray): +# if value.shape != (3,): +# msg = ( +# f"Invalid 'remanent_mag' with shape '{value.shape}'. " +# "It must be a (3,) array or None." +# ) +# raise ValueError(msg) +# elif value is not None: +# msg = ( +# f"Invalid 'remanent_mag' of type '{type(value).__name__}'. " +# "It must be a (3,) array or None." +# ) +# raise TypeError(msg) +# self._remanent_mag = value +# +# @property +# def rotation_matrix(self) -> npt.NDArray: +# """ +# Create a rotation matrix for the ellipsoid. +# +# Use this matrix to rotate from the local coordinate system (centered in the +# ellipsoid center) in to the global coordinate system +# (easting, northing, upward). +# +# Returns +# ------- +# rotation_matrix : (3, 3) array +# Rotation matrix that transforms coordinates from the local ellipsoid-aligned +# coordinate system to the global coordinate system. +# +# Notes +# ----- +# Generate the rotation matrix from Tait-Bryan intrinsic angles: yaw, pitch, and +# roll. The rotations are applied in the following order: (ZŶX). Yaw (Z) and roll +# (X) rotations are done using the right-hand rule. Rotations for the pitch (Ŷ) +# are carried out in the opposite direction, so positive pitch *lifts* the easting +# axis. +# """ +# return get_rotation_matrix(self.yaw, self.pitch, self.roll) +# +# def _check_positive_semiaxis(self, value: float, semiaxis: str): +# """ +# Raise error if passed semiaxis length value is not positive. +# +# Parameters +# ---------- +# value : float +# Value of the semiaxis length. +# semiaxis : {"a", "b", "c"} +# Semiaxis name. Used to generate the error message. +# """ +# if value <= 0: +# msg = f"Invalid value of '{semiaxis}' equal to '{value}'. It must be positive." +# raise ValueError(msg) +# +# def to_pyvista(self, **kwargs): +# """ +# Export ellipsoid to a :class:`pyvista.PolyData` object. +# +# .. important:: +# +# The :mod:`pyvista` optional dependency must be installed to use this method. +# +# Parameters +# ---------- +# kwargs : dict +# Keyword arguments passed to :func:`pyvista.ParametricEllipsoid`. +# +# Returns +# ------- +# ellipsoid : pyvista.PolyData +# A PyVista's parametric ellipsoid. +# """ +# if pyvista is None: +# msg = ( +# "Missing optional dependency 'pyvista' required for " +# "exporting ellipsoids to PyVista." +# ) +# raise ImportError(msg) +# ellipsoid = pyvista.ParametricEllipsoid( +# xradius=self.a, yradius=self.b, zradius=self.c, **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..f6dc78d2d 100644 --- a/harmonica/_forward/ellipsoids/gravity.py +++ b/harmonica/_forward/ellipsoids/gravity.py @@ -20,6 +20,7 @@ from .utils import ( calculate_lambda, get_elliptical_integrals, + get_permutation_matrix, is_almost_a_sphere, is_internal, ) @@ -97,19 +98,23 @@ 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 + # Get permutation matrix and order the semiaxes + permutation_matrix = get_permutation_matrix(ellipsoid) + a, b, c = permutation_matrix @ np.array([ellipsoid.a, ellipsoid.b, ellipsoid.c]) + + # Combine the rotation and the permutation matrices + transformation = ellipsoid.rotation_matrix.T @ permutation_matrix # Rotate observation points - x, y, z = r.T @ np.vstack(coords_shifted) + x, y, z = transformation @ 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 = transformation.T @ np.vstack(gravity_ellipsoid) # sum contributions from each ellipsoid ge += ge_i diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index 9d29b3d56..c1ca4792c 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -14,6 +14,101 @@ SEMIAXES_RTOL = 1e-5 +def get_permutation_matrix(ellipsoid): + """ + Get permutation matrix for the given ellipsoid. + + Build a permutation matrix that sorts the ellipsoid's semiaxes lengths + ``ellipsoid.a``, ``ellipsoid.b``, ``ellipsoid.c`` into ``a``, ``b``, ``c`` values + that satisfy that: + + - Triaxial: ``a > b > c`` + - Prolate: ``a > b == c`` + - Oblate: ``a < b == c`` + - Sphere: ``a == b == c`` + + + Parameters + ---------- + ellipsoid : Ellipsoid + Ellipsoid object. + + Returns + ------- + permutation_matrix : (3, 3) np.ndarray + Permutation matrix. + + Examples + -------- + >>> from harmonica import Ellipsoid + >>> ellipsoid = Ellipsoid( + ... a=43.0, b=50.0, c=83.0, yaw=0, pitch=0, roll=0, center=(0, 0, 0) + ... ) + >>> get_permutation_matrix(ellipsoid) + [[0 0 1] + [0 1 0] + [1 0 0]] + + >>> ellipsoid = Ellipsoid( + ... a=43.0, b=83.0, c=50.0, yaw=0, pitch=0, roll=0, center=(0, 0, 0) + ... ) + >>> get_permutation_matrix(ellipsoid) + [[0 1 0] + [0 0 1] + [1 0 0]] + """ + # Sphere + if ellipsoid.a == ellipsoid.b == ellipsoid.c: + a, b, c = ellipsoid.a, ellipsoid.b, ellipsoid.c + permutation_matrix = np.eye(3) + return a, b, c, permutation_matrix + + # Get permutation matrix + permutation_matrix = _get_simple_permutation_matrix( + ellipsoid.a, ellipsoid.b, ellipsoid.c + ) + + # Sort semiaxes so a >= b >= c + c, b, a = sorted((ellipsoid.a, ellipsoid.b, ellipsoid.c)) + + # Treat oblate ellipsoids + if a == b: + # Swap `a` with `c` so: a < b == c + a, c = c, a + permutation_matrix[0, :], permutation_matrix[2, :] = ( + permutation_matrix[2, :], + permutation_matrix[0, :], + ) + + return permutation_matrix + + +def _get_simple_permutation_matrix(a, b, c): + """ + Build permutation matrix for the given semiaxes. + + This permutation matrix sorts the semiaxes ``a``, ``b``, ``c`` in reverse order. + + Parameters + ---------- + a, b, c : float + Semiaxes lenghts. + + Returns + ------- + permutation_matrix : (3, 3) np.ndarray + Permutation matrix. + """ + # Get sort indices + semiaxes = np.array([a, b, c]) + argsort = np.argsort(semiaxes) + # Build permutation matrix + p_matrix = np.zeros((3, 3), dtype=np.int32) + sorted_indices = np.array([2, 1, 0]) + p_matrix[sorted_indices, argsort] = 1 + return p_matrix + + def is_internal(x, y, z, a, b, c): """ Check if a given point(s) is internal or external to the ellipsoid. From 963292a0b91f07f5c1a8ce7342a51bfef28e30db Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 27 Nov 2025 15:08:20 -0800 Subject: [PATCH 02/49] Rename variable holding rotation and permutation matrix --- harmonica/_forward/ellipsoids/gravity.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/harmonica/_forward/ellipsoids/gravity.py b/harmonica/_forward/ellipsoids/gravity.py index f6dc78d2d..0ce8e4658 100644 --- a/harmonica/_forward/ellipsoids/gravity.py +++ b/harmonica/_forward/ellipsoids/gravity.py @@ -103,10 +103,10 @@ def ellipsoid_gravity( a, b, c = permutation_matrix @ np.array([ellipsoid.a, ellipsoid.b, ellipsoid.c]) # Combine the rotation and the permutation matrices - transformation = ellipsoid.rotation_matrix.T @ permutation_matrix + rotation = ellipsoid.rotation_matrix.T @ permutation_matrix # Rotate observation points - x, y, z = transformation @ np.vstack(coords_shifted) + x, y, z = rotation @ np.vstack(coords_shifted) # Calculate gravity components on local coordinate system gravity_ellipsoid = _compute_gravity_ellipsoid( @@ -114,7 +114,7 @@ def ellipsoid_gravity( ) # project onto upward unit vector, axis U - ge_i, gn_i, gu_i = transformation.T @ np.vstack(gravity_ellipsoid) + ge_i, gn_i, gu_i = rotation.T @ np.vstack(gravity_ellipsoid) # sum contributions from each ellipsoid ge += ge_i From b1d369edb4e6fd5aa1bc6b0bb21373d1d7c16b49 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 27 Nov 2025 15:42:33 -0800 Subject: [PATCH 03/49] Fix get_permutation_matrix for spheres --- harmonica/_forward/ellipsoids/utils.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index c1ca4792c..911c470b0 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -59,9 +59,8 @@ def get_permutation_matrix(ellipsoid): """ # Sphere if ellipsoid.a == ellipsoid.b == ellipsoid.c: - a, b, c = ellipsoid.a, ellipsoid.b, ellipsoid.c permutation_matrix = np.eye(3) - return a, b, c, permutation_matrix + return permutation_matrix # Get permutation matrix permutation_matrix = _get_simple_permutation_matrix( From 22c8547c51d47a8b84ed24b00d669ed2f802086b Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 28 Nov 2025 09:51:48 -0800 Subject: [PATCH 04/49] Set angles and center as optional arguments --- harmonica/_forward/ellipsoids/ellipsoids.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/harmonica/_forward/ellipsoids/ellipsoids.py b/harmonica/_forward/ellipsoids/ellipsoids.py index 7163b9fb1..f8866bd35 100644 --- a/harmonica/_forward/ellipsoids/ellipsoids.py +++ b/harmonica/_forward/ellipsoids/ellipsoids.py @@ -30,16 +30,16 @@ class Ellipsoid: ---------- a, b, c : float Semi-axis lengths of the ellipsoid in meters. - yaw : float + yaw : float, optional Rotation angle about the upward axis, in degrees. - pitch : float + pitch : float, optional 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 + roll : float, optional Rotation angle about the easting axis (after yaw and pitch rotation), in degrees. - center : tuple of float + 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 @@ -71,10 +71,10 @@ def __init__( a: float, b: float, c: float, - yaw: float, - pitch: float, - roll: float, - center: tuple[float, float, 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, From abb937975271654cc5a9a5203e9fa69141ebe7f5 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 28 Nov 2025 10:15:08 -0800 Subject: [PATCH 05/49] Update elliptical integrals to oblates as `a == b > c` Update the elliptical integrals for oblates so they assume that `a == b > c`. Simplify the permutation matrix since now we don't need to take care of oblates as a special case. --- harmonica/_forward/ellipsoids/utils.py | 104 +++++++++---------------- 1 file changed, 38 insertions(+), 66 deletions(-) diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index 911c470b0..23f536585 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -24,7 +24,7 @@ def get_permutation_matrix(ellipsoid): - Triaxial: ``a > b > c`` - Prolate: ``a > b == c`` - - Oblate: ``a < b == c`` + - Oblate: ``a == b > c`` - Sphere: ``a == b == c`` @@ -57,55 +57,21 @@ def get_permutation_matrix(ellipsoid): [0 0 1] [1 0 0]] """ - # Sphere + # Return identity matrix for a sphere if ellipsoid.a == ellipsoid.b == ellipsoid.c: permutation_matrix = np.eye(3) return permutation_matrix - # Get permutation matrix - permutation_matrix = _get_simple_permutation_matrix( - ellipsoid.a, ellipsoid.b, ellipsoid.c - ) - - # Sort semiaxes so a >= b >= c - c, b, a = sorted((ellipsoid.a, ellipsoid.b, ellipsoid.c)) - - # Treat oblate ellipsoids - if a == b: - # Swap `a` with `c` so: a < b == c - a, c = c, a - permutation_matrix[0, :], permutation_matrix[2, :] = ( - permutation_matrix[2, :], - permutation_matrix[0, :], - ) - - return permutation_matrix - - -def _get_simple_permutation_matrix(a, b, c): - """ - Build permutation matrix for the given semiaxes. - - This permutation matrix sorts the semiaxes ``a``, ``b``, ``c`` in reverse order. - - Parameters - ---------- - a, b, c : float - Semiaxes lenghts. - - Returns - ------- - permutation_matrix : (3, 3) np.ndarray - Permutation matrix. - """ # Get sort indices - semiaxes = np.array([a, b, c]) + semiaxes = np.array([ellipsoid.a, ellipsoid.b, ellipsoid.c]) argsort = np.argsort(semiaxes) + # Build permutation matrix - p_matrix = np.zeros((3, 3), dtype=np.int32) + permutation_matrix = np.zeros((3, 3), dtype=np.int32) sorted_indices = np.array([2, 1, 0]) - p_matrix[sorted_indices, argsort] = 1 - return p_matrix + permutation_matrix[sorted_indices, argsort] = 1 + + return permutation_matrix def is_internal(x, y, z, a, b, c): @@ -296,8 +262,8 @@ def get_elliptical_integrals( g1, g2, g3 = _get_elliptical_integrals_triaxial(a, b, c, lambda_) elif a > b and 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 a == b and b > c: + g1, g2, g3 = _get_elliptical_integrals_oblate(b, c, lambda_) else: msg = f"Invalid semiaxis lenghts: a={a}, b={b}, c={c}." raise ValueError(msg) @@ -486,7 +452,7 @@ def _get_elliptical_integrals_prolate(a, b, lmbda): return g1, g2, g2 -def _get_elliptical_integrals_oblate(a, b, lmbda): +def _get_elliptical_integrals_oblate(b, c, lmbda): r""" Compute elliptical integrals for a oblate ellipsoid. @@ -512,46 +478,52 @@ 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 one 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))) + arctan = np.arctan(np.sqrt((b**2 - c**2) / (c**2 + lmbda))) 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 + lmbda))) / (b**2 + lmbda)) ) - return g1, g2, g2 + g3 = ( + 2 + / ((b**2 - c**2) ** (3 / 2)) + * ((np.sqrt((b**2 - c**2) / (c**2 + lmbda))) - arctan) + ) + return g1, g1, g3 def get_derivatives_of_elliptical_integrals( From 322edb40a80b7f0021f8e07867d9d183d47dcb71 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 28 Nov 2025 10:17:54 -0800 Subject: [PATCH 06/49] Update gravity tests --- harmonica/tests/ellipsoids/test_gravity.py | 201 +++++---------------- 1 file changed, 45 insertions(+), 156 deletions(-) diff --git a/harmonica/tests/ellipsoids/test_gravity.py b/harmonica/tests/ellipsoids/test_gravity.py index cb626dc84..990280e37 100644 --- a/harmonica/tests/ellipsoids/test_gravity.py +++ b/harmonica/tests/ellipsoids/test_gravity.py @@ -22,12 +22,7 @@ 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.errors import NoPhysicalPropertyWarning @@ -46,49 +41,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 +66,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 +93,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 +291,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 @@ -382,47 +343,33 @@ class TestEllipsoidVsSphere: @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, - ) + b = a + self.delta + c = b case "prolate": - ellipsoid = ProlateEllipsoid( - a=a, - b=a - self.delta, - yaw=yaw, - pitch=pitch, - center=self.center, - density=self.density, - ) + b = a - self.delta + c = b 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, - ) + b = a - self.delta + c = a - 2 * self.delta case _: raise ValueError() + ellipsoid = Ellipsoid(a, b, c, center=self.center, density=self.density) return ellipsoid @pytest.fixture @@ -467,8 +414,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,36 +428,19 @@ 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): @@ -558,23 +487,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 +538,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 +595,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 +609,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 +629,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 +650,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) From 8ffc88c2092fc1504975b235a46509cdfae8cfb7 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 28 Nov 2025 10:35:14 -0800 Subject: [PATCH 07/49] Update test_ellipsoids.py Remove quite a lot of tests since we now have a single ellipsoid class. --- harmonica/tests/ellipsoids/test_ellipsoids.py | 441 ++---------------- 1 file changed, 40 insertions(+), 401 deletions(-) diff --git a/harmonica/tests/ellipsoids/test_ellipsoids.py b/harmonica/tests/ellipsoids/test_ellipsoids.py index 1655a2243..51e7690cb 100644 --- a/harmonica/tests/ellipsoids/test_ellipsoids.py +++ b/harmonica/tests/ellipsoids/test_ellipsoids.py @@ -15,13 +15,7 @@ import numpy as np import pytest -from harmonica import ( - OblateEllipsoid, - ProlateEllipsoid, - Sphere, - TriaxialEllipsoid, - create_ellipsoid, -) +from harmonica import Ellipsoid try: import pyvista @@ -29,161 +23,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 +40,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 +54,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 +68,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 +132,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 +165,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 +178,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): From 2be3c559ae5a1b4e0a26142d5c641d9537339b79 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 28 Nov 2025 11:28:46 -0800 Subject: [PATCH 08/49] Update docstrings of get_permutation_matrix --- harmonica/_forward/ellipsoids/utils.py | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index 23f536585..71438f196 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -19,14 +19,8 @@ def get_permutation_matrix(ellipsoid): Get permutation matrix for the given ellipsoid. Build a permutation matrix that sorts the ellipsoid's semiaxes lengths - ``ellipsoid.a``, ``ellipsoid.b``, ``ellipsoid.c`` into ``a``, ``b``, ``c`` values - that satisfy that: - - - Triaxial: ``a > b > c`` - - Prolate: ``a > b == c`` - - Oblate: ``a == b > c`` - - Sphere: ``a == b == c`` - + ``ellipsoid.a``, ``ellipsoid.b``, ``ellipsoid.c`` in reverse order, into ``a``, + ``b``, ``c`` values that verify ``a >= b >= c``. Parameters ---------- @@ -44,7 +38,8 @@ def get_permutation_matrix(ellipsoid): >>> ellipsoid = Ellipsoid( ... a=43.0, b=50.0, c=83.0, yaw=0, pitch=0, roll=0, center=(0, 0, 0) ... ) - >>> get_permutation_matrix(ellipsoid) + >>> matrix = get_permutation_matrix(ellipsoid) + >>> print(matrix) [[0 0 1] [0 1 0] [1 0 0]] @@ -52,7 +47,8 @@ def get_permutation_matrix(ellipsoid): >>> ellipsoid = Ellipsoid( ... a=43.0, b=83.0, c=50.0, yaw=0, pitch=0, roll=0, center=(0, 0, 0) ... ) - >>> get_permutation_matrix(ellipsoid) + >>> matrix = get_permutation_matrix(ellipsoid) + >>> print(matrix) [[0 1 0] [0 0 1] [1 0 0]] From d430c7f72a0fa8d946fc1eed7daa502b2a1e7174 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 28 Nov 2025 11:29:25 -0800 Subject: [PATCH 09/49] Use `sorted` to sort semiaxes in gravity function It makes it more explicit about what's happening with the semiaxes. The permutation matrix should match that sorting, but we can cover that in tests. --- harmonica/_forward/ellipsoids/gravity.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/harmonica/_forward/ellipsoids/gravity.py b/harmonica/_forward/ellipsoids/gravity.py index 0ce8e4658..5fffc1eff 100644 --- a/harmonica/_forward/ellipsoids/gravity.py +++ b/harmonica/_forward/ellipsoids/gravity.py @@ -100,7 +100,7 @@ def ellipsoid_gravity( # Get permutation matrix and order the semiaxes permutation_matrix = get_permutation_matrix(ellipsoid) - a, b, c = permutation_matrix @ np.array([ellipsoid.a, ellipsoid.b, ellipsoid.c]) + a, b, c = sorted((ellipsoid.a, ellipsoid.b, ellipsoid.c), reverse=True) # Combine the rotation and the permutation matrices rotation = ellipsoid.rotation_matrix.T @ permutation_matrix From a1cdfc8dec7144fc8ac40313f78eceaa0d166add Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 28 Nov 2025 11:33:21 -0800 Subject: [PATCH 10/49] Update magnetic code to Ellipsoid class and oblates as `a == b > c` Update the magnetic code to use the single `Ellipsoid` class, make use of the permutation matrix, and update the internal demagnetization tensor equations for the oblate ellipsoid defined as `a == b > c`. --- harmonica/_forward/ellipsoids/magnetic.py | 49 +-- harmonica/tests/ellipsoids/test_magnetic.py | 314 +++++--------------- 2 files changed, 108 insertions(+), 255 deletions(-) diff --git a/harmonica/_forward/ellipsoids/magnetic.py b/harmonica/_forward/ellipsoids/magnetic.py index 1075f755d..b568d1102 100644 --- a/harmonica/_forward/ellipsoids/magnetic.py +++ b/harmonica/_forward/ellipsoids/magnetic.py @@ -25,6 +25,7 @@ calculate_lambda, get_derivatives_of_elliptical_integrals, get_elliptical_integrals, + get_permutation_matrix, is_almost_a_sphere, is_internal, ) @@ -153,28 +154,32 @@ def _single_ellipsoid_magnetic( easting, northing, upward = coordinates coords_shifted = (easting - origin_e, northing - origin_n, upward - origin_u) + # Get permutation matrix and order the semiaxes + permutation_matrix = get_permutation_matrix(ellipsoid) + a, b, c = sorted((ellipsoid.a, ellipsoid.b, ellipsoid.c), reverse=True) + + # Combine the rotation and the permutation matrices + rotation = ellipsoid.rotation_matrix.T @ permutation_matrix + # 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 +192,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 @@ -366,8 +371,8 @@ def get_demagnetization_tensor_internal(a: float, b: float, c: float): 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) + elif a == b and b > c: + n_diagonal = _demag_tensor_oblate_internal(b, c) else: msg = "Could not determine ellipsoid type for values given." raise ValueError(msg) @@ -434,26 +439,26 @@ def _demag_tensor_prolate_internal(a: float, b: float): 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: floats + Semi-axes lengths of the oblate ellipsoid (``a == b > c``). Returns ------- nxx, nyy, nzz : floats individual diagonal components of the x, y, z matrix. """ - m = a / b + m = c / b if not 0 < m < 1: - msg = f"Invalid aspect ratio for oblate ellipsoid: a={a}, b={b}, a/b={m}" + msg = f"Invalid aspect ratio for oblate ellipsoid: b={b}, c={c}, c/b={m}" raise ValueError(msg) - nxx = 1 / (1 - m**2) * (1 - (m / np.sqrt(1 - m**2)) * np.arccos(m)) - nyy = nzz = 0.5 * (1 - nxx) + 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 diff --git a/harmonica/tests/ellipsoids/test_magnetic.py b/harmonica/tests/ellipsoids/test_magnetic.py index 8fab168c5..1bc23785d 100644 --- a/harmonica/tests/ellipsoids/test_magnetic.py +++ b/harmonica/tests/ellipsoids/test_magnetic.py @@ -24,12 +24,7 @@ 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, @@ -47,22 +42,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 +55,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 +69,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 +91,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 +109,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 +129,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 +155,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 +213,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 +248,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 +288,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 +338,11 @@ class TestMagnetizationVersusSphere: Test if ellipsoid's magnetization approximates the one of the sphere. """ + # Difference between ellipsoid's semiaxes. + # It should be small compared to the sphere radius, so the ellipsoid approximates + # a sphere. + delta = 1e-2 + @pytest.fixture def radius(self): """ @@ -427,12 +359,13 @@ def ellipsoid_semiaxes(self, radius, request): ellipsoid_type = request.param match ellipsoid_type: case "oblate": - b = c = a + 1e-2 + b = a + c = a - self.delta case "prolate": - b = c = a - 1e-2 + b = c = a - self.delta case "triaxial": - b = a - 1e-3 - c = a - 1e-2 + b = a - self.delta + c = a - 2 * self.delta case _: raise ValueError() return a, b, c @@ -505,42 +438,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 = a - self.delta 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 = a - self.delta 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 = a - self.delta + c = a - 2 * self.delta case _: raise ValueError() + ellipsoid = Ellipsoid(a, b, c, susceptibility=self.susceptibility) return ellipsoid @pytest.mark.parametrize("ellipsoid_type", ["oblate", "prolate", "triaxial"]) @@ -548,8 +461,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) @@ -603,44 +516,22 @@ def get_ellipsoid(self, ellipsoid_type: str): ------- OblateEllipsoid, ProlateEllipsoid, Sphere or TriaxialEllipsoid """ - 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 = a - self.delta 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 = a - self.delta 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 = a - self.delta + c = a - 2 * self.delta 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): @@ -692,8 +583,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 +596,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 +668,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 +780,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( From 5a8211d4d3befc0e1092992c89f1e6a92853c623 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 28 Nov 2025 11:34:31 -0800 Subject: [PATCH 11/49] Remove unused import --- harmonica/tests/ellipsoids/test_ellipsoids.py | 1 - 1 file changed, 1 deletion(-) diff --git a/harmonica/tests/ellipsoids/test_ellipsoids.py b/harmonica/tests/ellipsoids/test_ellipsoids.py index 51e7690cb..5ed545faa 100644 --- a/harmonica/tests/ellipsoids/test_ellipsoids.py +++ b/harmonica/tests/ellipsoids/test_ellipsoids.py @@ -8,7 +8,6 @@ Test ellipsoid classes. """ -import itertools import re from unittest.mock import patch From 9ceb55245dfcce4e991a916e34f7a180bacfa8ee Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 28 Nov 2025 11:51:56 -0800 Subject: [PATCH 12/49] Get rid of mentions of the old ellipsoid classes --- doc/api/index.rst | 8 +- harmonica/__init__.py | 11 +- harmonica/_forward/ellipsoids/__init__.py | 9 +- harmonica/_forward/ellipsoids/ellipsoids.py | 729 -------------------- harmonica/_forward/ellipsoids/gravity.py | 4 +- harmonica/_forward/ellipsoids/magnetic.py | 4 +- harmonica/tests/ellipsoids/test_magnetic.py | 2 +- 7 files changed, 7 insertions(+), 760 deletions(-) diff --git a/doc/api/index.rst b/doc/api/index.rst index 481d6fa4c..7671ddea6 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: diff --git a/harmonica/__init__.py b/harmonica/__init__.py index 5f9e7a7d1..3acc064d8 100644 --- a/harmonica/__init__.py +++ b/harmonica/__init__.py @@ -16,16 +16,7 @@ from ._equivalent_sources.spherical import EquivalentSourcesSph from ._euler_deconvolution import EulerDeconvolution from ._forward.dipole import dipole_magnetic -from ._forward.ellipsoids import ( - Ellipsoid, - # 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 e8a37c7b1..9b3a6885c 100644 --- a/harmonica/_forward/ellipsoids/__init__.py +++ b/harmonica/_forward/ellipsoids/__init__.py @@ -8,13 +8,6 @@ Forward modelling of ellipsoids. """ -from .ellipsoids import ( - # OblateEllipsoid, - # ProlateEllipsoid, - # Sphere, - # TriaxialEllipsoid, - # create_ellipsoid, - 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 f8866bd35..c56e8a966 100644 --- a/harmonica/_forward/ellipsoids/ellipsoids.py +++ b/harmonica/_forward/ellipsoids/ellipsoids.py @@ -259,732 +259,3 @@ def to_pyvista(self, **kwargs): ellipsoid.rotate(rotation=self.rotation_matrix, inplace=True) ellipsoid.translate(self.center, inplace=True) return ellipsoid - - -# 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": -# """ -# 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. -# -# 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`. -# yaw : float, optional -# Rotation angle about the upward axis, in degrees. -# pitch : float, optional -# 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, 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, -# ... ) -# """ -# # 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) -# -# # Sort semiaxes so a >= b >= c -# c, b, a = sorted((a, b, c)) -# -# 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 -# ) -# -# return ellipsoid -# -# -# class BaseEllipsoid: -# """ -# Base class for ellipsoids. -# -# .. important:: -# -# This class is not meant to be instantiated. -# """ -# -# @property -# def density(self) -> float | None: -# """Density of the ellipsoid in :math:`kg/m^3`.""" -# return self._density -# -# @density.setter -# def density(self, value) -> None: -# if not isinstance(value, Real) and value is not None: -# msg = ( -# f"Invalid 'density' of type '{type(value).__name__}'. " -# "It must be a float or None." -# ) -# raise TypeError(msg) -# if isinstance(value, Real): -# value = float(value) -# self._density = value -# -# @property -# def susceptibility(self) -> float | npt.NDArray | None: -# """Magnetic susceptibility of the ellipsoid in SI units.""" -# return self._susceptibility -# -# @susceptibility.setter -# def susceptibility(self, value) -> None: -# if isinstance(value, np.ndarray): -# if value.shape != (3, 3): -# msg = ( -# f"Invalid 'susceptibility' as an array with shape '{value.shape}'. " -# "It must be a (3, 3) array, a single float or None." -# ) -# raise ValueError(msg) -# elif not isinstance(value, Real) and value is not None: -# msg = ( -# f"Invalid 'susceptibility' of type '{type(value).__name__}'. " -# "It must be a (3, 3) array, a single float or None." -# ) -# raise TypeError(msg) -# if isinstance(value, Real): -# value = float(value) -# self._susceptibility = value -# -# @property -# def remanent_mag(self) -> npt.NDArray | None: -# """Remanent magnetization of the ellipsoid in A/m.""" -# return self._remanent_mag -# -# @remanent_mag.setter -# def remanent_mag(self, value) -> None: -# if isinstance(value, Sequence): -# value = np.asarray(value) -# if isinstance(value, np.ndarray): -# if value.shape != (3,): -# msg = ( -# f"Invalid 'remanent_mag' with shape '{value.shape}'. " -# "It must be a (3,) array or None." -# ) -# raise ValueError(msg) -# elif value is not None: -# msg = ( -# f"Invalid 'remanent_mag' of type '{type(value).__name__}'. " -# "It must be a (3,) array or None." -# ) -# raise TypeError(msg) -# self._remanent_mag = value -# -# @property -# def rotation_matrix(self) -> npt.NDArray: -# """ -# Create a rotation matrix for the ellipsoid. -# -# Use this matrix to rotate from the local coordinate system (centered in the -# ellipsoid center) in to the global coordinate system -# (easting, northing, upward). -# -# Returns -# ------- -# rotation_matrix : (3, 3) array -# Rotation matrix that transforms coordinates from the local ellipsoid-aligned -# coordinate system to the global coordinate system. -# -# Notes -# ----- -# Generate the rotation matrix from Tait-Bryan intrinsic angles: yaw, pitch, and -# roll. The rotations are applied in the following order: (ZŶX). Yaw (Z) and roll -# (X) rotations are done using the right-hand rule. Rotations for the pitch (Ŷ) -# are carried out in the opposite direction, so positive pitch *lifts* the easting -# axis. -# """ -# return get_rotation_matrix(self.yaw, self.pitch, self.roll) -# -# def _check_positive_semiaxis(self, value: float, semiaxis: str): -# """ -# Raise error if passed semiaxis length value is not positive. -# -# Parameters -# ---------- -# value : float -# Value of the semiaxis length. -# semiaxis : {"a", "b", "c"} -# Semiaxis name. Used to generate the error message. -# """ -# if value <= 0: -# msg = f"Invalid value of '{semiaxis}' equal to '{value}'. It must be positive." -# raise ValueError(msg) -# -# def to_pyvista(self, **kwargs): -# """ -# Export ellipsoid to a :class:`pyvista.PolyData` object. -# -# .. important:: -# -# The :mod:`pyvista` optional dependency must be installed to use this method. -# -# Parameters -# ---------- -# kwargs : dict -# Keyword arguments passed to :func:`pyvista.ParametricEllipsoid`. -# -# Returns -# ------- -# ellipsoid : pyvista.PolyData -# A PyVista's parametric ellipsoid. -# """ -# if pyvista is None: -# msg = ( -# "Missing optional dependency 'pyvista' required for " -# "exporting ellipsoids to PyVista." -# ) -# raise ImportError(msg) -# ellipsoid = pyvista.ParametricEllipsoid( -# xradius=self.a, yradius=self.b, zradius=self.c, **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 5fffc1eff..ce01eb776 100644 --- a/harmonica/_forward/ellipsoids/gravity.py +++ b/harmonica/_forward/ellipsoids/gravity.py @@ -55,9 +55,7 @@ def ellipsoid_gravity( 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. + Ellipsoidal body represented by an instance of :class:`harmonica.Ellipsoid`. Returns ------- diff --git a/harmonica/_forward/ellipsoids/magnetic.py b/harmonica/_forward/ellipsoids/magnetic.py index b568d1102..f0b89768c 100644 --- a/harmonica/_forward/ellipsoids/magnetic.py +++ b/harmonica/_forward/ellipsoids/magnetic.py @@ -53,9 +53,7 @@ def ellipsoid_magnetic( 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. + Ellipsoidal body represented by an instance of :class:`harmonica.Ellipsoid`. external_field : tuple The uniform magnetic field (B) as and array with values of (magnitude, inclination, declination). The magnitude should be in nT, diff --git a/harmonica/tests/ellipsoids/test_magnetic.py b/harmonica/tests/ellipsoids/test_magnetic.py index 1bc23785d..460ea3770 100644 --- a/harmonica/tests/ellipsoids/test_magnetic.py +++ b/harmonica/tests/ellipsoids/test_magnetic.py @@ -514,7 +514,7 @@ def get_ellipsoid(self, ellipsoid_type: str): Returns ------- - OblateEllipsoid, ProlateEllipsoid, Sphere or TriaxialEllipsoid + Ellipsoid """ match ellipsoid_type: case "oblate": From 6084607c6852a3d46f5f032fca3ed107d39da74a Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 28 Nov 2025 16:20:07 -0800 Subject: [PATCH 13/49] Rename matrix to `transformation` Do this to differentiate it from just the rotation matrix, since it also includes the permutation one. --- harmonica/_forward/ellipsoids/gravity.py | 6 +++--- harmonica/_forward/ellipsoids/magnetic.py | 10 +++++----- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/harmonica/_forward/ellipsoids/gravity.py b/harmonica/_forward/ellipsoids/gravity.py index ce01eb776..5d138d32b 100644 --- a/harmonica/_forward/ellipsoids/gravity.py +++ b/harmonica/_forward/ellipsoids/gravity.py @@ -101,10 +101,10 @@ def ellipsoid_gravity( a, b, c = sorted((ellipsoid.a, ellipsoid.b, ellipsoid.c), reverse=True) # Combine the rotation and the permutation matrices - rotation = ellipsoid.rotation_matrix.T @ permutation_matrix + transformation = ellipsoid.rotation_matrix.T @ permutation_matrix # Rotate observation points - x, y, z = rotation @ np.vstack(coords_shifted) + x, y, z = transformation @ np.vstack(coords_shifted) # Calculate gravity components on local coordinate system gravity_ellipsoid = _compute_gravity_ellipsoid( @@ -112,7 +112,7 @@ def ellipsoid_gravity( ) # project onto upward unit vector, axis U - ge_i, gn_i, gu_i = rotation.T @ np.vstack(gravity_ellipsoid) + ge_i, gn_i, gu_i = transformation.T @ np.vstack(gravity_ellipsoid) # sum contributions from each ellipsoid ge += ge_i diff --git a/harmonica/_forward/ellipsoids/magnetic.py b/harmonica/_forward/ellipsoids/magnetic.py index f0b89768c..1e72d4062 100644 --- a/harmonica/_forward/ellipsoids/magnetic.py +++ b/harmonica/_forward/ellipsoids/magnetic.py @@ -157,10 +157,10 @@ def _single_ellipsoid_magnetic( a, b, c = sorted((ellipsoid.a, ellipsoid.b, ellipsoid.c), reverse=True) # Combine the rotation and the permutation matrices - rotation = ellipsoid.rotation_matrix.T @ permutation_matrix + transformation = ellipsoid.rotation_matrix.T @ permutation_matrix # Rotate observation points - x, y, z = rotation @ np.vstack(coords_shifted) + x, y, z = transformation @ np.vstack(coords_shifted) # Build internal demagnetization tensor n_tensor_internal = get_demagnetization_tensor_internal(a, b, c) @@ -170,8 +170,8 @@ def _single_ellipsoid_magnetic( remanent_mag = cast_remanent_magnetization(remanent_mag) # Rotate the external field and the remanent magnetization - h0_field_rotated = rotation @ h0_field - remnant_mag_rotated = rotation @ remanent_mag + h0_field_rotated = transformation @ h0_field + remnant_mag_rotated = transformation @ remanent_mag # Get magnetization of the ellipsoid magnetization = get_magnetisation( @@ -202,7 +202,7 @@ def _single_ellipsoid_magnetic( b_field[~internal, :] = mu_0 * (-n_tensors @ magnetization) # Rotate the b fields - be, bn, bu = rotation.T @ b_field.T + be, bn, bu = transformation.T @ b_field.T return be, bn, bu From acca7b7276fc380539bfb9a1762e4456f6e1a382 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 1 Dec 2025 16:49:34 -0800 Subject: [PATCH 14/49] Add test comparing gravity of triaxial vs oblate and prolate --- harmonica/tests/ellipsoids/test_gravity.py | 63 ++++++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/harmonica/tests/ellipsoids/test_gravity.py b/harmonica/tests/ellipsoids/test_gravity.py index 990280e37..203cc5c79 100644 --- a/harmonica/tests/ellipsoids/test_gravity.py +++ b/harmonica/tests/ellipsoids/test_gravity.py @@ -658,3 +658,66 @@ 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 + ) From 5383d85176fdab6a982f4b48cf75459f3c20f9a8 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 1 Dec 2025 16:55:41 -0800 Subject: [PATCH 15/49] Add test comparing mag field of triaxials vs oblate and prolate --- harmonica/tests/ellipsoids/test_magnetic.py | 66 +++++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/harmonica/tests/ellipsoids/test_magnetic.py b/harmonica/tests/ellipsoids/test_magnetic.py index 460ea3770..7d4e467e2 100644 --- a/harmonica/tests/ellipsoids/test_magnetic.py +++ b/harmonica/tests/ellipsoids/test_magnetic.py @@ -802,3 +802,69 @@ def test_warning(self): # Check the gravity acceleration 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 + ) From ff85ce99cc99b0cb8ec984a36b0622d3028be32a Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 1 Dec 2025 16:56:09 -0800 Subject: [PATCH 16/49] Fix comment in test --- harmonica/tests/ellipsoids/test_magnetic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/harmonica/tests/ellipsoids/test_magnetic.py b/harmonica/tests/ellipsoids/test_magnetic.py index 7d4e467e2..d4a0e417c 100644 --- a/harmonica/tests/ellipsoids/test_magnetic.py +++ b/harmonica/tests/ellipsoids/test_magnetic.py @@ -799,7 +799,7 @@ def test_warning(self): 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 From 83998ccd10dc8fa95ea9b18c5559553a073c1f34 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 1 Dec 2025 16:56:34 -0800 Subject: [PATCH 17/49] Fix comments in tests --- harmonica/tests/ellipsoids/test_gravity.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/harmonica/tests/ellipsoids/test_gravity.py b/harmonica/tests/ellipsoids/test_gravity.py index 203cc5c79..4f5af2ba4 100644 --- a/harmonica/tests/ellipsoids/test_gravity.py +++ b/harmonica/tests/ellipsoids/test_gravity.py @@ -445,7 +445,7 @@ def ellipsoid(self, request): 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. @@ -458,7 +458,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) From bd92efd87e7ec7fbc877ec81c345ff972b057eab Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 2 Dec 2025 15:35:34 -0800 Subject: [PATCH 18/49] Fix bug in permutation matrix The permutation matrix used wasn't right: we need to apply a proper rotation matrix (with 90 degrees) to ensure that the local coordinate system is a right-handed one. --- harmonica/_forward/ellipsoids/gravity.py | 16 ++-- harmonica/_forward/ellipsoids/magnetic.py | 20 ++--- harmonica/_forward/ellipsoids/utils.py | 90 ++++++++++++----------- 3 files changed, 68 insertions(+), 58 deletions(-) diff --git a/harmonica/_forward/ellipsoids/gravity.py b/harmonica/_forward/ellipsoids/gravity.py index 5d138d32b..d5b490ad1 100644 --- a/harmonica/_forward/ellipsoids/gravity.py +++ b/harmonica/_forward/ellipsoids/gravity.py @@ -20,7 +20,7 @@ from .utils import ( calculate_lambda, get_elliptical_integrals, - get_permutation_matrix, + get_semiaxes_rotation_matrix, is_almost_a_sphere, is_internal, ) @@ -96,15 +96,17 @@ def ellipsoid_gravity( origin_e, origin_n, origin_u = ellipsoid.center coords_shifted = (easting - origin_e, northing - origin_n, upward - origin_u) - # Get permutation matrix and order the semiaxes - permutation_matrix = get_permutation_matrix(ellipsoid) + # Sort the semiaxes (a >= b >= c) a, b, c = sorted((ellipsoid.a, ellipsoid.b, ellipsoid.c), reverse=True) - # Combine the rotation and the permutation matrices - transformation = ellipsoid.rotation_matrix.T @ permutation_matrix + # 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 = transformation @ np.vstack(coords_shifted) + x, y, z = rotation @ np.vstack(coords_shifted) # Calculate gravity components on local coordinate system gravity_ellipsoid = _compute_gravity_ellipsoid( @@ -112,7 +114,7 @@ def ellipsoid_gravity( ) # project onto upward unit vector, axis U - ge_i, gn_i, gu_i = transformation.T @ np.vstack(gravity_ellipsoid) + ge_i, gn_i, gu_i = rotation.T @ np.vstack(gravity_ellipsoid) # sum contributions from each ellipsoid ge += ge_i diff --git a/harmonica/_forward/ellipsoids/magnetic.py b/harmonica/_forward/ellipsoids/magnetic.py index 1e72d4062..4239d8f45 100644 --- a/harmonica/_forward/ellipsoids/magnetic.py +++ b/harmonica/_forward/ellipsoids/magnetic.py @@ -25,7 +25,7 @@ calculate_lambda, get_derivatives_of_elliptical_integrals, get_elliptical_integrals, - get_permutation_matrix, + get_semiaxes_rotation_matrix, is_almost_a_sphere, is_internal, ) @@ -152,15 +152,17 @@ def _single_ellipsoid_magnetic( easting, northing, upward = coordinates coords_shifted = (easting - origin_e, northing - origin_n, upward - origin_u) - # Get permutation matrix and order the semiaxes - permutation_matrix = get_permutation_matrix(ellipsoid) + # Sort the semiaxes (a >= b >= c) a, b, c = sorted((ellipsoid.a, ellipsoid.b, ellipsoid.c), reverse=True) - # Combine the rotation and the permutation matrices - transformation = ellipsoid.rotation_matrix.T @ permutation_matrix + # 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 = transformation @ np.vstack(coords_shifted) + x, y, z = rotation @ np.vstack(coords_shifted) # Build internal demagnetization tensor n_tensor_internal = get_demagnetization_tensor_internal(a, b, c) @@ -170,8 +172,8 @@ def _single_ellipsoid_magnetic( remanent_mag = cast_remanent_magnetization(remanent_mag) # Rotate the external field and the remanent magnetization - h0_field_rotated = transformation @ h0_field - remnant_mag_rotated = transformation @ remanent_mag + h0_field_rotated = rotation @ h0_field + remnant_mag_rotated = rotation @ remanent_mag # Get magnetization of the ellipsoid magnetization = get_magnetisation( @@ -202,7 +204,7 @@ def _single_ellipsoid_magnetic( b_field[~internal, :] = mu_0 * (-n_tensors @ magnetization) # Rotate the b fields - be, bn, bu = transformation.T @ b_field.T + be, bn, bu = rotation.T @ b_field.T return be, bn, bu diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index 71438f196..ec768d80a 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -8,19 +8,36 @@ 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 get_permutation_matrix(ellipsoid): +def get_semiaxes_rotation_matrix(ellipsoid): """ - Get permutation matrix for the given 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 - Build a permutation matrix that sorts the ellipsoid's semiaxes lengths - ``ellipsoid.a``, ``ellipsoid.b``, ``ellipsoid.c`` in reverse order, into ``a``, - ``b``, ``c`` values that verify ``a >= b >= c``. + a, b, c = sorted((ellipsoid.a, ellipsoid.b, ellipsoid.c), reverse=True) Parameters ---------- @@ -29,45 +46,34 @@ def get_permutation_matrix(ellipsoid): Returns ------- - permutation_matrix : (3, 3) np.ndarray - Permutation matrix. - - Examples - -------- - >>> from harmonica import Ellipsoid - >>> ellipsoid = Ellipsoid( - ... a=43.0, b=50.0, c=83.0, yaw=0, pitch=0, roll=0, center=(0, 0, 0) - ... ) - >>> matrix = get_permutation_matrix(ellipsoid) - >>> print(matrix) - [[0 0 1] - [0 1 0] - [1 0 0]] - - >>> ellipsoid = Ellipsoid( - ... a=43.0, b=83.0, c=50.0, yaw=0, pitch=0, roll=0, center=(0, 0, 0) - ... ) - >>> matrix = get_permutation_matrix(ellipsoid) - >>> print(matrix) - [[0 1 0] - [0 0 1] - [1 0 0]] - """ - # Return identity matrix for a sphere - if ellipsoid.a == ellipsoid.b == ellipsoid.c: - permutation_matrix = np.eye(3) - return permutation_matrix - - # Get sort indices - semiaxes = np.array([ellipsoid.a, ellipsoid.b, ellipsoid.c]) - argsort = np.argsort(semiaxes) + rotation_matrix : (3, 3) np.ndarray + Rotation matrix. - # Build permutation matrix - permutation_matrix = np.zeros((3, 3), dtype=np.int32) - sorted_indices = np.array([2, 1, 0]) - permutation_matrix[sorted_indices, argsort] = 1 + 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: + raise ValueError() - return permutation_matrix + matrix = get_rotation_matrix(yaw, pitch, roll).astype(int) + return matrix def is_internal(x, y, z, a, b, c): From 3efba9f65f3fd1a3efd0bd7f51ca498103e17646 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 2 Dec 2025 15:36:22 -0800 Subject: [PATCH 19/49] Test semiaxes in any order with equivalent ellipsoids Test ellipsoids defined with semiaxes passed in arbitrary orders against equivalent ellipsoids defined with ``a >= b >= c``, but with necessary rotations. --- harmonica/tests/ellipsoids/test_gravity.py | 93 ++++++++++++++ harmonica/tests/ellipsoids/test_magnetic.py | 130 ++++++++++++++++++++ 2 files changed, 223 insertions(+) diff --git a/harmonica/tests/ellipsoids/test_gravity.py b/harmonica/tests/ellipsoids/test_gravity.py index 4f5af2ba4..21ef08263 100644 --- a/harmonica/tests/ellipsoids/test_gravity.py +++ b/harmonica/tests/ellipsoids/test_gravity.py @@ -721,3 +721,96 @@ def test_triaxial_vs_oblate(self): 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) diff --git a/harmonica/tests/ellipsoids/test_magnetic.py b/harmonica/tests/ellipsoids/test_magnetic.py index d4a0e417c..fd97b7cfd 100644 --- a/harmonica/tests/ellipsoids/test_magnetic.py +++ b/harmonica/tests/ellipsoids/test_magnetic.py @@ -868,3 +868,133 @@ def test_triaxial_vs_oblate(self): 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) From 29923456263030c3a5edabc56233086759dee2a4 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 2 Dec 2025 15:40:29 -0800 Subject: [PATCH 20/49] Update tests for `__str__` method --- harmonica/tests/ellipsoids/test_ellipsoids.py | 105 +++--------------- 1 file changed, 16 insertions(+), 89 deletions(-) diff --git a/harmonica/tests/ellipsoids/test_ellipsoids.py b/harmonica/tests/ellipsoids/test_ellipsoids.py index ebd1fd181..b069bda5e 100644 --- a/harmonica/tests/ellipsoids/test_ellipsoids.py +++ b/harmonica/tests/ellipsoids/test_ellipsoids.py @@ -242,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" @@ -310,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 @@ -336,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 @@ -348,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 @@ -360,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 From 1a25659fac32f427d0ef4ce245b341fd0a4af063 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 2 Dec 2025 16:06:15 -0800 Subject: [PATCH 21/49] Add test to get_semiaxes_rotation_matrix function --- harmonica/tests/ellipsoids/test_utils.py | 49 +++++++++++++++++++++++- 1 file changed, 48 insertions(+), 1 deletion(-) diff --git a/harmonica/tests/ellipsoids/test_utils.py b/harmonica/tests/ellipsoids/test_utils.py index 914427083..63772d2f2 100644 --- a/harmonica/tests/ellipsoids/test_utils.py +++ b/harmonica/tests/ellipsoids/test_utils.py @@ -4,10 +4,16 @@ # # This code is part of the Fatiando a Terra project (https://www.fatiando.org) # +import itertools + 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, +) @pytest.mark.parametrize( @@ -94,3 +100,44 @@ 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) From 65d907ddd09359319d180497e85fc52b94048bca Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 2 Dec 2025 17:07:29 -0800 Subject: [PATCH 22/49] Fix is_almost_a_sphere function After changing the definition of the oblate to `a == b > c`, we need to change the condition for oblates in such function. --- harmonica/_forward/ellipsoids/utils.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index ec768d80a..4d899af9e 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -98,7 +98,7 @@ 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. Parameters @@ -114,10 +114,14 @@ def is_almost_a_sphere(a: float, b: float, c: float) -> bool: 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 From c723ca565a77f3a3e83cf84228260cbf428751f9 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 2 Dec 2025 17:15:48 -0800 Subject: [PATCH 23/49] Add test for is_almost_a_sphere --- harmonica/_forward/ellipsoids/utils.py | 7 +++++++ harmonica/tests/ellipsoids/test_utils.py | 18 ++++++++++++++++++ 2 files changed, 25 insertions(+) diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index 4d899af9e..80a6971fe 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -101,6 +101,10 @@ def is_almost_a_sphere(a: float, b: float, c: float) -> bool: 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 @@ -110,6 +114,9 @@ def is_almost_a_sphere(a: float, b: float, c: float) -> bool: ------- bool """ + if not (a >= b >= c): + raise ValueError() + # Exactly a sphere if a == b == c: return True diff --git a/harmonica/tests/ellipsoids/test_utils.py b/harmonica/tests/ellipsoids/test_utils.py index 63772d2f2..4ff233b69 100644 --- a/harmonica/tests/ellipsoids/test_utils.py +++ b/harmonica/tests/ellipsoids/test_utils.py @@ -13,6 +13,7 @@ from harmonica._forward.ellipsoids.utils import ( calculate_lambda, get_semiaxes_rotation_matrix, + is_almost_a_sphere, ) @@ -141,3 +142,20 @@ def test_example(self): 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 + ], +) +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 From 5ff0b35b8be2db559c1a348e0f640c9e8f4c1984 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 2 Dec 2025 17:21:51 -0800 Subject: [PATCH 24/49] Fix numerical instabilities in triaxial's gravity Fix the numerical instabilities of gravity fields produced by triaxial ellipsoids that approximate a prolate and an oblate ellipsoid. Fallback to the prolate and oblate solutions when two of the three semiaxes are close enough to each other. --- harmonica/_forward/ellipsoids/utils.py | 84 ++++++++++++++++++++++++-- 1 file changed, 80 insertions(+), 4 deletions(-) diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index 80a6971fe..187d58b71 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -140,6 +140,74 @@ 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 + Ellipsoid's semiaxes lenghts. + + Returns + ------- + bool + """ + if not (a >= b >= c): + raise ValueError() + + # 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 + Ellipsoid's semiaxes lenghts. + + Returns + ------- + bool + """ + if not (a >= b >= c): + raise ValueError() + + # 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. @@ -225,6 +293,11 @@ 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 @@ -271,12 +344,15 @@ 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: + if is_almost_a_sphere(a, b, c): + raise ValueError() + + if is_almost_prolate(a, b, c): g1, g2, g3 = _get_elliptical_integrals_prolate(a, b, lambda_) - elif a == b and b > c: + 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}." raise ValueError(msg) From 6549a179b643e16314efbedb396c441c18f9aaed Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 2 Dec 2025 17:23:15 -0800 Subject: [PATCH 25/49] Add test for numerical instabilities on gravity fields --- harmonica/tests/ellipsoids/test_gravity.py | 97 ++++++++++++++++++++++ 1 file changed, 97 insertions(+) diff --git a/harmonica/tests/ellipsoids/test_gravity.py b/harmonica/tests/ellipsoids/test_gravity.py index 21ef08263..20e0a0bec 100644 --- a/harmonica/tests/ellipsoids/test_gravity.py +++ b/harmonica/tests/ellipsoids/test_gravity.py @@ -814,3 +814,100 @@ def test_oblate(self, semiaxes, yaw, pitch, roll): 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) From 528ffebc266a958af5c55bd2f22504e4888ef3ce Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 2 Dec 2025 17:42:57 -0800 Subject: [PATCH 26/49] Fix numerical instabilities in internal demagnetization tensor Fallback to the prolate and oblate solutions for the internal demagnetization tensor when triaxial has two semiaxes close enough to each other. --- harmonica/_forward/ellipsoids/magnetic.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/harmonica/_forward/ellipsoids/magnetic.py b/harmonica/_forward/ellipsoids/magnetic.py index 4239d8f45..0ff194a70 100644 --- a/harmonica/_forward/ellipsoids/magnetic.py +++ b/harmonica/_forward/ellipsoids/magnetic.py @@ -27,6 +27,8 @@ get_elliptical_integrals, get_semiaxes_rotation_matrix, is_almost_a_sphere, + is_almost_oblate, + is_almost_prolate, is_internal, ) @@ -367,12 +369,12 @@ def get_demagnetization_tensor_internal(a: float, b: float, c: float): """ if is_almost_a_sphere(a, b, c): n_diagonal = 1 / 3 * np.ones(3) - elif a > b > c: - n_diagonal = _demag_tensor_triaxial_internal(a, b, c) - elif a > b and b == c: + elif is_almost_prolate(a, b, c): n_diagonal = _demag_tensor_prolate_internal(a, b) - elif a == b and b > c: + 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) else: msg = "Could not determine ellipsoid type for values given." raise ValueError(msg) From 384df38a6420d78078cb88f3eee152375699fd7e Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 2 Dec 2025 17:43:52 -0800 Subject: [PATCH 27/49] Add tests for the magnetic numerical instabilities --- harmonica/tests/ellipsoids/test_magnetic.py | 110 ++++++++++++++++++++ 1 file changed, 110 insertions(+) diff --git a/harmonica/tests/ellipsoids/test_magnetic.py b/harmonica/tests/ellipsoids/test_magnetic.py index fd97b7cfd..0f4decbb4 100644 --- a/harmonica/tests/ellipsoids/test_magnetic.py +++ b/harmonica/tests/ellipsoids/test_magnetic.py @@ -998,3 +998,113 @@ def test_oblate(self, semiaxes, yaw, pitch, roll): ) 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) From 202be31843d2d483a378ede32f1328d892362280 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Wed, 3 Dec 2025 15:43:10 -0800 Subject: [PATCH 28/49] Update gravity tests vs sphere Use the ratio instead of a static delta, so we are sure we are comparing the two different codes and not just triggering the numerical instability fix that falls back to the spherical solutions. --- harmonica/tests/ellipsoids/test_gravity.py | 27 +++++++++++----------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/harmonica/tests/ellipsoids/test_gravity.py b/harmonica/tests/ellipsoids/test_gravity.py index 20e0a0bec..68a55d6e0 100644 --- a/harmonica/tests/ellipsoids/test_gravity.py +++ b/harmonica/tests/ellipsoids/test_gravity.py @@ -23,6 +23,7 @@ from harmonica import ellipsoid_gravity, point_gravity from harmonica._forward.ellipsoids import Ellipsoid +from harmonica._forward.ellipsoids.utils import SEMIAXES_RTOL from harmonica.errors import NoPhysicalPropertyWarning @@ -335,10 +336,11 @@ 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): @@ -359,14 +361,15 @@ def ellipsoid(self, request): a = self.radius match request.param: case "oblate": - b = a + self.delta - c = b + a = b = self.radius + c = (1 - self.ratio) * b case "prolate": - b = a - self.delta - c = b + a = self.radius + b = c = (1 - self.ratio) * a case "triaxial": - b = a - self.delta - c = a - 2 * self.delta + 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) @@ -395,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: From d40d9fdca12d5237735bb2b9f85774cdb1e000d3 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Wed, 3 Dec 2025 15:51:00 -0800 Subject: [PATCH 29/49] Update magnetic tests vs sphere Use the ratio instead of a static delta, so we are sure we are comparing the two different codes and not just triggering the numerical instability fix that falls back to the spherical solutions. --- harmonica/tests/ellipsoids/test_magnetic.py | 63 +++++++++++---------- 1 file changed, 34 insertions(+), 29 deletions(-) diff --git a/harmonica/tests/ellipsoids/test_magnetic.py b/harmonica/tests/ellipsoids/test_magnetic.py index 0f4decbb4..6f714a3ad 100644 --- a/harmonica/tests/ellipsoids/test_magnetic.py +++ b/harmonica/tests/ellipsoids/test_magnetic.py @@ -29,6 +29,7 @@ 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 @@ -338,10 +339,11 @@ class TestMagnetizationVersusSphere: Test if ellipsoid's magnetization approximates the one of the sphere. """ - # Difference between ellipsoid's semiaxes. - # It should be small compared to the sphere radius, so the ellipsoid approximates - # a sphere. - delta = 1e-2 + # 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): @@ -355,17 +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 = a - c = a - self.delta + a = b = radius + c = (1 - self.ratio) * b case "prolate": - b = c = a - self.delta + a = radius + b = c = (1 - self.ratio) * a case "triaxial": - b = a - self.delta - c = a - 2 * self.delta + a = radius + b = (1 - self.ratio) * a + c = (1 - 2 * self.ratio) * a case _: raise ValueError() return a, b, c @@ -410,10 +413,11 @@ class TestMagneticFieldVersusSphere: 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) @@ -443,14 +447,14 @@ def get_ellipsoid(self, ellipsoid_type: str): match ellipsoid_type: case "oblate": a = b = self.radius - c = a - self.delta + c = (1 - self.ratio) * b case "prolate": a = self.radius - b = c = a - self.delta + b = c = (1 - self.ratio) * a case "triaxial": a = self.radius - b = a - self.delta - c = a - 2 * self.delta + b = (1 - self.ratio) * a + c = (1 - 2 * self.ratio) * a case _: raise ValueError() ellipsoid = Ellipsoid(a, b, c, susceptibility=self.susceptibility) @@ -472,7 +476,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) @@ -486,10 +490,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) @@ -519,14 +524,14 @@ def get_ellipsoid(self, ellipsoid_type: str): match ellipsoid_type: case "oblate": a = b = self.radius - c = a - self.delta + c = (1 - self.ratio) * b case "prolate": a = self.radius - b = c = a - self.delta + b = c = (1 - self.ratio) * a case "triaxial": a = self.radius - b = a - self.delta - c = a - 2 * self.delta + b = (1 - self.ratio) * a + c = (1 - 2 * self.ratio) * a case "sphere": a = b = c = self.radius case _: @@ -565,8 +570,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: From d000fe42eb04ff43a7edcfb7557ea2255b987cae Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 17:38:22 +0000 Subject: [PATCH 30/49] Fix grammar in docstring --- harmonica/_forward/ellipsoids/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index 187d58b71..4f517e423 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -597,7 +597,7 @@ def _get_elliptical_integrals_oblate(b, c, lmbda): .. important:: - These equations are modified versions of the one in Takahashi (2018), adapted to + These equations are modified versions of the ones in Takahashi (2018), adapted to any oblate ellipsoid defined as: ``a = b > c``. """ From 3e287075da94714cdc5a8c0bc458e931cbc5a91c Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 17:39:01 +0000 Subject: [PATCH 31/49] Use class in parameter description --- harmonica/_forward/ellipsoids/gravity.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/harmonica/_forward/ellipsoids/gravity.py b/harmonica/_forward/ellipsoids/gravity.py index d5b490ad1..15e8f0f92 100644 --- a/harmonica/_forward/ellipsoids/gravity.py +++ b/harmonica/_forward/ellipsoids/gravity.py @@ -54,7 +54,7 @@ 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 + ellipsoid : harmonica.Ellipsoid or list of harmonica.Ellipsoid Ellipsoidal body represented by an instance of :class:`harmonica.Ellipsoid`. Returns From d19ba94f546e8a5827c3002b269f54abb68bb69b Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 17:39:41 +0000 Subject: [PATCH 32/49] Fix docstyle in parameter descriptions --- harmonica/_forward/ellipsoids/utils.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index 4f517e423..85423bc65 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -107,7 +107,7 @@ def is_almost_a_sphere(a: float, b: float, c: float) -> bool: Parameters ---------- - a, b, c: float + a, b, c : float Ellipsoid's semiaxes lenghts. Returns @@ -153,7 +153,7 @@ def is_almost_prolate(a: float, b: float, c: float) -> bool: Parameters ---------- - a, b, c: float + a, b, c : float Ellipsoid's semiaxes lenghts. Returns @@ -187,7 +187,7 @@ def is_almost_oblate(a: float, b: float, c: float) -> bool: Parameters ---------- - a, b, c: float + a, b, c : float Ellipsoid's semiaxes lenghts. Returns From 74437dec13cfd92d5ce428a60383755a2d53a65b Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 17:39:59 +0000 Subject: [PATCH 33/49] Fix docstring --- harmonica/_forward/ellipsoids/gravity.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/harmonica/_forward/ellipsoids/gravity.py b/harmonica/_forward/ellipsoids/gravity.py index 15e8f0f92..7fe462ae0 100644 --- a/harmonica/_forward/ellipsoids/gravity.py +++ b/harmonica/_forward/ellipsoids/gravity.py @@ -55,7 +55,7 @@ def ellipsoid_gravity( coordinates of the computation points defined on a Cartesian coordinate system. All coordinates should be in meters. ellipsoid : harmonica.Ellipsoid or list of harmonica.Ellipsoid - Ellipsoidal body represented by an instance of :class:`harmonica.Ellipsoid`. + Ellipsoidal body represented by an instance of :class:`harmonica.Ellipsoid` or a list of them. Returns ------- From 97a68240940ed5bd2c9b174b4fbf76b43f65b0cf Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 09:41:15 -0800 Subject: [PATCH 34/49] Remove unused variable in magnetic test --- harmonica/tests/ellipsoids/test_magnetic.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/harmonica/tests/ellipsoids/test_magnetic.py b/harmonica/tests/ellipsoids/test_magnetic.py index 6f714a3ad..c94e8e4b5 100644 --- a/harmonica/tests/ellipsoids/test_magnetic.py +++ b/harmonica/tests/ellipsoids/test_magnetic.py @@ -408,9 +408,8 @@ 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 # Ratio between ellipsoid's semiaxes, defined as ratio = | a - b | / max(a, b). From ae1e1b3d6dbfe513be8782b687db53747a6acf98 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 09:42:35 -0800 Subject: [PATCH 35/49] Add ids to parametrized test function --- harmonica/tests/ellipsoids/test_utils.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/harmonica/tests/ellipsoids/test_utils.py b/harmonica/tests/ellipsoids/test_utils.py index 4ff233b69..16fe2d885 100644 --- a/harmonica/tests/ellipsoids/test_utils.py +++ b/harmonica/tests/ellipsoids/test_utils.py @@ -155,6 +155,15 @@ def test_example(self): (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.""" From d335125b875381bcfa9f7494fc86c720e971e1c1 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 09:44:03 -0800 Subject: [PATCH 36/49] Add error message to ValueError --- harmonica/_forward/ellipsoids/utils.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index 85423bc65..dde68eb7e 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -70,7 +70,8 @@ def get_semiaxes_rotation_matrix(ellipsoid): elif c >= a >= b: yaw, pitch, roll = 90, 90, 0 else: - raise ValueError() + msg = f"Invalid semiaxes: a={a}, b={b}, c={c}." + raise ValueError(msg) matrix = get_rotation_matrix(yaw, pitch, roll).astype(int) return matrix From 1a3e4183650ec95418521b4368e609b9a8bb38c3 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 09:44:59 -0800 Subject: [PATCH 37/49] Fix typo in error message --- harmonica/_forward/ellipsoids/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index dde68eb7e..2728613aa 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -355,7 +355,7 @@ def get_elliptical_integrals( 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 lenghts: a={a}, b={b}, c={c}." raise ValueError(msg) return g1, g2, g3 From 8b5495adb41b8c1adfa401c01e5db2a00f8566f7 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 09:48:13 -0800 Subject: [PATCH 38/49] Remove `typing.Ellipsoid` class Since we merged all classes to a single one, we don't actually need an `Ellipsoid` type class. It's rare that users will pass an object that is not an instance of `hm.Ellipsoid` or a child of it. --- doc/api/index.rst | 1 - harmonica/_forward/ellipsoids/gravity.py | 3 +- harmonica/_forward/ellipsoids/magnetic.py | 3 +- harmonica/typing.py | 60 ----------------------- 4 files changed, 4 insertions(+), 63 deletions(-) diff --git a/doc/api/index.rst b/doc/api/index.rst index 7671ddea6..65328119f 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -170,4 +170,3 @@ Type hints :toctree: generated/ typing.Coordinates - typing.Ellipsoid diff --git a/harmonica/_forward/ellipsoids/gravity.py b/harmonica/_forward/ellipsoids/gravity.py index 7fe462ae0..530abe095 100644 --- a/harmonica/_forward/ellipsoids/gravity.py +++ b/harmonica/_forward/ellipsoids/gravity.py @@ -16,7 +16,8 @@ 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, get_elliptical_integrals, diff --git a/harmonica/_forward/ellipsoids/magnetic.py b/harmonica/_forward/ellipsoids/magnetic.py index 0ff194a70..3fdbbfc13 100644 --- a/harmonica/_forward/ellipsoids/magnetic.py +++ b/harmonica/_forward/ellipsoids/magnetic.py @@ -20,7 +20,8 @@ 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, get_derivatives_of_elliptical_integrals, 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 From b09b7dcb21e2c409606965117672b1751e47ec94 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 11:00:50 -0800 Subject: [PATCH 39/49] Add messages to several ValueErrors --- harmonica/_forward/ellipsoids/utils.py | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index 2728613aa..0a878c740 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -116,7 +116,11 @@ def is_almost_a_sphere(a: float, b: float, c: float) -> bool: bool """ if not (a >= b >= c): - raise ValueError() + msg = ( + "Invalid semiaxes not properly sorted (a >= b >= c): " + f"a={a}, b={b}, c={c}." + ) + raise ValueError(msg) # Exactly a sphere if a == b == c: @@ -162,7 +166,11 @@ def is_almost_prolate(a: float, b: float, c: float) -> bool: bool """ if not (a >= b >= c): - raise ValueError() + msg = ( + "Invalid semiaxes not properly sorted (a >= b >= c): " + f"a={a}, b={b}, c={c}." + ) + raise ValueError(msg) # Exactly a prolate if a > b == c: @@ -196,7 +204,11 @@ def is_almost_oblate(a: float, b: float, c: float) -> bool: bool """ if not (a >= b >= c): - raise ValueError() + msg = ( + "Invalid semiaxes not properly sorted (a >= b >= c): " + f"a={a}, b={b}, c={c}." + ) + raise ValueError(msg) # Exactly an oblate if a == b > c: @@ -346,7 +358,11 @@ def get_elliptical_integrals( oblate and prolate). """ if is_almost_a_sphere(a, b, c): - raise ValueError() + 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_) From 915b831d710c5ca73b6927c4d4fb2519f3648425 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 11:04:15 -0800 Subject: [PATCH 40/49] Fix typo in lenghts --- harmonica/_forward/ellipsoids/utils.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index 0a878c740..aa834bf38 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -109,7 +109,7 @@ def is_almost_a_sphere(a: float, b: float, c: float) -> bool: Parameters ---------- a, b, c : float - Ellipsoid's semiaxes lenghts. + Ellipsoid's semiaxes lengths. Returns ------- @@ -159,7 +159,7 @@ def is_almost_prolate(a: float, b: float, c: float) -> bool: Parameters ---------- a, b, c : float - Ellipsoid's semiaxes lenghts. + Ellipsoid's semiaxes lengths. Returns ------- @@ -197,7 +197,7 @@ def is_almost_oblate(a: float, b: float, c: float) -> bool: Parameters ---------- a, b, c : float - Ellipsoid's semiaxes lenghts. + Ellipsoid's semiaxes lengths. Returns ------- @@ -371,7 +371,7 @@ def get_elliptical_integrals( elif a > b > c: g1, g2, g3 = _get_elliptical_integrals_triaxial(a, b, c, lambda_) else: - msg = f"Invalid semiaxes 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 From b3cc620f8372aa7da3e3dc7dbf2c4f1b6c8cc32e Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 11:05:35 -0800 Subject: [PATCH 41/49] Update ellipsoid's description in docstring --- harmonica/_forward/ellipsoids/gravity.py | 3 ++- harmonica/_forward/ellipsoids/magnetic.py | 5 +++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/harmonica/_forward/ellipsoids/gravity.py b/harmonica/_forward/ellipsoids/gravity.py index 530abe095..778ef224c 100644 --- a/harmonica/_forward/ellipsoids/gravity.py +++ b/harmonica/_forward/ellipsoids/gravity.py @@ -56,7 +56,8 @@ def ellipsoid_gravity( coordinates of the computation points defined on a Cartesian coordinate system. All coordinates should be in meters. ellipsoid : harmonica.Ellipsoid or list of harmonica.Ellipsoid - Ellipsoidal body represented by an instance of :class:`harmonica.Ellipsoid` or a list of them. + Ellipsoidal body represented by an instance of :class:`harmonica.Ellipsoid` or + a list of them. Returns ------- diff --git a/harmonica/_forward/ellipsoids/magnetic.py b/harmonica/_forward/ellipsoids/magnetic.py index 3fdbbfc13..b2503005f 100644 --- a/harmonica/_forward/ellipsoids/magnetic.py +++ b/harmonica/_forward/ellipsoids/magnetic.py @@ -55,8 +55,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.Ellipsoid`. + 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, From 59d69c7f361c8f2bd57b25223fdddfb96c62daa6 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 11:34:55 -0800 Subject: [PATCH 42/49] Check for sorted semiaxes when required Add sanity checks to some of the private functions in which the semiaxes are assumed to be sorted such as `a >= b >= c`. Most of these errors are never going to be accessed by code ran by users, but they might catch bugs when refactoring. Adding those checks there since it's hard to test the behaviour from outside. --- harmonica/_forward/ellipsoids/gravity.py | 6 +- harmonica/_forward/ellipsoids/magnetic.py | 21 +++-- harmonica/_forward/ellipsoids/utils.py | 98 +++++++++++++---------- 3 files changed, 76 insertions(+), 49 deletions(-) diff --git a/harmonica/_forward/ellipsoids/gravity.py b/harmonica/_forward/ellipsoids/gravity.py index 778ef224c..3e30d3351 100644 --- a/harmonica/_forward/ellipsoids/gravity.py +++ b/harmonica/_forward/ellipsoids/gravity.py @@ -20,6 +20,7 @@ from .ellipsoids import Ellipsoid from .utils import ( calculate_lambda, + check_semiaxes_sorted, get_elliptical_integrals, get_semiaxes_rotation_matrix, is_almost_a_sphere, @@ -150,8 +151,7 @@ def _compute_gravity_ellipsoid( x, y, z : arrays 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. + Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. density : float Density of the ellipsoidal body in kg/m³. @@ -161,6 +161,8 @@ def _compute_gravity_ellipsoid( 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 b2503005f..5690e0fe4 100644 --- a/harmonica/_forward/ellipsoids/magnetic.py +++ b/harmonica/_forward/ellipsoids/magnetic.py @@ -24,6 +24,7 @@ from .ellipsoids import Ellipsoid from .utils import ( calculate_lambda, + check_semiaxes_sorted, get_derivatives_of_elliptical_integrals, get_elliptical_integrals, get_semiaxes_rotation_matrix, @@ -238,7 +239,7 @@ def get_magnetisation( Parameters ---------- a, b, c : floats - Semi-axes lengths of the ellipsoid. + Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. susceptibility : (3, 3) array Susceptibility tensor. h0_field : (3,) array @@ -277,6 +278,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) @@ -349,7 +352,7 @@ def get_demagnetization_tensor_internal(a: float, b: float, c: float): Parameters ---------- a, b, c : floats - Semi-axes lengths of the given ellipsoid. + Semi-axes lengths of the given ellipsoid sorted such as ``a >= b >= c``. Returns ------- @@ -369,6 +372,8 @@ 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): @@ -392,13 +397,17 @@ 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). + Semi-axes lengths of the given ellipsoid sorted such as ``a > b > c``. Returns ------- nxx, nyy, nzz : floats 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) @@ -426,12 +435,12 @@ def _demag_tensor_prolate_internal(a: float, b: float): Parameters ---------- a, b: floats - Semi-axes lengths of the prolate ellipsoid (a > b = c). + Semi-axes lengths of the prolate ellipsoid (``a > b = c``). Returns ------- nxx, nyy, nzz : floats - individual diagonal components of the x, y, z matrix. + Diagonal components of the internal demagnetization tensor. """ m = a / b if not m > 1: @@ -455,7 +464,7 @@ def _demag_tensor_oblate_internal(b: float, c: float): Returns ------- nxx, nyy, nzz : floats - individual diagonal components of the x, y, z matrix. + Diagonal components of the internal demagnetization tensor. """ m = c / b if not 0 < m < 1: diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index aa834bf38..53e8c30e9 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -16,6 +16,27 @@ 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 = ( + "Invalid semiaxes not properly sorted (a >= b >= c): " + f"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. @@ -109,18 +130,13 @@ def is_almost_a_sphere(a: float, b: float, c: float) -> bool: Parameters ---------- a, b, c : float - Ellipsoid's semiaxes lengths. + Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. Returns ------- bool """ - if not (a >= b >= c): - msg = ( - "Invalid semiaxes not properly sorted (a >= b >= c): " - f"a={a}, b={b}, c={c}." - ) - raise ValueError(msg) + check_semiaxes_sorted(a, b, c) # Exactly a sphere if a == b == c: @@ -159,18 +175,13 @@ def is_almost_prolate(a: float, b: float, c: float) -> bool: Parameters ---------- a, b, c : float - Ellipsoid's semiaxes lengths. + Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. Returns ------- bool """ - if not (a >= b >= c): - msg = ( - "Invalid semiaxes not properly sorted (a >= b >= c): " - f"a={a}, b={b}, c={c}." - ) - raise ValueError(msg) + check_semiaxes_sorted(a, b, c) # Exactly a prolate if a > b == c: @@ -197,18 +208,13 @@ def is_almost_oblate(a: float, b: float, c: float) -> bool: Parameters ---------- a, b, c : float - Ellipsoid's semiaxes lengths. + Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. Returns ------- bool """ - if not (a >= b >= c): - msg = ( - "Invalid semiaxes not properly sorted (a >= b >= c): " - f"a={a}, b={b}, c={c}." - ) - raise ValueError(msg) + check_semiaxes_sorted(a, b, c) # Exactly an oblate if a == b > c: @@ -231,19 +237,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 ------- @@ -251,7 +248,11 @@ def calculate_lambda(x, y, z, a, b, c): The computed value(s) of the lambda parameter. """ + check_semiaxes_sorted(a, b, c) + # Solve lambda for prolate and oblate ellipsoids + # TODO: update this code. I think this is not right for the new definition of + # oblate. Maybe add a notes section explaining why splitting the computations. 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 @@ -281,7 +282,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) @@ -314,7 +315,7 @@ def get_elliptical_integrals( Parameters ---------- a, b, c : floats - Semi-axes lengths of the given ellipsoid. + 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. @@ -357,6 +358,8 @@ def get_elliptical_integrals( Expressions of the elliptic integrals vary for each type of ellipsoid (triaxial, oblate and prolate). """ + check_semiaxes_sorted(a, b, c) + if is_almost_a_sphere(a, b, c): msg = ( "Invalid semiaxes that create (almost) a spherical ellipsoid: " @@ -383,7 +386,7 @@ def _get_elliptical_integrals_triaxial(a, b, c, lmbda): Parameters ---------- a, b, c : floats - Semi-axes lengths of the given ellipsoid. + Semi-axes lengths of the ellipsoid sorted such as ``a > b > c``. lmbda : float The given lambda value for the point we are considering. @@ -448,6 +451,10 @@ 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)) phi = np.arcsin(int_arcsin) @@ -487,7 +494,7 @@ def _get_elliptical_integrals_prolate(a, b, lmbda): Parameters ---------- a, b : floats - Semi-axes lengths of the given ellipsoid. + Semi-axes lengths of the given ellipsoid, where ``a > b``. lmbda : float The given lambda value for the point we are considering. @@ -546,6 +553,10 @@ 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) @@ -564,8 +575,8 @@ def _get_elliptical_integrals_oblate(b, c, lmbda): Parameters ---------- - a, b : floats - Semi-axes lengths of the given ellipsoid. + b, c : floats + Semi-axes lengths of the given ellipsoid, where ``b > c``. lmbda : float The given lambda value for the point we are considering. @@ -618,6 +629,10 @@ def _get_elliptical_integrals_oblate(b, c, lmbda): any oblate ellipsoid defined as: ``a = b > c``. """ + 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 + lmbda))) g1 = ( 1 @@ -644,7 +659,7 @@ def get_derivatives_of_elliptical_integrals( Parameters ---------- a, b, c : floats - Semi-axes lengths of the given ellipsoid. + Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. lambda_ : float The given lambda value for the point we are considering. @@ -653,6 +668,7 @@ def get_derivatives_of_elliptical_integrals( hx, hy, hz : tuple of floats 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 From aaf3dbc3885879c83eb4678bd44a33a49c443f94 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 11:36:22 -0800 Subject: [PATCH 43/49] Update definition of oblate ellipsoid in lambda tests --- harmonica/tests/ellipsoids/test_utils.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/harmonica/tests/ellipsoids/test_utils.py b/harmonica/tests/ellipsoids/test_utils.py index 16fe2d885..46781c0d4 100644 --- a/harmonica/tests/ellipsoids/test_utils.py +++ b/harmonica/tests/ellipsoids/test_utils.py @@ -5,6 +5,7 @@ # This code is part of the Fatiando a Terra project (https://www.fatiando.org) # import itertools +import re import numpy as np import pytest @@ -19,7 +20,7 @@ @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): @@ -46,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): @@ -75,7 +87,7 @@ def test_zero_cases_for_lambda(a, b, c): @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): From 28fa366ecdfff38cd838b68412ff36049e2caf53 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 11:36:46 -0800 Subject: [PATCH 44/49] Fix style --- harmonica/_forward/ellipsoids/utils.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index 53e8c30e9..36c4150d1 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -31,8 +31,7 @@ def check_semiaxes_sorted(a, b, c): """ if not (a >= b >= c): msg = ( - "Invalid semiaxes not properly sorted (a >= b >= c): " - f"a={a}, b={b}, c={c}." + f"Invalid semiaxes not properly sorted (a >= b >= c): a={a}, b={b}, c={c}." ) raise ValueError(msg) From 9af1039d41ab8eea9a05b5076e966623a2ef4d7c Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 11:52:21 -0800 Subject: [PATCH 45/49] Update calculate_lambda for new definition of oblate Split the computation of lambda for oblates and prolates. Since we changed the definition of oblate ellipsoids to `a = b > c`, the solution of lambda through the second degree equation is different than the one for prolates. --- harmonica/_forward/ellipsoids/utils.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index 36c4150d1..e26633985 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -249,15 +249,19 @@ def calculate_lambda(x, y, z, a, b, c): """ check_semiaxes_sorted(a, b, c) - # Solve lambda for prolate and oblate ellipsoids - # TODO: update this code. I think this is not right for the new definition of - # oblate. Maybe add a notes section explaining why splitting the computations. + # 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 From ee4d04f24278ae867a6b2ba2cb105f14fa2503e3 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 12:04:24 -0800 Subject: [PATCH 46/49] Improve a few docstrings and rename lmbda for lambda_ --- harmonica/_forward/ellipsoids/magnetic.py | 6 +-- harmonica/_forward/ellipsoids/utils.py | 62 +++++++++++------------ harmonica/tests/ellipsoids/test_utils.py | 18 +++---- 3 files changed, 43 insertions(+), 43 deletions(-) diff --git a/harmonica/_forward/ellipsoids/magnetic.py b/harmonica/_forward/ellipsoids/magnetic.py index 5690e0fe4..d022a706b 100644 --- a/harmonica/_forward/ellipsoids/magnetic.py +++ b/harmonica/_forward/ellipsoids/magnetic.py @@ -569,17 +569,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 diff --git a/harmonica/_forward/ellipsoids/utils.py b/harmonica/_forward/ellipsoids/utils.py index e26633985..2aa02509e 100644 --- a/harmonica/_forward/ellipsoids/utils.py +++ b/harmonica/_forward/ellipsoids/utils.py @@ -103,9 +103,9 @@ def is_internal(x, y, z, a, b, c): 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 @@ -317,14 +317,14 @@ def get_elliptical_integrals( Parameters ---------- - a, b, c : floats + 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. @@ -382,22 +382,22 @@ def get_elliptical_integrals( 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 + a, b, c : float Semi-axes lengths of the ellipsoid sorted such as ``a > b > c``. - lmbda : float + 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 ----- @@ -459,7 +459,7 @@ def _get_elliptical_integrals_triaxial(a, b, c, lmbda): 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) @@ -474,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) @@ -483,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 + a, b : float Semi-axes lengths of the given ellipsoid, where ``a > b``. - lmbda : float + 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 ----- @@ -563,31 +563,31 @@ def _get_elliptical_integrals_prolate(a, b, lmbda): # 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(b, c, lmbda): +def _get_elliptical_integrals_oblate(b, c, lambda_): r""" Compute elliptical integrals for a oblate ellipsoid. Parameters ---------- - b, c : floats + b, c : float Semi-axes lengths of the given ellipsoid, where ``b > c``. - lmbda : float + 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 ----- @@ -636,16 +636,16 @@ def _get_elliptical_integrals_oblate(b, c, lmbda): 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 + lmbda))) + arctan = np.arctan(np.sqrt((b**2 - c**2) / (c**2 + lambda_))) g1 = ( 1 / ((b**2 - c**2) ** (3 / 2)) - * (arctan - (np.sqrt((b**2 - c**2) * (c**2 + lmbda))) / (b**2 + lmbda)) + * (arctan - (np.sqrt((b**2 - c**2) * (c**2 + lambda_))) / (b**2 + lambda_)) ) g3 = ( 2 / ((b**2 - c**2) ** (3 / 2)) - * ((np.sqrt((b**2 - c**2) / (c**2 + lmbda))) - arctan) + * ((np.sqrt((b**2 - c**2) / (c**2 + lambda_))) - arctan) ) return g1, g1, g3 @@ -661,14 +661,14 @@ def get_derivatives_of_elliptical_integrals( Parameters ---------- - a, b, c : floats + a, b, c : float Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. - lambda_ : float + 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) diff --git a/harmonica/tests/ellipsoids/test_utils.py b/harmonica/tests/ellipsoids/test_utils.py index 46781c0d4..3a97957cb 100644 --- a/harmonica/tests/ellipsoids/test_utils.py +++ b/harmonica/tests/ellipsoids/test_utils.py @@ -69,19 +69,19 @@ 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"]) From ec9f6d99d033ff5550292573df66e90d65bf20b7 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 14:23:54 -0800 Subject: [PATCH 47/49] Specify definition of susceptibility tensor --- harmonica/_forward/ellipsoids/ellipsoids.py | 3 ++- harmonica/_forward/ellipsoids/magnetic.py | 8 +++++--- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/harmonica/_forward/ellipsoids/ellipsoids.py b/harmonica/_forward/ellipsoids/ellipsoids.py index 8472f26cc..ea0ef79da 100644 --- a/harmonica/_forward/ellipsoids/ellipsoids.py +++ b/harmonica/_forward/ellipsoids/ellipsoids.py @@ -47,7 +47,8 @@ class Ellipsoid: 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. + 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 diff --git a/harmonica/_forward/ellipsoids/magnetic.py b/harmonica/_forward/ellipsoids/magnetic.py index d022a706b..6e1795cab 100644 --- a/harmonica/_forward/ellipsoids/magnetic.py +++ b/harmonica/_forward/ellipsoids/magnetic.py @@ -136,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 @@ -241,7 +242,7 @@ def get_magnetisation( a, b, c : floats 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 @@ -300,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 From 147c0374e252f103b41f7b7c70abec8e82dffbae Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 14:31:01 -0800 Subject: [PATCH 48/49] Improve error messages in private functions --- harmonica/_forward/ellipsoids/magnetic.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/harmonica/_forward/ellipsoids/magnetic.py b/harmonica/_forward/ellipsoids/magnetic.py index 6e1795cab..a53188a26 100644 --- a/harmonica/_forward/ellipsoids/magnetic.py +++ b/harmonica/_forward/ellipsoids/magnetic.py @@ -444,10 +444,10 @@ def _demag_tensor_prolate_internal(a: float, b: float): nxx, nyy, nzz : floats 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) @@ -468,10 +468,10 @@ def _demag_tensor_oblate_internal(b: float, c: float): nxx, nyy, nzz : floats Diagonal components of the internal demagnetization tensor. """ - m = c / b - if not 0 < m < 1: - msg = f"Invalid aspect ratio for oblate ellipsoid: b={b}, c={c}, c/b={m}" + if not b > c: + msg = f"Invalid semiaxes for oblate ellipsoid (not b > c): b={b}, c={c}." raise ValueError(msg) + 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 From 8e05ea95f69cd8f77076e893024e5605863958e1 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 5 Jan 2026 14:34:21 -0800 Subject: [PATCH 49/49] Fix a few more docstrings --- harmonica/_forward/ellipsoids/gravity.py | 6 +++--- harmonica/_forward/ellipsoids/magnetic.py | 24 +++++++++++------------ 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/harmonica/_forward/ellipsoids/gravity.py b/harmonica/_forward/ellipsoids/gravity.py index 3e30d3351..b3baaa986 100644 --- a/harmonica/_forward/ellipsoids/gravity.py +++ b/harmonica/_forward/ellipsoids/gravity.py @@ -148,16 +148,16 @@ 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 + 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). """ diff --git a/harmonica/_forward/ellipsoids/magnetic.py b/harmonica/_forward/ellipsoids/magnetic.py index a53188a26..9159ff6f0 100644 --- a/harmonica/_forward/ellipsoids/magnetic.py +++ b/harmonica/_forward/ellipsoids/magnetic.py @@ -239,7 +239,7 @@ def get_magnetisation( Parameters ---------- - a, b, c : floats + a, b, c : float Semi-axes lengths of the ellipsoid sorted such as ``a >= b >= c``. susceptibility : (3, 3) array Susceptibility tensor defined in the local coordinate system of the ellipsoid. @@ -353,7 +353,7 @@ def get_demagnetization_tensor_internal(a: float, b: float, c: float): Parameters ---------- - a, b, c : floats + a, b, c : float Semi-axes lengths of the given ellipsoid sorted such as ``a >= b >= c``. Returns @@ -398,12 +398,12 @@ def _demag_tensor_triaxial_internal(a: float, b: float, c: float): Parameters ---------- - a, b, c : floats + 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: @@ -436,12 +436,12 @@ def _demag_tensor_prolate_internal(a: float, b: float): Parameters ---------- - a, b: floats + a, b : float Semi-axes lengths of the prolate ellipsoid (``a > b = c``). Returns ------- - nxx, nyy, nzz : floats + nxx, nyy, nzz : float Diagonal components of the internal demagnetization tensor. """ if not a > b: @@ -460,12 +460,12 @@ def _demag_tensor_oblate_internal(b: float, c: float): Parameters ---------- - b, c: floats + b, c : float Semi-axes lengths of the oblate ellipsoid (``a == b > c``). Returns ------- - nxx, nyy, nzz : floats + nxx, nyy, nzz : float Diagonal components of the internal demagnetization tensor. """ if not b > c: @@ -494,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 @@ -601,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. """