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,   which is:
   N 36°7.2',   W 86°40.2'     (36.12,   -86.67)           -and-
  • Los Angeles International Airport (LAX)  in Los Angeles, CA, USA,   which is:
   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). 

Most of the examples below adopted Kaimbridge's recommended value of 6372.8 km for the earth radius. However, the derivation of this ellipsoidal quadratic mean radius is wrong (the averaging over azimuth is biased). When applying these examples in real applications, it is better to use the mean earth radius, 6371 km. This value is recommended by the International Union of Geodesy and Geophysics and it minimizes the RMS relative error between the great circle and geodesic distance.


11l

Translation of: Python
F haversine(=lat1, lon1, =lat2, lon2)
   V r = 6372.8
   V dLat = radians(lat2 - lat1)
   V dLon = radians(lon2 - lon1)
   lat1 = radians(lat1)
   lat2 = radians(lat2)
   V a = sin(dLat / 2) ^ 2 + cos(lat1) * cos(lat2) * sin(dLon / 2) ^ 2
   V c = 2 * asin(sqrt(a))
   R r * c

print(haversine(36.12, -86.67, 33.94, -118.40))
Output:
2887.26

ABAP

  DATA: X1 TYPE F, Y1 TYPE F,
        X2 TYPE F, Y2 TYPE F, YD TYPE F,
        PI TYPE F,
        PI_180 TYPE F,
        MINUS_1 TYPE F VALUE '-1'.

PI     = ACOS( MINUS_1 ).
PI_180 = PI / 180.

LATITUDE1 = 36,12 . LONGITUDE1 = -86,67 .
LATITUDE2 = 33,94 . LONGITUDE2 = -118,4 .

  X1 = LATITUDE1  * PI_180.
  Y1 = LONGITUDE1 * PI_180.
  X2 = LATITUDE2  * PI_180.
  Y2 = LONGITUDE2 * PI_180.
  YD = Y2 - Y1.

  DISTANCE = 20000 / PI *
    ACOS( SIN( X1 ) * SIN( X2 ) + COS( X1 ) * COS( X2 ) * COS( YD ) ).

WRITE : 'Distance between given points = ' , distance , 'km .' .
Output:
Distance between given points = 2.884,2687 km . 

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;

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
#!/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))
)
Output:
dist: 2887.3 km (1794.1 mi.)

ALGOL W

Translation of: ALGOL 68

Using the mean radius value suggested in the task.

begin % compute the distance between places using the Haversine formula %
    real procedure arcsin( real value x ) ; arctan( x / sqrt( 1 - ( x * x ) ) );
    real procedure distance ( real value th1Deg, ph1Deg, th2Deg, ph2Deg ) ;
    begin
        real ph1, th1, th2, toRad, dz, dx, dy;
        toRad := pi / 180;
        ph1   := ( ph1Deg - ph2Deg ) * toRad;
        th1   := th1Deg * toRad;
        th2   := th2Deg * toRad;
        dz    := sin( th1 ) - sin( th2 );
        dx    := cos( ph1 ) * cos( th1 ) - cos( th2 );
        dy    := sin( ph1 ) * cos( th1 );
        arcsin( sqrt( dx * dx + dy * dy + dz * dz ) / 2 ) * 2 * 6371
    end distance ;
    begin
        real d;
        integer mi, km;
        d  := distance( 36.12, -86.67, 33.94, -118.4 );
        mi := round( d );
        km := round( d / 1.609344 );
        writeon( i_w := 4, s_w := 0, "distance: ", mi, " km (", km, " mi.)" )
    end
end.
Output:
distance: 2886 km (1794 mi.)

Amazing Hopper

Translation of: Python
/* fórmula de Haversine para distancias en una
   superficie esférica */

#include <basico.h>

#define MIN      60
#define SEG      3600
#define RADIO     6372.8
#define UNAMILLA 1.609344

algoritmo
    
    números ( lat1, lon1, lat2, lon2, dlat, dlon, millas )
    
    si ' total argumentos es (9) '  // LAT1 M LON1 M LAT2 M LON2 M
        #basic{ lat1 = narg(2) + narg(3)/MIN 
                lon1 = narg(4) + narg(5)/MIN 
                lat2 = narg(6) + narg(7)/MIN 
                lon2 = narg(8) + narg(9)/MIN }
        
    sino si ' total argumentos es (13) ' // LAT1 M LON1 M S LAT2 M LON2 M S
        #basic{ lat1 = narg(2) + narg(3)/MIN + narg(4)/SEG 
                lon1 = narg(5) + narg(6)/MIN + narg(7)/SEG  
                lat2 = narg(8) + narg(9)/MIN + narg(10)/SEG  
                lon2 = narg(11) + narg(12)/MIN + narg(13)/SEG  }
    sino
        imprimir("Modo de uso:\ndist.bas La1 M [S] Lo1 M [S] La2 M [S] Lo2 M [S]\n")
        término prematuro
    fin si

    #basic{
        dlat = sin(radian(lat2 - lat1)/2)^2
        dlon = sin(radian(lon2 - lon1)/2)^2
        RADIO*(2*arc sin(sqrt(dlat + cos(radian(lat1)) * cos(radian(lat2)) * dlon )))
    }
    
    ---copiar en 'millas'---, " km. (",millas,dividido por (UNA MILLA)," mi.)\n"
    decimales '2', imprimir
    
terminar
Output:
$ hopper3 basica/dist.bas -x -o bin/dist
Generating binary ‘bin/dist’...Ok! 
Symbols: 27
Total size: 0.75 Kb
$ ./bin/dist 36 7.2 86 40.2 33 56.4  118 24
2887.26 km. (1794.06 mi.)

AMPL

set location;
set geo;

param coord{i in location, j in geo};
param dist{i in location, j in location};

data;

set location := BNA LAX;
set geo := LAT LON;

param coord:
               LAT      LON :=
      BNA    36.12   -86.67
      LAX    33.94   -118.4
;

let dist['BNA','LAX'] := 2 * 6372.8 * asin (sqrt(sin(atan(1)/45*(coord['LAX','LAT']-coord['BNA','LAT'])/2)^2 + cos(atan(1)/45*coord['BNA','LAT']) * cos(atan(1)/45*coord['LAX','LAT']) * sin(atan(1)/45*(coord['LAX','LON'] - coord
['BNA','LON'])/2)^2));

printf "The distance between the two points is approximately %f km.\n", dist['BNA','LAX'];
Output:
The distance between the two points is approximately 2887.259951 km.

APL

r6371
hf{(p q) ÷180  2×rׯ1(+/(2*1(p-q)÷2)×1(×/2○⊃¨p q))2}
36.12 ¯86.67 hf 33.94 ¯118.40
Output:
2886.44

AppleScript

AppleScript provides no trigonometric functions.

Here we reach through a foreign function interface to a temporary instance of a JavaScript interpreter.

use AppleScript version "2.4" -- Yosemite (10.10) or later
use framework "Foundation"
use framework "JavaScriptCore"
use scripting additions

property js : missing value


-- haversine :: (Num, Num) -> (Num, Num) -> Num
on haversine(latLong, latLong2)
    set {lat, lon} to latLong
    set {lat2, lon2} to latLong2
    
    set {rlat1, rlat2, rlon1, rlon2} to ¬
        map(my radians, {lat, lat2, lon, lon2})
    
    set dLat to rlat2 - rlat1
    set dLon to rlon2 - rlon1
    set radius to 6372.8 -- km
    
    set asin to math("asin")
    set sin to math("sin")
    set cos to math("cos")
    
    |round|((2 * radius * ¬
        (asin's |λ|((sqrt(((sin's |λ|(dLat / 2)) ^ 2) + ¬
            (((sin's |λ|(dLon / 2)) ^ 2) * ¬
                (cos's |λ|(rlat1)) * (cos's |λ|(rlat2)))))))) * 100) / 100
end haversine


-- math :: String -> Num -> Num
on math(f)
    script
        on |λ|(x)
            if missing value is js then ¬
                set js to current application's JSContext's new()
            (js's evaluateScript:("Math." & f & "(" & x & ")"))'s toDouble()
        end |λ|
    end script
end math


-------------------------- TEST ---------------------------
on run
    set distance to haversine({36.12, -86.67}, {33.94, -118.4})
    
    set js to missing value -- Clearing a c pointer.
    return distance
end run


-------------------- GENERIC FUNCTIONS --------------------

-- map :: (a -> b) -> [a] -> [b]
on map(f, xs)
    -- The list obtained by applying f
    -- to each element of xs.
    tell mReturn(f)
        set lng to length of xs
        set lst to {}
        repeat with i from 1 to lng
            set end of lst to |λ|(item i of xs, i, xs)
        end repeat
        return lst
    end tell
end map


-- mReturn :: First-class m => (a -> b) -> m (a -> b)
on mReturn(f)
    -- 2nd class handler function lifted into 1st class script wrapper. 
    if script is class of f then
        f
    else
        script
            property |λ| : f
        end script
    end if
end mReturn


-- radians :: Float x => Degrees x -> Radians x
on radians(x)
    (pi / 180) * x
end radians


-- round :: a -> Int
on |round|(n)
    round n
end |round|


-- sqrt :: Num -> Num
on sqrt(n)
    if n  0 then
        n ^ (1 / 2)
    else
        missing value
    end if
end sqrt
Output:
2887.26

Arturo

radians: function [x]-> x * pi // 180

haversine: function [src,tgt][
    dLat: radians tgt\0 - src\0
    dLon: radians tgt\1 - src\1
    lat1: radians src\0
    lat2: radians tgt\0

    a: add product @[cos lat1, cos lat2, sin dLon/2, sin dLon/2] (sin dLat/2) ^ 2
    c: 2 * asin sqrt a
    return 6372.8 * c
]

print haversine @[36.12 neg 86.67] @[33.94, neg 118.40]
Output:
2887.259950607111

ATS

#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"

#define R 6372.8
#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]
Output:
dist: 2887.3 km (1794.1 mi.)

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
}
Output:
2887.259951 km

AWK

# syntax: GAWK -f HAVERSINE_FORMULA.AWK
# 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.)
}
Output:
distance: 2887.2599 km

BASIC

Applesoft BASIC

Works with: Chipmunk Basic
Works with: GW-BASIC
Works with: Minimal BASIC
Works with: MSX BASIC
Translation of: Commodore BASIC
100 HOME : rem  100 CLS for GW-BASIC and MSX BASIC : DELETE for Minimal BASIC
110 LET P = ATN(1)*4
120 LET D = P/180
130 LET M = 36.12
140 LET K = -86.67
150 LET N = 33.94
160 LET L = -118.4
170 LET R = 6372.8
180 PRINT " DISTANCIA DE HAVERSINE ENTRE BNA Y LAX = ";
190 LET A = SIN((L-K)*D/2)
200 LET A = A*A
210 LET B = COS(M*D)*COS(N*D)
220 LET C = SIN((N-M)*D/2)
230 LET C = C*C
240 LET D = SQR(C+B*A)
250 LET E = D/SQR(1-D*D)
260 LET F = ATN(E)
270 PRINT 2*R*F;"KM"
280 END

BASIC256

global radioTierra    # radio de la tierra en km
radioTierra = 6372.8

function Haversine(lat1, long1, lat2, long2 , radio)
	d_long = radians(long1 - long2)
	theta1 = radians(lat1)
	theta2 = radians(lat2)
	dx = cos(d_long) * cos(theta1) - cos(theta2)
	dy = sin(d_long) * cos(theta1)
	dz = sin(theta1) - sin(theta2)
	return asin(sqr(dx*dx + dy*dy + dz*dz) / 2) * radio * 2
end function

print
print " Distancia de Haversine entre BNA y LAX = ";
print Haversine(36.12, -86.67, 33.94, -118.4, radioTierra); " km"
end
Output:
 Distancia de Haversine entre BNA y LAX = 2887.25994877 km.

Chipmunk Basic

