Add reference Sphere class

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)