Add reference Sphere class (#42)

Add a new Sphere class to represent reference ellipsoids with zero
flattening. The Sphere class is a subclass of Ellipsoid, inheriting its
methods and properties. Some methods were overridden to avoid division
by zero errors due to singularities caused by zero flattening. Override
normal_gravity method for computing the normal gravity generated by the
sphere as the self gravitational field plus the centrifugal term. Raise
warning when initializing an Ellipsoid with zero (or almost zero)
flattening. Update Ellipsoid class example: replace sphere with an
oblate ellipsoid. Add test functions for the overridden and inherited
methods.

boule/__init__.py

@@ -2,6 +2,7 @@
 # Import functions/classes to make the public API
 from . import version
 from .ellipsoid import Ellipsoid
+from .sphere import Sphere
 from .realizations import WGS84, GRS80, MARS
 
 

boule/ellipsoid.py

@@ -1,6 +1,7 @@
 """
 Module for defining and setting the reference ellipsoid.
 """
+from warnings import warn
 import attr
 import numpy as np
 
@@ -41,30 +42,30 @@ class Ellipsoid:
     Examples
     --------
 
-    We can define a reference unit sphere by using 0 as the flattening:
+    We can define an ellipsoid with flattening equal to 0.5 and unit semimajor axis:
 
-    >>> sphere = Ellipsoid(
-    ...     name="sphere",
-    ...     long_name="Unit sphere",
+    >>> ellipsoid = Ellipsoid(
+    ...     name="oblate-ellipsoid",
+    ...     long_name="Oblate Ellipsoid",
     ...     semimajor_axis=1,
-    ...     flattening=0,
+    ...     flattening=0.5,
     ...     geocentric_grav_const=1,
     ...     angular_velocity=0
     ... )
-    >>> print(sphere) # doctest: +ELLIPSIS
-    Ellipsoid(name='sphere', ...)
-    >>> print(sphere.long_name)
-    Unit sphere
-    >>> print("{:.2f}".format(sphere.semiminor_axis))
-    1.00
-    >>> print("{:.2f}".format(sphere.mean_radius))
-    1.00
-    >>> print("{:.2f}".format(sphere.linear_eccentricity))
-    0.00
-    >>> print("{:.2f}".format(sphere.first_eccentricity))
-    0.00
-    >>> print("{:.2f}".format(sphere.second_eccentricity))
-    0.00
+    >>> print(ellipsoid) # doctest: +ELLIPSIS
+    Ellipsoid(name='oblate-ellipsoid', ...)
+    >>> print(ellipsoid.long_name)
+    Oblate Ellipsoid
+    >>> print("{:.2f}".format(ellipsoid.semiminor_axis))
+    0.50
+    >>> print("{:.2f}".format(ellipsoid.mean_radius))
+    0.83
+    >>> print("{:.2f}".format(ellipsoid.linear_eccentricity))
+    0.87
+    >>> print("{:.2f}".format(ellipsoid.first_eccentricity))
+    0.87
+    >>> print("{:.2f}".format(ellipsoid.second_eccentricity))
+    1.73
 
     """
 
@@ -76,6 +77,19 @@ class Ellipsoid:
     long_name = attr.ib(default=None)
     reference = attr.ib(default=None)
 
+    @flattening.validator
+    def check_flattening(
+        self, flattening, value
+    ):  # pylint: disable=no-self-use,unused-argument
+        """
+        Check if flattening is not equal (or almost) zero
+        """
+        warn_msg = "Use boule.Sphere for representing ellipsoids with zero flattening."
+        if value == 0:
+            warn("Flattening equal to zero. " + warn_msg)
+        if value < 1e-7:
+            warn("Flattening '{}' too close to zero. ".format(value) + warn_msg)
+
     @property
     def semiminor_axis(self):
         "The small (polar) axis of the ellipsoid [meters]"

boule/sphere.py

@@ -0,0 +1,161 @@
+"""
+Module for defining and setting the reference sphere.
+"""
+import attr
+import numpy as np
+
+from . import Ellipsoid
+
+
+# Don't let ellipsoid parameters be changed to avoid messing up calculations
+# accidentally.
+@attr.s(frozen=True)
+class Sphere(Ellipsoid):
+    """
+    Reference spherical ellipsoid
+
+    Represents ellipsoids with zero flattening (spheres). Inherits methods and
+    properties of the :class:`boule.Ellipsoid`, guaranteeing no singularities
+    due to zero flattening (and thus zero eccentricity).
+
+    All parameters are in SI units.
+
+    Parameters
+    ----------
+    name : str
+        A short name for the ellipsoid, for example ``'MOON'``.
+    radius : float
+        The radius of the sphere [meters].
+    geocentric_grav_const : float
+        The geocentric gravitational constant (GM) [m^3 s^-2].
+    angular_velocity : float
+        The angular velocity of the rotating ellipsoid (omega) [rad s^-1].
+    long_name : str or None
+        A long name for the ellipsoid, for example ``"Moon Reference System"``
+        (optional).
+    reference : str or None
+        Citation for the ellipsoid parameter values (optional).
+
+    Examples
+    --------
+
+    We can define a unit sphere:
+
+    >>> sphere = Sphere(
+    ...     name="sphere",
+    ...     radius=1,
+    ...     geocentric_grav_const=1,
+    ...     angular_velocity=0,
+    ...     long_name="Spherical Ellipsoid",
+    ... )
+    >>> print(sphere) # doctest: +ELLIPSIS
+    Sphere(name='sphere', ...)
+    >>> print(sphere.long_name)
+    Spherical Ellipsoid
+    >>> print("{:.2f}".format(sphere.semiminor_axis))
+    1.00
+    >>> print("{:.2f}".format(sphere.mean_radius))
+    1.00
+    >>> print("{:.2f}".format(sphere.flattening))
+    0.00
+    >>> print("{:.2f}".format(sphere.linear_eccentricity))
+    0.00
+    >>> print("{:.2f}".format(sphere.first_eccentricity))
+    0.00
+    >>> print("{:.2f}".format(sphere.second_eccentricity))
+    0.00
+    >>> print(sphere.normal_gravity(latitude=0, height=1))
+    25000.0
+    >>> print(sphere.normal_gravity(latitude=90, height=1))
+    25000.0
+    """
+
+    name = attr.ib()
+    radius = attr.ib()
+    geocentric_grav_const = attr.ib()
+    angular_velocity = attr.ib()
+    long_name = attr.ib(default=None)
+    reference = attr.ib(default=None)
+    # semimajor_axis and flattening shouldn't be defined on initialization:
+    #   - semimajor_axis will be equal to radius
+    #   - flattening will be equal to zero
+    semimajor_axis = attr.ib(init=False)
+    flattening = attr.ib(init=False, default=0)
+
+    @semimajor_axis.default
+    def _set_semimajor_axis(self):
+        "The semimajor axis should be the radius"
+        return self.radius
+
+    def normal_gravity(self, latitude, height):
+        """
+        Calculate normal gravity at any latitude and height
+
+        Computes the magnitude of the gradient of the gravity potential
+        (gravitational + centrifugal) generated by the sphere at the given
+        latitude and height.
+
+        Parameters
+        ----------
+        latitude : float or array
+            The latitude where the normal gravity will be computed (in degrees).
+        height : float or array
+            The height (above the sphere) of computation point (in meters).
+
+        Returns
+        -------
+        gamma : float or array
+            The normal gravity in mGal.
+
+        References
+        ----------
+        [Heiskanen-Moritz]_
+        """
+        return self._gravity_sphere(height) + self._centrifugal_force(latitude, height)
+
+    def _gravity_sphere(self, height):
+        """
+        Calculate the gravity generated by a solid sphere (mGal)
+        """
+        return 1e5 * self.geocentric_grav_const / (self.radius + height) ** 2
+
+    def _centrifugal_force(self, latitude, height):
+        """
+        Calculate the centrifugal force due to the rotation of the sphere (mGal)
+        """
+        return 1e5 * (
+            (-1)
+            * self.angular_velocity ** 2
+            * (self.radius + height)
+            * np.cos(np.radians(latitude))
+        )
+
+    @property
+    def gravity_equator(self):
+        """
+        The norm of the gravity vector at the equator on the sphere [m/s^2]
+
+        Overrides the inherited method from :class:`boule.Ellipsoid` to avoid
+        singularities due to zero flattening.
+        """
+        return self._gravity_on_surface
+
+    @property
+    def gravity_pole(self):
+        """
+        The norm of the gravity vector at the poles on the sphere [m/s^2]
+
+        Overrides the inherited method from :class:`boule.Ellipsoid` to avoid
+        singularities due to zero flattening.
+        """
+        return self._gravity_on_surface
+
+    @property
+    def _gravity_on_surface(self):
+        """
+        Compute norm of the gravity vector on the surface of the sphere [m/s^2]
+
+        Due to rotational symmetry, the norm of the gravity vector is the same
+        on every point of the surface.
+        """
+        return self.geocentric_grav_const / self.radius ** 2

boule/tests/test_ellipsoid.py

@@ -2,6 +2,7 @@
 """
 Test the base Ellipsoid class.
 """
+import warnings
 import pytest
 import numpy as np
 import numpy.testing as npt
@@ -12,47 +13,30 @@
 ELLIPSOID_NAMES = [e.name for e in ELLIPSOIDS]
 
 
-@pytest.fixture
-def sphere():
-    "A spherical ellipsoid"
-    ellipsoid = Ellipsoid(
-        name="unit_sphere",
-        semimajor_axis=1.0,
-        flattening=0,
-        geocentric_grav_const=0,
-        angular_velocity=0,
-    )
-    return ellipsoid
-
-
-def test_geodetic_to_spherical_with_spherical_ellipsoid(sphere):
-    "Test geodetic to geocentric conversion if ellipsoid is a sphere."
-    rtol = 1e-10
-    size = 5
-    longitude = np.linspace(0, 180, size)
-    latitude = np.linspace(-90, 90, size)
-    height = np.linspace(-0.2, 0.2, size)
-    sph_longitude, sph_latitude, radius = sphere.geodetic_to_spherical(
-        longitude, latitude, height
-    )
-    npt.assert_allclose(sph_longitude, longitude, rtol=rtol)
-    npt.assert_allclose(sph_latitude, latitude, rtol=rtol)
-    npt.assert_allclose(radius, sphere.mean_radius + height, rtol=rtol)
-
-
-def test_spherical_to_geodetic_with_spherical_ellipsoid(sphere):
-    "Test spherical to geodetic conversion if ellipsoid is a sphere."
-    rtol = 1e-10
-    size = 5
-    spherical_longitude = np.linspace(0, 180, size)
-    spherical_latitude = np.linspace(-90, 90, size)
-    radius = np.linspace(0.8, 1.2, size)
-    longitude, latitude, height = sphere.spherical_to_geodetic(
-        spherical_longitude, spherical_latitude, radius
-    )
-    npt.assert_allclose(spherical_longitude, longitude, rtol=rtol)
-    npt.assert_allclose(spherical_latitude, latitude, rtol=rtol)
-    npt.assert_allclose(radius, sphere.mean_radius + height, rtol=rtol)
+def test_ellipsoid_zero_flattening():
+    """
+    Check if error is raised after passing zero flattening
+    """
+    # Test with zero flattening
+    with warnings.catch_warnings(record=True) as warn:
+        Ellipsoid(
+            name="zero-flattening",
+            semimajor_axis=1,
+            flattening=0,
+            geocentric_grav_const=1,
+            angular_velocity=0,
+        )
+        assert len(warn) >= 1
+    # Test with almost zero flattening
+    with warnings.catch_warnings(record=True) as warn:
+        Ellipsoid(
+            name="almost-zero-flattening",
+            semimajor_axis=1,
+            flattening=1e-8,
+            geocentric_grav_const=1,
+            angular_velocity=0,
+        )
+        assert len(warn) >= 1
 
 
 @pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES)