Adding greatCircleVolume function to approximate volume of lat lon voxel.

radixmath/tests/CMakeLists.txt

 INCLUDE(GoogleTest)
 
 ADD_GOOGLE_TEST(tstMatrix.cc NP 1)
+ADD_GOOGLE_TEST(tstUtil.cc NP 1)

radixmath/tests/tstUtil.cc

+#include <iostream>
+
+#include "gtest/gtest.h"
+#include "radixmath/util.hh"
+
+using namespace radix;
+TEST(radix, greatCircleDistance){
+    { // Quadrant 2 north west hemisphere
+        double delta = 0.05;
+        double lat1=37.1665;
+        double lon1=-116.15;
+        double lat2=lat1+delta;
+        double lon2=lon1+delta;
+
+        //
+        // south
+        double distance = greatCircleDistance(lat1
+                                              , lon1
+                                              , lat1
+                                              , lon2
+                                              , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(4430.4688, distance);
+
+        // north
+        distance = greatCircleDistance(lat2
+                                       , lon1
+                                       , lat2
+                                       , lon2
+                                       , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(4427.5361, distance);
+
+        // west
+        double west = greatCircleDistance(lat1
+                                          , lon1
+                                          , lat2
+                                          , lon1
+                                          , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(5559.7461, west);
+
+        // east (should equal west)
+        distance = greatCircleDistance(lat1
+                                       , lon2
+                                       , lat2
+                                       , lon2
+                                       , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(west, distance);
+    }
+    { // Quadrant 3 south west hemisphere
+        double delta = -0.05;
+        double lat1=-37.1665;
+        double lon1=-116.15;
+        double lat2=lat1+delta;
+        double lon2=lon1+delta;
+
+        //
+        // south
+        double distance = greatCircleDistance(lat1
+                                              , lon1
+                                              , lat1
+                                              , lon2
+                                              , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(4430.4688, distance);
+
+        // north
+        distance = greatCircleDistance(lat2
+                                       , lon1
+                                       , lat2
+                                       , lon2
+                                       , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(4427.5361, distance);
+
+        // west
+        double west = greatCircleDistance(lat1
+                                          , lon1
+                                          , lat2
+                                          , lon1
+                                          , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(5559.7461, west);
+
+        // east (should equal west)
+        distance = greatCircleDistance(lat1
+                                       , lon2
+                                       , lat2
+                                       , lon2
+                                       , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(west, distance);
+    }
+    { // Quadrant 1 north east hemisphere
+        double delta = 0.05;
+        double lat1=37.1665;
+        double lon1=116.15;
+        double lat2=lat1+delta;
+        double lon2=lon1+delta;
+
+        //
+        // south
+        double distance = greatCircleDistance(lat1
+                                              , lon1
+                                              , lat1
+                                              , lon2
+                                              , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(4430.4688, distance);
+
+        // north
+        distance = greatCircleDistance(lat2
+                                       , lon1
+                                       , lat2
+                                       , lon2
+                                       , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(4427.5361, distance);
+
+        // west
+        double west = greatCircleDistance(lat1
+                                          , lon1
+                                          , lat2
+                                          , lon1
+                                          , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(5559.7461, west);
+
+        // east (should equal west)
+        distance = greatCircleDistance(lat1
+                                       , lon2
+                                       , lat2
+                                       , lon2
+                                       , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(west, distance);
+    }
+    { // Quadrant 4 south east hemisphere
+        double delta = -0.05;
+        double lat1=-37.1665;
+        double lon1=116.15;
+        double lat2=lat1+delta;
+        double lon2=lon1+delta;
+
+        //
+        // south
+        double distance = greatCircleDistance(lat1
+                                              , lon1
+                                              , lat1
+                                              , lon2
+                                              , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(4430.4688, distance);
+
+        // north
+        distance = greatCircleDistance(lat2
+                                       , lon1
+                                       , lat2
+                                       , lon2
+                                       , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(4427.5361, distance);
+
+        // west
+        double west = greatCircleDistance(lat1
+                                          , lon1
+                                          , lat2
+                                          , lon1
+                                          , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(5559.7461, west);
+
+        // east (should equal west)
+        distance = greatCircleDistance(lat1
+                                       , lon2
+                                       , lat2
+                                       , lon2
+                                       , EARTH_RADIUS_MEAN);
+        EXPECT_FLOAT_EQ(west, distance);
+    }
+}
+
+TEST(radix, greatCircleVolume)
+{
+    double delta = 0.05;
+    double lat1=37.1665;
+    double lon1=-116.15;
+    double lat2=lat1+delta;
+    double lon2=lon1+delta;
+
+    double volume = greatCircleVolume(lat1
+                                      , lon1
+                                      , lat2
+                                      , lon2
+                                      , EARTH_RADIUS_MEAN
+                                      , EARTH_RADIUS_MEAN+1);
+    EXPECT_FLOAT_EQ(24624134, volume);
+}

radixmath/util.cc

@@ -112,4 +112,38 @@ Real greatCircleDistance(Real lat1, Real lon1, Real lat2, Real lon2, Real radius
     return radius * haversine(lat1, lon1, lat2, lon2);
 }
 
+Real greatCircleVolume(Real lat1, Real lon1, Real lat2, Real lon2, Real r1, Real r2)
+{
+    Real r11 = r1, r22 = r2;
+    r1 = std::min(r11, r22);
+    r2 = std::max(r11, r22);
+
+    double northa = radix::greatCircleDistance(lat1, lon1
+                                            , lat1, lon2
+                                            , r1);
+    double southa = radix::greatCircleDistance(lat2, lon1
+                                               , lat2, lon2
+                                               , r1);
+    double westa = radix::greatCircleDistance(lat1, lon1
+                                             , lat2, lon1
+                                             , r1);
+
+    double northb = radix::greatCircleDistance(lat1, lon1
+                                            , lat1, lon2
+                                            , r2);
+    double southb = radix::greatCircleDistance(lat2, lon1
+                                               , lat2, lon2
+                                               , r2);
+    double westb = radix::greatCircleDistance(lat1, lon1
+                                             , lat2, lon1
+                                             , r2);
+    //
+    // approximate with a square
+    double north =(northa + northb + southa + southb)/4.0;
+    double west = (westa + westb)/2.0;
+    double height = r2 - r1;
+
+    return (north*west*height);
+}
+
 } // namespace radix

radixmath/util.hh

@@ -28,15 +28,27 @@ Real haversine(Real lat1, Real lon1, Real lat2, Real lon2);
 
 /**
  * @brief greatCircleDistance Calculates the shortest distance between two points on the surface of a sphere
- * @param lat1
- * @param lon1
- * @param lat2
- * @param lon2
+ * @param lat1 Real point1 latitude
+ * @param lon1 Real point1 longitude
+ * @param lat2 Real point2 latitude
+ * @param lon2 Real point2 longitude
  * @param radius
  * @return
  */
 Real greatCircleDistance(Real lat1, Real lon1, Real lat2, Real lon2, Real radius);
 
+/**
+ * @brief greatCircleVolume Calculates the volume giving latitude, longitude and radius bounds.
+ * @param lat1 Real min latitude
+ * @param lon1 Real min longitude
+ * @param lat2 Real max latitude
+ * @param lon2 Real max longitude
+ * @param r1 Real min radius
+ * @param r2 Real max radius
+ * @return Volume in cubic meters
+ */
+Real greatCircleVolume(Real lat1, Real lon1, Real lat2, Real lon2, Real r1, Real r2);
+
 std::vector<Real> getEqualAreaRadii( const std::vector<Real>& radii, const std::vector<unsigned short>& subrings );
 
 inline Real toDegrees( Real radians )