Merge branch 'GEOMETRY-71__Matt'

Closes #57.
diff --git a/commons-geometry-core/src/main/java/org/apache/commons/geometry/core/Region.java b/commons-geometry-core/src/main/java/org/apache/commons/geometry/core/Region.java
index 233cd4b..2c15111 100644
--- a/commons-geometry-core/src/main/java/org/apache/commons/geometry/core/Region.java
+++ b/commons-geometry-core/src/main/java/org/apache/commons/geometry/core/Region.java
@@ -49,8 +49,8 @@
      */
     double getBoundarySize();
 
-    /** Get the barycenter of the region or null if none exists. A barycenter
-     * will not exist for empty or infinite regions.
+    /** Get the barycenter of the region or null if no single, unique barycenter exists.
+     * A barycenter will not exist for empty or infinite regions.
      * @return the barycenter of the region or null if none exists
      */
     P getBarycenter();
diff --git a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/RegionBSPTree1S.java b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/RegionBSPTree1S.java
index b056043..2baeeb2 100644
--- a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/RegionBSPTree1S.java
+++ b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/RegionBSPTree1S.java
@@ -22,7 +22,6 @@
 import java.util.List;
 import java.util.Objects;
 
-import org.apache.commons.numbers.angle.PlaneAngleRadians;
 import org.apache.commons.geometry.core.Transform;
 import org.apache.commons.geometry.core.partitioning.Hyperplane;
 import org.apache.commons.geometry.core.partitioning.HyperplaneLocation;
@@ -31,6 +30,8 @@
 import org.apache.commons.geometry.core.partitioning.bsp.AbstractBSPTree;
 import org.apache.commons.geometry.core.partitioning.bsp.AbstractRegionBSPTree;
 import org.apache.commons.geometry.core.precision.DoublePrecisionContext;
+import org.apache.commons.geometry.euclidean.twod.Vector2D;
+import org.apache.commons.numbers.angle.PlaneAngleRadians;
 
 /** BSP tree representing regions in 1D spherical space.
  */
@@ -337,10 +338,12 @@
     protected RegionSizeProperties<Point1S> computeRegionSizeProperties() {
         if (isFull()) {
             return new RegionSizeProperties<>(PlaneAngleRadians.TWO_PI, null);
+        } else if (isEmpty()) {
+            return new RegionSizeProperties<>(0, null);
         }
 
         double size = 0;
-        double scaledBarycenterSum = 0;
+        Vector2D scaledBarycenterSum = Vector2D.ZERO;
 
         double intervalSize;
 
@@ -348,12 +351,14 @@
             intervalSize = interval.getSize();
 
             size += intervalSize;
-            scaledBarycenterSum += intervalSize * interval.getBarycenter().getNormalizedAzimuth();
+            scaledBarycenterSum = scaledBarycenterSum.add(interval.getBarycenter().getVector().withNorm(intervalSize));
         }
 
-        final Point1S barycenter = size > 0 ?
-                Point1S.of(scaledBarycenterSum / size) :
-                null;
+        final DoublePrecisionContext precision = ((CutAngle) getRoot().getCutHyperplane()).getPrecision();
+
+        final Point1S barycenter = scaledBarycenterSum.eq(Vector2D.ZERO, precision) ?
+                 null :
+                 Point1S.from(scaledBarycenterSum);
 
         return new RegionSizeProperties<>(size, barycenter);
     }
diff --git a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
index c187cbb..4c56473 100644
--- a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
+++ b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
@@ -24,7 +24,6 @@
 import java.util.stream.Collectors;
 import java.util.stream.Stream;
 
-import org.apache.commons.numbers.angle.PlaneAngleRadians;
 import org.apache.commons.geometry.core.Transform;
 import org.apache.commons.geometry.core.partitioning.AbstractConvexHyperplaneBoundedRegion;
 import org.apache.commons.geometry.core.partitioning.ConvexSubHyperplane;
