Jump to content

Haversine formula

From Rosetta Code
Task
Haversine formula
You are encouraged to solve this task according to the task description, using any language you may know.
This page uses content from Wikipedia. The original article was at Haversine formula. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)

The haversine formula is an equation important in navigation,

giving great-circle distances between two points on a sphere from their longitudes and latitudes. It is a special case of a more general formula in spherical trigonometry, the law of haversines, relating the sides and angles of spherical "triangles".

Task: Implement a great-circle distance function, or use a library function, to show the great-circle distance between Nashville International Airport (BNA) in Nashville, TN, USA: N 36°7.2', W 86°40.2' (36.12, -86.67) and Los Angeles International Airport (LAX) in Los Angeles, CA, USA: N 33°56.4', W 118°24.0' (33.94, -118.40). 

User Kaimbridge clarified on the Talk page:

 -- 6371.0 km is the authalic radius based on/extracted from surface area;
 -- 6372.8 km is an approximation of the radius of the average circumference
    (i.e., the average great-elliptic or great-circle radius), where the
     boundaries are the meridian (6367.45 km) and the equator (6378.14 km).

Using either of these values results, of course, in differing distances:

 6371.0 km -> 2886.44444283798329974715782394574671655 km;
 6372.8 km -> 2887.25995060711033944886005029688505340 km;
 (results extended for accuracy check:  Given that the radii are only
  approximations anyways, .01' ≈ 1.0621333 km and .001" ≈ .00177 km,
  practical precision required is certainly no greater than about
  .0000001——i.e., .1 mm!)

As distances are segments of great circles/circumferences, it is
recommended that the latter value (r = 6372.8 km) be used (which
most of the given solutions have already adopted, anyways).

Ada

<lang ada>with Ada.Text_IO; use Ada.Text_IO; with Ada.Long_Float_Text_IO; use Ada.Long_Float_Text_IO; with Ada.Numerics.Generic_Elementary_Functions;

procedure Haversine_Formula is 

  package Math is new Ada.Numerics.Generic_Elementary_Functions (Long_Float); use Math;

  -- Compute great circle distance, given latitude and longitude of two points, in radians
  function Great_Circle_Distance (lat1, long1, lat2, long2 : Long_Float) return Long_Float is
     Earth_Radius : constant := 6371.0; -- in kilometers
     a : Long_Float := Sin (0.5 * (lat2 - lat1));
     b : Long_Float := Sin (0.5 * (long2 - long1));
  begin
     return 2.0 * Earth_Radius * ArcSin (Sqrt (a * a + Cos (lat1) * Cos (lat2) * b * b));
  end Great_Circle_Distance;

  -- convert degrees, minutes and seconds to radians
  function DMS_To_Radians (Deg, Min, Sec : Long_Float := 0.0) return Long_Float is
     Pi_Over_180 : constant := 0.017453_292519_943295_769236_907684_886127;
  begin
     return (Deg + Min/60.0 + Sec/3600.0) * Pi_Over_180;
  end DMS_To_Radians;

begin 

  Put_Line("Distance in kilometers between BNA and LAX");
  Put (Great_Circle_Distance (
        DMS_To_Radians (36.0, 7.2), DMS_To_Radians (86.0, 40.2),       -- Nashville International Airport (BNA)
        DMS_To_Radians (33.0, 56.4), DMS_To_Radians (118.0, 24.0)),    -- Los Angeles International Airport (LAX)
     Aft=>3, Exp=>0);

end Haversine_Formula;</lang>

ALGOL 68

Translation of: C
Works with: ALGOL 68 version Revision 1.
Works with: ALGOL 68G version Any - tested with release algol68g-2.3.5.

File: Haversine_formula.a68<lang algol68>#!/usr/local/bin/a68g --script #

REAL r = 20 000/pi + 6.6 # km #, 

    to rad = pi/180;

PROC dist = (REAL th1 deg, ph1 deg, th2 deg, ph2 deg)REAL: ( 

       REAL ph1 = (ph1 deg - ph2 deg) * to rad,
            th1 = th1 deg * to rad, th2 = th2 deg * to rad,

            dz = sin(th1) - sin(th2),
            dx = cos(ph1) * cos(th1) - cos(th2),
            dy = sin(ph1) * cos(th1);
       arc sin(sqrt(dx * dx + dy * dy + dz * dz) / 2) * 2 * r

);

main: ( 

       REAL d = dist(36.12, -86.67, 33.94, -118.4);
       # Americans don't know kilometers #
       printf(($"dist: "g(0,1)" km ("g(0,1)" mi.)"l$, d, d / 1.609344))

)</lang>

Output:
dist: 2887.3 km (1794.1 mi.)

ATS

<lang ATS>

  1. include

"share/atspre_staload.hats"

staload "libc/SATS/math.sats" staload _ = "libc/DATS/math.dats" staload "libc/SATS/stdio.sats" staload "libc/SATS/stdlib.sats"

  1. define R 6372.8
  2. define TO_RAD (3.1415926536 / 180)

typedef d = double

fun dist ( 

 th1: d, ph1: d, th2: d, ph2: d

) : d = let 

 val ph1 = ph1 - ph2
 val ph1 = TO_RAD * ph1
 val th1 = TO_RAD * th1
 val th2 = TO_RAD * th2
 val dz = sin(th1) - sin(th2)
 val dx = cos(ph1) * cos(th1) - cos(th2)
 val dy = sin(ph1) * cos(th1)

in 

 asin(sqrt(dx*dx + dy*dy + dz*dz)/2)*2*R

end // end of [dist]

implement main0((*void*)) = let 

 val d = dist(36.12, ~86.67, 33.94, ~118.4);
 /* Americans don't know kilometers */

in 

 $extfcall(void, "printf", "dist: %.1f km (%.1f mi.)\n", d, d / 1.609344)

end // end of [main0] </lang>

Output:
dist: 2887.3 km (1794.1 mi.)

AutoHotkey

<lang AutoHotkey>MsgBox, % GreatCircleDist(36.12, 33.94, -86.67, -118.40, 6372.8, "km")

GreatCircleDist(La1, La2, Lo1, Lo2, R, U) { return, 2 * R * ASin(Sqrt(Hs(Rad(La2 - La1)) + Cos(Rad(La1)) * Cos(Rad(La2)) * Hs(Rad(Lo2 - Lo1)))) A_Space U }

Hs(n) { return, (1 - Cos(n)) / 2 }

Rad(Deg) { return, Deg * 4 * ATan(1) / 180 }</lang>

Output:
2887.259951 km

AWK

<lang AWK>

  1. syntax: GAWK -f HAVERSINE_FORMULA.AWK
  2. converted from Python

BEGIN { 

   distance(36.12,-86.67,33.94,-118.40) # BNA to LAX
   exit(0)

} function distance(lat1,lon1,lat2,lon2, a,c,dlat,dlon) { 

   dlat = radians(lat2-lat1)
   dlon = radians(lon2-lon1)
   lat1 = radians(lat1)
   lat2 = radians(lat2)
   a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
   c = 2 * atan2(sqrt(a),sqrt(1-a))
   printf("distance: %.4f km\n",6372.8 * c)

} function radians(degree) { # degrees to radians 

   return degree * (3.1415926 / 180.)

} </lang>

Output:
distance: 2887.2599 km

BBC BASIC

Works with: BBC BASIC for Windows

Uses BBC BASIC's MOD(array()) function which calculates the square-root of the sum of the squares of the elements of an array. <lang bbcbasic> PRINT "Distance = " ; FNhaversine(36.12, -86.67, 33.94, -118.4) " km" 

     END
     
     DEF FNhaversine(n1, e1, n2, e2)
     LOCAL d() : DIM d(2)
     d() = COSRAD(e1-e2) * COSRAD(n1) - COSRAD(n2), \
     \     SINRAD(e1-e2) * COSRAD(n1), \
     \     SINRAD(n1) - SINRAD(n2)
     = ASN(MOD(d()) / 2) * 6372.8 * 2</lang>
Output:
Distance = 2887.25995 km

C

<lang c>#include <stdio.h>

  1. include <stdlib.h>
  2. include <math.h>
  1. define R 6371
  2. define TO_RAD (3.1415926536 / 180)

double dist(double th1, double ph1, double th2, double ph2) { double dx, dy, dz; ph1 -= ph2; ph1 *= TO_RAD, th1 *= TO_RAD, th2 *= TO_RAD;

dz = sin(th1) - sin(th2); dx = cos(ph1) * cos(th1) - cos(th2); dy = sin(ph1) * cos(th1); return asin(sqrt(dx * dx + dy * dy + dz * dz) / 2) * 2 * R; }

int main() { double d = dist(36.12, -86.67, 33.94, -118.4); /* Americans don't know kilometers */ printf("dist: %.1f km (%.1f mi.)\n", d, d / 1.609344);

return 0; }</lang>

C#

Translation of: Groovy

<lang csharp>public static class Haversine { 

 public static double calculate(double lat1, double lon1, double lat2, double lon2) {
   var R = 6372.8; // In kilometers
   var dLat = toRadians(lat2 - lat1);
   var dLon = toRadians(lon2 - lon1);
   lat1 = toRadians(lat1);
   lat2 = toRadians(lat2);
  
   var a = Math.Sin(dLat / 2) * Math.Sin(dLat / 2) + Math.Sin(dLon / 2) * Math.Sin(dLon / 2) * Math.Cos(lat1) * Math.Cos(lat2);
   var c = 2 * Math.Asin(Math.Sqrt(a));
   return R * 2 * Math.Asin(Math.Sqrt(a));
 }
 
 public static double toRadians(double angle) {
   return Math.PI * angle / 180.0;
 }

}

void Main() { 

 Console.WriteLine(String.Format("The distance between coordinates {0},{1} and {2},{3} is: {4}", 36.12, -86.67, 33.94, -118.40, Haversine.calculate(36.12, -86.67, 33.94, -118.40)));

}

// Returns: The distance between coordinates 36.12,-86.67 and 33.94,-118.4 is: 2887.25995060711 </lang>

clojure

Translation of: Java

<lang clojure> (defn haversine 

 [{lon1 :longitude lat1 :latitude} {lon2 :longitude lat2 :latitude}]
 (let [R 6372.8 ; kilometers
       dlat (Math/toRadians (- lat2 lat1))
       dlon (Math/toRadians (- lon2 lon1))
       lat1 (Math/toRadians lat1)
       lat2 (Math/toRadians lat2)
       a (+ (* (Math/sin (/ dlat 2)) (Math/sin (/ dlat 2))) (* (Math/sin (/ dlon 2)) (Math/sin (/ dlon 2)) (Math/cos lat1) (Math/cos lat2)))]
   (* R 2 (Math/asin (Math/sqrt a)))))

(haversine {:latitude 36.12 :longitude -86.67} {:latitude 33.94 :longitude -118.40})

=> 2887.2599506071106

</lang>

CoffeeScript

Translation of: JavaScript

<lang coffee>haversine = (args...) -> 

 R = 6372.8; # km
 radians = args.map (deg) -> deg/180.0 * Math.PI
 lat1 = radians[0]; lon1 = radians[1]; lat2 = radians[2]; lon2 = radians[3]
 dLat = lat2 - lat1
 dLon = lon2 - lon1
 a = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.sin(dLon / 2) * Math.sin(dLon / 2) * Math.cos(lat1) * Math.cos(lat2)
 R * 2 * Math.asin(Math.sqrt(a))

console.log haversine(36.12, -86.67, 33.94, -118.40)</lang>

Output:
2887.2599506071124

Common Lisp

<lang lisp>(defparameter *earth-radius* 6372.8)

(defparameter *rad-conv* (/ pi 180))

(defun deg->rad (x) 

 (* x *rad-conv*))

(defun haversine (x) 

 (expt (sin (/ x 2)) 2))

(defun dist-rad (lat1 lng1 lat2 lng2) 

 (let* ((hlat (haversine (- lat2 lat1)))
        (hlng (haversine (- lng2 lng1)))
        (root (sqrt (+ hlat (* (cos lat1) (cos lat2) hlng)))))
   (* 2 *earth-radius* (asin root))))

(defun dist-deg (lat1 lng1 lat2 lng2) 

 (dist-rad (deg->rad lat1)
           (deg->rad lng1)
           (deg->rad lat2)
           (deg->rad lng2)))</lang>
Output:
CL-USER> (format t "~%The distance between BNA and LAX is about ~$ km.~%" 
		 (dist-deg 36.12 -86.67 33.94 -118.40))

The distance between BNA and LAX is about 2887.26 km.

D

<lang d>import std.stdio, std.math;

real haversineDistance(in real dth1, in real dph1, 

                      in real dth2, in real dph2)

pure nothrow @nogc { 

   enum real R = 6371;
   enum real TO_RAD = PI / 180;

   alias imr = immutable real;
   imr ph1d = dph1 - dph2;
   imr ph1 = ph1d * TO_RAD;
   imr th1 = dth1 * TO_RAD;
   imr th2 = dth2 * TO_RAD;

   imr dz = th1.sin - th2.sin;
   imr dx = ph1.cos * th1.cos - th2.cos;
   imr dy = ph1.sin * th1.cos;
   return asin(sqrt(dx ^^ 2 + dy ^^ 2 + dz ^^ 2) / 2) * 2 * R;

}

void main() { 

   writefln("Haversine distance: %.1f km",
            haversineDistance(36.12, -86.67, 33.94, -118.4));

}</lang>

Output:
Haversine distance: 2887.3 km

Alternative Version

An alternate direct implementation of the haversine formula as shown at wikipedia. The same length, but perhaps a little more clear about what is being done.

<lang d>import std.stdio, std.math;

real toRad(in real degrees) pure nothrow @safe @nogc { 

   return degrees * PI / 180;

}

real haversin(in real theta) pure nothrow @safe @nogc { 

   return (1 - theta.cos) / 2;

}

real greatCircleDistance(in real lat1, in real lng1, 

                        in real lat2, in real lng2,
                        in real radius)

pure nothrow @safe @nogc { 

   immutable h = haversin(lat2.toRad - lat1.toRad) +
                 lat1.toRad.cos * lat2.toRad.cos *
                 haversin(lng2.toRad - lng1.toRad);
   return 2 * radius * h.sqrt.asin;

}

void main() { 

   enum real earthRadius = 6372.8L; // Average earth radius.

   writefln("Great circle distance: %.1f km",
            greatCircleDistance(36.12, -86.67, 33.94, -118.4,
                                earthRadius));

}</lang>

Output:
Great circle distance: 2887.3 km


Delphi

<lang delphi> program HaversineDemo; uses Math;

function HaversineDist(th1, ph1, th2, ph2:double):double; const diameter = 2 * 6372.8; var dx, dy, dz:double; begin 

 ph1    := degtorad(ph1 - ph2);
 th1    := degtorad(th1);
 th2    := degtorad(th2);

 dz     := sin(th1) - sin(th2);
 dx     := cos(ph1) * cos(th1) - cos(th2);
 dy     := sin(ph1) * cos(th1);
 Result := arcsin(sqrt(sqr(dx) + sqr(dy) + sqr(dz)) / 2) * diameter;

end;

begin 

 Writeln('Haversine distance: ', HaversineDist(36.12, -86.67, 33.94, -118.4):7:2, ' km.');

end. </lang>

Output:
Haversine distance: 2887.26 km.

Erlang

<lang erlang>% Implementer by Arjun Sunel -module(haversine). -export([main/0]).

main() -> haversine(36.12, -86.67, 33.94, -118.40).

haversine(Lat1, Long1, Lat2, Long2) -> V = math:pi()/180, R = 6372.8, % In kilometers Diff_Lat = (Lat2 - Lat1)*V , Diff_Long = (Long2 - Long1)*V, NLat = Lat1*V, NLong = Lat2*V, A = math:sin(Diff_Lat/2) * math:sin(Diff_Lat/2) + math:sin(Diff_Long/2) * math:sin(Diff_Long/2) * math:cos(NLat) * math:cos(NLong), C = 2 * math:asin(math:sqrt(A)), R*C. </lang>

Output:
2887.2599506071106

ERRE

<lang ERRE>% Implemented by Claudio Larini

PROGRAM HAVERSINE_DEMO

!$DOUBLE

CONST DIAMETER=12745.6

FUNCTION DEG2RAD(X) 

   DEG2RAD=X*π/180

END FUNCTION

FUNCTION RAD2DEG(X) 

   RAD2DEG=X*180/π

END FUNCTION

PROCEDURE HAVERSINE_DIST(TH1,PH1,TH2,PH2->RES) 

   LOCAL DX,DY,DZ
   PH1=DEG2RAD(PH1-PH2)
   TH1=DEG2RAD(TH1)
   TH2=DEG2RAD(TH2)
   DZ=SIN(TH1)-SIN(TH2)
   DX=COS(PH1)*COS(TH1)-COS(TH2)
   DY=SIN(PH1)*COS(TH1)
   RES=ASN(SQR(DX^2+DY^2+DZ^2)/2)*DIAMETER

END PROCEDURE

BEGIN 

   HAVERSINE_DIST(36.12,-86.67,33.94,-118.4->RES)
   PRINT("HAVERSINE DISTANCE: ";RES;" KM.")

END PROGRAM </lang> Using double-precision variables output is 2887.260209071741 km, while using single-precision variable output is 2887.261 Km.

Euler Math Toolbox

Euler has a package for spherical geometry, which is used in the following code. The distances are then computed with the average radius between the two positions. Overwriting the rearth function with the given value yields the known result. 

>load spherical
 Spherical functions for Euler. 
>TNA=[rad(36,7.2),-rad(86,40.2)];
>LAX=[rad(33,56.4),-rad(118,24)];
>esdist(TNA,LAX)->km
 2886.48817482
>type esdist
 function esdist (frompos: vector, topos: vector)
     r1=rearth(frompos[1]); 
     r2=rearth(topos[1]);
     xfrom=spoint(frompos)*r1; 
     xto=spoint(topos)*r2;
     delta=xto-xfrom;
     return asin(norm(delta)/(r1+r2))*(r1+r2);
 endfunction
>function overwrite rearth (x) := 6372.8*km$
>esdist(TNA,LAX)->km
 2887.25995061


F#

Translation of: Go

using units of measure

<lang fsharp>open System

[<Measure>] type deg [<Measure>] type rad [<Measure>] type km

let haversine (θ: float<rad>) = 0.5 * (1.0 - Math.Cos(θ/1.0<rad>))

let radPerDeg = (Math.PI / 180.0) * 1.0<rad/deg>