Works with: Chipmunk Basic version 3.6.4
100 cls
110 pi = arctan(1)*4 : rem define pi = 3.1415...
120 deg2rad = pi/180 : rem define grados a radianes 0.01745..
130 lat1 = 36.12
140 long1 = -86.67
150 lat2 = 33.94
160 long2 = -118.4
170 radio = 6372.8
180 print " Distancia de Haversine entre BNA y LAX = ";
190 d_long = deg2rad*(long1-long2)
200 theta1 = deg2rad*(lat1)
210 theta2 = deg2rad*(lat2)
220 dx = cos(d_long)*cos(theta1)-cos(theta2)
230 dy = sin(d_long)*cos(theta1)
240 dz = sin(theta1)-sin(theta2)
250 print (asin(sqr(dx*dx+dy*dy+dz*dz)/2)*radio*2);"km"
260 end

Gambas

Public deg2rad As Float = Pi / 180    ' define grados a radianes 0.01745..
Public radioTierra As Float = 6372.8  ' radio de la tierra en km

Public Sub Main() 
  
  Print "\n Distancia de Haversine entre BNA y LAX = "; Haversine(36.12, -86.67, 33.94, -118.4, radioTierra); " km"
  
End

Function Haversine(lat1 As Float, long1 As Float, lat2 As Float, long2 As Float, radio As Float) As Float 
  
  Dim d_long As Float = deg2rad * (long1 - long2) 
  Dim theta1 As Float = deg2rad * lat1 
  Dim theta2 As Float = deg2rad * lat2 
  Dim dx As Float = Cos(d_long) * Cos(theta1) - Cos(theta2) 
  Dim dy As Float = Sin(d_long) * Cos(theta1) 
  Dim dz As Float = Sin(theta1) - Sin(theta2) 

  Return ASin(Sqr(dx * dx + dy * dy + dz * dz) / 2) * radio * 2 

End Function
Output:
Same as FreeBASIC entry.

GW-BASIC

Works with: Applesoft BASIC
Works with: BASICA
Works with: Chipmunk Basic
Works with: MSX BASIC
Works with: Minimal BASIC
Works with: PC-BASIC version any
Translation of: Commodore BASIC
100 CLS : rem  100 HOME for Applesoft BASIC : DELETE for Minimal BASIC
110 LET P = ATN(1)*4
120 LET D = P/180
130 LET M = 36.12
140 LET K = -86.67
150 LET N = 33.94
160 LET L = -118.4
170 LET R = 6372.8
180 PRINT " DISTANCIA DE HAVERSINE ENTRE BNA Y LAX = ";
190 LET A = SIN((L-K)*D/2)
200 LET A = A*A
210 LET B = COS(M*D)*COS(N*D)
220 LET C = SIN((N-M)*D/2)
230 LET C = C*C
240 LET D = SQR(C+B*A)
250 LET E = D/SQR(1-D*D)
260 LET F = ATN(E)
270 PRINT 2*R*F;"KM"
280 END

Minimal BASIC

Works with: Applesoft BASIC
Works with: Chipmunk Basic
Works with: GW-BASIC
Works with: MSX BASIC
Translation of: Commodore BASIC
110 LET P = ATN(1)*4
120 LET D = P/180
130 LET M = 36.12
140 LET K = -86.67
150 LET N = 33.94
160 LET L = -118.4
170 LET R = 6372.8
180 PRINT " DISTANCIA DE HAVERSINE ENTRE BNA Y LAX = ";
190 LET A = SIN((L-K)*D/2)
200 LET A = A*A
210 LET B = COS(M*D)*COS(N*D)
220 LET C = SIN((N-M)*D/2)
230 LET C = C*C
240 LET D = SQR(C+B*A)
250 LET E = D/SQR(1-D*D)
260 LET F = ATN(E)
270 PRINT 2*R*F;"KM"
280 END

MSX Basic

The GW-BASIC solution works without any changes.

QBasic

Works with: QBasic version 1.1
Works with: QuickBasic version 4.5
CONST pi = 3.141593     ' define pi
CONST radio = 6372.8    ' radio de la tierra en km

PRINT : PRINT " Distancia de Haversine:";
PRINT Haversine!(36.12, -86.67, 33.94, -118.4); "km"
END

FUNCTION Haversine! (lat1!, long1!, lat2!, long2!)
    deg2rad! = pi / 180    ' define grados a radianes 0.01745..

    dLong! = deg2rad! * (long1! - long2!)
    theta1! = deg2rad! * lat1!
    theta2! = deg2rad! * lat2!
    dx! = COS(dLong!) * COS(theta1!) - COS(theta2!)
    dy! = SIN(dLong!) * COS(theta1!)
    dz! = SIN(theta1!) - SIN(theta2!)
   
    Haversine! = (SQR(dx! * dx! + dy! * dy! + dz! * dz!) / 2) * radio * 2
END FUNCTION
Output:
 Distancia de Haversine: 2862.63 km

Quite BASIC

100 CLS
110 LET p = atan(1)*4
120 LET d = p/180
130 LET k = 36.12
140 LET m = -86.67
150 LET l = 33.94
160 LET n = -118.4
170 LET r = 6372.8
180 PRINT " Distancia de Haversine entre BNA y LAX = ";
190 LET g = d*(m-n)
200 LET t = d*(k)
210 LET s = d*(l)
220 LET x = COS(g)*COS(t)-COS(s)
230 LET y = SIN(g)*COS(t)
240 LET z = SIN(t)-SIN(s)
250 PRINT (ASIN(SQR(x*x+y*y+z*z)/2)*r*2);"km"
260 END

True BASIC

DEF Haversine (lat1, long1, lat2, long2)
    OPTION ANGLE RADIANS
    LET R = 6372.8                !radio terrestre en km.
    LET dLat = RAD(lat2-lat1)
    LET dLong = RAD(long2-long1)
    LET lat1 = RAD(lat1)
    LET lat2 = RAD(lat2)
    LET Haversine = R *2 * ASIN(SQR(SIN(dLat/2)^2 + SIN(dLong/2)^2 *COS(lat1) * COS(lat2)))
END DEF
PRINT

PRINT "Distancia de Haversine:"; Haversine(36.12, -86.67, 33.94, -118.4); "km"
END
Output:
Distancia de Haversine: 2887.26 km

Yabasic

Translation of: FreeBASIC
//pi está predefinido en Yabasic
deg2rad = pi / 180    // define grados a radianes 0.01745..
radioTierra = 6372.8  // radio de la tierra en km

sub Haversine(lat1, long1, lat2, long2 , radio)   
    d_long = deg2rad * (long1 - long2)
    theta1 = deg2rad * lat1
    theta2 = deg2rad * lat2
    dx = cos(d_long) * cos(theta1) - cos(theta2)
    dy = sin(d_long) * cos(theta1)
    dz = sin(theta1) - sin(theta2)
    return asin(sqr(dx*dx + dy*dy + dz*dz) / 2) * radio * 2
end sub

print " Distancia de Haversine entre BNA y LAX = ", Haversine(36.12, -86.67, 33.94, -118.4, radioTierra), " km"
end
Output:
 Distancia de Haversine entre BNA y LAX = 259.478 km

BBC BASIC

Uses BBC BASIC's MOD(array()) function which calculates the square-root of the sum of the squares of the elements of an array.

      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
Output:
Distance = 2887.25995 km

bc

Works with: GNU bc(1) version 1.07.1-2
Works with: dc(1)-based MirBSD bc(1) version as of 2021-06-11
Works with: bc(1) version POSIX

… supplied with a small POSIX shell wrapper to feed the arguments to bc. (see also)

#!/bin/sh
#-
# © 2021 mirabilos Ⓕ CC0; implementation of Haversine GCD from public sources
#
# now developed online:
# https://evolvis.org/plugins/scmgit/cgi-bin/gitweb.cgi?p=useful-scripts/mirkarte.git;a=blob;f=geo.sh;hb=HEAD

if test "$#" -ne 4; then
	echo >&2 "E: syntax: $0 lat1 lon1 lat2 lon2"
	exit 1
fi

set -e

# make GNU bc use POSIX mode and shut up
BC_ENV_ARGS=-qs
export BC_ENV_ARGS

# assignment of constants, variables and functions
# p: multiply with to convert from degrees to radians (π/180)
# r: earth radius in metres
# d: distance
# h: haversine intermediate
# i,j: (lat,lon) point 1
# x,y: (lat,lon) point 2
# k: delta lat
# l: delta lon
# m: sin(k/2) (square root of hav(k))
# n: sin(l/2) (  partial haversine  )
# n(x): arcsin(x)
# r(x,n): round x to n decimal digits
# v(x): sign (Vorzeichen)
# w(x): min(1, sqrt(x)) (Wurzel)

bc -l <<-EOF
scale=64
define n(x) {
	if (x == -1) return (-2 * a(1))
	if (x == 1) return (2 * a(1))
	return (a(x / sqrt(1 - x*x)))
}
define v(x) {
	if (x < 0) return (-1)
	if (x > 0) return (1)
	return (0)
}
define r(x, n) {
	auto o
	o = scale
	if (scale < (n + 1)) scale = (n + 1)
	x += v(x) * 0.5 * A^-n
	scale = n
	x /= 1
	scale = o
	return (x)
}
define w(x) {
	if (x >= 1) return (1)
	return (sqrt(x))
}
/* WGS84 reference ellipsoid: große Halbachse (metres), Abplattung */
i = 6378137.000
x = 1/298.257223563
/* other axis */
j = i * (1 - x)
/* mean radius resulting */
r = (2 * i + j) / 3
/* coordinates */
p = (4 * a(1) / 180)
i = (p * $1)
j = (p * $2)
x = (p * $3)
y = (p * $4)
/* calculation */
k = (x - i)
l = (y - j)
m = s(k / 2)
n = s(l / 2)
h = ((m * m) + (c(i) * c(x) * n * n))
d = 2 * r * n(w(h))
r(d, 3)
EOF

# output is in metres, rounded to millimetres, error < ¼% in WGS84
Output:
$ sh dist.sh 36.12 -86.67 33.94 -118.4
2886448.430

Note I used a more precise earth radius; this matches the other implementations when choosing the same.

C

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define R 6371
#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;
}

C#

Translation of: Groovy
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

C++

#define _USE_MATH_DEFINES

#include <math.h>
#include <iostream>

const static double EarthRadiusKm = 6372.8;

inline double DegreeToRadian(double angle)
{
	return M_PI * angle / 180.0;
}

class Coordinate
{
public:
	Coordinate(double latitude ,double longitude):myLatitude(latitude), myLongitude(longitude)
	{}

	double Latitude() const
	{
		return myLatitude;
	}

	double Longitude() const
	{
		return myLongitude;
	}

private:

	double myLatitude;
	double myLongitude;
};

double HaversineDistance(const Coordinate& p1, const Coordinate& p2)
{
	double latRad1 = DegreeToRadian(p1.Latitude());
	double latRad2 = DegreeToRadian(p2.Latitude());
	double lonRad1 = DegreeToRadian(p1.Longitude());
	double lonRad2 = DegreeToRadian(p2.Longitude());

	double diffLa = latRad2 - latRad1;
	double doffLo = lonRad2 - lonRad1;

	double computation = asin(sqrt(sin(diffLa / 2) * sin(diffLa / 2) + cos(latRad1) * cos(latRad2) * sin(doffLo / 2) * sin(doffLo / 2)));
	return 2 * EarthRadiusKm * computation;
}

int main()
{
	Coordinate c1(36.12, -86.67);
	Coordinate c2(33.94, -118.4);

	std::cout << "Distance = " << HaversineDistance(c1, c2) << std::endl;
	return 0;
}

clojure

Translation of: Java
(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

CoffeeScript

Translation of: JavaScript
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)
Output:
2887.2599506071124

Commodore BASIC

PETSCII has the pi symbol π in place of the ASCII tilde ~; Commodore BASIC interprets this symbol as the mathematical constant.

10 REM================================
15 REM HAVERSINE FORMULA
20 REM
25 REM 2021-09-24
30 REM EN.WIKIPEDIA.ORG/WIKI/HAVERSINE_FORMULA
35 REM
40 REM C64 HAS PI CONSTANT
45 REM X1 LONGITUDE 1
50 REM Y1 LATITUDE 1
55 REM X2 LONGITUDE 2
60 REM Y2 LATITUDE 2
65 REM
70 REM V1, 2021-10-02, ALVALONGO
75 REM ===============================
100 REM MAIN
105 DR=π/180:REM DEGREES TO RADIANS
110 PRINT CHR$(147);CHR$(5);"HAVERSINE FORMULA"
120 PRINT "GREAT-CIRCLE DISTANCE"
130 R=6372.8:REM AVERAGE EARTH RADIUS IN KILOMETERS
200 REM GET DATA
210 PRINT
220 INPUT "LONGITUDE 1=";X1
230 INPUT "LATITUDE  1=";Y1
240 PRINT
250 INPUT "LONGITUDE 2=";X2
260 INPUT "LATITUDE  2=";Y2
270 GOSUB 500
280 PRINT
290 PRINT "DISTANCE=";D;"KM"
300 GET K$:IF K$="" THEN 300
310 GOTO 210
490 END
500 REM HAVERSINE FORMULA ------------
520 A=SIN((X2-X1)*DR/2)
530 A=A*A
540 B=COS(Y1*DR)*COS(Y2*DR)
550 C=SIN((Y2-Y1)*DR/2)
560 C=C*C
570 D=SQR(C+B*A)
580 E=D/SQR(1-D*D)
590 F=ATN(E)
600 D=2*R*F
610 RETURN
Output:
HAVERSINE FORMULA
GREAT-CIRCLE DISTANCE

LONGITUDE 1=? -86.67
LATITUDE  1=? 36.12

LONGITUDE 2=? -118.40
LATITUDE  2=? 33.94

DISTANCE= 2887.25995 KM

Common 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)))
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.

Crystal

Translation of: Python
include Math
 
def haversine(lat1, lon1, lat2, lon2)
    r = 6372.8        # Earth radius in kilometers
    deg2rad = PI/180  # convert degress to radians
 
    dLat = (lat2 - lat1) * deg2rad
    dLon = (lon2 - lon1) * deg2rad
    lat1 = lat1 * deg2rad
    lat2 = lat2 * deg2rad
 
    a = sin(dLat / 2)**2 + cos(lat1) * cos(lat2) * sin(dLon / 2)**2
    c = 2 * asin(sqrt(a))
    r * c
end

puts "distance is #{haversine(36.12, -86.67, 33.94, -118.40)} km "
Output:
distance is 2887.2599506071106 km 

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));
}
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.

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));
}
Output:
Great circle distance: 2887.3 km

Dart

Translation of: Java
import 'dart:math';

class Haversine {
  static final R = 6372.8; // In kilometers

  static double haversine(double lat1, lon1, lat2, lon2) {
    double dLat = _toRadians(lat2 - lat1);
    double dLon = _toRadians(lon2 - lon1);
    lat1 = _toRadians(lat1);
    lat2 = _toRadians(lat2);
    double a = pow(sin(dLat / 2), 2) + pow(sin(dLon / 2), 2) * cos(lat1) * cos(lat2);
    double c = 2 * asin(sqrt(a));
    return R * c;
  }

  static double _toRadians(double degree) {
    return degree * pi / 180;
  }

  static void main() {
    print(haversine(36.12, -86.67, 33.94, -118.40));
  }
}
Output:
2887.2599506071106

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.
Output:
Haversine distance: 2887.26 km.

EasyLang

Translation of: C
func dist th1 ph1 th2 ph2 .
   r = 6371
   ph1 -= ph2
   dz = sin th1 - sin th2
   dx = cos ph1 * cos th1 - cos th2
   dy = sin ph1 * cos th1
   return 2 * r * pi / 180 * asin (sqrt (dx * dx + dy * dy + dz * dz) / 2)
.
print dist 36.12 -86.67 33.94 -118.4

Elena

ELENA 4.x:

import extensions;
import system'math;
 
Haversine(lat1,lon1,lat2,lon2)
{
    var R := 6372.8r;
    var dLat := (lat2 - lat1).Radian;
    var dLon := (lon2 - lon1).Radian;
 
    var dLat1 := lat1.Radian;
    var dLat2 := lat2.Radian;
 
    var a := (dLat / 2).sin() * (dLat / 2).sin() + (dLon / 2).sin() * (dLon / 2).sin() * dLat1.cos() * dLat2.cos();
 
    ^ R * 2 * a.sqrt().arcsin()
}
 
public program()
{
    console.printLineFormatted("The distance between coordinates {0},{1} and {2},{3} is: {4}", 36.12r, -86.67r, 33.94r, -118.40r, 
        Haversine(36.12r, -86.67r, 33.94r, -118.40r))
}
Output:
The distance between coordinates 36.12,-86.67 and 33.94,-118.4 is: 2887.259950607

Elixir

defmodule Haversine do
  @v  :math.pi / 180
  @r  6372.8            # km for the earth radius
  def distance({lat1, long1}, {lat2, long2}) do
    dlat  = :math.sin((lat2 - lat1) * @v / 2)
    dlong = :math.sin((long2 - long1) * @v / 2)
    a = dlat * dlat + dlong * dlong * :math.cos(lat1 * @v) * :math.cos(lat2 * @v)
    @r * 2 * :math.asin(:math.sqrt(a))
  end
end

bna = {36.12,  -86.67}
lax = {33.94, -118.40}
IO.puts Haversine.distance(bna, lax)
Output:
2887.2599506071106

Elm

haversine : ( Float, Float ) -> ( Float, Float ) -> Float
haversine ( lat1, lon1 ) ( lat2, lon2 ) =
    let
        r =
            6372.8

        dLat =
            degrees (lat2 - lat1)

        dLon =
            degrees (lon2 - lon1)

        a =
            (sin (dLat / 2))
                ^ 2
                + (sin (dLon / 2))
                ^ 2
                * cos (degrees lat1)
                * cos (degrees lat2)
    in
        r * 2 * asin (sqrt a)

view =
    Html.div []
      [ Html.text (toString (haversine ( 36.12, -86.67 ) ( 33.94, -118.4 )))
      ]
Output:
2887.2599506071106

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.
Output:
2887.2599506071106

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

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

Excel

LAMBDA

Binding the name HAVERSINE to the following lambda expression in the Name Manager of the Excel workbook:

(See LAMBDA: The ultimate Excel worksheet function)

HAVERSINE
=LAMBDA(lla,
    LAMBDA(llb,
        LET(
            REM, "Approximate radius of Earth in km.",
            earthRadius, 6372.8,
            
            sinHalfDeltaSquared, LAMBDA(x, SIN(x / 2) ^ 2)(
                RADIANS(llb - lla)
            ),
    
            2 * earthRadius * ASIN(
                SQRT(
                    INDEX(sinHalfDeltaSquared, 1) + (
                        PRODUCT(COS(RADIANS(
                            CHOOSE({1,2}, 
                                INDEX(lla, 1), 
                                INDEX(llb, 1)
                            )
                        )))
                    ) * INDEX(sinHalfDeltaSquared, 2)
                )
            )
        )
    )
)

Each of the two arguments in the example below is an Excel dynamic array of two adjacent values. The # character yields a reference to the array with the given top-left grid address.

Cell B2 is formatted to display only two decimal places.

Output:
fx =HAVERSINE(E2#)(H2#)
A B C D E F G H I
1 Distance BNA LAX
2 2887.26 km 36.12 -86.67 33.94 -118.4

F#

Translation of: Go
using units of measure
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
Output:
2887.259951

Factor

Translation of: J
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 * ;
( scratchpad ) { 36.12 -86.67 } { 33.94 -118.4 } haversineDist .
2887.259950607113

FBSL

Based on the Fortran and Groovy versions.

#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
Output:
Distance = 2887.25995060711 km
Press any key to continue...

FOCAL

Translation of: ZX Spectrum Basic
1.01 S BA = 36.12; S LA = -86.67
1.02 S BB = 33.94; S LB = -118.4
1.03 S DR = 3.1415926536 / 180; S D = 2 * 6372.8
1.04 S TA = (LB - LA) * DR
1.05 S TB = DR * BA
1.06 S TC = DR * BB
1.07 S DZ = FSIN(TB) - FSIN(TC)
1.08 S DX = FCOS(TA) * FCOS(TB) - FCOS(TC)
1.09 S DY = FSIN(TA) * FCOS(TB)
1.10 S AS = DX * DX + DY * DY + DZ * DZ
1.11 S AS = FSQT(AS) / 2
1.12 S HDIST = D * FATN(AS / FSQT(1 - AS^2))
1.13 T %6.2,"Haversine distance ",HDIST,!
Output:
Haversine distance = 2887.26

Note that FOCAL lacks a built-in arcsine function, but appendix D of the FOCAL manual shows how to compute it using arctangent and square root instead.

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.
Output:
2887.25995060711

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

Free Pascal

Here is a Free Pascal version, works in most Pascal dialects, but also note the Delphi entry that also works in Free Pascal.

program HaversineDemo;
uses 
  Math;

function HaversineDistance(const lat1, lon1, lat2, lon2:double):double;inline;
const 
  rads = pi / 180;
  dia  = 2 * 6372.8;
begin
  HaversineDistance := dia * arcsin(sqrt(sqr(cos(rads * (lon1 - lon2)) * cos(rads * lat1)  
                     - cos(rads * lat2)) + sqr(sin(rads * (lon1 - lon2)) 
                     * cos(rads * lat1)) + sqr(sin(rads * lat1) - sin(rads * lat2))) / 2);
end;

begin
  Writeln('Haversine distance between BNA and LAX: ', HaversineDistance(36.12, -86.67, 33.94, -118.4):7:2, ' km.');
end.

FreeBASIC

' version 09-10-2016
' compile with: fbc -s console

' Nashville International Airport (BNA) in Nashville, TN, USA,
' N 36°07.2',  W  86°40.2' (36.12,  -86.67)
' Los Angeles International Airport (LAX) in Los Angeles, CA, USA,
' N 33°56.4', W 118°24.0'  (33.94, -118.40).
' 6372.8 km is an approximation of the radius of the average circumference

#Define Pi Atn(1) * 4        ' define Pi = 3.1415..
#Define deg2rad Pi / 180     ' define deg to rad 0.01745..
#Define earth_radius 6372.8  ' earth radius in km.

Function Haversine(lat1 As Double, long1 As Double, lat2 As Double, _
                                long2 As Double , radius As Double) As Double

  Dim As Double d_long = deg2rad * (long1 - long2)
  Dim As Double theta1 = deg2rad * lat1
  Dim As Double theta2 = deg2rad * lat2
  Dim As Double dx = Cos(d_long) * Cos(theta1) - Cos(theta2)
  Dim As Double dy = Sin(d_long) * Cos(theta1)
  Dim As Double dz = Sin(theta1) - Sin(theta2)
  Return Asin(Sqr(dx*dx + dy*dy + dz*dz) / 2) * radius * 2

End Function

Print
Print " Haversine distance between BNA and LAX = "; _
      Haversine(36.12, -86.67, 33.94, -118.4, earth_radius); " km."


' empty keyboard buffer
While Inkey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End
Output:
 Haversine distance between BNA and LAX =  2887.259950607111 km.

Frink

Frink has built-in constants for the radius of the earth, whether it is the mean radius earthradius, the equatorial radius earthradius_equatorial, or the polar radius earthradius_polar. Below calculates the distance between the points using the haversine formula on a sphere using the mean radius, but we can do better:

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"]
Output:
2886.4489734366999158 km

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.

The calculations are due to:

"Direct and Inverse Solutions of Geodesics on the Ellipsoid with Application of Nested Equations", T.Vincenty, Survey Review XXII, 176, April 1975. http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf

There is also a slightly higher-accuracy algorithm (closer to nanometers instead of fractions of millimeters LOL):

"Algorithms for geodesics", Charles F. F. Karney, Journal of Geodesy, January 2013, Volume 87, Issue 1, pp 43-55 http://link.springer.com/article/10.1007%2Fs00190-012-0578-z

use navigation.frink

d = earthDistance[36.12 deg North, 86.67 deg West, 33.94 deg North, 118.40 deg West]
println[d-> "km"]
Output:
2892.7769573807044975 km

Which should be treated as the closest-to-right answer for the actual distance on the earth's geoid, based on the WGS84 geoid datum.

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)) )
Output:
2887.259950607111

FutureBasic

Note: The Haversine function returns an approximate theoretical value of the Great Circle Distance between two points because it does not factor the ellipsoidal shape of Earth -- fat in the middle from centrifugal force, and squashed at the ends. Navigators once relied on trigonometric functions like versine (versed sine) where angle A is 1-cos(A), and haversine (half versine) or ( 1-cos(A) ) / 2. Also, the radius of the Earth varies, at least depending on who you talk to. Here's NASA's take on it: http://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html

