Add a new create_ellipsoid function (#606)

Add a new `create_ellipsoid` function that takes values of the three
semiaxes lenghts (`a`, `b`, and `c`) in any particular order and returns
an ellipsoid object of the corresponding type. It takes the `center` as
required argument as well, but defaults yaw, pitch, and roll angles to
zero. It can also take physical properties of the ellipsoids that are
passed through ``kwargs``. Add tests for the new function, and
reorganize the ellipsoid code and tests under new `ellipsoids`
submodule.

**Relevant issues/PRs:**

Part of #594

Author: santisoler

doc/api/index.rst

@@ -91,6 +91,7 @@ Ellipsoids:
 .. autosummary::
     :toctree: generated/
 
+    create_ellipsoid
     OblateEllipsoid
     ProlateEllipsoid
     Sphere

harmonica/__init__.py

@@ -16,13 +16,14 @@
 from ._equivalent_sources.spherical import EquivalentSourcesSph
 from ._euler_deconvolution import EulerDeconvolution
 from ._forward.dipole import dipole_magnetic
-from ._forward.ellipsoid_gravity import ellipsoid_gravity
-from ._forward.ellipsoid_magnetics import ellipsoid_magnetic
 from ._forward.ellipsoids import (
     OblateEllipsoid,
     ProlateEllipsoid,
     Sphere,
     TriaxialEllipsoid,
+    create_ellipsoid,
+    ellipsoid_gravity,
+    ellipsoid_magnetic,
 )
 from ._forward.point import point_gravity
 from ._forward.prism_gravity import prism_gravity

harmonica/_forward/ellipsoids/__init__.py

@@ -0,0 +1,19 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Forward modelling of ellipsoids.
+"""
+
+from .ellipsoids import (
+    OblateEllipsoid,
+    ProlateEllipsoid,
+    Sphere,
+    TriaxialEllipsoid,
+    create_ellipsoid,
+)
+from .gravity import ellipsoid_gravity
+from .magnetic import ellipsoid_magnetic

harmonica/_forward/ellipsoids.py …monica/_forward/ellipsoids/ellipsoids.pyharmonica/_forward/ellipsoids.py renamed to harmonica/_forward/ellipsoids/ellipsoids.py harmonica/_forward/ellipsoids.py renamed to harmonica/_forward/ellipsoids/ellipsoids.py

@@ -14,7 +14,134 @@
 import numpy as np
 import numpy.typing as npt
 
-from .utils import get_rotation_matrix
+from ..utils import get_rotation_matrix
+
+
+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:

harmonica/_forward/ellipsoid_gravity.py harmonica/_forward/ellipsoids/gravity.pyharmonica/_forward/ellipsoid_gravity.py renamed to harmonica/_forward/ellipsoids/gravity.py

@@ -15,9 +15,9 @@
 import numpy.typing as npt
 from choclo.constants import GRAVITATIONAL_CONST
 
-from ..errors import NoPhysicalPropertyWarning
-from ..typing import Coordinates, Ellipsoid
-from .utils_ellipsoids import calculate_lambda, get_elliptical_integrals, is_internal
+from ...errors import NoPhysicalPropertyWarning
+from ...typing import Coordinates, Ellipsoid
+from .utils import calculate_lambda, get_elliptical_integrals, is_internal
 
 
 def ellipsoid_gravity(

harmonica/_forward/ellipsoid_magnetics.py harmonica/_forward/ellipsoids/magnetic.pyharmonica/_forward/ellipsoid_magnetics.py renamed to harmonica/_forward/ellipsoids/magnetic.py

@@ -18,10 +18,10 @@
 from scipy.constants import mu_0
 from scipy.special import ellipeinc, ellipkinc
 
-from .._utils import magnetic_angles_to_vec
-from ..errors import NoPhysicalPropertyWarning
-from ..typing import Coordinates, Ellipsoid
-from .utils_ellipsoids import (
+from ..._utils import magnetic_angles_to_vec
+from ...errors import NoPhysicalPropertyWarning
+from ...typing import Coordinates, Ellipsoid
+from .utils import (
     calculate_lambda,
     get_derivatives_of_elliptical_integrals,
     get_elliptical_integrals,

harmonica/_forward/utils_ellipsoids.py harmonica/_forward/ellipsoids/utils.pyharmonica/_forward/utils_ellipsoids.py renamed to harmonica/_forward/ellipsoids/utils.py

@@ -0,0 +1,6 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#

harmonica/tests/ellipsoids/__init__.py

@@ -8,12 +8,19 @@
 Test ellipsoid classes.
 """
 
+import itertools
 import re
 
 import numpy as np
 import pytest
 
-from harmonica import OblateEllipsoid, ProlateEllipsoid, Sphere, TriaxialEllipsoid
+from harmonica import (
+    OblateEllipsoid,
+    ProlateEllipsoid,
+    Sphere,
+    TriaxialEllipsoid,
+    create_ellipsoid,
+)
 
 
 class TestProlateEllipsoid:
@@ -434,3 +441,106 @@ def test_invalid_remanent_mag_shape(self, ellipsoid_class, ellipsoid_args):
         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)

harmonica/tests/test_ellipsoids.py …nica/tests/ellipsoids/test_ellipsoids.pyharmonica/tests/test_ellipsoids.py renamed to harmonica/tests/ellipsoids/test_ellipsoids.py harmonica/tests/test_ellipsoids.py renamed to harmonica/tests/ellipsoids/test_ellipsoids.py

@@ -21,18 +21,14 @@
 import pytest
 import verde as vd
 
-from harmonica import point_gravity
-from harmonica.errors import NoPhysicalPropertyWarning
-
-from .._forward.ellipsoid_gravity import (
-    ellipsoid_gravity,
-)
-from .._forward.ellipsoids import (
+from harmonica import ellipsoid_gravity, point_gravity
+from harmonica._forward.ellipsoids import (
     OblateEllipsoid,
     ProlateEllipsoid,
     Sphere,
     TriaxialEllipsoid,
 )
+from harmonica.errors import NoPhysicalPropertyWarning
 
 
 def build_ellipsoid(ellipsoid_type, *, center=(0, 0, 0), density=None):