Adding greatCircleDistance and haversine formula.

radixmath/constants.hh

-/*
- * File:   constants.hh
- * Author: Jordan P. Lefebvre
- *
- * Created on August 13, 2012, 8:26 PM
- */
-
-#ifndef RADIX_RADIXMATH_CONSTANTS_H
-#define	RADIX_RADIXMATH_CONSTANTS_H
-#include <cmath>
-#include <limits>
-namespace radix
-{
-    /**
-     * Type def double precision
-     */
-    typedef unsigned short index_t;
-    typedef double Real;
-    const Real kEpsilon = 1.0E-12;
-    const Real lookAheadEpsilon = 1.0E-9;
-    const Real halfKEpsilon = 0.5 * kEpsilon;
-    const Real kHugeValue = 1.0E10;
-    const Real maxReal = std::numeric_limits<Real>::max();
-    const Real minReal = std::numeric_limits<Real>::min();
-    const Real PI = 2 * std::acos( Real( 0 ) );//3.1415926535897932384;
-    const Real TWO_PI = PI * 2;//6.2831853071795864769;
-    const Real PI_ON_180 = PI / 180;//0.0174532925199432957;
-    const Real invPI = 1 / PI;//0.3183098861837906715;
-    const Real invTWO_PI = 1 / TWO_PI;//0.1591549430918953358;
-    const Real infinity = std::numeric_limits<Real>::infinity();
-    const Real sqrtHALF = std::sqrt(Real(0.5));
-    const Real sqrtTWO = std::sqrt(Real(2.0));
-    const Real sqrtTHREE = std::sqrt(Real(3.0));
-    /**
-     * Type def unsigned long to Identifier
-     */
-    typedef unsigned long Identifier;
-
-    /**
-     * @function clamp
-     * @param x 1d point in space
-     * @param min
-     * @param max
-     * @return
-     */
-    inline Real clamp(const Real x, const Real min, const Real max) {
-        return (x < min ? min : (x > max ? max : x));
-    }
-    inline bool isWithin(Real x, Real y, Real within=kEpsilon){
-        return std::abs(x-y) <= within;
-    }
-    /**
-     * Determine if the two values are within kepsilon of one-another
-     * @param x
-     * @param y
-     * @return true, iff abs(x-y) <= kEpsilon
-     */
-    inline bool isWithinKEpsilon(Real x, Real y){
-        return isWithin(x,y,kEpsilon);
-    }
-
-} // namespace radix
-#endif	/* RADIX_RADIXMATH_CONSTANTS_H */
-
+/*
+ * File:   constants.hh
+ * Author: Jordan P. Lefebvre
+ *
+ * Created on August 13, 2012, 8:26 PM
+ */
+
+#ifndef RADIX_RADIXMATH_CONSTANTS_H
+#define	RADIX_RADIXMATH_CONSTANTS_H
+#include <cmath>
+#include <limits>
+namespace radix
+{
+    /**
+     * Type def double precision
+     */
+    typedef unsigned short index_t;
+    typedef double Real;
+    const Real EARTH_RADIUS_MEAN = 6371000.0;
+    const Real wgs84_minor = 6356752.314245;
+    const Real kEpsilon = 1.0E-12;
+    const Real lookAheadEpsilon = 1.0E-9;
+    const Real halfKEpsilon = 0.5 * kEpsilon;
+    const Real kHugeValue = 1.0E10;
+    const Real maxReal = std::numeric_limits<Real>::max();
+    const Real minReal = std::numeric_limits<Real>::min();
+    const Real PI = 2 * std::acos( Real( 0 ) );//3.1415926535897932384;
+    const Real TWO_PI = PI * 2;//6.2831853071795864769;
+    const Real PI_ON_180 = PI / 180;//0.0174532925199432957;
+    const Real invPI = 1 / PI;//0.3183098861837906715;
+    const Real invTWO_PI = 1 / TWO_PI;//0.1591549430918953358;
+    const Real infinity = std::numeric_limits<Real>::infinity();
+    const Real sqrtHALF = std::sqrt(Real(0.5));
+    const Real sqrtTWO = std::sqrt(Real(2.0));
+    const Real sqrtTHREE = std::sqrt(Real(3.0));
+    /**
+     * Type def unsigned long to Identifier
+     */
+    typedef unsigned long Identifier;
+
+    /**
+     * @function clamp
+     * @param x 1d point in space
+     * @param min
+     * @param max
+     * @return
+     */
+    inline Real clamp(const Real x, const Real min, const Real max) {
+        return (x < min ? min : (x > max ? max : x));
+    }
+    inline bool isWithin(Real x, Real y, Real within=kEpsilon){
+        return std::abs(x-y) <= within;
+    }
+    /**
+     * Determine if the two values are within kepsilon of one-another
+     * @param x
+     * @param y
+     * @return true, iff abs(x-y) <= kEpsilon
+     */
+    inline bool isWithinKEpsilon(Real x, Real y){
+        return isWithin(x,y,kEpsilon);
+    }
+
+} // namespace radix
+#endif	/* RADIX_RADIXMATH_CONSTANTS_H */
+

radixmath/util.cc

@@ -89,4 +89,27 @@ Real lineIntersect(const Point3D &p, const Vector3D &v
     t *= denominator < 0 ? -1 : 1;
     return t;
 }
+
+Real haversine(Real lat1, Real lon1, Real lat2, Real lon2)
+{
+    // Haversine formula
+    Real dlat = toRadians(lat2 - lat1);
+    Real dlon = toRadians(lon2 - lon1);
+    Real haversine_dlat = std::sin(dlat / 2.0);
+    haversine_dlat *= haversine_dlat;
+    Real haversine_dlon = std::sin(dlon / 2.0);
+    haversine_dlon *= haversine_dlon;
+    Real delta_sigma = haversine_dlat
+            + std::cos(toRadians(lat1))
+            * std::cos(toRadians(lat2))
+            * haversine_dlon;
+    delta_sigma = 2.0 * std::asin(std::sqrt(delta_sigma));
+    return delta_sigma;
+}
+
+Real greatCircleDistance(Real lat1, Real lon1, Real lat2, Real lon2, Real radius)
+{
+    return radius * haversine(lat1, lon2, lat2, lon2);
+}
+
 } // namespace radix

radixmath/util.hh

@@ -15,6 +15,28 @@
 
 namespace radix
 {
+
+/**
+ * @brief haversine The Haversine formula for calculating great-circle
+ * @param lat1
+ * @param lon1
+ * @param lat2
+ * @param lon2
+ * @return angle(radians) between two points
+ */
+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 radius
+ * @return
+ */
+Real greatCircleDistance(Real lat1, Real lon1, Real lat2, Real lon2, Real radius);
+
 std::vector<Real> getEqualAreaRadii( const std::vector<Real>& radii, const std::vector<unsigned short>& subrings );
 
 inline Real toDegrees( Real radians )