Add normal gravity and gravitational potential to sphere class.

boule/_ellipsoid.py

@@ -529,7 +529,7 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height):
             The coordinate u, which is the semiminor axis of the ellipsoid that
             passes through the input coordinates.
         """
-        if np.array(height).all() == 0:
+        if (np.array(height) == 0).all():
             # Use simple formula that relates geodetic and reduced latitude
             beta = np.degrees(
                 np.arctan(
@@ -538,7 +538,7 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height):
                     * np.tan(np.radians(latitude))
                 )
             )
-            u = self.semiminor_axis
+            u = np.full_like(height, fill_value=self.semiminor_axis)
 
             return longitude, beta, u
 
@@ -569,12 +569,13 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height):
             # cos(reduced latitude) squared of the computation point
             cosbeta_p2 = 0.5 + big_r / 2 - np.sqrt(0.25 + big_r**2 / 4 - big_d / 2)
 
-            # Semiminor axis of ellipsoid passing through the computation point
+            # Semiminor axis of the ellipsoid passing through the computation
+            # point. This is the coordinate u
             b_p = np.sqrt(r_p2 + z_p2 - self.linear_eccentricity**2 * cosbeta_p2)
 
             # Note that the sign of beta_p needs to be fixed as it is defined
-            # from -90 to 90 degrees.
-            beta_p = np.sign(latitude) * np.degrees(np.arccos(np.sqrt(cosbeta_p2)))
+            # from -90 to 90 degrees, but arccos is from 0 to 180.
+            beta_p = np.copysign(np.degrees(np.arccos(np.sqrt(cosbeta_p2))), latitude)
 
             return longitude, beta_p, b_p
 
@@ -940,7 +941,7 @@ def centrifugal_potential(self, latitude, height):
                 \omega^2 \left(N(\phi) + h\right)^2 \cos^2(\phi)
 
         in which :math:`N(\phi)` is the prime vertical radius of curvature of
-        the ellipsoid.
+        the ellipsoid and :math:`\omega` is the angular velocity.
         """
         # Pre-compute to avoid repeated calculations
         sinlat = np.sin(np.radians(latitude))

boule/_sphere.py

@@ -204,8 +204,8 @@ def normal_gravity(self, latitude, height, si_units=False):
 
         .. caution::
 
-            The current implementation is only valid for heights on or above
-            the surface of the sphere.
+            These expressions are only valid for heights on or above the
+            surface of the sphere.
 
         Parameters
         ----------
@@ -292,10 +292,9 @@ def normal_gravity(self, latitude, height, si_units=False):
             }
 
         It's worth noting that a sphere under rotation is not in hydrostatic
