Add area and area_equivalent_radius to all three ellipsoid classes (#178

)

Add the properties `area` and `area_equivalent_radius` to the three 
ellipsoid classes. The "area_equivalent_radius" is what WGS84
calls `R_2`.

The equations for the area of bi-axial and tri-axial ellipsoids are
taken from https://en.wikipedia.org/wiki/Ellipsoid#Surface_area. The
equation for a bi-axial ellispoid is relatively simple, and only
requires to use `np.arctanh`. The equation for a triaxial ellipsoid,
however, requires computing an elliptic integral. To do this, I have
used the elliptic integral function from scipy: `scipy.special.elliprg`. 
This required adding scipy>=1.8.0 as a dependency.

I have verified that the area equivalent radius for WGS84 in boule is
what is given in the literature. I have verified that the triaxial
ellipsoid routine gives the correct value when a=b for WGS84.

---------

Co-authored-by: Leonardo Uieda <leo@uieda.com>

AUTHORS.md

@@ -10,3 +10,4 @@ order by last name) and are considered "The Boule Developers":
 * [Agustina Pesce](https://github.com/aguspesce) - Unaffiliated (ORCID: [0000-0002-5538-8845](https://orcid.org//0000-0002-5538-8845))
 * [Santiago Soler](https://github.com/santisoler) - CONICET, Argentina; Instituto Geofísico Sismológico Volponi, Universidad Nacional de San Juan, Argentina (ORCID: [0000-0001-9202-5317](https://www.orcid.org/0000-0001-9202-5317))
 * [Leonardo Uieda](https://github.com/leouieda) - Universidade de São Paulo, Brazil (ORCID: 0000-0001-6123-9515)
+* [Mark Wieczorek](https://github.com/MarkWieczorek) - Institut de physique du globe de Paris, France (ORCID: [0000-0001-7007-4222](https://orcid.org/0000-0001-7007-4222))

boule/_ellipsoid.py

@@ -113,6 +113,10 @@ class Ellipsoid:
     5513 kg/m³
     >>> print(f"{ellipsoid.volume * 1e-9:.5e} km³")
     1.08321e+12 km³
+    >>> print(f"{ellipsoid.area:.10e} m²")
+    5.1006562172e+14 m²
+    >>> print(f"{ellipsoid.area_equivalent_radius:0.4f} m")
+    6371007.1809 m
     >>> print(f"{ellipsoid.gravity_equator:.10f} m/s²")
     9.7803253359 m/s²
     >>> print(f"{ellipsoid.gravity_pole:.10f} m/s²")
@@ -255,6 +259,27 @@ def semiaxes_mean_radius(self):
         """
         return 1 / 3 * (2 * self.semimajor_axis + self.semiminor_axis)
 
+    @property
+    def area(self):
+        r"""
+        The area of the ellipsoid.
+        Definition: :math:`A = 2 \pi a^2 \left(1 + \dfrac{b^2}{e a^2}
+        \text{arctanh}\,e \right)`.
+        Units: :math:`m^2`.
+        """
+        # see https://en.wikipedia.org/wiki/Ellipsoid#Surface_area
+        return (
+            2
+            * np.pi
+            * self.semimajor_axis**2
+            * (
+                1
+                + (self.semiminor_axis / self.semimajor_axis) ** 2
+                / self.first_eccentricity
+                * np.arctanh(self.first_eccentricity)
+            )
+        )
+
     @property
     def volume(self):
         r"""
@@ -264,6 +289,15 @@ def volume(self):
         """
         return (4 / 3 * np.pi) * self.semimajor_axis**2 * self.semiminor_axis
 
+    @property
+    def area_equivalent_radius(self):
+        r"""
+        The area equivalent radius of the ellipsoid.
+        Definition: :math:`R_2 = \sqrt{A / (4 \pi)}`.
+        Units: :math:`m`.
+        """
+        return np.sqrt(self.area / (4 * np.pi))
+
     @property
     def mass(self):
         r"""
@@ -289,7 +323,7 @@ def volume_equivalent_radius(self):
         Definition: :math:`R_3 = \left(\dfrac{3}{4 \pi} V \right)^{1/3}`.
         Units: :math:`m`.
         """
-        return (self.volume * 3 / 4 / np.pi) ** (1 / 3)
+        return (self.volume * 3 / (4 * np.pi)) ** (1 / 3)
 
     @property
     def _emm(self):

boule/_sphere.py

@@ -107,6 +107,10 @@ class Sphere:
     1.0 m
     >>> print(f"{sphere.volume:.10f} m³")
     4.1887902048 m³
+    >>> print(f"{sphere.area:.10f} m²")
+    12.5663706144 m²
+    >>> print(sphere.area_equivalent_radius)
+    1
     >>> print(f"{sphere.mass:.12e} kg")
     2.996568928577e+10 kg
     >>> print(f"{sphere.mean_density:.0f} kg/m³")
@@ -188,6 +192,15 @@ def first_eccentricity(self):
         """
         return 0
 
+    @property
+    def area(self):
+        r"""
+        The area of the sphere.
+        Definition: :math:`A = 4 \pi r^2`.
+        Units: :math:`m^2`.
+        """
+        return 4 * np.pi * self.radius**2
+
     @property
     def mean_radius(self):
         """
@@ -217,6 +230,15 @@ def volume(self):
         """
         return (4 / 3 * np.pi) * self.radius**3
 
+    @property
+    def area_equivalent_radius(self):
+        """
+        The area equivalent radius of the sphere is equal to its radius.
+        Definition: :math:`R_2 = R`.
+        Units: :math:`m`.
+        """
+        return self.radius
+
     @property
     def mass(self):
         r"""

boule/_triaxialellipsoid.py

@@ -11,6 +11,7 @@
 
 import attr
 import numpy as np
+import scipy.special
 
 from ._constants import G
 
@@ -102,6 +103,10 @@ class TriaxialEllipsoid:
     259813 m
     >>> print(f"{ellipsoid.semiaxes_mean_radius:.0f} m")
     262700 m
+    >>> print(f"{ellipsoid.area:.10e} m²")
+    8.6562393883e+11 m²
+    >>> print(f"{ellipsoid.area_equivalent_radius:0.0f} m")
+    262458 m
     >>> print(f"{ellipsoid.volume_equivalent_radius:.0f} m")
     261115 m
     >>> print(f"{ellipsoid.mass:.10e} kg")
@@ -222,6 +227,26 @@ def semiaxes_mean_radius(self):
         """
         return (self.semimajor_axis + self.semimedium_axis + self.semiminor_axis) / 3
 
+    @property
+    def area(self):
+        r"""
+        The area of the ellipsoid.
+        Definition: :math:`A = 3 V R_G(a^{-2}, b^{-2}, c^{-2})`, in which
+        :math:`R_G` is the completely-symmetric elliptic integral of the second
+        kind.
+        Units: :math:`m^2`.
+        """
+        # see https://en.wikipedia.org/wiki/Ellipsoid#Surface_area
+        return (
+            3
+            * self.volume
+            * scipy.special.elliprg(
+                1 / self.semimajor_axis**2,
+                1 / self.semimedium_axis**2,
+                1 / self.semiminor_axis**2,
+            )
+        )
+
     @property
     def volume(self):
         r"""
@@ -236,6 +261,15 @@ def volume(self):
             * self.semiminor_axis
         )
 
+    @property
+    def area_equivalent_radius(self):
+        r"""
+        The area equivalent radius of the ellipsoid.
+        Definition: :math:`R_2 = \sqrt{A / (4 \pi)}`.
+        Units: :math:`m`.
+        """
+        return np.sqrt(self.area / (4 * np.pi))
+
     @property
     def mass(self):
         r"""
@@ -261,7 +295,7 @@ def volume_equivalent_radius(self):
         Definition: :math:`R_3 = \left(\dfrac{3}{4 \pi} V \right)^{1/3}`.
         Units: :math:`m`.
         """
-        return (self.volume * 3 / 4 / np.pi) ** (1 / 3)
+        return (self.volume * 3 / (4 * np.pi)) ** (1 / 3)
 
     @property
     def equatorial_flattening(self):

environment.yml

@@ -9,6 +9,7 @@ dependencies:
   # Run
   - numpy
   - attrs
+  - scipy
   # Build
   - build
   - twine

pyproject.toml

@@ -31,6 +31,7 @@ requires-python = ">=3.8"
 dependencies = [
     "numpy>=1.19",
     "attrs>=19.3.0",
+    "scipy>=1.8.0",
 ]
 
 [project.urls]