Since it was trivial, this functions returns the distance in miles and kilometers.

window 1

local fn Haversine( lat1 as double, lon1 as double, lat2 as double, lon2 as double, miles as ^double, kilometers as ^double )
  double deg2rad, dLat, dLon, a, c, earth_radius_miles, earth_radius_kilometers
  
  earth_radius_miles = 3959.0 // Radius of the Earth in miles
  earth_radius_kilometers = 6372.8 // Radius of the Earth in kilometers
  deg2rad = Pi / 180 // Pi is predefined in FutureBasic
  
  dLat = deg2rad * ( lat2  - lat1 )
  dLon = deg2rad * ( lon2 - lon1 )
  a = sin( dLat / 2 ) * sin( dLat / 2 ) + cos( deg2rad * lat1 ) * cos( deg2rad * lat2 ) * sin( dLon / 2 ) * sin( dLon / 2 )
  c = 2 * asin( sqr(a) )
  
  miles.nil# =  earth_radius_miles * c
  kilometers.nil# = earth_radius_kilometers * c
end fn

double miles, kilometers

fn Haversine( 36.12, -86.67, 33.94, -118.4, @miles, @kilometers )

print "Distance in miles between BNA and LAX: "; using "####.####"; miles; " miles."
print "Distance in kilometers between BNA LAX: "; using "####.####"; kilometers; " km."

HandleEvents

Output:

Distance in miles between BNA and LAX: 1793.6640 miles.
Distance in kilometers between BNA LAX: 2887.2600 km.

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)))
}
Output:
2887.2599506071097

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

Haskell

import Control.Monad (join)
import Data.Bifunctor (bimap)
import Text.Printf (printf)

-------------------- HAVERSINE FORMULA -------------------

-- The haversine of an angle.
haversine :: Float -> Float
haversine = (^ 2) . sin . (/ 2)

-- The approximate distance, in kilometers,
-- between two points on Earth.
-- The latitude and longtitude are assumed to be in degrees.
greatCircleDistance ::
  (Float, Float) ->
  (Float, Float) ->
  Float
greatCircleDistance = distDeg 6371
  where
    distDeg radius p1 p2 =
      distRad
        radius
        (deg2rad p1)
        (deg2rad p2)
    distRad radius (lat1, lng1) (lat2, lng2) =
      (2 * radius)
        * asin
          ( min
              1.0
              ( sqrt $
                  haversine (lat2 - lat1)
                    + ( (cos lat1 * cos lat2)
                          * haversine (lng2 - lng1)
                      )
              )
          )
    deg2rad = join bimap ((/ 180) . (pi *))

--------------------------- TEST -------------------------
main :: IO ()
main =
  printf
    "The distance between BNA and LAX is about %0.f km.\n"
    (greatCircleDistance bna lax)
  where
    bna = (36.12, -86.67)
    lax = (33.94, -118.40)
Output:
The distance between BNA and LAX is about 2886 km.

Icon and Unicon

Translation of: C
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

printf.icn provides formatting

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

Idris

Translation of: Haskell
module Main
 
-- The haversine of an angle.
hsin : Double -> Double
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 : Double -> (Double, Double) -> (Double, Double) -> Double
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 : Double -> (Double, Double) -> (Double, Double) -> Double 
distDeg radius p1 p2 = distRad radius (deg2rad p1) (deg2rad p2)
  where 
        d2r : Double -> Double
        d2r t = t * pi / 180 
        deg2rad (t, u) = (d2r t, d2r u)
 
-- The approximate distance, in kilometers, between two points on Earth.  
-- The latitude and longtitude are assumed to be in degrees.
earthDist : (Double, Double) -> (Double, Double) -> Double
earthDist = distDeg 6372.8

main : IO () 
main = putStrLn $ "The distance between BNA and LAX is about " ++ show (floor dst) ++ " km."
 where 
      bna : (Double, Double)
      bna = (36.12,  -86.67)

      lax : (Double, Double)
      lax = (33.94, -118.40)

      dst : Double
      dst = earthDist bna lax
Output:
The distance between BNA and LAX is about 2887 km.

IS-BASIC

100 PROGRAM "Haversine.bas"
110 PRINT "Haversine distance:";HAVERSINE(36.12,-86.67,33.94,-118.4);"km"
120 DEF HAVERSINE(LAT1,LON1,LAT2,LON2)
130   OPTION ANGLE RADIANS
140   LET R=6372.8
150   LET DLAT=RAD(LAT2-LAT1):LET DLON=RAD(LON2-LON1)
160   LET LAT1=RAD(LAT1):LET LAT2=RAD(LAT2)
170   LET HAVERSINE=R*2*ASIN(SQR(SIN(DLAT/2)^2+SIN(DLON/2)^2*COS(LAT1)*COS(LAT2)))
190 END DEF

J

Solution:

require 'trig'
haversin=: 0.5 * 1 - cos
Rearth=: 6372.8
haversineDist=: Rearth * haversin^:_1@((1 , *&(cos@{.)) +/ .* [: haversin -)&rfd

Note: J derives the inverse haversin ( haversin^:_1 ) from the definition of haversin.

Example Use:

   36.12 _86.67 haversineDist 33.94 _118.4
2887.26

Java

Translation of: Groovy
public class Haversine {
    public static final double R = 6372.8; // In kilometers

    public static double haversine(double lat1, double lon1, double lat2, double lon2) {
        lat1 = Math.toRadians(lat1);
        lat2 = Math.toRadians(lat2);
        double dLat = lat2 - lat1;
        double dLon = Math.toRadians(lon2 - lon1);

        double a = Math.pow(Math.sin(dLat / 2), 2) + Math.pow(Math.sin(dLon / 2), 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));
    }
}
Output:
2887.2599506071106

JavaScript

ES5

Translation of: Java
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));
Output:
2887.2599506071124

ES6

((x, y) => {
    'use strict';

    // haversine :: (Num, Num) -> (Num, Num) -> Num
    const haversine = ([lat1, lon1], [lat2, lon2]) => {
        // Math lib function names
        const [pi, asin, sin, cos, sqrt, pow, round] = [
            'PI', 'asin', 'sin', 'cos', 'sqrt', 'pow', 'round'
        ]
        .map(k => Math[k]),

            // degrees as radians
            [rlat1, rlat2, rlon1, rlon2] = [lat1, lat2, lon1, lon2]
            .map(x => x / 180 * pi),

            dLat = rlat2 - rlat1,
            dLon = rlon2 - rlon1,
            radius = 6372.8; // km

        // km
        return round(
            radius * 2 * asin(
                sqrt(
                    pow(sin(dLat / 2), 2) +
                    pow(sin(dLon / 2), 2) *
                    cos(rlat1) * cos(rlat2)
                )
            ) * 100
        ) / 100;
    };

    // TEST
    return haversine(x, y);

    // --> 2887.26

})([36.12, -86.67], [33.94, -118.40]);
Output:
2887.26

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) ;

Example:

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

Jsish

From Javascript, ES5, except the arguments value is an Array in jsish, not an Object.

/* Haversine formula, in Jsish */
function haversine() {
       var radians = arguments.map(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;
}

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

/*
=!EXPECTSTART!=
haversine(36.12, -86.67, 33.94, -118.40) ==> 2887.259950607112
=!EXPECTEND!=
*/
Output:
prompt$ jsish -u haversineFormula.jsi
[PASS] haversineFormula.jsi

Julia

Works with: Julia version 0.6
haversine(lat1, lon1, lat2, lon2) =
    2 * 6372.8 * asin(sqrt(sind((lat2 - lat1) / 2) ^ 2 +
    cosd(lat1) * cosd(lat2) * sind((lon2 - lon1) / 2) ^ 2))

@show haversine(36.12, -86.67, 33.94, -118.4)
Output:
haversine(36.12, -86.67, 33.94, -118.4) = 2887.2599506071106

Kotlin

Translation of: Groovy

Use Unicode characters.

import java.lang.Math.*

const val R = 6372.8 // in kilometers

fun haversine(lat1: Double, lon1: Double, lat2: Double, lon2: Double): Double {
    val λ1 = toRadians(lat1)
    val λ2 = toRadians(lat2)
    val Δλ = toRadians(lat2 - lat1)
    val Δφ = toRadians(lon2 - lon1)
    return 2 * R * asin(sqrt(pow(sin(Δλ / 2), 2.0) + pow(sin(Δφ / 2), 2.0) * cos(λ1) * cos(λ2)))
}

fun main(args: Array<String>) = println("result: " + haversine(36.12, -86.67, 33.94, -118.40))

Lambdatalk

Translation of: Python
{def haversine
 {def diameter {* 6372.8 2}}
 {def radians {lambda {:a} {* {/ {PI} 180} :a}}}
 {lambda {:lat1 :lon1 :lat2 :lon2}
  {let { {:dLat {radians {- :lat2 :lat1}}}
         {:dLon {radians {- :lon2 :lon1}}}
         {:lat1 {radians :lat1}}
         {:lat2 {radians :lat2}}
       } {* {diameter}
            {asin {sqrt {+ {pow {sin {/ :dLat 2}} 2}
                           {* {cos :lat1}
                              {cos :lat2} 
                              {pow {sin {/ :dLon 2}} 2} }}}}}}}}
-> haversine

{haversine 36.12 -86.67 33.94 -118.40}
-> 2887.2599506071106

or, using 

{def deg2dec
 {lambda {:s :w}
  {let { {:s {if {or {W.equal? :s W}
                     {W.equal? :s S}} then - else +}}
         {:dm {S.replace ° by space in
              {S.replace ' by in :w}}}
         } :s{S.get 0 :dm}.{round {* {/ 100 60} {S.get 1 :dm}}}}}}
-> deg2dec

we can just write

{haversine
 {deg2dec N 36°7.2'}
 {deg2dec W 86°40.2'}
 {deg2dec N 33°56.4'}
 {deg2dec W 118°24.0'}} 
-> 2887.2599506071106

Liberty BASIC

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
Haversine distance: 2887.25995060711  km.

LiveCode

function radians n
    return n * (3.1415926 / 180)
end radians

function haversine lat1, lng1, lat2, lng2
    local radiusEarth 
    local lat3, lng3
    local lat1Rad, lat2Rad, lat3Rad
    local lngRad1, lngRad2, lngRad3
    local haver
    put 6372.8 into radiusEarth
    put (lat2 - lat1) into lat3
    put (lng2 - lng1) into lng3
    put radians(lat1) into lat1Rad
    put radians(lat2) into lat2Rad
    put radians(lat3) into lat3Rad
    put radians(lng1) into lngRad1
    put radians(lng2) into lngRad2
    put radians(lng3) into lngRad3
    
    put (sin(lat3Rad/2.0)^2) + (cos(lat1Rad)) \
          * (cos(lat2Rad)) \
          * (sin(lngRad3/2.0)^2) \
          into haver 
    return (radiusEarth * (2.0 * asin(sqrt(haver))))
    
end haversine

Test

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

Lua

local function haversine(x1, y1, x2, y2)
r=0.017453292519943295769236907684886127;
x1= x1*r; x2= x2*r; y1= y1*r; y2= y2*r; dy = y2-y1; dx = x2-x1;
a = math.pow(math.sin(dx/2),2) + math.cos(x1) * math.cos(x2) * math.pow(math.sin(dy/2),2); c = 2 * math.asin(math.sqrt(a)); d = 6372.8 * c;
return d;
end

Usage:

print(haversine(36.12, -86.67, 33.94, -118.4));

Output:

2887.2599506071

Maple

Inputs assumed to be in radians.

distance := (theta1, phi1, theta2, phi2)->2*6378.14*arcsin( sqrt((1-cos(theta2-theta1))/2 + cos(theta1)*cos(theta2)*(1-cos(phi2-phi1))/2) );
If you prefer, you can define a haversine function to clarify the definition:
haversin := theta->(1-cos(theta))/2;
distance := (theta1, phi1, theta2, phi2)->2*6378.14*arcsin( sqrt(haversin(theta2-theta1) + cos(theta1)*cos(theta2)*haversin(phi2-phi1)) );

Usage:

distance(0.6304129261, -1.512676863, 0.5923647483, -2.066469834)
Output:
2889.679287

Mathematica / Wolfram Language

Inputs assumed in degrees. Sin and Haversine expect arguments in radians; the built-in variable 'Degree' converts from degrees to radians.

distance[{theta1_, phi1_}, {theta2_, phi2_}] := 
 2*6378.14 ArcSin@
   Sqrt[Haversine[(theta2 - theta1) Degree] + 
     Cos[theta1*Degree] Cos[theta2*Degree] Haversine[(phi2 - phi1) Degree]]

Usage:

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

MATLAB / Octave

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
Output:
distance: 2887.2600 km

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 */

МК-61/52

П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	*	С/П

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.

MySQL

DELIMITER $$

CREATE FUNCTION haversine (
		lat1 FLOAT, lon1 FLOAT,
		lat2 FLOAT, lon2 FLOAT
	) RETURNS FLOAT
	NO SQL DETERMINISTIC
BEGIN
	DECLARE r FLOAT unsigned DEFAULT 6372.8;
	DECLARE dLat FLOAT unsigned;
	DECLARE dLon FLOAT unsigned;
	DECLARE a FLOAT unsigned;
	DECLARE c FLOAT unsigned;
	
	SET dLat = ABS(RADIANS(lat2 - lat1));
	SET dLon = ABS(RADIANS(lon2 - lon1));
	SET lat1 = RADIANS(lat1);
	SET lat2 = RADIANS(lat2);
	
	SET a = POW(SIN(dLat / 2), 2) + COS(lat1) * COS(lat2) * POW(SIN(dLon / 2), 2);
	SET c = 2 * ASIN(SQRT(a));
	
	RETURN (r * c);
END$$

DELIMITER ;

Usage:

SELECT haversine(36.12, -86.67, 33.94, -118.4);
Output:
2887.260009765625

Nim

import std/math

proc haversine(lat1, lon1, lat2, lon2: float): float =
  const r = 6372.8 # Earth radius in kilometers
  let
    dLat = degToRad(lat2 - lat1)
    dLon = degToRad(lon2 - lon1)
    lat1 = degToRad(lat1)
    lat2 = degToRad(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)
Output:
2887.259950607111

Oberon-2

Works with oo2c version2

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.

Output:

2887.2602975600

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));
    }
  }
}
Output:
distance: 2886.44

Objective-C

+ (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;
}

OCaml

The core calculation is fairly straightforward, but with an eye toward generality and reuse, this is how I might start:

(* 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;;

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

# earth_dist bna lax;;
- : float = 2887.25995060711102

Oforth

import: math

: 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
Output:
2887.25995060711

ooRexx

Translation of: REXX

The rxmath library provides the required functions.

/*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
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

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)
Output:
%1 = 2886.44444

Pascal

Works with: Free_Pascal
Library: Math
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.
Output:
Haversine distance: 2887.26 km.

Perl

Low-Level

Library: ntheory
use ntheory qw/Pi/;

sub asin { my $x = shift; atan2($x, sqrt(1-$x*$x)); }

sub surfacedist {
  my($lat1, $lon1, $lat2, $lon2) = @_;
  my $radius = 6372.8;
  my $radians = Pi() / 180;;
  my $dlat = ($lat2 - $lat1) * $radians;
  my $dlon = ($lon2 - $lon1) * $radians;
  $lat1 *= $radians;
  $lat2 *= $radians;
  my $a = sin($dlat/2)**2 + cos($lat1) * cos($lat2) * sin($dlon/2)**2;
  my $c = 2 * asin(sqrt($a));
  return $radius * $c;
}
my @BNA = (36.12, -86.67);
my @LAX = (33.94, -118.4);
printf "Distance: %.3f km\n", surfacedist(@BNA, @LAX);
Output:
Distance: 2887.260 km

Idiomatic

Contrary to ntheory, Math::Trig is part of the Perl core distribution. It comes with a great circle distance built-in.

use Math::Trig qw(great_circle_distance deg2rad);
 
# Notice the 90 - latitude: phi zero is at the North Pole.
# Parameter order is: LON, LAT
my @BNA = (deg2rad(-86.67), deg2rad(90 - 36.12));
my @LAX = (deg2rad(-118.4), deg2rad(90 - 33.94));

print "Distance: ", great_circle_distance(@BNA, @LAX, 6372.8), " km\n";
Output:
Distance: 2887.25995060711 km

Phix

function haversine(atom lat1, long1, lat2, long2)
    constant MER = 6371, -- mean earth radius(km)
             DEG_TO_RAD = PI/180
    lat1 *= DEG_TO_RAD
    lat2 *= DEG_TO_RAD
    long1 *= DEG_TO_RAD
    long2 *= DEG_TO_RAD
    return MER*arccos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(long2-long1))
end function
 
atom d = haversine(36.12,-86.67,33.94,-118.4)
printf(1,"Distance is %f km (%f miles)\n",{d,d/1.609344})
Output:
Distance is 2886.444443 km (1793.553425 miles)

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 = 6371; // Earth's radius in kilometers.

        $diffLatitude = $other->getLatitude() - $this->latitude;
        $diffLongitude = $other->getLongitude() - $this->longitude;

        $a = sin($diffLatitude / 2) ** 2 +
             cos($this->latitude) *
             cos($other->getLatitude()) *
             sin($diffLongitude / 2) ** 2;

        $c = 2 * asin(sqrt($a));
        $distance = $radiusOfEarth * $c;

        return $distance;
    }
}

Test:

$bna = new POI(36.12, -86.67); // Nashville International Airport
$lax = new POI(33.94, -118.40); // Los Angeles International Airport
printf('%.2f km', $bna->getDistanceInMetersTo($lax));
Output:
2886.44 km

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 ) ) ) ) )

Test:

(prinl
   "Haversine distance: "
   (round (haversine 36.12 -86.67 33.94 -118.4))
   " km" )
Output:
Haversine distance: 2,886.444 km

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;
Output:
distance:   2887.260 km

PowerShell

Works with: PowerShell version 3
Add-Type -AssemblyName System.Device
 
$BNA = New-Object System.Device.Location.GeoCoordinate 36.12, -86.67
$LAX = New-Object System.Device.Location.GeoCoordinate 33.94, -118.40
 
$BNA.GetDistanceTo( $LAX ) / 1000
Output:
2888.93627213254
Works with: PowerShell version 2
function Get-GreatCircleDistance ( $Coord1, $Coord2 )
    {
    #  Convert decimal degrees to radians
    $Lat1  = $Coord1[0] / 180 * [math]::Pi
    $Long1 = $Coord1[1] / 180 * [math]::Pi
    $Lat2  = $Coord2[0] / 180 * [math]::Pi
    $Long2 = $Coord2[1] / 180 * [math]::Pi
 
    #  Mean Earth radius (km)
    $R = 6371
   
    #  Haversine formula
    $ArcLength = 2 * $R *
                    [math]::Asin(
                        [math]::Sqrt(
                            [math]::Sin( ( $Lat1 - $Lat2 ) / 2 ) *
                            [math]::Sin( ( $Lat1 - $Lat2 ) / 2 ) +
                            [math]::Cos( $Lat1 ) *
                            [math]::Cos( $Lat2 ) *
                            [math]::Sin( ( $Long1 - $Long2 ) / 2 ) *
                            [math]::Sin( ( $Long1 - $Long2 ) / 2 ) ) )
    return $ArcLength
    }
 
$BNA = 36.12,  -86.67
$LAX = 33.94, -118.40
 
Get-GreatCircleDistance $BNA $LAX
Output:
2886.44444283799

Pure Data

Up until now there is no 64bit float in Pure Data, so the result of the calculation might not be completely accurate.

#N canvas 527 1078 450 686 10;
#X obj 28 427 atan2;
#X obj 28 406 sqrt;
#X obj 62 405 sqrt;
#X obj 28 447 * 2;
#X obj 62 384 -;
#X msg 62 362 1 \$1;
#X obj 28 339 t f f;
#X obj 28 210 sin;
#X obj 83 207 sin;
#X obj 138 206 cos;
#X obj 193 206 cos;
#X obj 28 179 / 2;
#X obj 83 182 / 2;
#X obj 28 74 unpack f f;
#X obj 28 98 t f f;
#X obj 28 301 expr $f1 + ($f2 * $f3 * $f4);
#X obj 28 148 deg2rad;
#X obj 83 149 deg2rad;
#X obj 138 148 deg2rad;
#X obj 193 149 deg2rad;
#X obj 28 232 t f f;
#X obj 28 257 *;
#X obj 83 232 t f f;
#X obj 83 257 *;
#X obj 83 98 t f b;
#X obj 28 542 * 6372.8;
#X obj 193 120 f 33.94;
#X obj 28 125 - 33.94;
#X msg 28 45 36.12 -86.67;
#X obj 83 123 - -118.4;
#X floatatom 28 577 8 0 0 0 - - -, f 8;
#X connect 0 0 3 0;
#X connect 1 0 0 0;
#X connect 2 0 0 1;
#X connect 3 0 25 0;
#X connect 4 0 2 0;
#X connect 5 0 4 0;
#X connect 6 0 1 0;
#X connect 6 1 5 0;
#X connect 7 0 20 0;
#X connect 8 0 22 0;
#X connect 9 0 15 2;
#X connect 10 0 15 3;
#X connect 11 0 7 0;
#X connect 12 0 8 0;
#X connect 13 0 14 0;
#X connect 13 1 24 0;
#X connect 14 0 27 0;
#X connect 14 1 18 0;
#X connect 15 0 6 0;
#X connect 16 0 11 0;
#X connect 17 0 12 0;
#X connect 18 0 9 0;
#X connect 19 0 10 0;
#X connect 20 0 21 0;
#X connect 20 1 21 1;
#X connect 21 0 15 0;
#X connect 22 0 23 0;
#X connect 22 1 23 1;
#X connect 23 0 15 1;
#X connect 24 0 29 0;
#X connect 24 1 26 0;
#X connect 25 0 30 0;
#X connect 26 0 19 0;
#X connect 27 0 16 0;
#X connect 28 0 13 0;
#X connect 29 0 17 0;

PureBasic

Translation of: Pascal
#DIA=2*6372.8

Procedure.d Haversine(th1.d,ph1.d,th2.d,ph2.d)
  Define dx.d,
         dy.d,
         dz.d
  
  ph1=Radian(ph1-ph2)
  th1=Radian(th1)
  th2=Radian(th2)
  
  dz=Sin(th1)-Sin(th2)
  dx=Cos(ph1)*Cos(th1)-Cos(th2)
  dy=Sin(ph1)*Cos(th1)
  ProcedureReturn ASin(Sqr(Pow(dx,2)+Pow(dy,2)+Pow(dz,2))/2)*#DIA
EndProcedure

OpenConsole("Haversine distance")
Print("Haversine distance: ")
Print(StrD(Haversine(36.12,-86.67,33.94,-118.4),7)+" km.")
Input()
Output:
Haversine distance: 2887.2599506 km.

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
>>>

QB64

Translation of: BASIC
SCREEN _NEWIMAGE(800, 100, 32)

'*** Units: K=kilometers  M=miles  N=nautical miles
DIM UNIT AS STRING
DIM Distance AS STRING
DIM Result AS DOUBLE
DIM ANSWER AS DOUBLE

'*** Change the To/From Latittude/Logitudes for your run

'*** LAT/LON for Nashville International Airport (BNA)
lat1 = 36.12
Lon1 = -86.67

'*** LAT/LONG for Los Angeles International Airport (LAX)
Lat2 = 33.94
Lon2 = -118.40

'*** Initialize Values
UNIT = "K"
Distance = ""
'Radius = 6378.137
Radius = 6372.8

'*** Calculate distance using Haversine Function
lat1 = (lat1 * _PI / 180)
Lon1 = (Lon1 * _PI / 180)
Lat2 = (Lat2 * _PI / 180)
Lon2 = (Lon2 * _PI / 180)
DLon = Lon1 - Lon2

ANSWER = _ACOS(SIN(lat1) * SIN(Lat2) + COS(lat1) * COS(Lat2) * COS(DLon)) * Radius

'*** Adjust Answer based on Distance Unit (kilometers, miles, nautical miles)
SELECT CASE UNIT
       CASE "M"
            Result = ANSWER * 0.621371192
            Distance = "miles"
       CASE "N"
            Result = ANSWER * 0.539956803
            Distance = "nautical miles"
       CASE ELSE
            Result = ANSWER
            Distance = "kilometers"
END SELECT

'*** Change PRINT statement with your labels for FROM/TO locations
PRINT "The distance from Nashville International to Los Angeles International in "; Distance;
PRINT USING " is: ##,###.##"; Result;
PRINT "."

END

R

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

# Volumetric mean radius is 6371 km, see http://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html
# 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))
}

# Coordinates are found here:
#     http://www.airport-data.com/airport/BNA/
#     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)

# Output:  2886.327

Racket

Almost the same as the Scheme version.

#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))
Output:
2886.444442837984

Raku

(formerly Perl 6)

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

Raven

Translation of: Groovy
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
Output:
haversine: 2887.25995060711

REXX

The use of normalization for angles isn't required for the Haversine formula, but those normalization functions were included
herein anyway   (to support normalization of input arguments to the trigonometric functions for the general case).

/*REXX program  calculates  the  distance between  Nashville  and  Los Angles  airports.*/
call pi;    numeric digits length(pi) % 2        /*use half of PI dec. digits for output*/
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º"
@using_radius= 'using the mean radius of the earth as '            /*a literal for  SAY.*/
radii.=.;    radii.1=6372.8;     radii.2=6371    /*mean radii of the earth in kilometers*/
say;                         m=1/0.621371192237  /*M:   one statute mile  in      "     */
    do radius=1  while radii.radius\==.          /*calc. distance using specific radii. */
    d= surfaceDist( 36.12,    -86.67,    33.94,   -118.4,    radii.radius);         say
    say center(@using_radius     radii.radius         ' kilometers', 75, '─')
    say ' Distance between:  '   format(d/1            ,,2)    " kilometers,"
    say '               or   '   format(d/m            ,,2)    " statute miles,"
    say '               or   '   format(d/m*5280/6076.1,,2)    " nautical (or air miles)."
    end   /*radius*/                             /*show──┘   2 dec. digs past dec. point*/
exit                                             /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
surfaceDist: parse arg th1,ph1,th2,ph2,r         /*use  haversine  formula for distance.*/
      numeric digits digits() * 2                /*double number of decimal digits used.*/
              ph1 = d2r(ph1 - ph2)               /*convert degrees ──► radians & reduce.*/
              th1 = d2r(th1)                     /*   "       "           "    "    "   */
              th2 = d2r(th2)                     /*   "       "           "    "    "   */
      cosTH1= cos(th1)                           /*compute a shortcut (it's used twice).*/
                x = cos(ph1) * cosTH1 - cos(th2) /*   "    X   coordinate.              */
                y = sin(ph1) * cosTH1            /*   "    Y       "                    */
                z = sin(th1)          - sin(th2) /*   "    Z       "                    */
      return Asin(sqrt(x*x + y*y + z*z)*.5) *r*2 /*compute the arcsin and return value. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
Acos: return pi() * .5   -   aSin( arg(1) )      /*calculate the ArcCos of an argument. */
d2d:  return arg(1)               //  360        /*normalize degrees to a  unit circle. */
d2r:  return r2r(  arg(1) * pi()  /   180)       /*normalize and convert deg ──► radians*/
r2d:  return d2d( (arg(1) * 180   /   pi()))     /*normalize and convert rad ──► degrees*/
r2r:  return arg(1)               // (pi() * 2)  /*normalize radians to a  unit circle. */
pi:   pi= 3.141592653589793238462643383279502884197169399375105820975;           return pi
/*──────────────────────────────────────────────────────────────────────────────────────*/
Asin: procedure; parse arg x 1 z 1 o 1 p;     a= abs(x);               aa= a * a
      if a >= sqrt(2) * .5  then return sign(x) * Acos( sqrt(1 - aa) )
        do j=2  by 2  until p=z;    p= z;     o= o * aa * (j-1) / j;   z= z  +  o / (j+1)
        end   /*j*/;                return z      /* [↑]  compute until no more noise.  */
/*──────────────────────────────────────────────────────────────────────────────────────*/
cos:  procedure; parse arg x;       x= r2r(x);       a= abs(x);              Hpi= pi * .5
      numeric fuzz min(6, digits() - 3)  ;     if a=pi    then return -1
      if a=Hpi | a=Hpi*3  then return 0  ;     if a=pi/3  then return .5
      if a=pi* 2/3        then return -.5;     q= x*x;    p= 1;     z= 1;     _= 1
        do k=2  by 2;  _= -_*q / (k*(k-1)); z= z+_; if z=p  then leave; p=z; end; return z
/*──────────────────────────────────────────────────────────────────────────────────────*/
sin:  procedure; parse arg x;  x= r2r(x);                numeric fuzz min(5, digits() - 3)
      if abs(x)=pi  then  return 0;            q= x*x;    p= x;     z= x;      _= x
        do k=2  by 2; _= -_*q / (k*(k+1));  z= z+_; if z=p  then leave; p=z; end; return z
/*──────────────────────────────────────────────────────────────────────────────────────*/
sqrt: procedure; parse arg x; if x=0  then return 0; d=digits(); m.=9; numeric form; h=d+6
      numeric digits;  parse value format(x,2,1,,0) 'E0' with g "E" _ .;  g=g * .5'e'_ % 2
        do j=0  while h>9;      m.j= h;              h= h%2 + 1;    end  /*j*/
        do k=j+5  to 0  by -1;  numeric digits m.k;  g= (g+x/g)*.5; end  /*k*/;   return g

REXX doesn't have most of the higher math functions, so they are included here (above) as subroutines (functions).

      ╔════════════════════════════════════════════════════════════════════════╗
      ║ A note on built─in functions:  REXX doesn't have a lot of mathematical ║
      ║ or  (particularly) trigonometric 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.                                        ║
      ║                                                                        ║
      ║ Programming note:  the  "general 1─liner"  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    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.  ║
      ╚════════════════════════════════════════════════════════════════════════╝ 
output   when using the in-line defaults:
       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º


─────────using the mean radius of the earth as  6372.8  kilometers─────────
 Distance between:   2887.26  kilometers,
               or    1794.06  statute miles,
               or    1559.00  nautical (or air miles).

──────────using the mean radius of the earth as  6371  kilometers──────────
 Distance between:   2886.44  kilometers,
               or    1793.55  statute miles,
               or    1558.56  nautical (or air miles).

Ring

decimals(8)
see haversine(36.12, -86.67, 33.94, -118.4) + nl

func haversine x1, y1, x2, y2
     r=0.01745
     x1= x1*r
     x2= x2*r
     y1= y1*r
     y2= y2*r
     dy = y2-y1
     dx = x2-x1
     a = pow(sin(dx/2),2) + cos(x1) * cos(x2) * pow(sin(dy/2),2)
     c = 2 * asin(sqrt(a))
     d = 6372.8 * c
     return d

RPL

Works with: Halcyon Calc version 4.2.7
Code Comments
≪ 
  ROT - 2 / DEG SIN SQ OVER COS * 3 PICK COS * 
  ROT ROT - 2 / SIN SQ + √ RAD ASIN 6372.8 * 2 *
≫ 'AHAV' STO
 ( lat1 lon1 lat2 lon2 -- distance )
 Start by the end of the formula, in degree mode
 Switch to radian mode to compute Arcsin

The following line of command delivers what is required:

36.12 -86.67 33.94 -118.4 AHAV

Due to the uncertainty in values of Earth radius and airports coordinates, the result shall be announced as 2887 ± 1 km even if the calculation provides many digits after the decimal point

Output:
1:       2887.25995061

Ruby

include Math

Radius = 6372.8  # 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)
Output:
2887.3

Alternatively:

Translation of: Python
include Math
 
def haversine(lat1, lon1, lat2, lon2)
    r = 6372.8        # Earth radius in kilometers
    deg2rad = PI/180  # convert degress to radians
 
    dLat = (lat2 - lat1) * deg2rad
    dLon = (lon2 - lon1) * deg2rad
    lat1 = lat1 * deg2rad
    lat2 = lat2 * deg2rad
 
    a = sin(dLat / 2)**2 + cos(lat1) * cos(lat2) * sin(dLon / 2)**2
    c = 2 * asin(sqrt(a))
    r * c
end

puts "distance is #{haversine(36.12, -86.67, 33.94, -118.40)} km "
Output:
distance is 2887.2599506071106 km 

Run BASIC

    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.
 '
 '      http://maps.google.com
 '      Search   36.12,-86.67
 '      Earth.
 '      Center the pin, zoom airport.
 '      Directions (destination).
 '      36.12.-86.66999
 '      Distance is 35.37 inches.
Output
Haversine distance: 2887.2599506071104 km.

Rust

struct Point {
    lat: f64,
    lon: f64,
}

fn haversine(origin: Point, destination: Point) -> f64 {
    const R: f64 = 6372.8;

    let lat1 = origin.lat.to_radians();
    let lat2 = destination.lat.to_radians();
    let d_lat = lat2 - lat1;
    let d_lon = (destination.lon - origin.lon).to_radians();

    let a = (d_lat / 2.0).sin().powi(2) + (d_lon / 2.0).sin().powi(2) * lat1.cos() * lat2.cos();
    let c = 2.0 * a.sqrt().asin();
    R * c
}

#[cfg(test)]
mod test {
    use super::{Point, haversine};

    #[test]
    fn test_haversine() {
        let origin: Point = Point {
            lat: 36.12,
            lon: -86.67
        };
        let destination: Point = Point {
            lat: 33.94,
            lon: -118.4
        };
        let d: f64 = haversine(origin, destination);
        println!("Distance: {} km ({} mi)", d, d / 1.609344);
        assert_eq!(d, 2887.2599506071106);
    }

}
Output
Distance: 2887.2599506071106 km (1794.060157807846 mi)

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);
Output:
Distance is 2887.2599506 K

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))
  }
}
Output:
2887.2599506071106

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

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;
Output:
2887.26