@@ -32,6 +31,7 @@
 import org.apache.commons.geometry.core.partitioning.Split;
 import org.apache.commons.geometry.core.precision.DoublePrecisionContext;
 import org.apache.commons.geometry.euclidean.threed.Vector3D;
+import org.apache.commons.numbers.angle.PlaneAngleRadians;
 
 /** Class representing a convex area in 2D spherical space. The boundaries of this
  * area, if any, are composed of convex great circle arcs.
@@ -128,16 +128,19 @@
     /** {@inheritDoc} */
     @Override
     public Point2S getBarycenter() {
-        Vector3D weighted = getWeightedBarycenter();
+        Vector3D weighted = getWeightedBarycenterVector();
         return weighted == null ? null : Point2S.from(weighted);
     }
 
-    /**
-     * Returns the weighted vector for barycenter.
-     *
-     * @return Weighted barycenter
+    /** Returns the weighted vector for the barycenter. This vector is computed by scaling the
+     * pole vector of the great circle of each boundary arc by the size of the arc and summing
+     * the results. By combining the weighted barycenter vectors of multiple areas, a single
+     * barycenter can be computed for the whole group.
+     * @return weighted barycenter vector.
+     * @see <a href="https://archive.org/details/centroidinertiat00broc">
+     *  <em>The Centroid and Inertia Tensor for a Spherical Triangle</em> - John E. Brock</a>
      */