type pos(latitude: float<deg>, longitude: float<deg>) = 

   member this.φ = latitude * radPerDeg
   member this.ψ = longitude * radPerDeg

let rEarth = 6372.8<km>

let hsDist (p1: pos) (p2: pos) = 

   2.0 * rEarth *
       Math.Asin(Math.Sqrt(haversine(p2.φ - p1.φ)+
                   Math.Cos(p1.φ/1.0<rad>)*Math.Cos(p2.φ/1.0<rad>)*haversine(p2.ψ - p1.ψ)))

[<EntryPoint>] let main argv = 

   printfn "%A" (hsDist (pos(36.12<deg>, -86.67<deg>)) (pos(33.94<deg>, -118.40<deg>)))
   0</lang>
Output:
2887.259951


Factor

Translation of: J

<lang factor>USING: arrays kernel math math.constants math.functions math.vectors sequences ;

haversin ( x -- y ) cos 1 swap - 2 / ;
haversininv ( y -- x ) 2 * 1 swap - acos ;
haversineDist ( as bs -- d )

[ [ 180 / pi * ] map ] bi@ 

 [ [ swap - haversin ] 2map ]
 [ [ first cos ] bi@ * 1 swap 2array ]
 2bi

v. haversininv R_earth * ;</lang> <lang factor>( scratchpad ) { 36.12 -86.67 } { 33.94 -118.4 } haversineDist . 2887.259950607113</lang>

FBSL

Based on the Fortran and Groovy versions. <lang qbasic>#APPTYPE CONSOLE

PRINT "Distance = ", Haversine(36.12, -86.67, 33.94, -118.4), " km" PAUSE

FUNCTION Haversine(DegLat1 AS DOUBLE, DegLon1 AS DOUBLE, DegLat2 AS DOUBLE, DegLon2 AS DOUBLE) AS DOUBLE 

   CONST radius = 6372.8
   DIM dLat AS DOUBLE = D2R(DegLat2 - DegLat1)
   DIM dLon AS DOUBLE = D2R(DegLon2 - DegLon1)
   DIM lat1 AS DOUBLE = D2R(DegLat1)
   DIM lat2 AS DOUBLE = D2R(DegLat2)
   DIM a AS DOUBLE = SIN(dLat / 2) * SIN(dLat / 2) + SIN(dLon / 2) * SIN(dLon / 2) * COS(lat1) * COS(lat2)
   DIM c AS DOUBLE = 2 * ASIN(SQRT(a))
   RETURN radius * c

END FUNCTION </lang>

Output:
Distance = 2887.25995060711 km
Press any key to continue...

Forth

<lang forth>: s>f s>d d>f ;

deg>rad 174532925199433e-16 f* ;
difference f- deg>rad 2 s>f f/ fsin fdup f* ;
haversine ( lat1 lon1 lat2 lon2 -- haversine)
 frot difference                      ( lat1 lat2 dLon^2)
 frot frot fover fover                ( dLon^2 lat1 lat2 lat1 lat2)
 fswap difference                     ( dLon^2 lat1 lat2 dLat^2)
 fswap deg>rad fcos                   ( dLon^2 lat1 dLat^2 lat2)
 frot  deg>rad fcos f*                ( dLon^2 dLat2 lat1*lat2)
 frot  f* f+                          ( lat1*lat2*dLon^2+dLat^2)
 fsqrt fasin 127456 s>f f* 10 s>f f/  ( haversine)

36.12e -86.67e 33.94e -118.40e haversine cr f.</lang>

Output:
2887.25995060711

Fortran

<lang Fortran> program example implicit none real :: d

d = haversine(36.12,-86.67,33.94,-118.40) ! BNA to LAX print '(A,F9.4,A)', 'distance: ',d,' km' ! distance: 2887.2600 km

contains 

     function to_radian(degree) result(rad)
         ! degrees to radians
         real,intent(in) :: degree
         real, parameter :: deg_to_rad = atan(1.0)/45 ! exploit intrinsic atan to generate pi/180 runtime constant
         real :: rad

         rad = degree*deg_to_rad
     end function to_radian

     function haversine(deglat1,deglon1,deglat2,deglon2) result (dist)
         ! great circle distance -- adapted from Matlab 
         real,intent(in) :: deglat1,deglon1,deglat2,deglon2
         real :: a,c,dist,dlat,dlon,lat1,lat2
         real,parameter :: radius = 6372.8 

         dlat = to_radian(deglat2-deglat1)
         dlon = to_radian(deglon2-deglon1)
         lat1 = to_radian(deglat1)
         lat2 = to_radian(deglat2)
         a = (sin(dlat/2))**2 + cos(lat1)*cos(lat2)*(sin(dlon/2))**2
         c = 2*asin(sqrt(a))
         dist = radius*c
     end function haversine

end program example </lang>

Frink

<lang frink> haversine[theta] := (1-cos[theta])/2

dist[lat1, long1, lat2, long2] := 2 earthradius arcsin[sqrt[haversine[lat2-lat1] + cos[lat1] cos[lat2] haversine[long2-long1]]]

d = dist[36.12 deg, -86.67 deg, 33.94 deg, -118.40 deg] println[d-> "km"] </lang>

Note that physical constants like degrees, kilometers, and the average radius of the earth (as well as the polar and equatorial radii) are already known to Frink. Also note that units of measure are tracked throughout all calculations, and results can be displayed in a huge number of units of distance (miles, km, furlongs, chains, feet, statutemiles, etc.) by changing the final "km" to something like "miles".

However, Frink's library/sample program navigation.frink (included in larger distributions) contains a much higher-precision calculation that uses ellipsoidal (not spherical) calculations to determine the distance on earth's geoid with far greater accuracy:

<lang frink> use navigation.frink

d = earthDistance[36.12 deg North, 86.67 deg West, 33.94 deg North, 118.40 deg West] println[d-> "km"] </lang>

FunL

<lang funl>import math.*

def haversin( theta ) = (1 - cos( theta ))/2

def radians( deg ) = deg Pi/180

def haversine( (lat1, lon1), (lat2, lon2) ) = 

 R = 6372.8
 h = haversin( radians(lat2 - lat1) ) + cos( radians(lat1) ) cos( radians(lat2) ) haversin( radians(lon2 - lon1) )
 2R asin( sqrt(h) )

println( haversine((36.12, -86.67), (33.94, -118.40)) )</lang>

Output:
2887.259950607111

Go

<lang go>package main

import ( 

   "fmt"
   "math"

)

func haversine(θ float64) float64 { 

   return .5 * (1 - math.Cos(θ))

}

type pos struct { 

   φ float64 // latitude, radians
   ψ float64 // longitude, radians

}

func degPos(lat, lon float64) pos { 

   return pos{lat * math.Pi / 180, lon * math.Pi / 180}

}

const rEarth = 6372.8 // km

func hsDist(p1, p2 pos) float64 { 

   return 2 * rEarth * math.Asin(math.Sqrt(haversine(p2.φ-p1.φ)+
       math.Cos(p1.φ)*math.Cos(p2.φ)*haversine(p2.ψ-p1.ψ)))

}

func main() { 

   fmt.Println(hsDist(degPos(36.12, -86.67), degPos(33.94, -118.40)))

}</lang>

Output:
2887.2599506071097

Groovy

<lang Groovy>def haversine(lat1, lon1, lat2, lon2) { 

 def R = 6372.8
 // In kilometers
 def dLat = Math.toRadians(lat2 - lat1)
 def dLon = Math.toRadians(lon2 - lon1)
 lat1 = Math.toRadians(lat1)
 lat2 = Math.toRadians(lat2)

 def a = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.sin(dLon / 2) * Math.sin(dLon / 2) * Math.cos(lat1) * Math.cos(lat2)
 def c = 2 * Math.asin(Math.sqrt(a))
 R * c

}

haversine(36.12, -86.67, 33.94, -118.40)

> 2887.25995060711</lang>

Haskell

<lang Haskell>import Text.Printf

-- The haversine of an angle. hsin t = let u = sin (t/2) in u*u

-- The distance between two points, given by latitude and longtitude, on a -- circle. The points are specified in radians. distRad radius (lat1, lng1) (lat2, lng2) = 

 let hlat = hsin (lat2 - lat1)
     hlng = hsin (lng2 - lng1)
     root = sqrt (hlat + cos lat1 * cos lat2 * hlng)
 in 2 * radius * asin (min 1.0 root)

-- The distance between two points, given by latitude and longtitude, on a -- circle. The points are specified in degrees. distDeg radius p1 p2 = distRad radius (deg2rad p1) (deg2rad p2) 

 where deg2rad (t, u) = (d2r t, d2r u)
       d2r t = t * pi / 180

-- The approximate distance, in kilometers, between two points on Earth. -- The latitude and longtitude are assumed to be in degrees. earthDist = distDeg 6372.8

main = do 

 let bna = (36.12,  -86.67)
     lax = (33.94, -118.40)
     dst = earthDist bna lax :: Double
 printf "The distance between BNA and LAX is about %0.f km.\n" dst</lang>
Output:
The distance between BNA and LAX is about 2887 km.

Icon and Unicon

Translation of: C

<lang Icon>link printf

procedure main() #: Haversine formula 

  printf("BNA to LAX is %d km (%d miles)\n",
     d := gcdistance([36.12, -86.67],[33.94, -118.40]),d*3280/5280)  # with cute km2mi conversion