Sidef

Translation of: Raku
class EarthPoint(lat, lon) {

    const earth_radius = 6371       # mean earth radius
    const radian_ratio = Num.pi/180

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

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

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

        earth_radius * a.sqrt.asin * 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.444442837983299747157823945746716...

smart BASIC

Translation of: BASIC
'*** LAT/LONG for Nashville International Airport (BNA)
lat1=36.12
Lon1=-86.67

'*** LAT/LONG for Los Angeles International Airport (LAX)
Lat2=33.94
Lon2=-118.40

'*** Units: K=kilometers  M=miles  N=nautical miles
Unit$ = "K"	

Result=HAVERSINE(Lat1,Lon1,Lat2,Lon2,Unit$)
R$=STR$(Result,"#,###.##")

PRINT "The distance between Nashville International Airport and Los Angeles International Airport in kilometers is: "&R$

STOP

DEF HAVERSINE(Lat1,Lon1,Lat2,Lon2,Unit$)
'---------------------------------------------------------------
'*** Haversine Formula - Calculate distances by LAT/LONG
'

'*** Pass to it the LAT/LONG of the two locations, and then unit of measure
'*** Usage: X=HAVERSINE(Lat1,Lon1,Lat2,Lon2,Unit$)

    PI=3.14159265358979323846
    Radius=6372.8
    Lat1=(Lat1*PI/180)
    Lon1=(Lon1*PI/180)
    Lat2=(Lat2*PI/180)
    Lon2=(Lon2*PI/180)
    DLon=Lon1-Lon2
    Answer=ACOS(SIN(Lat1)*SIN(Lat2)+COS(Lat1)*COS(Lat2)*COS(DLon))*Radius

    IF UNIT$="M" THEN Answer=Answer*0.621371192
    IF UNIT$="N" THEN Answer=Answer*0.539956803

  RETURN Answer

ENDDEF
Output:
The distance between Nashville International Airport and Los Angeles International Airport in kilometers is: 2,887.26

Stata

First, a program to add a distance variable to a dataset, given variables for LAT/LON of two points.

program spheredist
	version 15.0
	syntax varlist(min=4 max=4 numeric), GENerate(namelist max=1) ///
		[Radius(real 6371) ALTitude(real 0) LABel(string)]
	confirm new variable `generate'
	local lat1 : word 1 of `varlist'
	local lon1 : word 2 of `varlist'
	local lat2 : word 3 of `varlist'
	local lon2 : word 4 of `varlist'
	local r=2*(`radius'+`altitude'/1000)
	local k=_pi/180
	gen `generate'=`r'*asin(sqrt(sin((`lat2'-`lat1')*`k'/2)^2+ ///
		cos(`lat1'*`k')*cos(`lat2'*`k')*sin((`lon2'-`lon1')*`k'/2)^2))
	if `"`label'"' != "" {
		label variable `generate' `"`label'"'
	}