-    Vector3D getWeightedBarycenter() {
+    Vector3D getWeightedBarycenterVector() {
         List<GreatArc> arcs = getBoundaries();
         switch (arcs.size()) {
         case 0:
diff --git a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java
index 17ddd45..b304186 100644
--- a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java
+++ b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java
@@ -161,19 +161,18 @@
             return new RegionSizeProperties<>(0, null);
         }
 
-        List<ConvexArea2S> areas = toConvex();
-        DoublePrecisionContext precision = ((GreatArc) getRoot().getCut()).getPrecision();
+        final List<ConvexArea2S> areas = toConvex();
+        final DoublePrecisionContext precision = ((GreatArc) getRoot().getCut()).getPrecision();
 
         double sizeSum = 0;
         Vector3D barycenterVector = Vector3D.ZERO;
 
-        double size;
         for (ConvexArea2S area : areas) {
             sizeSum += area.getSize();
-            barycenterVector = barycenterVector.add(area.getWeightedBarycenter());
+            barycenterVector = barycenterVector.add(area.getWeightedBarycenterVector());
         }
 
-        Point2S barycenter = barycenterVector.eq(Vector3D.ZERO, precision) ?
+        final Point2S barycenter = barycenterVector.eq(Vector3D.ZERO, precision) ?
                 null :
                 Point2S.from(barycenterVector);
 
diff --git a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/RegionBSPTree1STest.java b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/RegionBSPTree1STest.java
index db37af8..d5fe5e1 100644
--- a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/RegionBSPTree1STest.java
+++ b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/RegionBSPTree1STest.java
@@ -24,6 +24,7 @@
 import org.apache.commons.geometry.core.partitioning.SplitLocation;
 import org.apache.commons.geometry.core.precision.DoublePrecisionContext;
 import org.apache.commons.geometry.core.precision.EpsilonDoublePrecisionContext;
+import org.apache.commons.geometry.euclidean.twod.Vector2D;
 import org.apache.commons.numbers.angle.PlaneAngleRadians;
 import org.junit.Assert;
 import org.junit.Test;
@@ -713,7 +714,23 @@
         // act/assert
         Assert.assertEquals(0.6, tree.getSize(), TEST_EPS);
         Assert.assertEquals(0, tree.getBoundarySize(), TEST_EPS);
-        Assert.assertEquals(2.2 / 6, tree.getBarycenter().getAzimuth(), TEST_EPS);
+
+        Vector2D barycenterVector = Point1S.of(0.1).getVector().withNorm(0.2)
+                .add(Point1S.of(0.5).getVector().withNorm(0.4));
+        Assert.assertEquals(Point1S.from(barycenterVector).getAzimuth(), tree.getBarycenter().getAzimuth(), TEST_EPS);
+    }
+
+    @Test
+    public void testRegionProperties_equalAndOppositeIntervals() {
+        // arrange
+        RegionBSPTree1S tree = RegionBSPTree1S.empty();
+        tree.add(AngularInterval.of(-1, 1, TEST_PRECISION));
+        tree.add(AngularInterval.of(Math.PI - 1, Math.PI + 1, TEST_PRECISION));
+
+        // act/assert
+        Assert.assertEquals(4, tree.getSize(), TEST_EPS);
+        Assert.assertEquals(0, tree.getBoundarySize(), TEST_EPS);
+        Assert.assertNull(tree.getBarycenter()); // no unique barycenter exists
     }
 
     @Test
diff --git a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java
index 0307804..759e9db 100644
--- a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java
+++ b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java
@@ -42,10 +42,6 @@
     private static final DoublePrecisionContext TEST_PRECISION =
             new EpsilonDoublePrecisionContext(TEST_EPS);
 
-    // epsilon value for use when comparing computed barycenter locations;
-    // this must currently be set much higher than the other epsilon
-    private static final double BARYCENTER_EPS = 1e-2;
-
     @Test
     public void testFull() {
         // act
@@ -257,7 +253,7 @@
         Assert.assertEquals(size, area.getSize(), TEST_EPS);
 
         Point2S expectedBarycenter = triangleBarycenter(p1, p2, p3);
-        SphericalTestUtils.assertPointsEq(expectedBarycenter, area.getBarycenter(), BARYCENTER_EPS);
+        SphericalTestUtils.assertPointsEq(expectedBarycenter, area.getBarycenter(), TEST_EPS);
 
         checkBarycenterConsistency(area);
 
@@ -759,11 +755,18 @@
     }
 
     private static Point2S triangleBarycenter(Point2S p1, Point2S p2, Point2S p3) {
-        // compute the barycenter using intersection mid point arcs
-        GreatCircle c1 = GreatCircle.fromPoints(p1, p2.slerp(p3, 0.5), TEST_PRECISION);
-        GreatCircle c2 = GreatCircle.fromPoints(p2, p1.slerp(p3, 0.5), TEST_PRECISION);
+        // compute the barycenter as the sum of the cross product of each point pair weighted by
+        // the angle between the points
+        Vector3D v1 = p1.getVector();
+        Vector3D v2 = p2.getVector();
+        Vector3D v3 = p3.getVector();
 
-        return c1.intersection(c2);
+        Vector3D sum = Vector3D.ZERO;
+        sum = sum.add(v1.cross(v2).withNorm(v1.angle(v2)));
+        sum = sum.add(v2.cross(v3).withNorm(v2.angle(v3)));
+        sum = sum.add(v3.cross(v1).withNorm(v3.angle(v1)));
+
+        return Point2S.from(sum);
     }
 
     private static void checkArc(GreatArc arc, Point2S start, Point2S end) {
@@ -802,21 +805,15 @@
 
             ConvexArea2S minus = split.getMinus();
             double minusSize = minus.getSize();
-            Point2S minusBc = minus.getBarycenter();
-
-            Vector3D weightedMinus = minusBc.getVector()
-                    .multiply(minus.getSize());
 
             ConvexArea2S plus = split.getPlus();
             double plusSize = plus.getSize();
-            Point2S plusBc = plus.getBarycenter();
 
-            Vector3D weightedPlus = plusBc.getVector()
-                    .multiply(plus.getSize());
-            Point2S computedBarycenter = Point2S.from(weightedMinus.add(weightedPlus));
+            Point2S computedBarycenter = Point2S.from(minus.getWeightedBarycenterVector()
+                    .add(plus.getWeightedBarycenterVector()));
 
             Assert.assertEquals(size, minusSize + plusSize, TEST_EPS);
-            SphericalTestUtils.assertPointsEq(barycenter, computedBarycenter, BARYCENTER_EPS);
+            SphericalTestUtils.assertPointsEq(barycenter, computedBarycenter, TEST_EPS);
         }
     }
 }