end

procedure gcdistance(a,b) a[2] -:= b[2] 

  every (x := a|b)[i := 1 to 2] := dtor(x[i])

dz := sin(a[1]) - sin(b[1]) dx := cos(a[2]) * cos(a[1]) - cos(b[1]) dy := sin(a[2]) * cos(a[1]) return asin(sqrt(dx * dx + dy * dy + dz * dz) / 2) * 2 * 6371 end</lang>

Library: Icon Programming Library

printf.icn provides formatting

Output:
BNA to LAX is 2886 km (1793 miles)

J

Solution: <lang j>require 'trig' haversin=: 0.5 * 1 - cos Rearth=: 6372.8 haversineDist=: Rearth * haversin^:_1@((1 , *&(cos@{.)) +/ .* [: haversin -)&rfd </lang> Note: J derives the inverse haversin ( haversin^:_1 ) from the definition of haversin.

Example Use: <lang j> 36.12 _86.67 haversineDist 33.94 _118.4 2887.26</lang>

Java

Translation of: Groovy

<lang java>public class Haversine { 

   public static final double R = 6372.8; // In kilometers
   public static double haversine(double lat1, double lon1, double lat2, double lon2) {
       double dLat = Math.toRadians(lat2 - lat1);
       double dLon = Math.toRadians(lon2 - lon1);
       lat1 = Math.toRadians(lat1);
       lat2 = Math.toRadians(lat2);

       double a = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.sin(dLon / 2) * Math.sin(dLon / 2) * Math.cos(lat1) * Math.cos(lat2);
       double c = 2 * Math.asin(Math.sqrt(a));
       return R * c;
   }
   public static void main(String[] args) {
       System.out.println(haversine(36.12, -86.67, 33.94, -118.40));
   }

}</lang>

Output:
2887.2599506071106

JavaScript

Translation of: Java

<lang javascript>function haversine() { 

      var radians = Array.prototype.map.call(arguments, function(deg) { return deg/180.0 * Math.PI; });
      var lat1 = radians[0], lon1 = radians[1], lat2 = radians[2], lon2 = radians[3];
      var R = 6372.8; // km
      var dLat = lat2 - lat1;
      var dLon = lon2 - lon1;
      var a = Math.sin(dLat / 2) * Math.sin(dLat /2) + Math.sin(dLon / 2) * Math.sin(dLon /2) * Math.cos(lat1) * Math.cos(lat2);
      var c = 2 * Math.asin(Math.sqrt(a));
      return R * c;

} console.log(haversine(36.12, -86.67, 33.94, -118.40));</lang>

Output:
2887.2599506071124

jq

<lang jq>def haversine(lat1;lon1; lat2;lon2): 

 def radians: . * (1|atan)/45;
 def sind: radians|sin;
 def cosd: radians|cos;
 def sq: . * .;

   (((lat2 - lat1)/2) | sind | sq) as $dlat
 | (((lon2 - lon1)/2) | sind | sq) as $dlon
 | 2 * 6372.8 * (( $dlat + (lat1|cosd) * (lat2|cosd) * $dlon ) | sqrt | asin) ;</lang>

Example: 

haversine(36.12; -86.67; 33.94; -118.4)
# 2887.2599506071106

Julia

<lang julia>julia> haversine(lat1,lon1,lat2,lon2) = 2 * 6372.8 * asin(sqrt(sind((lat2-lat1)/2)^2 + cosd(lat1) * cosd(lat2) * sind((lon2 - lon1)/2)^2))

  1. method added to generic function haversine

julia> haversine(36.12,-86.67,33.94,-118.4) 2887.2599506071106</lang>

Liberty BASIC

<lang lb>print "Haversine distance: "; using( "####.###########", havDist( 36.12, -86.67, 33.94, -118.4)); " km." end function havDist( th1, ph1, th2, ph2) 

 degtorad   = acs(-1)/180
 diameter   = 2 * 6372.8
   LgD      = degtorad  * (ph1 - ph2)
   th1      = degtorad  * th1
   th2      = degtorad  * th2
   dz       = sin( th1) - sin( th2)
   dx       = cos( LgD) * cos( th1) - cos( th2)
   dy       = sin( LgD) * cos( th1)
   havDist  = asn( ( dx^2 +dy^2 +dz^2)^0.5 /2) *diameter

end function</lang> 

Haversine distance: 2887.25995060711  km.

Mathematica

Inputs assumed in degrees. Sin and Haversine expect arguments in radians; the built-in variable 'Degree' converts from degrees to radians. <lang Mathematica> distance[{theta1_, phi1_}, {theta2_, phi2_}] := 

2*6378.14 ArcSin@
  Sqrt[Haversine[(theta2 - theta1) Degree] + 
    Cos[theta1*Degree] Cos[theta2*Degree] Haversine[(phi2 - phi1) Degree]]

</lang> Usage: 

distance[{36.12, -86.67}, {33.94, -118.4}]
Output:
2889.68

MATLAB / Octave

<lang Matlab>function rad = radians(degree) % degrees to radians 

   rad = degree .* pi / 180;

end;

function [a,c,dlat,dlon]=haversine(lat1,lon1,lat2,lon2) % HAVERSINE_FORMULA.AWK - converted from AWK 

   dlat = radians(lat2-lat1);
   dlon = radians(lon2-lon1);
   lat1 = radians(lat1);
   lat2 = radians(lat2);
   a = (sin(dlat./2)).^2 + cos(lat1) .* cos(lat2) .* (sin(dlon./2)).^2;
   c = 2 .* asin(sqrt(a));
   arrayfun(@(x) printf("distance: %.4f km\n",6372.8 * x), c);

end;

[a,c,dlat,dlon] = haversine(36.12,-86.67,33.94,-118.40); % BNA to LAX</lang>

Output:
distance: 2887.2600 km

Maxima

<lang maxima>dms(d, m, s) := (d + m/60 + s/3600)*%pi/180$

great_circle_distance(lat1, long1, lat2, long2) := 

  12742*asin(sqrt(sin((lat2 - lat1)/2)^2 + cos(lat1)*cos(lat2)*sin((long2 - long1)/2)^2))$

/* Coordinates are found here: 

     http://www.airport-data.com/airport/BNA/
     http://www.airport-data.com/airport/LAX/   */

great_circle_distance(dms( 36, 7, 28.10), -dms( 86, 40, 41.50), 

                     dms( 33, 56, 32.98), -dms(118, 24, 29.05)), numer;

/* 2886.326609413624 */</lang>

МК-61/52

<lang>П3 -> П2 -> П1 -> П0 пи 1 8 0 / П4 ИП1 МГ ИП3 МГ - ИП4 * П1 ИП0 МГ ИП4 * П0 ИП2 МГ ИП4 * П2 ИП0 sin ИП2 sin - П8 ИП1 cos ИП0 cos * ИП2 cos - П6 ИП1 sin ИП0 cos * П7 ИП6 x^2 ИП7 x^2 ИП8 x^2 + + КвКор 2 / arcsin 2 * ИП5 * С/П</lang>

Input: 6371,1 as a radius of the Earth, taken as the ball, or 6367,554 as an average radius of the Earth, or 6367,562 as an approximation of the radius of the average circumference (by Krasovsky's ellipsoid) to Р5; В/О lat1 С/П long1 С/П lat2 С/П long2 С/П; the coordinates must be entered as degrees,minutes (example: 46°50' as 46,5).

Test:

  • N 36°7.2', W 86°40.2' - N 33°56.4', W 118°24.0' (Nashville - Los Angeles):
Input: 6371,1 П5 36,072 С/П -86,402 С/П 33,564 С/П -118,24 С/П
Output: 2886,4897.
  • N 54°43', E 20°3' - N 43°07', E 131°54' (Kaliningrad - Vladivostok):
Input: 6371,1 П5 54,43 С/П 20,3 С/П 43,07 С/П 131,54 С/П
Output: 7357,4526.

Nim

<lang nim>import math

proc radians(x): float = x * Pi / 180

proc haversine(lat1, lon1, lat2, lon2): float = 

 const r = 6372.8 # Earth radius in kilometers
 let
   dLat = radians(lat2 - lat1)
   dLon = radians(lon2 - lon1)
   lat1 = radians(lat1)
   lat2 = radians(lat2)

   a = sin(dLat/2)*sin(dLat/2) + cos(lat1)*cos(lat2)*sin(dLon/2)*sin(dLon/2)
   c = 2*arcsin(sqrt(a))

 result = r * c

echo haversine(36.12, -86.67, 33.94, -118.40)</lang>

Output:
2.8872599506071115e+03

Oberon-2

Works with oo2c version2 <lang oberon2> MODULE Haversines; IMPORT 

 LRealMath,
 Out;
 
 PROCEDURE Distance(lat1,lon1,lat2,lon2: LONGREAL): LONGREAL;
 CONST
   r = 6372.8D0; (* Earth radius as LONGREAL *)
   to_radians = LRealMath.pi / 180.0D0;
 VAR
   d,ph1,th1,th2: LONGREAL;
   dz,dx,dy: LONGREAL;
 BEGIN
   d := lon1 - lon2;
   ph1 := d * to_radians;
   th1 := lat1 * to_radians;
   th2 := lat2 * to_radians;
   
   dz := LRealMath.sin(th1) - LRealMath.sin(th2);
   dx := LRealMath.cos(ph1) * LRealMath.cos(th1) - LRealMath.cos(th2);
   dy := LRealMath.sin(ph1) * LRealMath.cos(th1);
   
   RETURN LRealMath.arcsin(LRealMath.sqrt(LRealMath.power(dx,2.0) + LRealMath.power(dy,2.0) + LRealMath.power(dz,2.0)) / 2.0) * 2.0 * r;
 END Distance;

BEGIN 

 Out.LongRealFix(Distance(36.12,-86.67,33.94,-118.4),6,10);Out.Ln

END Haversines. </lang> Output: 

2887.2602975600

Objeck

<lang objeck> bundle Default { 

 class Haversine {
   function : Dist(th1 : Float, ph1 : Float, th2 : Float, ph2 : Float) ~ Float {
     ph1 -= ph2;
     ph1 := ph1->ToRadians();
     th1 := th1->ToRadians();
     th2 := th2->ToRadians();

     dz := th1->Sin()- th2->Sin();
     dx := ph1->Cos() * th1->Cos() - th2->Cos();
     dy := ph1->Sin() * th1->Cos();

     return ((dx * dx + dy * dy + dz * dz)->SquareRoot() / 2.0)->ArcSin() * 2 * 6371.0;
   }

   function : Main(args : String[]) ~ Nil {
     IO.Console->Print("distance: ")->PrintLine(Dist(36.12, -86.67, 33.94, -118.4));
   }
 }

} </lang>

Output:
distance: 2886.44

Objective-C

<lang objc>+ (double) distanceBetweenLat1:(double)lat1 lon1:(double)lon1 

                         lat2:(double)lat2 lon2:(double)lon2 {
   //degrees to radians
   double lat1rad = lat1 * M_PI/180; 
   double lon1rad = lon1 * M_PI/180;
   double lat2rad = lat2 * M_PI/180;
   double lon2rad = lon2 * M_PI/180;
   
   //deltas
   double dLat = lat2rad - lat1rad;
   double dLon = lon2rad - lon1rad;

   double a = sin(dLat/2) * sin(dLat/2) + sin(dLon/2) * sin(dLon/2) * cos(lat1rad) * cos(lat2rad);
   double c = 2 * asin(sqrt(a));
   double R = 6372.8;
   return R * c;

}</lang>

OCaml

The core calculation is fairly straightforward, but with an eye toward generality and reuse, this is how I might start: <lang ocaml>(* Preamble -- some math, and an "angle" type which might be part of a common library. *) let pi = 4. *. atan 1. let radians_of_degrees = ( *. ) (pi /. 180.) let haversin theta = 0.5 *. (1. -. cos theta)

(* The angle type can track radians or degrees, which I'll use for automatic conversion. *) type angle = Deg of float | Rad of float let as_radians = function 

 | Deg d -> radians_of_degrees d
 | Rad r -> r

(* Demonstrating use of a module, and record type. *) module LatLong = struct 

 type t = { lat: float; lng: float }
 let of_angles lat lng = { lat = as_radians lat; lng = as_radians lng }
 let sub a b = { lat = a.lat-.b.lat; lng = a.lng-.b.lng }

 let dist radius a b =
   let d = sub b a in
   let h = haversin d.lat +. haversin d.lng *. cos a.lat *. cos b.lat in
   2. *. radius *. asin (sqrt h)

end

(* Now we can use the LatLong module to construct coordinates and calculate 

* great-circle distances.
* NOTE radius and resulting distance are in the same measure, and units could
* be tracked for this too... but who uses miles? ;) *)

let earth_dist = LatLong.dist 6372.8 and bna = LatLong.of_angles (Deg 36.12) (Deg (-86.67)) and lax = LatLong.of_angles (Deg 33.94) (Deg (-118.4)) in earth_dist bna lax;;</lang>

If the above is fed to the REPL, the last line will produce this: 

# earth_dist bna lax;;
- : float = 2887.25995060711102

Oforth

<lang Oforth>func: haversine(lat1, lon1, lat2, lon2) { | lat lon | 

  lat2 lat1 - asRadian ->lat
  lon2 lon1 - asRadian ->lon

  lon 2 / sin sq lat1 asRadian cos * lat2 asRadian cos * 
  lat 2 / sin sq + sqrt asin 2 * 6372.8 *

}

haversine(36.12, -86.67, 33.94, -118.40) println</lang>

Output:
2887.25995060711

ooRexx

Translation of: REXX

The rxmath library provides the required functions. <lang oorexx>/*REXX pgm calculates distance between Nashville & Los Angles airports. */ say " Nashville: north 36º 7.2', west 86º 40.2' = 36.12º, -86.67º" say "Los Angles: north 33º 56.4', west 118º 24.0' = 33.94º, -118.40º" say dist=surfaceDistance(36.12, -86.67, 33.94, -118.4) kdist=format(dist/1 ,,2) /*show 2 digs past decimal point.*/ mdist=format(dist/1.609344,,2) /* " " " " " " */ ndist=format(mdist*5280/6076.1,,2) /* " " " " " " */ say ' distance between= ' kdist " kilometers," say ' or ' mdist " statute miles," say ' or ' ndist " nautical or air miles." exit /*stick a fork in it, we're done.*/ /*----------------------------------SURFACEDISTANCE subroutine----------*/ surfaceDistance: arg th1,ph1,th2,ph2 /*use haversine formula for dist.*/ 

 radius = 6372.8                      /*earth's mean radius in km      */
 ph1 = ph1-ph2
 x = cos(ph1) * cos(th1) - cos(th2)
 y = sin(ph1) * cos(th1)
 z = sin(th1) - sin(th2)
 return radius * 2 * aSin(sqrt(x**2+y**2+z**2)/2 )

cos: Return RxCalcCos(arg(1)) sin: Return RxCalcSin(arg(1)) asin: Return RxCalcArcSin(arg(1),,'R') sqrt: Return RxCalcSqrt(arg(1))

requires rxMath library</lang>
Output:
 Nashville:  north 36º  7.2', west  86º 40.2'   =   36.12º,  -86.67º
Los Angles:  north 33º 56.4', west 118º 24.0'   =   33.94º, -118.40º

 distance between=   2887.26  kilometers,
               or    1794.06  statute miles,
               or    1559.00  nautical or air miles.

PARI/GP

<lang parigp>dist(th1, th2, ph)={ 

 my(v=[cos(ph)*cos(th1)-cos(th2),sin(ph)*cos(th1),sin(th1)-sin(th2)]);
 asin(sqrt(norml2(v))/2)

}; distEarth(th1, ph1, th2, ph2)={ 

 my(d=12742, deg=Pi/180); \\ Authalic diameter of the Earth
 d*dist(th1*deg, th2*deg, (ph1-ph2)*deg)

}; distEarth(36.12, -86.67, 33.94, -118.4)</lang>

Output:
%1 = 2886.44444

Pascal

Works with: Free_Pascal
Library: Math

<lang pascal>Program HaversineDemo(output);

uses 

 Math;

function haversineDist(th1, ph1, th2, ph2: double): double; 

 const
  diameter = 2 * 6372.8;
 var
   dx, dy, dz: double;
 begin
   ph1 := degtorad(ph1 - ph2);
   th1 := degtorad(th1);
   th2 := degtorad(th2);

   dz := sin(th1) - sin(th2);
   dx := cos(ph1) * cos(th1) - cos(th2);
   dy := sin(ph1) * cos(th1);
   haversineDist := arcsin(sqrt(dx**2 + dy**2 + dz**2) / 2) * diameter;
 end;

begin 

 writeln ('Haversine distance: ', haversineDist(36.12, -86.67, 33.94, -118.4):7:2, ' km.');

end.</lang>

Output:
Haversine distance: 2887.26 km.

Perl 6

<lang perl6>class EarthPoint { 

       has $.lat; # latitude
       has $.lon; # longitude

       has $earth_radius = 6371; # mean earth radius
       has $radian_ratio = pi / 180;

       # accessors for radians
       method latR { $.lat * $radian_ratio }
       method lonR { $.lon * $radian_ratio }

       method haversine-dist(EarthPoint $p) {

               my EarthPoint $arc .= new(
                       lat => $!lat - $p.lat,
                       lon => $!lon - $p.lon );

               my $a = sin($arc.latR/2) ** 2 + sin($arc.lonR/2) ** 2
                       * cos($.latR) * cos($p.latR);
               my $c = 2 * asin( sqrt($a) );

               return $earth_radius * $c;
       }

}

my EarthPoint $BNA .= new(lat => 36.12, lon => -86.67); my EarthPoint $LAX .= new(lat => 33.94, lon => -118.4);

say $BNA.haversine-dist($LAX); # 2886.44444099822</lang>

PHP

<lang php>class POI { 

   private $latitude;
   private $longitude;
   public function __construct($latitude, $longitude) {
       $this->latitude = deg2rad($latitude);
       $this->longitude = deg2rad($longitude);
   }
   public function getLatitude() return $this->latitude;
   public function getLongitude() return $this->longitude;
   public function getDistanceInMetersTo(POI $other) {
       $radiusOfEarth = 6371000;// Earth's radius in meters.
       $diffLatitude = $other->getLatitude() - $this->latitude;
       $diffLongitude = $other->getLongitude() - $this->longitude;
       $a = sin($diffLatitude / 2) * sin($diffLatitude / 2) +
           cos($this->latitude) * cos($other->getLatitude()) *
           sin($diffLongitude / 2) * sin($diffLongitude / 2);
       $c = 2 * asin(sqrt($a));
       $distance = $radiusOfEarth * $c;
       return $distance;
   }

}</lang> Test: <lang php>$user = new POI($_GET["latitude"], $_GET["longitude"]); $poi = new POI(19,69276, -98,84350); // Piramide del Sol, Mexico echo $user->getDistanceInMetersTo($poi);</lang>

PicoLisp

<lang PicoLisp>(scl 12) (load "@lib/math.l")

(de haversine (Th1 Ph1 Th2 Ph2) 

  (setq
     Ph1 (*/ (- Ph1 Ph2) pi 180.0)
     Th1 (*/ Th1 pi 180.0)
     Th2 (*/ Th2 pi 180.0) )
  (let
     (DX (- (*/ (cos Ph1) (cos Th1) 1.0) (cos Th2))
        DY (*/ (sin Ph1) (cos Th1) 1.0)
        DZ (- (sin Th1) (sin Th2)) )
     (* `(* 2 6371)
        (asin
           (/
              (sqrt (+ (* DX DX) (* DY DY) (* DZ DZ)))
              2 ) ) ) ) )</lang>

Test: <lang PicoLisp>(prinl 

  "Haversine distance: "
  (round (haversine 36.12 -86.67 33.94 -118.4))
  " km" )</lang>
Output:
Haversine distance: 2,886.444 km

PL/I

<lang PL/I>test: procedure options (main); /* 12 January 2014. Derived from Fortran version */ 

  declare d float;

  d = haversine(36.12, -86.67, 33.94, -118.40);  /* BNA to LAX */
  put edit ( 'distance: ', d, ' km') (A, F(10,3)); /* distance: 2887.2600 km */


degrees_to_radians: procedure (degree) returns (float); 

  declare degree float nonassignable;
  declare pi float (15) initial ( (4*atan(1.0d0)) );

  return ( degree*pi/180 );

end degrees_to_radians;

haversine: procedure (deglat1, deglon1, deglat2, deglon2) returns (float); 

  declare (deglat1, deglon1, deglat2, deglon2) float nonassignable;
  declare (a, c, dlat, dlon, lat1, lat2) float;
  declare radius float value (6372.8);

  dlat = degrees_to_radians(deglat2-deglat1);
  dlon = degrees_to_radians(deglon2-deglon1);
  lat1 = degrees_to_radians(deglat1);
  lat2 = degrees_to_radians(deglat2);
  a = (sin(dlat/2))**2 + cos(lat1)*cos(lat2)*(sin(dlon/2))**2;
  c = 2*asin(sqrt(a));
  return ( radius*c );

end haversine;

end test;</lang>

Output:
distance:   2887.260 km

Python

<lang python>from math import radians, sin, cos, sqrt, asin

def haversine(lat1, lon1, lat2, lon2): 

 R = 6372.8 # Earth radius in kilometers

 dLat = radians(lat2 - lat1)
 dLon = radians(lon2 - lon1)
 lat1 = radians(lat1)
 lat2 = radians(lat2)

 a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2
 c = 2*asin(sqrt(a))

 return R * c

>>> haversine(36.12, -86.67, 33.94, -118.40) 2887.2599506071106 >>> </lang>

R

<lang r>dms_to_rad <- function(d, m, s) (d + m / 60 + s / 3600) * pi / 180

  1. Volumetric mean radius is 6371 km, see http://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html
  2. The diameter is thus 12742 km

great_circle_distance <- function(lat1, long1, lat2, long2) { 

  a <- sin(0.5 * (lat2 - lat1))
  b <- sin(0.5 * (long2 - long1))
  12742 * asin(sqrt(a * a + cos(lat1) * cos(lat2) * b * b))

}

  1. Coordinates are found here:
  2. http://www.airport-data.com/airport/BNA/
  3. http://www.airport-data.com/airport/LAX/

great_circle_distance( 

  dms_to_rad(36,  7, 28.10), dms_to_rad( 86, 40, 41.50),   # Nashville International Airport (BNA)
  dms_to_rad(33, 56, 32.98), dms_to_rad(118, 24, 29.05))  # Los Angeles International Airport (LAX)
  1. Output: 2886.327</lang>

Racket

Almost the same as the Scheme version. <lang racket>

  1. lang racket

(require math) (define earth-radius 6371)

(define (distance lat1 long1 lat2 long2) 

 (define (h a b) (sqr (sin (/ (- b a) 2))))
 (* 2 earth-radius 
    (asin (sqrt (+ (h lat1 lat2) 
                   (* (cos lat1) (cos lat2) (h long1 long2)))))))

(define (deg-to-rad d m s) 

 (* (/ pi 180) (+ d (/ m 60) (/ s 3600))))

(distance (deg-to-rad 36 7.2 0) (deg-to-rad 86 40.2 0) 

         (deg-to-rad 33 56.4 0) (deg-to-rad 118 24.0 0))

</lang>

Output:
2886.444442837984

Raven

Translation of: Groovy

<lang Raven>define PI 

 -1 acos

define toRadians use $degree 

 $degree PI * 180 /

define haversine use $lat1, $lon1, $lat2, $lon2 

 6372.8 as $R
 # In kilometers
 $lat2 $lat1 - toRadians   as $dLat
 $lon2 $lon1 - toRadians   as $dLon
 $lat1 toRadians  as $lat1
 $lat2 toRadians  as $lat2

 $dLat 2 /  sin 
 $dLat 2 /  sin *
 $dLon 2 /  sin
 $dLon 2 /  sin *
 $lat1 cos * 
 $lat2 cos * +        as $a
 $a sqrt  asin  2 *   as $c
 $R $c *

}

-118.40 33.94 -86.67 36.12 haversine "haversine: %.15g\n" print</lang>

Output:
haversine: 2887.25995060711

REXX

REXX doen't have most of the higher math functions, so they are included here as subroutines. <lang rexx>/*REXX pgm calculates distance between Nashville & Los Angles airports. */ say " Nashville: north 36º 7.2', west 86º 40.2' = 36.12º, -86.67º" say "Los Angles: north 33º 56.4', west 118º 24.0' = 33.94º, -118.40º" say $.= /*set defaults for subroutines. */ dist=surfaceDistance(36.12, -86.67, 33.94, -118.4) kdist=format(dist/1 ,,2) /*show 2 digs past decimal point.*/ mdist=format(dist/1.609344,,2) /* " " " " " " */ ndist=format(mdist*5280/6076.1,,2) /* " " " " " " */ say ' distance between= ' kdist " kilometers," say ' or ' mdist " statute miles," say ' or ' ndist " nautical or air miles." exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────SURFACEDISTANCE subroutine──────────*/ surfaceDistance: arg th1,ph1,th2,ph2 /*use haversine formula for dist.*/ numeric digits digits()*2 /*double the number of digits. */ radius = 6372.8 /*earth's mean radius in km */ 

  ph1 = d2r(ph1-ph2)                  /*convert degs──►radians & reduce*/
  ph2 = d2r(ph2)                      /*   "      "       "    "    "  */
  th1 = d2r(th1)                      /*   "      "       "    "    "  */
  th2 = d2r(th2)
                  x = cos(ph1) * cos(th1) - cos(th2)
                  y = sin(ph1) * cos(th1)
                  z = sin(th1) - sin(th2)

return radius * 2 * aSin(sqrt(x**2+y**2+z**2)/2 ) /*═════════════════════════════general 1-line subs══════════════════════*/ d2d: return arg(1) // 360 /*normalize degrees. */ d2r: return r2r(arg(1)*pi() / 180) /*normalize and convert deg──►rad*/ r2d: return d2d((arg(1)*180 / pi())) /*normalize and convert rad──►deg*/ r2r: return arg(1) // (2*pi()) /*normalize radians. */ p: return word(arg(1),1) /*pick the first of two words. */ pi: if $.pipi== then $.pipi=$pi(); return $.pipi /*return π.*/

aCos: procedure expose $.; arg x; if x<-1|x>1 then call $81r -1,1,x,"ACOS" 

        return .5*pi()-aSin(x)        /*$81R  says arg is out of range,*/
                                      /*    and it isn't included here.*/

aSin: procedure expose $.; parse arg x 

        if x<-1 | x>1   then call $81r -1,1,x,"ASIN";    s=x*x
        if abs(x)>=.7   then return sign(x)*aCos(sqrt(1-s),'-ASIN')
        z=x;  o=x;  p=z;     do j=2 by 2;  o=o*s*(j-1)/j;  z=z+o/(j+1)
        if z=p then leave;   p=z;  end;                  return z

cos: procedure expose $.; parse arg x; x=r2r(x); a=abs(x) 

        numeric fuzz min(9,digits()-9);        if a=pi()   then return -1
        if a=pi()/2 | a=2*pi() then return 0;  if a=pi()/3 then return .5
        if a=2*pi()/3  then return -.5;        return .sinCos(1,1,-1)

sin: procedure expose $.; parse arg x; x=r2r(x); 

        numeric fuzz  min(5,digits()-3)
        if abs(x)=pi()  then return 0;         return .sinCos(x,x,1)

.sinCos: parse arg z,_,i; x=x*x; p=z; do k=2 by 2; _=-_*x/(k*(k+i)); z=z+_ 

        if z=p then leave; p=z;  end;  return z    /*used by SIN & COS.*/

sqrt: procedure; parse arg x; if x=0 then return 0; d=digits() 

 numeric digits 11; g=.sqrtGuess(); do j=0 while p>9; m.j=p; p=p%2+1; end
 do k=j+5 to 0 by -1; if m.k>11 then numeric digits m.k; g=.5*(g+x/g);end
 numeric digits d;  return g/1

.sqrtGuess: numeric form; m.=11; p=d+d%4+2 

 parse value format(x,2,1,,0) 'E0' with g 'E' _ .;    return  g*.5'E'_%2

$pi: return , '3.1415926535897932384626433832795028841971693993751058209749445923078'||, '164062862089986280348253421170679821480865132823066470938446095505822'||, '3172535940812848111745028410270193852110555964462294895493038196'</lang> <lang> ┌───────────────────────────────────────────────────────────────────────┐ 

 │ A note on built-in functions.  REXX doesn't have a lot of mathmatical │
 │ or  (particularly) trigomentric  functions,  so REXX programmers have │
 │ to write their own.  Usually, this is done once, or most likely,  one │
 │ is borrowed from another program.  Knowing this, the one that is used │
 │ has a lot of boilerplate in it.  Once coded and throughly debugged, I │
 │ put those commonly-used subroutines into the  "1-line sub"  section.  │
 │                                                                       │
 │ Programming note:  the  "general 1-line"  subroutines are taken from  │
 │ other programs that I wrote, but I broke up their one line of source  │
 │ so it can be viewed without shifting the viewing window.              │
 │                                                                       │
 │ The  "er 81"  [which won't happen here]  just shows an error telling  │
 │ the legal range for  ARCxxx  functions   (in this case:  -1 ──► +1).  │
 │                                                                       │
 │ Similarly,  the   SQRT   function checks for a negative argument      │
 │ [which again, won't happen here].                                     │
 │                                                                       │
 │ The pi constant (as used here) is actually a much more robust function│
 │ and will return up to one million digits in the real version.         │
 │                                                                       │
 │ One bad side effect is that, like a automobile without a hood, you see│
 │ all the dirty stuff going on.    Also, don't visit a sausage factory. │
 └───────────────────────────────────────────────────────────────────────┘</lang>
Output:
 Nashville:  north 36º  7.2', west  86º 40.2'   =   36.12º,  -86.67º
Los Angles:  north 33º 56.4', west 118º 24.0'   =   33.94º, -118.40º

 distance between=   2887.26  kilometers,
               or    1794.06  statute miles,
               or    1559.00  nautical or air miles.

Ruby

<lang ruby>include Math

Radius = 6371 # rough radius of the Earth, in kilometers

def spherical_distance(start_coords, end_coords) 

 lat1, long1 = deg2rad *start_coords
 lat2, long2 = deg2rad *end_coords
 2 * Radius * asin(sqrt(sin((lat2-lat1)/2)**2 + cos(lat1) * cos(lat2) * sin((long2 - long1)/2)**2))

end

def deg2rad(lat, long) 

 [lat * PI / 180, long * PI / 180]

end

bna = [36.12, -86.67] lax = [33.94, -118.4]

puts "%.1f" % spherical_distance(bna, lax)</lang>

Output:
2886.4

Run BASIC

<lang runbasic> D2R = atn(1)/45 

   diam  = 2 * 6372.8

Lg1m2 = ((-86.67)-(-118.4)) * D2R Lt1 = 36.12 * D2R ' degrees to rad Lt2 = 33.94 * D2R 

   dz    = sin(Lt1) - sin(Lt2)
   dx    = cos(Lg1m2) * cos(Lt1) - cos(Lt2)
   dy    = sin(Lg1m2) * cos(Lt1)
   hDist = asn((dx^2 + dy^2 + dz^2)^0.5 /2) * diam

print "Haversine distance: ";using("####.#############",hDist);" km." 

'Tips:  ( 36 deg 7 min 12 sec ) = print 36+(7/60)+(12/3600).  Produces: 36.12 deg
'
'       Put  "36.12,-86.67"  into http://maps.google.com  (no quotes).  Click map,
'       satellite, center the pin "A", zoom in, and see airport.  Extra: in "Get
'       Directions" enter  36.12,-86.66999 and see pin "B" about one meter away.
'       (.00089846878 km., or 35.37 in.)
'

' This code also works in Liberty BASIC.</lang>Output 

Haversine distance: 2887.2599506071104 km.

SAS

<lang SAS> options minoperator;

%macro haver(lat1, long1, lat2, long2, type=D, dist=K);

%if %upcase(&type) in (D DEG DEGREE DEGREES) %then %do; %let convert = constant('PI')/180; %end; %else %if %upcase(&type) in (R RAD RADIAN RADIANS) %then %do; %let convert = 1; %end; %else %do; %put ERROR - Enter RADIANS or DEGREES for type.; %goto exit; %end;

%if %upcase(&dist) in (M MILE MILES) %then %do; %let distrat = 1.609344; %end; %else %if %upcase(&dist) in (K KM KILOMETER KILOMETERS) %then %do; %let distrat = 1; %end; %else %do; %put ERROR - Enter M on KM for dist; %goto exit; %end;

data _null_; convert = &convert; lat1 = &lat1 * convert; lat2 = &lat2 * convert; long1 = &long1 * convert; long2 = &long2 * convert;

diff1 = lat2 - lat1; diff2 = long2 - long1;

part1 = sin(diff1/2)**2; part2 = cos(lat1)*cos(lat2); part3 = sin(diff2/2)**2;

root = sqrt(part1 + part2*part3);

dist = 2 * 6372.8 / &distrat * arsin(root);

put "Distance is " dist "%upcase(&dist)"; run;

%exit: %mend;

%haver(36.12, -86.67, 33.94, -118.40); </lang>

Output:
Distance is 2887.2599506 K

Scala

<lang scala>import math._

object Haversine { 

  val R = 6372.8  //radius in km

  def haversine(lat1:Double, lon1:Double, lat2:Double, lon2:Double)={
     val dLat=(lat2 - lat1).toRadians
     val dLon=(lon2 - lon1).toRadians

     val a = pow(sin(dLat/2),2) + pow(sin(dLon/2),2) * cos(lat1.toRadians) * cos(lat2.toRadians)
     val c = 2 * asin(sqrt(a))
     R * c
  }

  def main(args: Array[String]): Unit = {
     println(haversine(36.12, -86.67, 33.94, -118.40))
 }

}</lang>

Output:
2887.2599506071106

Scheme

<lang scheme>(define earth-radius 6371) (define pi (acos -1))

(define (distance lat1 long1 lat2 long2) (define (h a b) (expt (sin (/ (- b a) 2)) 2)) (* 2 earth-radius (asin (sqrt (+ (h lat1 lat2) (* (cos lat1) (cos lat2) (h long1 long2)))))))

(define (deg-to-rad d m s) (* (/ pi 180) (+ d (/ m 60) (/ s 3600))))

(distance (deg-to-rad 36 7.2 0) (deg-to-rad 86 40.2 0) 

         (deg-to-rad 33 56.4 0) (deg-to-rad 118 24.0 0))
2886.444442837984</lang>

Seed7

<lang seed7>$ include "seed7_05.s7i"; 

 include "float.s7i";
 include "math.s7i";

const func float: greatCircleDistance (in float: latitude1, in float: longitude1, 

   in float: latitude2, in float: longitude2) is func
 result
   var float: distance is 0.0;
 local
   const float: EarthRadius is 6372.8;  # Average great-elliptic or great-circle radius in kilometers
 begin
   distance := 2.0 * EarthRadius * asin(sqrt(sin(0.5 * (latitude2 - latitude1)) ** 2 +
                                             cos(latitude1) * cos(latitude2) *
                                             sin(0.5 * (longitude2 - longitude1)) ** 2));
 end func;

const func float: degToRad (in float: degrees) is 

 return degrees * 0.017453292519943295769236907684886127;

const proc: main is func 

 begin
   writeln("Distance in kilometers between BNA and LAX");
   writeln(greatCircleDistance(degToRad(36.12), degToRad(-86.67),  # Nashville International Airport (BNA)
                               degToRad(33.94), degToRad(-118.4))  # Los Angeles International Airport (LAX)
           digits 2);
 end func;</lang>
Output:
2887.26

Sidef

Translation of: Perl 6

<lang ruby>class EarthPoint(lat, lon) { 

   static earth_radius = 6371;       # mean earth radius
   static radian_ratio = Math.pi/180;

   # accessors for radians
   method latR { self.lat * radian_ratio };
   method lonR { self.lon * radian_ratio };

   method haversine_dist(p is EarthPoint) {
       var arc = __CLASS__.new(
           self.lat - p.lat,
           self.lon - p.lon,
       );

       var a = [ Math.pow(Math.sin(arc.latR / 2), 2),
                 Math.pow(Math.sin(arc.lonR / 2), 2)
                   * Math.cos(self.latR) * Math.cos(p.latR),
               ].sum;

       earth_radius * Math.asin(Math.sqrt(a)) * 2;
   }

}

var BNA = EarthPoint.new(lat: 36.12, lon: -86.67); var LAX = EarthPoint.new(lat: 33.94, lon: -118.4);

say BNA.haversine_dist(LAX); 

   # => 2886.4444428379832997471578239457467165513238</lang>

Swift

Translation of: Objective-C

<lang Swift>import Foundation

func haversine(lat1:Double, lon1:Double, lat2:Double, lon2:Double) -> Double { 

   let lat1rad = lat1 * M_PI/180
   let lon1rad = lon1 * M_PI/180
   let lat2rad = lat2 * M_PI/180
   let lon2rad = lon2 * M_PI/180
   
   let dLat = lat2rad - lat1rad
   let dLon = lon2rad - lon1rad
   let a = sin(dLat/2) * sin(dLat/2) + sin(dLon/2) * sin(dLon/2) * cos(lat1rad) * cos(lat2rad)
   let c = 2 * asin(sqrt(a))
   let R = 6372.8
   
   return R * c

}

println(haversine(36.12, -86.67, 33.94, -118.40))</lang>

Output:
2887.25995060711

Tcl

Translation of: Groovy

<lang tcl>package require Tcl 8.5 proc haversineFormula {lat1 lon1 lat2 lon2} { 

   set rads [expr atan2(0,-1)/180]
   set R 6372.8    ;# In kilometers

   set dLat [expr {($lat2-$lat1) * $rads}]
   set dLon [expr {($lon2-$lon1) * $rads}]
   set lat1 [expr {$lat1 * $rads}]
   set lat2 [expr {$lat2 * $rads}]

   set a [expr {sin($dLat/2)**2 + sin($dLon/2)**2*cos($lat1)*cos($lat2)}]
   set c [expr {2*asin(sqrt($a))}]
   return [expr {$R * $c}]

}

  1. Don't bother with too much inappropriate accuracy!

puts [format "distance=%.1f km" [haversineFormula 36.12 -86.67 33.94 -118.40]]</lang>

Output:
distance=2887.3 km

UBASIC

<lang basic> 

  10  Point 7    'Sets decimal display to 32 places (0+.1^56)
  20  Rf=#pi/180 'Degree -> Radian Conversion
 100 ?Using(,7),.DxH(36+7.2/60,-(86+40.2/60),33+56.4/60,-(118+24/60));" km"
 999  End
1000 '*** Haversine Distance Function ***
1010 .DxH(Lat_s,Long_s,Lat_f,Long_f)
1020  L_s=Lat_s*rf:L_f=Lat_f*rf:LD=L_f-L_s:MD=(Long_f-Long_s)*rf
1030  Return(12745.6*asin( (sin(.5*LD)^2+cos(L_s)*cos(L_f)*sin(.5*MD)^2)^.5))
 

Run
 2887.2599506 km
OK

</lang>

X86 Assembly

Assemble with tasm /m /l; tlink /t <lang asm>0000 .model tiny 0000 .code 

                                    .486
                                    org     100h            ;.com files start here

0100 9B DB E3 start: finit ;initialize floating-point unit (FPU) 

                            ;Great circle distance =
                            ; 2.0*Radius * ASin( sqrt( Haversine(Lat2-Lat1) +
                            ;                          Haversine(Lon2-Lon1)*Cos(Lat1)*Cos(Lat2) ) )

0103 D9 06 0191r fld Lat2 ;push real onto FPU stack 0107 D8 26 018Dr fsub Lat1 ;subtract real from top of stack (st(0) = st) 010B E8 0070 call Haversine  ;(1.0-cos(st)) / 2.0 010E D9 06 0199r fld Lon2 ;repeat for longitudes 0112 D8 26 0195r fsub Lon1 0116 E8 0065 call Haversine ;st(1)=Lats; st=Lons 0119 D9 06 018Dr fld Lat1 011D D9 FF fcos ;replace st with its cosine 011F D9 06 0191r fld Lat2 0123 D9 FF fcos ;st=cos(Lat2); st(1)=cos(Lat1); st(2)=Lats; st(3)=Lons 0125 DE C9 fmul ;st=cos(Lat2)*cos(Lat1); st(1)=Lats; st(2)=Lons 0127 DE C9 fmul ;st=cos(Lat2)*cos(Lat1)*Lats; st(1)=Lons 0129 DE C1 fadd ;st=cos(Lat2)*cos(Lat1)*Lats + Lons 012B D9 FA fsqrt ;replace st with its square root 

                            ;asin(x) = atan(x/sqrt(1-x^2))

012D D9 C0 fld st ;duplicate tos 012F D8 C8 fmul st, st ;x^2 0131 D9 E8 fld1 ;get 1.0 0133 DE E1 fsubr ;1 - x^2 0135 D9 FA fsqrt ;sqrt(1-x^2) 0137 D9 F3 fpatan ;take atan(st(1)/st) 0139 D8 0E 019Dr fmul Radius2  ;*2.0*Radius 

                            ;Display value in FPU's top of stack (st)
     =0004                  before  equ     4               ;places before
     =0002                  after   equ     2               ; and after decimal point
     =0001                  scaler  =       1               ;"=" allows scaler to be redefined, unlike equ
                                    rept    after           ;repeat block "after" times
                            scaler  =       scaler*10
                                    endm                    ;scaler now = 10^after

013D 66| 6A 64 push dword ptr scaler;use stack for convenient memory location 0140 67| DA 0C 24 fimul dword ptr [esp] ;st:= st*scaler 0144 67| DB 1C 24 fistp dword ptr [esp] ;round st to nearest integer 0148 66| 58 pop eax  ; and put it into eax

014A 66| BB 0000000A mov ebx, 10 ;set up for idiv instruction 0150 B9 0006 mov cx, before+after;set up loop counter 0153 66| 99 ro10: cdq ;convert double to quad; i.e: edx:= 0 0155 66| F7 FB idiv ebx ;eax:= edx:eax/ebx; remainder in edx 0158 52 push dx ;save least significant digit on stack 0159 E2 F8 loop ro10 ;cx--; loop back if not zero

015B B1 06 mov cl, before+after;(ch=0) 015D B3 00 mov bl, 0 ;used to suppress leading zeros 015F 58 ro20: pop ax ;get digit 0160 0A D8 or bl, al ;turn off suppression if not a zero 0162 80 F9 03 cmp cl, after+1 ;is digit immediately to left of decimal point? 0165 75 01 jne ro30 ;skip if not 0167 43 inc bx ;turn off leading zero suppression 0168 04 30 ro30: add al, '0' ;if leading zero then ' ' else add 0 016A 84 DB test bl, bl 016C 75 02 jne ro40 016E B0 20 mov al, ' ' 0170 CD 29 ro40: int 29h ;display character in al register 0172 80 F9 03 cmp cl, after+1 ;is digit immediately to left of decimal point? 0175 75 04 jne ro50 ;skip if not 0177 B0 2E mov al, '.' ;display decimal point 0179 CD 29 int 29h 017B E2 E2 ro50: loop ro20 ;loop until all digits displayed 017D C3 ret ;return to OS

017E Haversine: ;return (1.0-Cos(Ang)) / 2.0 in st 017E D9 FF fcos 0180 D9 E8 fld1 0182 DE E1 fsubr 0184 D8 36 0189r fdiv N2 0188 C3 ret

0189 40000000 N2 dd 2.0 018D 3F21628D Lat1 dd 0.63041 ;36.12*pi/180 0191 3F17A4E8 Lat2 dd 0.59236 ;33.94*pi/180 0195 BFC19F80 Lon1 dd -1.51268  ;-86.67*pi/180 0199 C004410B Lon2 dd -2.06647  ;-118.40*pi/180 019D 46472666 Radius2 dd 12745.6 ;6372.8 average radius of Earth (km) times 2 

                            ;(TASM isn't smart enough to do floating point constant calculations)
                                    end     start

</lang>

Output:
2887.25

XPL0

<lang XPL0>include c:\cxpl\codes; \intrinsic 'code' declarations

func real Haversine(Ang); real Ang; return (1.0-Cos(Ang)) / 2.0;

func real Dist(Lat1, Lat2, Lon1, Lon2); \Great circle distance real Lat1, Lat2, Lon1, Lon2; def R = 6372.8; \average radius of Earth (km) return 2.0*R * ASin( sqrt( Haversine(Lat2-Lat1) + 

      Cos(Lat1)*Cos(Lat2)*Haversine(Lon2-Lon1) ));

def D2R = 3.141592654/180.0; \degrees to radians RlOut(0, Dist(36.12*D2R, 33.94*D2R, -86.67*D2R, -118.40*D2R ));</lang>

Output:
 2887.25995

zkl

Translation of: Erlang

<lang zkl>haversine(36.12, -86.67, 33.94, -118.40).println();

fcn haversine(Lat1, Long1, Lat2, Long2){ 

  const R = 6372.8; 	// In kilometers;
  Diff_Lat  := (Lat2  - Lat1) .toRad();
  Diff_Long := (Long2 - Long1).toRad();
  NLat      := Lat1.toRad();
  NLong     := Lat2.toRad();
  A 	     := (Diff_Lat/2) .sin().pow(2) + 
               (Diff_Long/2).sin().pow(2) *

NLat.cos() * NLong.cos(); 

  C 	     := 2.0 * A.sqrt().asin();
  R*C;

}</lang>

Output:
2887.26
Retrieved from "https://rosettacode.org/wiki/Haversine_formula?oldid=208805"
Cookies help us deliver our services. By using our services, you agree to our use of cookies.
More information