end

Illustration with a sample dataset.

import delimited airports.csv, clear
format %9.4f l*
list

     +----------------------------------------------------------------------------------------------------+
     | iata                                   airport          city         country       lat         lon |
     |----------------------------------------------------------------------------------------------------|
  1. |  AMS                Amsterdam Airport Schiphol     Amsterdam     Netherlands   52.3086      4.7639 |
  2. |  BNA           Nashville International Airport     Nashville   United States   36.1245    -86.6782 |
  3. |  CDG   Charles de Gaulle International Airport         Paris          France   49.0128      2.5500 |
  4. |  CGN                      Cologne Bonn Airport       Cologne         Germany   50.8659      7.1427 |
  5. |  LAX         Los Angeles International Airport   Los Angeles   United States   33.9425   -118.4080 |
     |----------------------------------------------------------------------------------------------------|
  6. |  MEM             Memphis International Airport       Memphis   United States   35.0424    -89.9767 |
     +----------------------------------------------------------------------------------------------------+

MEM/CGN joins two Fedex Express hubs. The line AMS/LAX is operated by KLM Royal Dutch Airlines. We will compute the distance between each pair of airports, both at sea level and at typical cruising flight level (35000 ft).

Bear in mind that the actual route of an airliner is usually not a piece of great circle, so this will only give an idea. For instance, according to FlightAware, the route of a Fedex flight from Memphis to Paris is 7852 km long, at FL300 altitude (9150 m). The program given here would yield 7328.33 km instead.

