Add distance calculation based on haversine

src/grl-coordinates.adb

@@ -1,3 +1,5 @@
+with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;
+
 package body Grl.Coordinates is
 
    function To_DD (DM_Base : Coordinate) return Coordinate is
@@ -37,4 +39,39 @@ package body Grl.Coordinates is
       return DD_Target;
    end To_DD;
 
+   function Plane_Distance (Origin, Destination : Coordinate) return H is
+      Heuristic      : H;
+      A              : H;
+      DD_Origin      : constant Coordinate := To_DD (Origin);
+      DD_Destination : constant Coordinate := To_DD (Destination);
+      Earth_Radius   : constant H          := 6_371_000.0; --  M
+      Delta_Phi      : constant H          := --  Latitude difference
+        H (DD_Destination.DD.Latitude - DD_Origin.DD.Latitude);
+      Delta_Lambda   : constant H          := -- Longitude difference
+        H (DD_Destination.DD.Longitude - DD_Origin.DD.Longitude);
+
+      --  @name Hav
+      --  @return Float
+      --  @parameter Theta:Float
+      --  @variable Fraction : Float
+      --  @description Compute Haversine function.
+      function Hav (Theta : H) return Float is
+         Fraction : Float;
+      begin
+         Fraction := Float (Theta) / 2.0;
+         return Sin (Fraction)**2;
+      end Hav;
+   begin
+
+      A :=
+        H (Hav (Delta_Phi) +
+           Cos (Float (DD_Origin.DD.Latitude)) *
+             Cos (Float (DD_Destination.DD.Latitude)) * Hav (Delta_Lambda));
+
+      --  Haversine formula
+      Heuristic := Earth_Radius * H (2.0 * Arcsin (Sqrt (Float (A))));
+
+      return Heuristic;
+   end Plane_Distance;
+
 end Grl.Coordinates;

src/grl-coordinates.ads

@@ -118,4 +118,12 @@ is
       (DM_Base.Form = DD) or else (DM_Base.Form = DMS)
       or else (DM_Base.Form = DDM);
 
+   --  @name Plane_Distance
+   --  @return H
+   --  @parameter Origin:Coordinate
+   --  @parameter Destination:Coordinate
+   --  @variable Heuristic : H
+   --  @description Calculate heuristic between coordinates.
+   function Plane_Distance (Origin, Destination : Coordinate) return H;
+
 end Grl.Coordinates;