From 0518e680c3bddcae46108ca1d523da74ab37fb77 Mon Sep 17 00:00:00 2001
From: Klaus Basan <klaus.basan@klausbasan.de>
Date: Tue, 10 Dec 2013 19:32:02 +0000
Subject: [PATCH] add greatCircleDistance

refs #81
---
 src/blackmisc/coordinategeodetic.cpp | 26 ++++++++++++++++++++++++++
 src/blackmisc/coordinategeodetic.h   | 10 +++++++++-
 2 files changed, 35 insertions(+), 1 deletion(-)

diff --git a/src/blackmisc/coordinategeodetic.cpp b/src/blackmisc/coordinategeodetic.cpp
index d72217b2a..71471e6ae 100644
--- a/src/blackmisc/coordinategeodetic.cpp
+++ b/src/blackmisc/coordinategeodetic.cpp
@@ -68,5 +68,31 @@ namespace BlackMisc
         }
 
 
+        /*
+         * Great circle distance
+         */
+        PhysicalQuantities::CLength greatCircleDistance(const ICoordinateGeodetic &coordinate1, const ICoordinateGeodetic &coordinate2)
+        {
+            // same coordinate results in 0 distance
+            if (coordinate1.latitude() == coordinate2.latitude() && coordinate1.longitude() == coordinate2.longitude())
+            {
+                return CLength(0, CLengthUnit::NM());
+            }
+
+            // first, prelimary distance calculation
+            // http://www.geodatasource.com/developers/c
+            double dist;
+            double lon1rad = coordinate1.longitude().value(CAngleUnit::rad());
+            double lon2rad = coordinate2.longitude().value(CAngleUnit::rad());
+            double lat1rad = coordinate1.latitude().value(CAngleUnit::rad());
+            double lat2rad = coordinate2.latitude().value(CAngleUnit::rad());
+            double theta = lon1rad - lon2rad;
+            dist = qSin(lat1rad) * qSin(lat2rad) + qCos(lat1rad) * qCos(lat2rad) * cos(theta);
+            dist = qAcos(dist);
+            dist = CMath::rad2deg(dist);
+            dist = dist * 60; // dist in NM
+            return CLength(qAbs(dist), CLengthUnit::NM());
+        }
+
     } // namespace
 } // namespace
diff --git a/src/blackmisc/coordinategeodetic.h b/src/blackmisc/coordinategeodetic.h
index 839d72e15..a2ac465c5 100644
--- a/src/blackmisc/coordinategeodetic.h
+++ b/src/blackmisc/coordinategeodetic.h
@@ -55,10 +55,18 @@ namespace BlackMisc
 
         };
 
+        /*!
+         * \brief Great circle distance between points
+         * \param coordinate1
+         * \param coordinate2
+         * \return
+         */
+        BlackMisc::PhysicalQuantities::CLength greatCircleDistance(const ICoordinateGeodetic &coordinate1, const ICoordinateGeodetic &coordinate2);
+
         /*!
          * \brief Geodetic coordinate
          */
-        class CCoordinateGeodetic : public CValueObject
+        class CCoordinateGeodetic : public CValueObject, public ICoordinateGeodetic
         {
         private:
             BlackMisc::Geo::CLatitude m_latitude; //!< Latitude