keep iata lat lon
rename (iata lat lon) =2
gen k=0
tempfile tmp
save "`tmp'"
rename *2 *1
joinby k using `tmp'
drop if iata1>=iata2
drop k
list

     +-----------------------------------------------------------+
     | iata1      lat1        lon1   iata2      lat2        lon2 |
     |-----------------------------------------------------------|
  1. |   AMS   52.3086      4.7639     BNA   36.1245    -86.6782 |
  2. |   AMS   52.3086      4.7639     CGN   50.8659      7.1427 |
  3. |   AMS   52.3086      4.7639     LAX   33.9425   -118.4080 |
  4. |   AMS   52.3086      4.7639     CDG   49.0128      2.5500 |
  5. |   AMS   52.3086      4.7639     MEM   35.0424    -89.9767 |
     |-----------------------------------------------------------|
  6. |   BNA   36.1245    -86.6782     CGN   50.8659      7.1427 |
  7. |   BNA   36.1245    -86.6782     CDG   49.0128      2.5500 |
  8. |   BNA   36.1245    -86.6782     LAX   33.9425   -118.4080 |
  9. |   BNA   36.1245    -86.6782     MEM   35.0424    -89.9767 |
 10. |   CDG   49.0128      2.5500     LAX   33.9425   -118.4080 |
     |-----------------------------------------------------------|
 11. |   CDG   49.0128      2.5500     MEM   35.0424    -89.9767 |
 12. |   CDG   49.0128      2.5500     CGN   50.8659      7.1427 |
 13. |   CGN   50.8659      7.1427     LAX   33.9425   -118.4080 |
 14. |   CGN   50.8659      7.1427     MEM   35.0424    -89.9767 |
 15. |   LAX   33.9425   -118.4080     MEM   35.0424    -89.9767 |
     +-----------------------------------------------------------+

Now compute the distances and print the result.

spheredist lat1 lon1 lat2 lon2, gen(dist) lab(Distance at sea level)
spheredist lat1 lon1 lat2 lon2, gen(fl350) alt(10680) lab(Distance at FL350 altitude)
format %9.2f dist fl350
list iata* dist fl350

     +-----------------------------------+
     | iata1   iata2      dist     fl350 |
     |-----------------------------------|
  1. |   AMS     CGN    229.64    230.03 |
  2. |   AMS     CDG    398.27    398.94 |
  3. |   AMS     MEM   7295.19   7307.56 |
  4. |   AMS     BNA   7004.61   7016.48 |
  5. |   AMS     LAX   8955.95   8971.13 |
     |-----------------------------------|
  6. |   BNA     LAX   2886.32   2891.21 |
  7. |   BNA     CGN   7222.75   7234.99 |
  8. |   BNA     CDG   7018.39   7030.29 |
  9. |   BNA     MEM    321.62    322.16 |
 10. |   CDG     LAX   9102.51   9117.94 |
     |-----------------------------------|
 11. |   CDG     CGN    387.82    388.48 |
 12. |   CDG     MEM   7317.82   7330.23 |
 13. |   CGN     LAX   9185.47   9201.04 |
 14. |   CGN     MEM   7514.96   7527.70 |
 15. |   LAX     MEM   2599.71   2604.12 |
     +-----------------------------------+

Notice that the distance from Nashville to Los Angeles is given as 2886.32 km, which is slightly different from the task description. The coordinates come from OpenFlights and are supposably more accurate. Using the data in the task description, one gets 2886.44 as expected.

Swift

Translation of: Objective-C
import Foundation

func haversine(lat1:Double, lon1:Double, lat2:Double, lon2:Double) -> Double {
    let lat1rad = lat1 * Double.pi/180
    let lon1rad = lon1 * Double.pi/180
    let lat2rad = lat2 * Double.pi/180
    let lon2rad = lon2 * Double.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
}

print(haversine(lat1:36.12, lon1:-86.67, lat2:33.94, lon2:-118.40))
Output:
2887.25995060711

Symsyn

lat1 : 36.12
lon1 : -86.67
lat2 : 33.94
lon2 : -118.4

dx : 0.
dy : 0.
dz : 0.
kms : 0.

 {degtorad(lon2 - lon1)} lon1
 {degtorad lat1} lat1
 {degtorad lat2} lat2

 {sin lat1 - sin lat2} dz
 {cos lon1 * cos lat1 - cos lat2} dx
 {sin lon1 * cos lat1} dy

 {arcsin(sqrt(dx^2 + dy^2 + dz^2)/2) * 12745.6} kms

 "'Haversine distance: ' kms ' kms'" []
Output:
Haversine distance:        2887.259951 kms

tbas

option angle radians ' the default
sub haversine(lat1, lon1, lat2, lon2)
	dim EarthRadiusKm = 6372.8        ' Earth radius in kilometers
	dim latRad1 = RAD(lat1)
	dim latRad2 = RAD(lat2)
	dim lonRad1 = RAD(lon1)
	dim lonRad2 = RAD(lon2)
	dim _diffLa = latRad2 - latRad1
	dim _doffLo = lonRad2 - lonRad1
	dim sinLaSqrd = sin(_diffLa / 2) ^ 2
	dim sinLoSqrd = sin(_doffLo / 2) ^ 2
	dim computation = asin(sqrt(sinLaSqrd + cos(latRad1) * cos(latRad2) * sinLoSqrd))
	return 2 * EarthRadiusKm * computation
end sub

print using "Nashville International Airport to Los Angeles International Airport ####.########### km", haversine(36.12, -86.67, 33.94, -118.40)
print using "Perth, WA Australia to Baja California, Mexico #####.########### km", haversine(-31.95, 115.86, 31.95, -115.86)
Nashville International Airport to Los Angeles International Airport  2887.25995060712 km
Perth, WA Australia to Baja California, Mexico 15188.70229560390 km

Tcl

Translation of: Groovy
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}]
}

# Don't bother with too much inappropriate accuracy!
puts [format "distance=%.1f km" [haversineFormula 36.12 -86.67 33.94 -118.40]]
Output:
distance=2887.3 km

TechBASIC

{{trans|BASIC}}
FUNCTION HAVERSINE
!---------------------------------------------------------------
!*** Haversine Formula - Calculate distances by LAT/LONG
!

!*** LAT/LON of the two locations and Unit of measure are GLOBAL
!*** as they are defined in the main logic of the program, so they
!*** available for use in the Function.
!*** Usage: X=HAVERSINE

    
    Radius=6378.137
    Lat1=(Lat1*MATH.PI/180)
    Lon1=(Lon1*MATH.PI/180)
    Lat2=(Lat2*MATH.PI/180)
    Lon2=(Lon2*MATH.PI/180)
    DLon=Lon1-Lon2
    ANSWER=ACOS(SIN(Lat1)*SIN(Lat2)+COS(Lat1)*COS(Lat2)*COS(DLon))*Radius

    DISTANCE="kilometers"
    SELECT CASE UNIT
           CASE "M"
                HAVERSINE=ANSWER*0.621371192
                Distance="miles"
           CASE "N"
                HAVERSINE=ANSWER*0.539956803
                Distance="nautical miles"
    END SELECT       

END FUNCTION


The following is the main code that invokes the function. It takes your location and determines how far away you are from Tampa, Florida. You can change UNIT to either M=Miles, N=Nautical Miles, or K (or leave blank) as default is in Kilometers:

!*** In techBASIC, all variables defined in the main program act as GLOBAL
!*** variables and are available to all SUBROUTINES and FUNCTIONS. So in the
!*** HAVERSINE Function being used, no paramaters need to be passed to it, so
!*** it acts as a variable when I use it as Result=HAVERSINE. The way that
!*** the Function is setup, it returns its value back as HAVERSINE.

BASE 1

!*** Get the GPS LAT/LONG of current location
location = sensors.location(30)
Lat1=location(1) 
Lon1=location(2) 

!*** LAT/LONG For Tampa, FL
Lat2=27.9506
Lon2=-82.4572

!*** Units: K=kilometers  M=miles  N=nautical miles
DIM UNIT      AS STRING 
DIM Distance  AS STRING
DIM Result    AS SINGLE
UNIT = "M"	

!*** Calculate distance using Haversine Function
Result=HAVERSINE

PRINT "The distance from your current location to Tampa, FL in ";Distance;" is: ";
PRINT USING "#,###.##";Result;"."

STOP

OUTPUT: *** NOTE: When I run this, I am in my house in Venice, Florida, and that distance is correct (as the crow flies). ***

The distance from your current location to Tampa, FL in miles is:    57.94

Teradata Stored Procedure

# syntax: call SP_HAVERSINE(36.12,33.94,-86.67,-118.40,x);

CREATE PROCEDURE SP_HAVERSINE
(
IN lat1 FLOAT,
IN lat2 FLOAT,
IN lon1 FLOAT,
IN lon2 FLOAT,
OUT distance FLOAT)
 
BEGIN 
    DECLARE dLat FLOAT;
    DECLARE dLon FLOAT;
    DECLARE c FLOAT;
    DECLARE a FLOAT;    
    DECLARE km FLOAT;

    SET dLat = RADIANS(lat2-lat1);
    SET dLon = RADIANS(lon2-lon1);

    SET a = SIN(dLat / 2) * SIN(dLat / 2) + SIN(dLon / 2) * SIN(dLon / 2) * COS(RADIANS(lat1)) * COS(RADIANS(lat2));
    SET c = 2 * ASIN(SQRT(a));
    SET km = 6372.8 * c;
    
    select km into distance;
END;
Output:
distance: 2887.2599 km

Transact-SQL

Translation of: C#
CREATE FUNCTION [dbo].[Haversine](@Lat1 AS DECIMAL(9,7), @Lon1 AS DECIMAL(10,7), @Lat2 AS DECIMAL(9,7), @Lon2 AS DECIMAL(10,7))
RETURNS DECIMAL(12,7)
AS
BEGIN
	DECLARE @R	DECIMAL(11,7);
	DECLARE @dLat	DECIMAL(9,7);
	DECLARE @dLon	DECIMAL(10,7);
	DECLARE @a	DECIMAL(10,7);
	DECLARE @c	DECIMAL(10,7);

	SET @R		= 6372.8;
	SET @dLat	= RADIANS(@Lat2 - @Lat1);
	SET @dLon	= RADIANS(@Lon2 - @Lon1);
	SET @Lat1	= RADIANS(@Lat1);
	SET @Lat2	= RADIANS(@Lat2);
	SET @a		= SIN(@dLat / 2) * SIN(@dLat / 2) + SIN(@dLon / 2) * SIN(@dLon / 2) * COS(@Lat1) * COS(@Lat2);
	SET @c		= 2 * ASIN(SQRT(@a));

	RETURN @R * @c;
END
GO

SELECT dbo.Haversine(36.12,-86.67,33.94,-118.4)
Output:
 2887.2594934

TypeScript

Translation of: Matlab
let radians = function (degree: number) {

    // degrees to radians
    let rad: number = degree * Math.PI / 180;

    return rad;
}

export const haversine = (lat1: number, lon1: number, lat2: number, lon2: number) => {

    // var dlat: number, dlon: number, a: number, c: number, R: number;
    let dlat, dlon, a, c, R: number;

    R = 6372.8; // km
    dlat = radians(lat2 - lat1);
    dlon = radians(lon2 - lon1);
    lat1 = radians(lat1);
    lat2 = radians(lat2);
    a = Math.sin(dlat / 2) * Math.sin(dlat / 2) + Math.sin(dlon / 2) * Math.sin(dlon / 2) * Math.cos(lat1) * Math.cos(lat2)
    c = 2 * Math.asin(Math.sqrt(a));
    return R * c;
}

console.log("Distance:" + haversine(36.12, -86.67, 33.94, -118.40));
Output:
Distance: 2887.2599506071106

UBASIC

   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

VBA

Translation of: Phix
Const MER = 6371         '-- mean earth radius(km)
Public DEG_TO_RAD As Double
 
Function haversine(lat1 As Double, long1 As Double, lat2 As Double, long2 As Double) As Double
    lat1 = lat1 * DEG_TO_RAD
    lat2 = lat2 * DEG_TO_RAD
    long1 = long1 * DEG_TO_RAD
    long2 = long2 * DEG_TO_RAD
    haversine = MER * WorksheetFunction.Acos(Sin(lat1) * Sin(lat2) + Cos(lat1) * Cos(lat2) * Cos(long2 - long1))
End Function
 
Public Sub main()
    DEG_TO_RAD = WorksheetFunction.Pi / 180
    d = haversine(36.12, -86.67, 33.94, -118.4)
    Debug.Print "Distance is "; Format(d, "#.######"); " km ("; Format(d / 1.609344, "#.######"); " miles)."
End Sub
Output:
Distance is 2886,444443 km (1793,553425 miles).

Visual Basic .NET

Translation of: C#
If you read the fine print in the Wikipedia article, you will find that the Haversine method of finding distances may have an error of up to 0.5%. This could lead one to believe that discussion about whether to use 6371.0 km or 6372.8 km for an approximation of the Earth's radius is moot.
Imports System.Math

Module Module1

  Const deg2rad As Double = PI / 180

  Structure AP_Loc
    Public IATA_Code As String, Lat As Double, Lon As Double

    Public Sub New(ByVal iata_code As String, ByVal lat As Double, ByVal lon As Double)
      Me.IATA_Code = iata_code : Me.Lat = lat * deg2rad : Me.Lon = lon * deg2rad
    End Sub

    Public Overrides Function ToString() As String
      Return String.Format("{0}: ({1}, {2})", IATA_Code, Lat / deg2rad, Lon / deg2rad)
    End Function
  End Structure

  Function Sin2(ByVal x As Double) As Double
    Return Pow(Sin(x / 2), 2)
  End Function

  Function calculate(ByVal one As AP_Loc, ByVal two As AP_Loc) As Double
    Dim R As Double = 6371, ' In kilometers, (as recommended by the International Union of Geodesy and Geophysics)
        a As Double = Sin2(two.Lat - one.Lat) + Sin2(two.Lon - one.Lon) * Cos(one.Lat) * Cos(two.Lat)
    Return R * 2 * Asin(Sqrt(a))
  End Function

  Sub ShowOne(pntA As AP_Loc, pntB as AP_Loc)
    Dim adst As Double = calculate(pntA, pntB), sfx As String = "km"
    If adst < 1000 Then adst *= 1000 : sfx = "m"
    Console.WriteLine("The approximate distance between airports {0} and {1} is {2:n2} {3}.", pntA, pntB, adst, sfx)
    Console.WriteLine("The uncertainty is under 0.5%, or {0:n1} {1}." & vbLf, adst / 200, sfx)
  End Sub

' Airport coordinate data excerpted from the data base at http://www.partow.net/miscellaneous/airportdatabase/

' The four additional airports are the furthest and closest pairs, according to the "Fun Facts..." section.

' KBNA, BNA, NASHVILLE INTERNATIONAL, NASHVILLE, USA, 036, 007, 028, N, 086, 040, 041, W, 00183, 36.124, -86.678
' KLAX, LAX, LOS ANGELES INTERNATIONAL, LOS ANGELES, USA, 033, 056, 033, N, 118, 024, 029, W, 00039, 33.942, -118.408
' SKNV, NVA, BENITO SALAS, NEIVA, COLOMBIA, 002, 057, 000, N, 075, 017, 038, W, 00439, 2.950, -75.294
' WIPP, PLM, SULTAN MAHMUD BADARUDDIN II, PALEMBANG, INDONESIA, 002, 053, 052, S, 104, 042, 004, E, 00012, -2.898, 104.701 
' LOWL, LNZ, HORSCHING INTERNATIONAL AIRPORT (AUS - AFB), LINZ, AUSTRIA, 048, 014, 000, N, 014, 011, 000, E, 00096, 48.233, 14.183
' LOXL, N/A, LINZ, LINZ, AUSTRIA, 048, 013, 059, N, 014, 011, 015, E, 00299, 48.233, 14.188

  Sub Main()
    ShowOne(New AP_Loc("BNA", 36.124, -86.678),  New AP_Loc("LAX", 33.942, -118.408))
    ShowOne(New AP_Loc("NVA",  2.95,  -75.294),  New AP_Loc("PLM", -2.898,  104.701))
    ShowOne(New AP_Loc("LNZ", 48.233,  14.183),  New AP_Loc("N/A", 48.233,   14.188))
  End Sub
End Module
Output:
The approximate distance between airports BNA: (36.124, -86.678) and LAX: (33.942, -118.408) is 2,886.36 km.
The uncertainty is under 0.5%, or 14.4 km.

The approximate distance between airports NVA: (2.95, -75.294) and PLM: (-2.898, 104.701) is 20,009.28 km.
The uncertainty is under 0.5%, or 100.0 km.

The approximate distance between airports LNZ: (48.233, 14.183) and N/A: (48.233, 14.188) is 370.34 m.
The uncertainty is under 0.5%, or 1.9 m.
Looking at the altitude difference between the last two airports, (299 - 96 = 203), the reported distance of 370 meters ought to be around 422 meters if you actually went there and saw it for yourself.

V (Vlang)

Translation of: go
import math

fn haversine(h f64) f64 {
    return .5 * (1 - math.cos(h))
}
 
struct Pos {
    lat f64 // latitude, radians
    long f64 // longitude, radians
}
 
fn deg_pos(lat f64, lon f64) Pos {
    return Pos{lat * math.pi / 180, lon * math.pi / 180}
}
 
const r_earth = 6372.8 // km
 
fn hs_dist(p1 Pos, p2 Pos) f64 {
    return 2 * r_earth * math.asin(math.sqrt(haversine(p2.lat-p1.lat)+
        math.cos(p1.lat)*math.cos(p2.lat)*haversine(p2.long-p1.long)))
}
 
fn main() {
    println(hs_dist(deg_pos(36.12, -86.67), deg_pos(33.94, -118.40)))
}
Output:
2887.2599506071

Wren

Translation of: Julia
var R = 6372.8  // Earth's approximate radius in kilometers.

/*  Class containing trig methods which work with degrees rather than radians. */
class D {
    static deg2Rad(deg) { (deg*Num.pi/180 + 2*Num.pi) % (2*Num.pi) }
    static sin(d) { deg2Rad(d).sin }
    static cos(d) { deg2Rad(d).cos }
}

var haversine = Fn.new { |lat1, lon1, lat2, lon2|
    var dlat = lat2 - lat1
    var dlon = lon2 - lon1
    return 2 * R * (D.sin(dlat/2).pow(2) + D.cos(lat1) * D.cos(lat2) * D.sin(dlon/2).pow(2)).sqrt.asin
}

System.print(haversine.call(36.12, -86.67, 33.94, -118.4))
Output:
2887.2599506071

X86 Assembly

Assemble with tasm /m /l; tlink /t

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
Output:
2887.25

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 ));
Output:
 2887.25995

XQuery

declare namespace xsd = "http://www.w3.org/2001/XMLSchema";
declare namespace math = "http://www.w3.org/2005/xpath-functions/math";

declare function local:haversine($lat1 as xsd:float, $lon1 as xsd:float, $lat2 as xsd:float, $lon2 as xsd:float)
    as xsd:float
{
    let $dlat  := ($lat2 - $lat1) * math:pi() div 180
    let $dlon  := ($lon2 - $lon1) * math:pi() div 180
    let $rlat1 := $lat1 * math:pi() div 180
    let $rlat2 := $lat2 * math:pi() div 180
    let $a     := math:sin($dlat div 2) * math:sin($dlat div 2) + math:sin($dlon div 2) * math:sin($dlon div 2) * math:cos($rlat1) * math:cos($rlat2)
    let $c     := 2 * math:atan2(math:sqrt($a), math:sqrt(1-$a))
    return xsd:float($c * 6371.0)
};

local:haversine(36.12, -86.67, 33.94, -118.4)
Output:
 2886.444

Zig

Translation of: R

When a Zig struct type can be inferred then anonymous structs .{} can be used for initialisation. This can be seen on the lines where the constants bna and lax are instantiated.

A Zig struct can have methods, the same as an enum and or a union. They are only namespaced functions that can be called with dot syntax.

const std = @import("std");
const math = std.math; // Save some typing, reduce clutter. Otherwise math.sin() would be std.math.sin() etc.

pub fn main() !void {
    // Coordinates are found here:
    //     http://www.airport-data.com/airport/BNA/
    //     http://www.airport-data.com/airport/LAX/

    const bna = LatLong{
        .lat = .{ .d = 36, .m = 7, .s = 28.10 },
        .long = .{ .d = 86, .m = 40, .s = 41.50 },
    };

    const lax = LatLong{
        .lat = .{ .d = 33, .m = 56, .s = 32.98 },
        .long = .{ .d = 118, .m = 24, .s = 29.05 },
    };

    const distance = calcGreatCircleDistance(bna, lax);

    std.debug.print("Output: {d:.6} km\n", .{distance});

    // Output: 2886.326609 km
}

const LatLong = struct { lat: DMS, long: DMS };

/// degrees, minutes, decimal seconds
const DMS = struct {
    d: f64,
    m: f64,
    s: f64,

    fn toRadians(self: DMS) f64 {
        return (self.d + self.m / 60 + self.s / 3600) * math.pi / 180;
    }
};

// Volumetric mean radius is 6371 km, see http://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html
// The diameter is thus 12742 km
fn calcGreatCircleDistance(lat_long1: LatLong, lat_long2: LatLong) f64 {
    const lat1 = lat_long1.lat.toRadians();
    const lat2 = lat_long2.lat.toRadians();
    const long1 = lat_long1.long.toRadians();
    const long2 = lat_long2.long.toRadians();

    const a = math.sin(0.5 * (lat2 - lat1));
    const b = math.sin(0.5 * (long2 - long1));

    return 12742 * math.asin(math.sqrt(a * a + math.cos(lat1) * math.cos(lat2) * b * b));
}

zkl

Translation of: Erlang
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;
}
Output:
2887.26

ZX Spectrum Basic

Translation of: Run_BASIC
10 LET diam=2*6372.8
20 LET Lg1m2=FN r((-86.67)-(-118.4))
30 LET Lt1=FN r(36.12)
40 LET Lt2=FN r(33.94)
50 LET dz=SIN (Lt1)-SIN (Lt2)
60 LET dx=COS (Lg1m2)*COS (Lt1)-COS (Lt2)
70 LET dy=SIN (Lg1m2)*COS (Lt1)
80 LET hDist=ASN ((dx*dx+dy*dy+dz*dz)^0.5/2)*diam
90 PRINT "Haversine distance: ";hDist;" km."
100 STOP 
1000 DEF FN r(a)=a*0.017453293: REM convert degree to radians