-        equilibrium. Therefore unlike the oblate ellipsoid, it is not it's own
-        equipotential gravity surface (as is the case for the ellipsoid), the
-        gravity potential is not constant at the surface, and  the normal
-        gravity vector is not normal to the surface of the sphere.
+        equilibrium. Therefore, unlike the oblate ellipsoid, the gravity
+        potential is not constant at the surface, and the normal gravity vector
+        is not normal to the surface of the sphere.
 
         """
         if np.any(height < 0):
@@ -380,9 +379,113 @@ def normal_gravitation(self, height, si_units=False):
 
         return gamma
 
+    def normal_gravity_potential(self, latitude, height):
+        r"""
+        Normal gravity potential of the sphere at the given latitude and
+        height.
+
+        Computes the normal gravity potential (gravitational + centrifugal)
+        generated by the sphere at the given spherical latitude :math:`\theta`
+        and height above the surface of the sphere :math:`h`:
+
+        .. math::
+
+            U(\theta, h) = V(h) + \Phi(\theta, h) = \dfrac{GM}{(R + h)}
+            + \dfrac{1}{2} \omega^2 \left(R + h\right)^2 \cos^2(\theta)
+
+        in which :math:`U = V + \Phi` is the gravity potential of the sphere,
+        :math:`V` is the gravitational potential of the sphere, and
+        :math:`\Phi` is the centrifugal potential.
+
+        .. caution::
+
+            These expressions are only valid for heights on or above the
+            surface of the sphere.
+
+        Parameters
+        ----------
+        latitude : float or array
+            The spherical latitude where the normal gravity will be computed
+            (in degrees).
+        height : float or array
+            The height above the surface of the sphere of the computation point
+            (in meters).
+
+        Returns
+        -------
+        U : float or array
+            The normal gravity potential in m²/s².
+
+        Notes
+        -----
+        A sphere under rotation is not in hydrostatic equilibrium. Therefore,
+        unlike the oblate ellipsoid, it is not its own equipotential gravity
+        surface (as is the case for the ellipsoid), the
+        gravity potential is not constant at the surface, and  the normal
+        gravity vector is not normal to the surface of the sphere.
+
+        """
+        if np.any(height < 0):
+            warn(
+                "Formulas used are valid for points outside the sphere. "
+                "Height must be greater than or equal to zero."
+            )
+
+        radial_distance = self.radius + height
+        big_u = self.geocentric_grav_const / radial_distance
+        big_phi = (
+            0.5
+            * (
+                self.angular_velocity
+                * (self.radius + height)
+                * np.cos(np.radians(latitude))
+            )
+            ** 2
+        )
+
+        return big_u + big_phi
+
+    def normal_gravitational_potential(self, height):
+        r"""
+        Normal gravitational potential at a given height above a sphere.
+
+        .. math::
+
+            V(h) = \dfrac{GM}{(R + h)}
+
+        in which :math:`R` is the sphere radius, :math:`h` is the height above
+        the sphere, :math:`G` is the gravitational constant, and :math:`M` is
+        the mass of the sphere.
+
+        .. caution::
+
+            These expressions are only valid for heights on or above the
+            surface of the sphere.
+
+        Parameters
+        ----------
+        height : float or array
+            The height above the surface of the sphere of the computation point
+            (in meters).
+
+        Returns
+        -------
+        V : float or array
+            The normal gravitational potential in m²/s².
+
+        """
+        if np.any(height < 0):
+            warn(
+                "Formulas used are valid for points outside the sphere. "
+                "Height must be greater than or equal to zero."
+            )
+        radial_distance = self.radius + height
+        return self.geocentric_grav_const / radial_distance
+
     def centrifugal_potential(self, latitude, height):
         r"""
-        Centrifugal potential at the given latitude and height.
+        Centrifugal potential at the given latitude and height above the
+        sphere.
 
         Parameters
         ----------
@@ -400,18 +503,23 @@ def centrifugal_potential(self, latitude, height):
         Notes
         -----
 
-        The centrifugal potential :math:`\Phi` at latitude :math:`\phi` and
+        The centrifugal potential :math:`\Phi` at latitude :math:`\theta` and
         height above the sphere :math:`h` is
 
         .. math::
 
-            \Phi(\phi, h) = \dfrac{1}{2}
-                \omega^2 \left(R + h\right)^2 \cos^2(\phi)
+            \Phi(\theta, h) = \dfrac{1}{2}
+                \omega^2 \left(R + h\right)^2 \cos^2(\theta)
 
