Add method for computing prime vertical radius (#35)

Add `boule.Ellipsoid.prime_vertical_radius` (usually represented by N) instead
of calculating it in the coordinate conversion method. This is useful in other
areas (e.g. fatiando/harmonica#154). Method uses pre-computed sine of latitude
to avoid repeated computations of sines in other methods using this. The
definition is from Vermeille (2002; https://doi.org/10.1007/s00190-002-0273-6)
and Vajda (2004).

boule/ellipsoid.py

@@ -146,6 +146,35 @@ def gravity_pole(self):
  )
  return result
 
+ def prime_vertical_radius(self, sinlat):
+ r"""
+ Calculate the prime vertical radius for a given geodetic latitude
+
+ The prime vertical radius is defined as:
+
+ .. math::
+
+ N(\phi) = \frac{a}{\sqrt{1 - e^2 \sin^2(\phi)}}
+
+ Where :math:`a` is the semimajor axis and :math:`e` is the first eccentricity.
+
+ This function receives the sine of the latitude as input to avoid repeated
+ computations of trigonometric functions.
+
+ Parameters
+ ----------
+ sinlat : float or array-like
+ Sine of the latitude angle.
+
+ Returns
+ -------
+ prime_vertical_radius : float or array-like
+ Prime vertical radius given in the same units as the semimajor axis
+ """
+ return self.semimajor_axis / np.sqrt(
+ 1 - self.first_eccentricity ** 2 * sinlat ** 2
+ )
+
  def geodetic_to_spherical(self, longitude, latitude, height):
  """
  Convert from geodetic to geocentric spherical coordinates.
@@ -176,15 +205,14 @@ def geodetic_to_spherical(self, longitude, latitude, height):
 
  """
  latitude_rad = np.radians(latitude)
- prime_vertical_radius = self.semimajor_axis / np.sqrt(
- 1 - self.first_eccentricity ** 2 * np.sin(latitude_rad) ** 2
- )
+ coslat, sinlat = np.cos(latitude_rad), np.sin(latitude_rad)
+ prime_vertical_radius = self.prime_vertical_radius(sinlat)
  # Instead of computing X and Y, we only compute the projection on the
  # XY plane: xy_projection = sqrt( X**2 + Y**2 )
- xy_projection = (height + prime_vertical_radius) * np.cos(latitude_rad)
+ xy_projection = (height + prime_vertical_radius) * coslat
  z_cartesian = (
  height + (1 - self.first_eccentricity ** 2) * prime_vertical_radius
- ) * np.sin(latitude_rad)
+ ) * sinlat
  radius = np.sqrt(xy_projection ** 2 + z_cartesian ** 2)
  spherical_latitude = np.degrees(np.arcsin(z_cartesian / radius))
  return longitude, spherical_latitude, radius

boule/tests/test_ellipsoid.py

@@ -171,3 +171,40 @@ def test_normal_gravity_non_zero_height(ellipsoid):
  assert gamma_pole < ellipsoid.normal_gravity(90, -1000)
  assert gamma_pole < ellipsoid.normal_gravity(-90, -1000)
  assert gamma_eq < ellipsoid.normal_gravity(0, -1000)
+
+
+@pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES)
+def test_prime_vertical_radius(ellipsoid):
+ r"""
+ Check prime vertical radius on equator, poles and 45 degrees
+
+ The prime vertical radius can be also expressed in terms of the semi-major and
+ semi-minor axis:
+
+ .. math:
+
+ N(\phi) = \frac{a^2}{\sqrt{a^2 \cos^2 \phi + b^2 \sin^2 \phi}}
+ """
+ # Compute prime vertical radius on the equator and the poles
+ latitudes = np.array([0, 90, -90, 45])
+ prime_vertical_radii = ellipsoid.prime_vertical_radius(
+ np.sin(np.radians(latitudes))
+ )
+ # Computed expected values
+ prime_vertical_radius_equator = ellipsoid.semimajor_axis
+ prime_vertical_radius_pole = (
+ ellipsoid.semimajor_axis ** 2 / ellipsoid.semiminor_axis
+ )
+ prime_vertical_radius_45 = ellipsoid.semimajor_axis ** 2 / np.sqrt(
+ 0.5 * ellipsoid.semimajor_axis ** 2 + 0.5 * ellipsoid.semiminor_axis ** 2
+ )
+ expected_pvr = np.array(
+ [
+ prime_vertical_radius_equator,
+ prime_vertical_radius_pole,
+ prime_vertical_radius_pole,
+ prime_vertical_radius_45,
+ ]
+ )
+ # Compare calculated vs expected values
+ npt.assert_allclose(prime_vertical_radii, expected_pvr)