/*  File:    haversine.pl
    Author:  Jan Wielemaker,,,
    Created: Apr 28 2011
    Purpose: Compute distance between points on the globe
*/

:- module(haversine,
	  [ haversine/3			% +Point1, +Point2, -Distance
	  ]).

%%     haversine(Point1, Point2, Distance:meters) is det.
%
%      Implements the Haversine formula to compute the distance between
%      Point1 and Point2 in meters.
%
%      @see http://en.wikipedia.org/wiki/Haversine_formula

haversine(Point1, Point2, Distance) :-
       arc_in_radians(Point1, Point2, Arc),
       earth_radius(Meters),
       Distance is Meters*Arc.

earth_radius(6372797.560856).
deg_to_rad(0.017453292519943295769236907684886).

arc_in_radians(point(Lat1, Lon1), point(Lat2, Lon2), Arc) :-
       deg_to_rad(DegToRad),
       LatArc is (Lat1-Lat2) * DegToRad,
       LonArc is (Lon1-Lon2) * DegToRad,
       LatH is sin(LatArc*0.5)**2,
       LonH is sin(LonArc*0.5)**2,
       Tmp is cos(Lat1*DegToRad) * cos(Lat2*DegToRad),
       Arc is 2.0 * asin(sqrt(LatH+Tmp*LonH)).