WIP: Add option to use geocentric latitude with Ellipsoid.normal_gravity

boule/_ellipsoid.py

@@ -281,8 +281,8 @@ def geocentric_radius(self, latitude, geodetic=True):
         r"""
         Radial distance from the center of the ellipsoid to its surface.
 
-        Can be calculated from either geocentric geodetic or geocentric
-        spherical latitudes.
+        Can be calculated from either geodetic or geocentric spherical
+        latitudes.
 
         Parameters
         ----------
@@ -501,14 +501,17 @@ def spherical_to_geodetic(self, longitude, spherical_latitude, radius):
         height = (k + self.eccentricity**2 - 1) / k * hypot_dz
         return longitude, latitude, height
 
-    def normal_gravity(self, latitude, height, si_units=False):
+    def normal_gravity(self, latitude, height, geodetic=True, si_units=False):
         r"""
-        Normal gravity of the ellipsoid at the given latitude and height.
+        Normal gravity of the ellipsoid at the given latitude and geodetic
+        height.
 
         Computes the magnitude of the gradient of the gravity potential
         (gravitational + centrifugal; see [HofmannWellenhofMoritz2006]_)
-        generated by the ellipsoid at the given geodetic latitude :math:`\phi`
-        and height above the ellipsoid :math:`h` (geometric height).
+        generated by the ellipsoid at the given latitude :math:`\phi` and
+        geodetic height above the ellipsoid :math:`h`. The optional parameter
+        geodetic allows to specify whether the input latitude is geodetic
+        latitude (default) or geocentric latitude.
 
         .. math::
 
@@ -529,14 +532,24 @@ def normal_gravity(self, latitude, height, si_units=False):
             These expressions are only valid for heights on or above the
             surface of the ellipsoid.
 
+        .. note::
+
+           The normal gravity computation requires the use of geodetic
+           latitude. If the geodetic latitude is known, it should be input to
+           avoid unnecessary conversions from geocentric to geodetic latitude.
+
         Parameters
         ----------
         latitude : float or array
-            The geodetic latitude where the normal gravity will be computed (in
-            degrees).
+            The geodetic (default) or geocentric latitude of the point on the
+            ellipsoid where the normal gravity will be computed (in degrees).
         height : float or array
-            The ellipsoidal (geometric) height of computation the point (in
-            meters).
+            The geodetic height of the computation point in meters (normal to
+            the ellipsoid surface).
+        geodetic : bool
+            If True (default), will assume that latitudes are geodetic
+            latitudes. Otherwise, will assume that they are geocentric
+            spherical latitudes.
         si_units : bool
             Return the value in mGal (False, default) or m/s² (True)
 
@@ -553,6 +566,25 @@ def normal_gravity(self, latitude, height, si_units=False):
                 "Height must be greater than or equal to zero."
             )
 
+        # Convert from geocentric to geodetic latitude on the ellipsoid surface
+        if geodetic is False:
+            radius = self.geocentric_radius(latitude, geodetic=False)
+            sinlat = np.sin(np.radians(latitude))
+            coslat = np.sqrt(1 - sinlat**2)
+            big_z = radius * sinlat
+            p_0 = radius**2 * coslat**2 / self.semimajor_axis**2
+            q_0 = (1 - self.eccentricity**2) / self.semimajor_axis**2 * big_z**2
+            r_0 = (p_0 + q_0 - self.eccentricity**4) / 6
+            s_0 = self.eccentricity**4 * p_0 * q_0 / 4 / r_0**3
+            t_0 = np.cbrt(1 + s_0 + np.sqrt(2 * s_0 + s_0**2))
+            u_0 = r_0 * (1 + t_0 + 1 / t_0)
+            v_0 = np.sqrt(u_0**2 + q_0 * self.eccentricity**4)
+            w_0 = self.eccentricity**2 * (u_0 + v_0 - q_0) / 2 / v_0
+            k = np.sqrt(u_0 + v_0 + w_0**2) - w_0
+            big_d = k * radius * coslat / (k + self.eccentricity**2)
+            hypot_dz = np.hypot(big_d, big_z)
+            latitude = np.degrees(2 * np.arctan2(big_z, (big_d + hypot_dz)))
+
         # Pre-compute to avoid repeated calculations
         sinlat = np.sin(np.radians(latitude))
         coslat = np.sqrt(1 - sinlat**2)

boule/tests/test_ellipsoid.py

@@ -228,6 +228,20 @@ def test_normal_gravity_non_zero_height(ellipsoid):
     assert gamma_eq < ellipsoid.normal_gravity(0, -1000)
 
 
+@pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES)
+def test_normal_gravity_geocentric_latitude(ellipsoid):
+    "Check equality of normal gravity using geodetic and geocentric latitude"
+    geodetic_lat = np.linspace(-90, 90, 9)
+    normal_grav = ellipsoid.normal_gravity(geodetic_lat, height=0)
+    longitude, geocentric_lat, radius = ellipsoid.geodetic_to_spherical(
+        0.0, geodetic_lat, 0.0
+    )
+    normal_grav_geocentric = ellipsoid.normal_gravity(
+        geocentric_lat, height=0, geodetic=False
+    )
+    npt.assert_allclose(normal_grav, normal_grav_geocentric)
+
+
 @pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES)
 def test_prime_vertical_radius(ellipsoid):
     r"""