diff --git a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java
index 3e0e290..a7040e9 100644
--- a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java
+++ b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java
@@ -39,13 +39,13 @@
 
     private static final double TEST_EPS = 1e-10;
 
+    // alternative epsilon value for checking the barycenters of complex
+    // or very small regions
+    private static final double BARYCENTER_EPS = 1e-5;
+
     private static final DoublePrecisionContext TEST_PRECISION =
             new EpsilonDoublePrecisionContext(TEST_EPS);
 
-    // epsilon value for use when comparing computed barycenter locations;
-    // this must currently be set much higher than the other epsilon
-    private static final double BARYCENTER_EPS = 1e-2;
-
     private static final GreatCircle EQUATOR = GreatCircle.fromPoleAndU(
             Vector3D.Unit.PLUS_Z, Vector3D.Unit.PLUS_X, TEST_PRECISION);
 
@@ -339,7 +339,7 @@
     }
 
     @Test
-    public void testToConvex_doubleLune_comlement() {
+    public void testToConvex_doubleLune_complement() {
         // arrange
         RegionBSPTree2S tree = GreatArcPath.builder(TEST_PRECISION)
                 .append(EQUATOR.arc(0,  PlaneAngleRadians.PI))
@@ -523,12 +523,12 @@
         Assert.assertEquals(3.5 * PlaneAngleRadians.PI, tree.getSize(), TEST_EPS);
         Assert.assertEquals(1.5 * PlaneAngleRadians.PI, tree.getBoundarySize(), TEST_EPS);
 
-//        Point2S center = Point2S.from(Point2S.MINUS_K.getVector()
-//                .add(Point2S.PLUS_I.getVector())
-//                .add(Point2S.MINUS_J.getVector()));
-//        SphericalTestUtils.assertPointsEq(center.antipodal(), tree.getBarycenter(), TEST_EPS);
-//
-//        checkBarycenterConsistency(tree);
+        Point2S center = Point2S.from(Point2S.MINUS_K.getVector()
+                .add(Point2S.PLUS_I.getVector())
+                .add(Point2S.MINUS_J.getVector()));
+        SphericalTestUtils.assertPointsEq(center.antipodal(), tree.getBarycenter(), TEST_EPS);
+
+        checkBarycenterConsistency(tree);
 
         List<GreatArcPath> paths = tree.getBoundaryPaths();
         Assert.assertEquals(1, paths.size());
@@ -543,6 +543,105 @@
     }
 
     @Test
+    public void testGeometricProperties_polygonWithHole() {
+        // arrange
+        Point2S center = Point2S.of(0.5, 2);
+
+        double outerRadius = 1;
+        double innerRadius = 0.5;
+
+        RegionBSPTree2S outer = buildDiamond(center, outerRadius);
+        RegionBSPTree2S inner = buildDiamond(center, innerRadius);
+
+        // rotate the inner diamond a quarter turn to become a square
+        inner.transform(Transform2S.createRotation(center, 0.25 * Math.PI));
+
+        // act
+        RegionBSPTree2S tree = RegionBSPTree2S.empty();
+        tree.difference(outer, inner);
+
+        // assert
+        double area = 4 * (rightTriangleArea(outerRadius, outerRadius) - rightTriangleArea(innerRadius, innerRadius));
+        Assert.assertEquals(area, tree.getSize(), TEST_EPS);
+
+        double outerSideLength = sphericalHypot(outerRadius, outerRadius);
+        double innerSideLength = sphericalHypot(innerRadius, innerRadius);
+        double boundarySize = 4 * (outerSideLength + innerSideLength);
+        Assert.assertEquals(boundarySize, tree.getBoundarySize(), TEST_EPS);
+
+        SphericalTestUtils.assertPointsEq(center, tree.getBarycenter(), TEST_EPS);
+        checkBarycenterConsistency(tree);
+
+        SphericalTestUtils.checkClassify(tree, RegionLocation.OUTSIDE, center);
+    }
+
+    @Test
+    public void testGeometricProperties_polygonWithHole_complex() {
+        // arrange
+        Point2S center = Point2S.of(0.5, 2);
+
+        double outerRadius = 2;
+        double midRadius = 1;
+        double innerRadius = 0.5;
+
+        RegionBSPTree2S outer = buildDiamond(center, outerRadius);
+        RegionBSPTree2S mid = buildDiamond(center, midRadius);
+        RegionBSPTree2S inner = buildDiamond(center, innerRadius);
+
+        // rotate the middle diamond a quarter turn to become a square
+        mid.transform(Transform2S.createRotation(center, 0.25 * Math.PI));
+
+        // act
+        RegionBSPTree2S tree = RegionBSPTree2S.empty();
+        tree.difference(outer, mid);
+        tree.union(inner);
+        tree.complement();
+
+        // assert
+        // compute the area, adjusting the first computation for the fact that the triangles comprising the
+        // outer diamond have lengths greater than pi/2
+        double nonComplementedArea = 4 * ((PlaneAngleRadians.PI - rightTriangleArea(outerRadius, outerRadius) -
+                rightTriangleArea(midRadius, midRadius) + rightTriangleArea(innerRadius, innerRadius)));
+        double area = (4 * PlaneAngleRadians.PI) - nonComplementedArea;
+        Assert.assertEquals(area, tree.getSize(), TEST_EPS);
+
+        double outerSideLength = sphericalHypot(outerRadius, outerRadius);
+        double midSideLength = sphericalHypot(midRadius, midRadius);
+        double innerSideLength = sphericalHypot(innerRadius, innerRadius);
+        double boundarySize = 4 * (outerSideLength + midSideLength + innerSideLength);
+        Assert.assertEquals(boundarySize, tree.getBoundarySize(), TEST_EPS);
+
+        SphericalTestUtils.assertPointsEq(center.antipodal(), tree.getBarycenter(), TEST_EPS);
+        checkBarycenterConsistency(tree);
+
+        SphericalTestUtils.checkClassify(tree, RegionLocation.OUTSIDE, center);
+    }
+
+    @Test
+    public void testGeometricProperties_equalAndOppositeRegions() {
+        // arrange
+        Point2S center = Point2S.PLUS_I;
+        double radius = 0.25 * Math.PI;
+
+        RegionBSPTree2S a = buildDiamond(center, radius);
+        RegionBSPTree2S b = buildDiamond(center.antipodal(), radius);
+
+        // act
+        RegionBSPTree2S tree = RegionBSPTree2S.empty();
+        tree.union(a, b);
+
+        // assert
+        double area = 8 * rightTriangleArea(radius, radius);
+        Assert.assertEquals(area, tree.getSize(), TEST_EPS);
+
+        double boundarySize = 8 * sphericalHypot(radius, radius);
+        Assert.assertEquals(boundarySize, tree.getBoundarySize(), TEST_EPS);
+
+        // should be null since no unique barycenter exists
+        Assert.assertNull(tree.getBarycenter());
+    }
+
+    @Test
     public void testSplit_both() {
         // arrange
         GreatCircle c1 = GreatCircle.fromPole(Vector3D.Unit.MINUS_X, TEST_PRECISION);
@@ -700,6 +799,7 @@
         // assert
         Assert.assertEquals(0.6316801448267251, france.getBoundarySize(), TEST_EPS);
         Assert.assertEquals(0.013964220234478741, france.getSize(), TEST_EPS);
+
         SphericalTestUtils.assertPointsEq(Point2S.of(0.04368552749392928, 0.7590839905197961),
                 france.getBarycenter(), BARYCENTER_EPS);
 
@@ -809,8 +909,6 @@
         Point2S barycenter = region.getBarycenter();
         double size = region.getSize();
 
-        SphericalTestUtils.checkClassify(region, RegionLocation.INSIDE, barycenter);
-
         GreatCircle circle = GreatCircle.fromPole(barycenter.getVector(), TEST_PRECISION);
         for (double az = 0; az <= PlaneAngleRadians.TWO_PI; az += 0.2) {
             Point2S pt = circle.toSpace(Point1S.of(az));
@@ -822,24 +920,77 @@
 
             RegionBSPTree2S minus = split.getMinus();
             double minusSize = minus.getSize();
-            Point2S minusBc = minus.getBarycenter();
-
-            Vector3D weightedMinus = minusBc.getVector()
-                    .multiply(minus.getSize());
 
             RegionBSPTree2S plus = split.getPlus();
             double plusSize = plus.getSize();
-            Point2S plusBc = plus.getBarycenter();
 
-            Vector3D weightedPlus = plusBc.getVector()
-                    .multiply(plus.getSize());
-            Point2S computedBarycenter = Point2S.from(weightedMinus.add(weightedPlus));
+            Point2S computedBarycenter = Point2S.from(weightedBarycenterVector(minus)
+                    .add(weightedBarycenterVector(plus)));
 
             Assert.assertEquals(size, minusSize + plusSize, TEST_EPS);
-            SphericalTestUtils.assertPointsEq(barycenter, computedBarycenter, BARYCENTER_EPS);
+            SphericalTestUtils.assertPointsEq(barycenter, computedBarycenter, TEST_EPS);
         }
     }
 
+    private static Vector3D weightedBarycenterVector(RegionBSPTree2S tree) {
+        Vector3D sum = Vector3D.ZERO;
+        for (ConvexArea2S convex : tree.toConvex()) {
+            sum = sum.add(convex.getWeightedBarycenterVector());
+        }
+
+        return sum;
+    }
+
+    private static RegionBSPTree2S buildDiamond(Point2S center, double radius) {
+        Vector3D u = center.getVector();
+        Vector3D w = u.orthogonal(Vector3D.Unit.PLUS_Z);
+        Vector3D v = w.cross(u);
+
+        Transform2S rotV = Transform2S.createRotation(v, radius);
+        Transform2S rotW = Transform2S.createRotation(w, radius);
+
+        Point2S top = rotV.inverse().apply(center);
+        Point2S bottom = rotV.apply(center);
+
+        Point2S right = rotW.apply(center);
+        Point2S left = rotW.inverse().apply(center);
+
+        return GreatArcPath.fromVertexLoop(Arrays.asList(top, left, bottom, right), TEST_PRECISION)
+                .toTree();
+    }
+
+    /** Solve for the hypotenuse of a spherical right triangle, given the lengths of the
+     * other two side. The sides must have lengths less than pi/2.
+     * @param a first side; must be less than pi/2
+     * @param b second side; must be less than pi/2
+     * @return the hypotenuse of the spherical right triangle with sides of the given lengths
+     */
+    private static double sphericalHypot(double a, double b) {
+        // use the spherical law of cosines and the fact that cos(pi/2) = 0
+        // https://en.wikipedia.org/wiki/Spherical_trigonometry#Cosine_rules
+        return Math.acos(Math.cos(a) * Math.cos(b));
+    }
+
+    /**
+     * Compute the area of the spherical right triangle with the given sides. The sides must have lengths
+     * less than pi/2.
+     * @param a first side; must be less than pi/2
+     * @param b second side; must be less than pi/2
+     * @return the area of the spherical right triangle
+     */
+    private static double rightTriangleArea(double a, double b) {
+        double c = sphericalHypot(a, b);
+
+        // use the spherical law of sines to determine the interior angles
+        // https://en.wikipedia.org/wiki/Spherical_trigonometry#Sine_rules
+        double sinC = Math.sin(c);
+        double angleA = Math.asin(Math.sin(a) / sinC);
+        double angleB = Math.asin(Math.sin(b) / sinC);
+
+        // use Girard's theorem
+        return angleA + angleB  - (0.5 * PlaneAngleRadians.PI);
+    }
+
     private static RegionBSPTree2S circleToPolygon(Point2S center, double radius, int numPts, boolean clockwise) {
         List<Point2S> pts = new ArrayList<>(numPts);