-        in which :math:`R` is the sphere radius.
+        in which :math:`R` is the sphere radius and :math:`\omega` is the
+        angular velocity.
         """
-        return (1 / 2) * (
-            self.angular_velocity
-            * (self.radius + height)
-            * np.cos(np.radians(latitude))
-        ) ** 2
+        return (
+            0.5
+            * (
+                self.angular_velocity
+                * (self.radius + height)
+                * np.cos(np.radians(latitude))
+            )
+            ** 2
+        )

boule/tests/test_ellipsoid.py

@@ -170,28 +170,29 @@ def test_spherical_to_geodetic_on_poles(ellipsoid):
 
 @pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES)
 def test_geodetic_to_ellipsoidal_conversions(ellipsoid):
-    "Test the geodetic to ellipsoidal-harmonic coordinate conversions by
-    going from geodetic to ellipsoidal to geodetic."
-    rtol = 1e-6  #  The conversion is not too accurate for large height
+    """
+    Test the geodetic to ellipsoidal-harmonic coordinate conversions by
+    going from geodetic to ellipsoidal and back.
+    """
     size = 5
     geodetic_latitude_in = np.linspace(-90, 90, size)
     height_in = np.zeros(size)
-
     longitude, reduced_latitude, u = ellipsoid.geodetic_to_ellipsoidal_harmonic(
         None, geodetic_latitude_in, height_in
     )
-    longitude, geodetic_latitude_out, height_out = ellipsoid.ellipsoidal_harmonic_to_geodetic(
-        None, reduced_latitude, u
+    longitude, geodetic_latitude_out, height_out = (
+        ellipsoid.ellipsoidal_harmonic_to_geodetic(None, reduced_latitude, u)
     )
     npt.assert_allclose(geodetic_latitude_in, geodetic_latitude_out)
     npt.assert_allclose(height_in, height_out)
 
+    rtol = 1e-5  # The conversion is not too accurate for large heights
     height_in = np.array(size * [1000])
     longitude, reduced_latitude, u = ellipsoid.geodetic_to_ellipsoidal_harmonic(
         None, geodetic_latitude_in, height_in
     )
-    longitude, geodetic_latitude_out, height_out = ellipsoid.ellipsoidal_harmonic_to_geodetic(
-        None, reduced_latitude, u
+    longitude, geodetic_latitude_out, height_out = (
+        ellipsoid.ellipsoidal_harmonic_to_geodetic(None, reduced_latitude, u)
     )
     npt.assert_allclose(geodetic_latitude_in, geodetic_latitude_out, rtol=rtol)
     npt.assert_allclose(height_in, height_out, rtol=rtol)
@@ -379,15 +380,12 @@ def test_normal_gravity_gravitational_centrifugal_potential(ellipsoid):
     Test that the normal gravity potential is equal to the sum of the normal
     gravitational potential and centrifugal potential.
     """
-    latitude = np.array([-90, -45, 0, 45, 90])
-    big_u0 = ellipsoid.normal_gravity_potential(latitude, height=0)
-    big_v0 = ellipsoid.normal_gravitational_potential(latitude, height=0)
-    big_phi0 = ellipsoid.centrifugal_potential(latitude, height=0)
-    height = 1000
+    size = 5
+    latitude = np.array([np.linspace(-90, 90, size)] * 2)
+    height = np.array([[0] * size, [1000] * size])
     big_u = ellipsoid.normal_gravity_potential(latitude, height)
     big_v = ellipsoid.normal_gravitational_potential(latitude, height)
     big_phi = ellipsoid.centrifugal_potential(latitude, height)
-    npt.assert_allclose(big_u0, big_v0 + big_phi0)
     npt.assert_allclose(big_u, big_v + big_phi)
 
 

boule/tests/test_sphere.py

@@ -56,6 +56,12 @@ def test_normal_gravity_computed_on_internal_point(sphere):
     with pytest.warns(UserWarning) as warn:
         sphere.normal_gravity(latitude, height=-10)
         assert len(warn) >= 1
+    with warnings.catch_warnings(record=True) as warn:
+        sphere.normal_gravity_potential(latitude, height=-10)
+        assert len(warn) >= 1
+    with warnings.catch_warnings(record=True) as warn:
+        sphere.normal_gravitational_potential(height=-10)
+        assert len(warn) >= 1
 
 
 def test_check_geocentric_grav_const():
@@ -146,3 +152,17 @@ def test_normal_gravity_only_rotation():
             expected_value,
             sphere.normal_gravity(latitude=45, height=height),
         )
+
+
+def test_normal_gravity_gravitational_centrifugal_potential(sphere):
+    """
+    Test that the normal gravity potential is equal to the sum of the normal
+    gravitational potential and centrifugal potential.
+    """
+    size = 5
+    latitude = np.array([np.linspace(-90, 90, size)] * 2)
+    height = np.array([[0] * size, [1000] * size])
+    big_u = sphere.normal_gravity_potential(latitude, height)
+    big_v = sphere.normal_gravitational_potential(height)
+    big_phi = sphere.centrifugal_potential(latitude, height)
+    npt.assert_allclose(big_u, big_v + big_phi)