presto-geospatial/src/main/java/io/prestosql/plugin/geospatial/GeoFunctions.java/GeoFunctions/computeSphericalExcess
class GeoFunctions:    private static double computeSphericalExcess(Polygon polygon, int start, int end)
    {
        // Our calculations rely on not processing the same point twice
        if (polygon.getPoint(end - 1).equals(polygon.getPoint(start))) {
            end = end - 1;
        }

        if (end - start < 3) {
            // A path with less than 3 distinct points is not valid for calculating an area
            throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "Polygon is not valid: a loop contains less then 3 vertices.");
        }

        Point point = new Point();
        // Initialize the calculator with the last point
        polygon.getPoint(end - 1, point);

        double sphericalExcess = 0;
        double courseDelta = 0;
        boolean firstPoint = true;
        double firstInitialBearing = 0;
        double previousFinalBearing = 0;
        double previousPhi = toRadians(point.getY());
        double previousCos = Math.cos(previousPhi);
        double previousSin = Math.sin(previousPhi);
        double previousTan = Math.tan(previousPhi / 2);
        double previousLongitude = toRadians(point.getX());

        for (int i = start; i < end; i++) {
            polygon.getPoint(i, point);
            double phi = toRadians(point.getY());
            double tan = Math.tan(phi / 2);
            double longitude = toRadians(point.getX());

            // We need to check for that specifically
            // Otherwise calculating the bearing is not deterministic
            if (longitude == previousLongitude && phi == previousPhi) {
                throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "Polygon is not valid: it has two identical consecutive vertices");
            }

            double deltaLongitude = longitude - previousLongitude;
            sphericalExcess += 2 * Math.atan2(Math.tan(deltaLongitude / 2) * (previousTan + tan), 1 + previousTan * tan);

            double cos = Math.cos(phi);
            double sin = Math.sin(phi);
            double sinOfDeltaLongitude = Math.sin(deltaLongitude);
            double cosOfDeltaLongitude = Math.cos(deltaLongitude);

            // Initial bearing from previous to current
            double y = sinOfDeltaLongitude * cos;
            double x = previousCos * sin - previousSin * cos * cosOfDeltaLongitude;
            double initialBearing = (Math.atan2(y, x) + 2 * Math.PI) % (2 * Math.PI);

            // Final bearing from previous to current = opposite of bearing from current to previous
            double finalY = -sinOfDeltaLongitude * previousCos;
            double finalX = previousSin * cos - previousCos * sin * cosOfDeltaLongitude;
            double finalBearing = (Math.atan2(finalY, finalX) + PI) % (2 * Math.PI);

            // When processing our first point we don't yet have a previousFinalBearing
            if (firstPoint) {
                // So keep our initial bearing around, and we'll use it at the end
                // with the last final bearing
                firstInitialBearing = initialBearing;
                firstPoint = false;
            }
            else {
                courseDelta += (initialBearing - previousFinalBearing + 3 * Math.PI) % (2 * Math.PI) - PI;
            }

            courseDelta += (finalBearing - initialBearing + 3 * Math.PI) % (2 * Math.PI) - PI;

            previousFinalBearing = finalBearing;
            previousCos = cos;
            previousSin = sin;
            previousPhi = phi;
            previousTan = tan;
            previousLongitude = longitude;
        }

        // Now that we have our last final bearing, we can calculate the remaining course delta
        courseDelta += (firstInitialBearing - previousFinalBearing + 3 * Math.PI) % (2 * Math.PI) - PI;

        // The courseDelta should be 2Pi or - 2Pi, unless a pole is enclosed (and then it should be ~ 0)
        // In which case we need to correct the spherical excess by 2Pi
        if (Math.abs(courseDelta) < PI / 4) {
            sphericalExcess = Math.abs(sphericalExcess) - 2 * Math.PI;
        }

        return sphericalExcess;
    }