From Wikipedia, the free encyclopedia
(Redirected from )
This article needs additional citations for . Please help
by . Unsourced material may be challenged and removed. (March 2014)
The versine or versed sine, versin(θ), is a
equal to 1 - cos(θ), or 2sin2( 1/2 θ). The function appeared in some of the earliest trigonometric tables. There are several related functions, most notably the haversine, half the versine, known in the
of navigation.
It is also written as vers(θ) or ver(θ). In , it is known as the sinus versus (flipped sine), versinus, versus or the sagitta (arrow).
There are several other related functions:
The versed cosine, or vercosine, written
The coversed sine, coversine, cosinus versus, or coversinus, written
and sometimes abbreviated to
The coversed cosine, or covercosine, written
The haversed sine, haversine, or semiversus, written
or , most famous from the
used historically in
The haversed cosine, or havercosine, written
The hacoversed sine, also called hacoversine or cohaversine and written
The hacoversed cosine, also called hacovercosine or cohavercosine and written
The , written
The excosecant, written
The trigonometric functions can be constructed geometrically in terms of a unit circle centered at O. This figure also illustrates the reason why the versine was sometimes called the sagitta, Latin for . If the arc ADB is viewed as a "" and the chord AB as its "string", then the versine CD is clearly the "arrow shaft".
Historically, the versed sine was considered one of the most important trigonometric functions. As θ goes to zero, versin(θ) is the difference between two nearly equal quantities, so a user of a trigonometric table for the cosine alone would need a very high accuracy to obtain the versine in order to avoid , making separate tables for the latter convenient. Even with a calculator or computer,
make it advisable to use the sin2 formula for small θ. Another historical advantage of the versine is that it is always non-negative, so its
is defined everywhere except for the single angle (θ = 0, 2π,...) where it is zero—thus, one could use logarithmic tables for multiplications in formulas involving versines.
The haversine, in particular, was important in
because it appears in the , which is used to reasonably accurately compute distances on a sphere (see issues with the Earth`s radius vs. sphere) given angular positions (e.g.,
and ). One could also use sin2(θ/2) directly, but having a table of the haversine removed the need to compute squares and square roots. The term haversine was, apparently, coined in a navigation text for just such an application.
In fact, the earliest surviving table of
(half-) values (as opposed to the
and other Greek authors), calculated from the
of India dated back to the 3rd century BC, was a table of values for the sine and versed sine (in 3.75° increments from 0 to 90°). The versine appears as an intermediate step in the application of the half-angle formula sin2(θ/2) = versin(θ)/2, derived by , that was used to construct such tables.
Sine, cosine, and versine of θ in terms of a unit circle, centered at O
The ordinary sine function () was sometimes historically called the sinus rectus ("vertical sine"), to contrast it with the versed sine (sinus versus). The meaning of these terms is apparent if one looks at the functions in the original context for their definition, a unit circle, shown at right. For a vertical chord AB of the unit circle, the sine of the angle θ (half the subtended angle) is the distance AC (half of the chord). On the other hand, the versed sine of θ is the distance CD from the center of the chord to the center of the arc. Thus, the sum of cos(θ) = OC and versin(θ) = CD is the radius OD = 1. Illustrated this way, the sine is vertical (rectus, lit. "straight") while the versine is horizontal (versus, lit. "turned against, out-of-place"); both are distances from C to the circle.
This figure also illustrates the reason why the versine was sometimes called the sagitta, Latin for , from the Arabic usage sahem of the same meaning. This itself comes from the Indian word 'sara' (arrow) that was commonly used to refer to "". If the arc ADB is viewed as a "" and the chord AB as its "string", then the versine CD is clearly the "arrow shaft".
In further keeping with the interpretation of the sine as "vertical" and the versed sine as "horizontal", sagitta is also an obsolete synonym for the
(the horizontal axis of a graph).
One period (0 & θ & π/2) of a versine or, more commonly, a haversine waveform is also commonly used in
as the shape of a
or a , because it smoothly ( in value and ) "turns on" from
(for haversine) and back to zero. In these applications, it is given yet another name:
Comparison of the versine function with three approximations to the versine functions, for angles ranging from 0 to 2 Pi
Comparison of the versine function with three approximations to the versine functions, for angles ranging from 0 to Pi/2
When the versine v is small in comparison to the radius r, it may be approximated from the half-chord length L (the distance AC shown above) by the formula
Alternatively, if the versine is small and the versine, radius, and half-chord length are known, they may be used to estimate the arc length s (AD in the figure above) by the formula
This formula was known to the Chinese mathematician , and a more accurate formula also involving the sagitta was developed two centuries later by .
A more accurate approximation used in engineering is
The term versine is also sometimes used to describe deviations from straightness in an arbitrary planar curve, of which the above circle is a special case. Given a chord between two points in a curve, the perpendicular distance v from the chord to the curve (usually at the chord midpoint) is called a versine measurement. For a straight line, the versine of any chord is zero, so this measurement characterizes the straightness of the curve. In the
as the chord length L goes to zero, the ratio 8v/L2 goes to the instantaneous . This usage is especially common in , where it describes measurements of the straightness of the
and it is the basis of the
for rail surveying. The term '' (often abbreviated sag) is used similarly in , for describing the surfaces of
(3rd ed.). Oxford University Press. September 2005.
Boyer, Carl B. (1991). A History of Mathematics (2nd ed.). New York: .
Miller, J. .
Calvert, James B. .
(3rd ed.). Oxford University Press. September 2005. Cites coinage by Prof. Jas. Inman, D. D., in his Navigation and Nautical Astronomy, 3rd ed. (1835).
Woodward, Ernest (1978), , Research & Education Assoc., p. 359,  .
Needham, Joseph (1959), , Cambridge University Press, p. 39,  .
Boardman, Harry (1930), Table For Use in Computing Arcs, Chords and Versines, Chicago Bridge and Iron Company, p. 32.
Nair, Bhaskaran (1972). "Track measurement systems—concepts and techniques". Rail International 3 (3): 159–166.  .
: Hidden categories:From Rosetta Code
Haversine formula
You are encouraged to
according to the task description, using any language you may know.
This page uses content from . The original article was at . The list of authors can be seen in the . , the text of Wikipedia
under the . (See links for details on variance)
The haversine formula is an equation important in navigation,
giving great-circle distances between two points on a sphere
from their longitudes and latitudes.
It is a special case of a more general formula in spherical trigonometry,
the law of haversines, relating the sides and angles of spherical "triangles".
Task: Implement a great-circle distance function, or use a library function,
to show the great-circle distance between Nashville International Airport (BNA)
in Nashville, TN, USA: N 36°7.2', W 86°40.2' (36.12, -86.67)
and Los Angeles International Airport (LAX) in Los Angeles, CA, USA: N 33°56.4', W 118°24.0' (33.94, -118.40).
User Kaimbridge clarified on the Talk page:
-- 6371.0 km is the authalic radius based on/extract
-- 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 -& 1655
6372.8 km -& 5340
(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).
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 M 
-- 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));
return 2.0 * Earth_Radius * ArcSin (Sqrt (a * a + Cos (lat1) * Cos (lat2) * b * b));
end Great_Circle_D 
-- 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_295_684_886127;
return (Deg + Min/60.0 + Sec/3600.0) * Pi_Over_180;
end DMS_To_R 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_F
Translation of:
Works with:
version Revision 1.
Works with:
version Any - tested with release .
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)))
dist: 2887.3 km (1794.1 mi.)
 #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. / 180) typedef d = double fundist(
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*Rend // end of [dist] implementmain0((*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] 
dist: 2887.3 km (1794.1 mi.)
, % GreatCircleDist(36.12, 33.94, -86.67, -118.40, 6372.8, &km&) GreatCircleDist(La1, La2, Lo1, Lo2, R, U) { return, 2 * R * ((Hs(Rad(La2 - La1)) + (Rad(La1)) * (Rad(La2)) * Hs(Rad(Lo2 - Lo1))))
U} Hs(n) { return, (1 - (n)) / 2} Rad(Deg) { return, Deg * 4 * (1) / 180}
 # syntax: GAWK -f HAVERSINE_FORMULA.AWK# converted from PythonBEGIN {
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.)} 
Works with:
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&
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
Distance =
#include &stdio.h&#include &stdlib.h&#include &math.h& #define R 6371#define TO_RAD (3. / 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 = (th1) - (th2); dx = (ph1) * (th1) - (th2); dy = (ph1) * (th1); return ((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 */ (&dist: %.1f km (%.1f mi.)\n&, d, d / 1.609344);  return 0;}
Translation of:
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: 0711 
Translation of:
 (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});=& 071106 
Translation of:
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)
(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)))
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.
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));}
Haversine distance: 2887.3 km
An alternate direct implementation of the haversine formula as shown at . 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));}
Great circle distance: 2887.3 km
% 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
:pi()/180, R
% In kilometers Diff_Lat
(Lat2 - Lat1)*V ,
(Long2 - Long1)*V,
Lat1*V, NLong
:sin(Diff_Lat/2) * :sin(Diff_Lat/2) + :sin(Diff_Long/2) * :sin(Diff_Long/2) * :cos(NLat) * :cos(NLong), C
2 * :asin(:sqrt(A)), R*C. 
% Implemented by Claudio Larini PROGRAM HAVERSINE_DEMO !$DOUBLE CONST DIAMETER=12745.6 FUNCTION DEG2RAD(X)
DEG2RAD=X*π/180END 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)*DIAMETEREND 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 71741 km, while using single-precision variable output is
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
&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-
return asin(norm(delta)/(r1+r2))*(r1+r2);
endfunction
&function overwrite rearth (x) := 6372.8*km$
&esdist(TNA,LAX)-&km
Translation of:
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&)))
Translation of:
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 ]
2biv.haversininv R_earth * ;
( scratchpad ) { 36.12 -86.67 } { 33.94 -118.4 } haversineDist .07113
Based on the Fortran and Groovy versions.
#APPTYPE CONSOLE  &Distance = &, Haversine(36.12, -86.67, 33.94, -118.4), & km&PAUSE FUNCTION Haversine(DegLat1
radius = 6372.8
= D2R(DegLat2 - DegLat1)
= D2R(DegLon2 - DegLon1)
= D2R(DegLat1)
= D2R(DegLat2)
= (dLat / 2) * (dLat / 2) + (dLon / 2) * (dLon / 2) * (lat1) * (lat2)
= 2 * ASIN(SQRT(a))
RETURN radius * c FUNCTION 
Distance = 0711 km
Press any key to continue...
: s&f s&d d&f ;: deg&rad 433e-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)
deg&rad fcos f*
( dLon^2 dLat2 lat1*lat2)
( 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.
 program exampleimplicit nonereal :: d d = haversine(36.12,-86.67,33.94,-118.40) ! BNA to LAXprint '(A,F9.4,A)', 'distance: ',d,' km' ! distance:
km contains 
function to_radian(degree) result(rad)
! degrees to radians
real,intent(in) :: degree
real :: rad,pi 
pi = 4*atan(1.0)
! exploit intrinsic atan to generate pi
rad = degree*pi/180
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 
 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&] 
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
(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:
 use navigation.frink d = earthDistance[36.12 deg North, 86.67 deg West, 33.94 deg North, 118.40 deg West]println[d-& &km&] 
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)) )
package main import (
&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)))}
haversine(lat1, lon1, lat2, lon2) {
R = 6372.8
// In kilometers
dLat = .toRadians(lat2 - lat1)
dLon = .toRadians(lon2 - lon1)
lat1 = .toRadians(lat1)
lat2 = .toRadians(lat2) 
a = .sin(dLat / 2) * .sin(dLat / 2) + .sin(dLon / 2) * .sin(dLon / 2) * .cos(lat1) * .cos(lat2)
c = 2 * .asin(.sqrt(a))
R * c} haversine(36.12, -86.67, 33.94, -118.40) & 0711
import Text.Printf -- The haversine of an angle.hsin t = let u =
(t/2) in u*u -- The distance between two points, given by latitude and longtitude, on a-- circle.
The points are specified in radians.distRad radius (lat1, lng1) (lat2, lng2) =
let hlat = hsin (lat2 - lat1)
hlng = hsin (lng2 - lng1)
(hlat +
lat2 * hlng)
in 2 * radius *
( 1.0 root) -- The distance between two points, given by latitude and longtitude, on a-- circle.
The points are specified in degrees.distDeg radius p1 p2 = distRad radius (deg2rad p1) (deg2rad p2)
where deg2rad (t, u) = (d2r t, d2r u)
d2r t = t *
/ 180 -- The approximate distance, in kilometers, between two points on Earth.
-- The latitude and longtitude are assumed to be in degrees.earthDist = distDeg 6372.8 main = do
let bna = (36.12,
-86.67)
lax = (33.94, -118.40)
dst = earthDist bna lax ::
printf &The distance between BNA and LAX is about %0.f km.\n& dst
The distance between BNA and LAX is about 2887 km.
Translation of:
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 conversionend 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 * 6371end
BNA to LAX is 2886 km (1793 miles)
require 'trig'haversin=: 0.5 * 1 - cosRearth=: 6372.8haversineDist=: 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.42887.26
Translation of:
public class Haversine {
public static final double R = 6372.8; // In kilometers
public static double haversine(double lat1, double lon1, double lat2, double lon2) {
double dLat = .toRadians(lat2 - lat1);
double dLon = .toRadians(lon2 - lon1);
lat1 = .toRadians(lat1);
lat2 = .toRadians(lat2); 
double a = .sin(dLat / 2) * .sin(dLat / 2) + .sin(dLon / 2) * .sin(dLon / 2) * .cos(lat1) * .cos(lat2);
double c = 2 * .asin(.sqrt(a));
return R * c;
public static void main([] args) {
.out.println(haversine(36.12, -86.67, 33.94, -118.40));
}}
Translation of:
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));
def haversine(lat1;lon1; lat2;lon2):
def radians: . * (1|atan)/45;
def sind: radians|
def cosd: radians|
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) ;
haversine(36.12; -86.67; 33.94; -118.4)
julia& haversine(lat1,lon1,lat2,lon2) = 2 * 6372.8 * asin(sqrt(sind((lat2-lat1)/2)^2 + cosd(lat1) * cosd(lat2) * sind((lon2 - lon1)/2)^2))# method added to generic function haversine julia& haversine(36.12,-86.67,33.94,-118.4)071106
print &Haversine distance: &; using( &####.###########&, havDist( 36.12, -86.67, 33.94, -118.4)); & km.&endfunction havDist( th1, ph1, th2, ph2)
= acs(-1)/180
= 2 * 6372.8
= degtorad
* (ph1 - ph2)
= degtorad
= degtorad
= sin( th1) - sin( th2)
= cos( LgD) * cos( th1) - cos( th2)
= sin( LgD) * cos( th1)
= asn( ( dx^2 +dy^2 +dz^2)^0.5 /2) *diameterend function
Haversine distance: 0711
Inputs assumed in degrees. Sin and Haversine expect
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]] 
distance[{36.12, -86.67}, {33.94, -118.4}]
function rad = radians(degree) % degrees to radians
rad = degree .*
/ 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 = ((dlat./2)).^2 + (lat1) .* (lat2) .* ((dlon./2)).^2;
c = 2 .* ((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
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/BNA/
http://www./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)),/* 13624 */
П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
as an average radius of the Earth, or
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).
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 С/П
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 С/П
import math proc radians(x): float = x * Pi / 180 proc haversine(lat1, lon1, lat2, lon2): float =
const 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)*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)
2.1115e+03
Works with oo2c version2
 MODULE HIMPORT
LRealMath,
PROCEDURE Distance(lat1,lon1,lat2,lon2: LONGREAL): LONGREAL;
r = ; (* Earth radius as LONGREAL *)
to_radians = LRealMath.pi / 180.0D0;
d,ph1,th1,th2: LONGREAL;
dz,dx,dy: LONGREAL;
d := lon1 - lon2;
ph1 := d * to_
th1 := lat1 * to_
th2 := lat2 * to_ 
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 *
END DBEGIN
Out.LongRealFix(Distance(36.12,-86.67,33.94,-118.4),6,10);Out.LnEND Haversines. 
 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));
}} 
distance: 2886.44
+ (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; 
double dLat = lat2rad - lat1
double dLon = lon2rad - lon1 
double a = (dLat/2) * (dLat/2) + (dLon/2) * (dLon/2) * (lat1rad) * (lat2rad);
double c = 2 * ((a));
double R = 6372.8;
return R *}
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. *.
1.let radians_of_degrees = ( *. ) (pi /. 180.)let haversin theta = 0.5 *. (1. -.
theta) (* The angle type can track radians or degrees, which I'll use for automatic conversion. *)type angle = Deg of
| Rad of 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: ; lng:
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 *.
2. *. radius *.
( 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.8and bna = LatLong.of_angles (Deg 36.12) (Deg (-86.67))and lax = LatLong.of_angles (Deg 33.94) (Deg (-118.4))inearth_dist bna lax;;
If the above is fed to the REPL, the last line will produce this:
- : float = 0711102
func: haversine(lat1, lon1, lat2, lon2){| lat lon | 
lat2 lat1 - asRadian -&lat
lon2 lon1 - asRadian -&lon 
lon 2 / sin sq lat1 asRadian cos * lat2 asRadian cos *
lat 2 / sin sq + sqrt asin 2 * 6372.8 *} haversine(36.12, -86.67, 33.94, -118.40) println
Translation of:
The rxmath library provides the required functions.
/*REXX pgm calculates distance between Nashville & Los Angles airports. */say & Nashville:
7.2', west
-86.67?&say &Los Angles:
north 33? 56.4', west 118? 24.0'
33.94?, -118.40?&saydist=surfaceDistance(36.12,
-118.4)kdist=format(dist/1
/*show 2 digs past decimal point.*/mdist=format(dist/1.609344,,2)
*/ndist=format(mdist*5280/6076.1,,2)
*/say ' distance between=
& kilometers,&say '
& statute miles,&say '
& 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
Nashville:
7.2', west
Los Angles:
north 33? 56.4', west 118? 24.0'
33.94?, -118.40?
distance between=
kilometers,
statute miles,
nautical or air miles.
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)
Works with:
Program HaversineDemo(output); uses
Math; function haversineDist(th1, ph1, th2, ph2: double): double;
diameter = 2 * 6372.8;
dx, dy, dz: double;
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.
Haversine distance: 2887.26 km.
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); # 9822
class POI {
private $latitude;
private $longitude;
public function __construct($latitude, $longitude) {
$this-&latitude = ($latitude);
$this-&longitude = ($longitude);
public function getLatitude() return $this-&latitude;
public function getLongitude() return $this-&longitude;
public function getDistanceInMetersTo(POI $other) {
$radiusOfEarth = 6371000;// Earth's radius in meters.
$diffLatitude = $other-&getLatitude() - $this-&latitude;
$diffLongitude = $other-&getLongitude() - $this-&longitude;
$a = ($diffLatitude / 2) * ($diffLatitude / 2) +
($this-&latitude) * ($other-&getLatitude()) *
($diffLongitude / 2) * ($diffLongitude / 2);
$c = 2 * (($a));
$distance = $radiusOfEarth * $c;
return $distance;
}}
$user = new POI($_GET[&latitude&], $_GET[&longitude&]);$poi = new POI(19,69276, -98,84350); // Piramide del Sol, Mexicoecho $user-&getDistanceInMetersTo($poi);
(scl 12)(load &@lib/math.l&) (de haversine (Th1 Ph1 Th2 Ph2)
Ph1 (*/ (- Ph1 Ph2) pi 180.0)
Th1 (*/ Th1 pi 180.0)
Th2 (*/ Th2 pi 180.0) )
(DX (- (*/ (cos Ph1) (cos Th1) 1.0) (cos Th2))
DY (*/ (sin Ph1) (cos Th1) 1.0)
DZ (- (sin Th1) (sin Th2)) )
(* `(* 2 6371)
(sqrt (+ (* DX DX) (* DY DY) (* DZ DZ)))
2 ) ) ) ) )
&Haversine distance: &
(round (haversine 36.12 -86.67 33.94 -118.4))
Haversine distance: 2,886.444 km
test: procedure options (main); /* 12 January 2014.
Derived from Fortran version */ 
d = haversine(36.12, -86.67, 33.94, -118.40);
/* BNA to LAX */
put edit ( 'distance: ', d, ' km') (A, F(10,3)); /* distance:
km */  degrees_to_radians: procedure (degree) returns (float);
declare degree
declare pi float (15) initial ( (4*atan(1.0d0)) ); 
return ( degree*pi/180 );end degrees_to_ haversine: procedure (deglat1, deglon1, deglat2, deglon2) returns (float);
declare (deglat1, deglon1, deglat2, deglon2)
declare (a, c, dlat, dlon, lat1, lat2)
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 ); 
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)071106&&&
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/BNA/#
http://www./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:
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)) 
Translation of:
-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
$lon2 $lon1 - toRadians
$lat1 toRadians
$lat2 toRadians
as $lat2 
$lat1 cos *
$lat2 cos * +
$R $c *} -118.40 33.94 -86.67 36.12 haversine &haversine: %.15g\n& print
haversine: 0711
REXX doen't have most of the higher math functions, so they are included here as subroutines.
/*REXX pgm calculates distance between Nashville & Los Angles airports. */say & Nashville:
7.2', west
-86.67?&say &Los Angles:
north 33? 56.4', west 118? 24.0'
33.94?, -118.40?&say$.=
/*set defaults for subroutines.
*/dist=surfaceDistance(36.12,
-118.4)kdist=format(dist/1
/*show 2 digs past decimal point.*/mdist=format(dist/1.609344,,2)
*/ndist=format(mdist*5280/6076.1,,2)
*/say ' distance between=
& kilometers,&say '
& statute miles,&say '
& nautical or air miles.&exit
/*stick a fork in it, we're done.*//*──────────────────────────────────SURFACEDISTANCE subroutine──────────*/surfaceDistance: arg th1,ph1,th2,ph2
/*use haversine formula for dist.*/numeric digits digits()*2
/*double the number of digits.
*/radius = 6372.8
/*earth's mean radius in km
ph1 = d2r(ph1-ph2)
/*convert degs──?radians & reduce*/
ph2 = d2r(ph2)
th1 = d2r(th1)
th2 = d2r(th2)
x = cos(ph1) * cos(th1) - cos(th2)
y = sin(ph1) * cos(th1)
z = sin(th1) - sin(th2)return radius * 2 * aSin(sqrt(x**2+y**2+z**2)/2 )/*═════════════════════════════general 1-line subs══════════════════════*/d2d: return arg(1) // 360
/*normalize degrees.
*/d2r: return r2r(arg(1)*pi() / 180)
/*normalize and convert deg──?rad*/r2d: return d2d((arg(1)*180 / pi()))
/*normalize and convert rad──?deg*/r2r: return arg(1) // (2*pi())
/*normalize radians.
return word(arg(1),1)
/*pick the first of two words.
if $.pipi==''
then $.pipi=$pi();
return $.pipi
/*return π.*/ aCos: procedure expose $.; arg if x&-1|x&1 then call $81r -1,1,x,&ACOS&
return .5*pi()-aSin(x)
says arg is out of range,*/
and it isn't included here.*/aSin: procedure expose $.;
parse arg x
if x&-1 | x&1
then call $81r -1,1,x,&ASIN&;
if abs(x)&=.7
then return sign(x)*aCos(sqrt(1-s),'-ASIN')
do j=2 by 2;
o=o*s*(j-1)/j;
z=z+o/(j+1)
if z=p then leave;
return z cos:
procedure expose $.;
x=r2r(x);
a=abs(x)
numeric fuzz min(9,digits()-9);
if a=pi()
then return -1
if a=pi()/2 | a=2*pi() then return 0;
if a=pi()/3 then return .5
if a=2*pi()/3
then return -.5;
return .sinCos(1,1,-1) sin:
procedure expose $.;
x=r2r(x);
numeric fuzz
min(5,digits()-3)
if abs(x)=pi()
then return 0;
return .sinCos(x,x,1) .sinCos: parse arg z,_,i; x=x*x; p=z; do k=2 by 2; _=-_*x/(k*(k+i)); z=z+_
if z=p then leave; p=z;
/*used by SIN & COS.*/ sqrt: procedure;
if x=0 then return 0;
d=digits()
numeric digits 11; g=.sqrtGuess(); do j=0 while p&9; m.j=p; p=p%2+1; end
do k=j+5 to 0 by -1; if m.k&11 then numeric digits m.k; g=.5*(g+x/g);end
numeric digits
return g/1.sqrtGuess:
numeric form;
parse value format(x,2,1,,0) 'E0' with g 'E' _ .;
g*.5'E'_%2 $pi: return ,'3.5923078'||,''||,'8196'
┌───────────────────────────────────────────────────────────────────────┐
│ A note on built-in functions.
REXX doesn't have a lot of mathmatical │
(particularly) trigomentric
functions,
so REXX programmers have │
│ to write their own.
Usually, this is done once, or most likely,
│ is borrowed from another program.
Knowing this, the one that is used │
│ has a lot of boilerplate in it.
Once coded and throughly debugged, I │
│ put those commonly-used subroutines into the
&1-line sub&
│ Programming note:
&general 1-line&
subroutines are taken from
│ other programs that I wrote, but I broke up their one line of source
│ so it can be viewed without shifting the viewing window.
[which won't happen here]
just shows an error telling
│ the legal range for
(in this case:
-1 ──? +1).
│ Similarly,
function checks for a negative argument
│ [which again, won't happen here].
│ The pi constant (as used here) is actually a much more robust function│
│ and will return up to one million digits in the real version.
│ One bad side effect is that, like a automobile without a hood, you see│
│ all the dirty stuff going on.
Also, don't visit a sausage factory. │
└───────────────────────────────────────────────────────────────────────┘
Nashville:
7.2', west
Los Angles:
north 33? 56.4', west 118? 24.0'
33.94?, -118.40?
distance between=
kilometers,
statute miles,
nautical or air miles.
include Math Radius = 6371
# rough radius of the Earth, in kilometers def spherical_distance(start_coords, end_coords)
lat1, long1 = deg2rad *start_coords
lat2, long2 = deg2rad *end_coords
2 * Radius * asin(sqrt(sin((lat2-lat1)/2)**2 + cos(lat1) * cos(lat2) * sin((long2 - long1)/2)**2))end def deg2rad(lat, long)
[lat * PI / 180, long * PI / 180]end bna = [36.12, -86.67]lax = [33.94, -118.4] puts &%.1f& % spherical_distance(bna, lax)
D2R = atn(1)/45
= 2 * 6372.8Lg1m2
= ((-86.67)-(-118.4)) * D2RLt1
= 36.12 * D2R ' degrees to radLt2
= 33.94 * D2R
= sin(Lt1) - sin(Lt2)
= cos(Lg1m2) * cos(Lt1) - cos(Lt2)
= sin(Lg1m2) * cos(Lt1)
hDist = asn((dx^2 + dy^2 + dz^2)^0.5 /2) * diamprint &Haversine distance: &;using(&####.#############&,hDist);& km.&  'Tips:
( 36 deg 7 min 12 sec ) = print 36+(7/60)+(12/3600).
Produces: 36.12 deg ' '
&36.12,-86.67&
(no quotes).
Click map, '
satellite, center the pin &A&, zoom in, and see airport.
Extra: in &Get '
Directions& enter
36.12,-86.66999 and see pin &B& about one meter away. '
(. km., or 35.37 in.) ' '
This code also works in Liberty BASIC.Output Haversine distance: 071104 km.
 options %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.;
%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
%end; 
data _null_;
convert = &convert;
lat1 = &lat1 *
lat2 = &lat2 *
long1 = &long1 *
long2 = &long2 * 
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);  
Distance is
math._  Haversine {
R = 6372.8
//radius in km 
haversine(lat1:Double, lon1:Double, lat2:Double, lon2:Double)={
dLat=(lat2 - lat1).toRadians
dLon=(lon2 - lon1).toRadians 
a = pow(sin(dLat/2),2) + pow(sin(dLon/2),2) * cos(lat1.toRadians) * cos(lat2.toRadians)
c = 2 * asin(sqrt(a))
} 
main(args: Array[String]): Unit = {
println(haversine(36.12, -86.67, 33.94, -118.40))
}}
(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)); 37984
$ 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
var float: distance is 0.0;
const float: EarthRadius is 6372.8;
# Average great-elliptic or great-circle radius in kilometers
distance := 2.0 * EarthRadius * asin(sqrt(sin(0.5 * (latitude2 - latitude1)) ** 2 +
cos(latitude1) * cos(latitude2) *
sin(0.5 * (longitude2 - longitude1)) ** 2)); const func float: degToRad (in float: degrees) is
return degrees * 0.; const proc: main is func
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);
Translation of:
class EarthPoint(lat, lon) { 
static earth_radius = 6371;
# mean earth radius
static radian_ratio = Math.pi/180; 
# accessors for radians
method latR { self.lat * radian_ratio };
method lonR { self.lon * radian_ratio }; 
method haversine_dist(p is EarthPoint) {
var arc = __CLASS__.new(
self.lat - p.lat,
self.lon - p.lon,
); 
var a = [ Math.pow(Math.sin(arc.latR / 2), 2),
Math.pow(Math.sin(arc.lonR / 2), 2)
* Math.cos(self.latR) * Math.cos(p.latR),
].sum; 
earth_radius * Math.asin(Math.sqrt(a)) * 2;
}} var BNA = EarthPoint.new(lat: 36.12, lon: -86.67);var LAX = EarthPoint.new(lat: 33.94, lon: -118.4); say BNA.haversine_dist(LAX);
Translation of:
import Foundation func haversine(lat1:Double, lon1:Double, lat2:Double, lon2:Double) -& Double {
let lat1rad = lat1 * M_PI/180
let lon1rad = lon1 * M_PI/180
let lat2rad = lat2 * M_PI/180
let lon2rad = lon2 * M_PI/180 
let dLat = lat2rad - lat1rad
let dLon = lon2rad - lon1rad
let a = sin(dLat/2) * sin(dLat/2) + sin(dLon/2) * sin(dLon/2) * cos(lat1rad) * cos(lat2rad)
let c = 2 * asin(sqrt(a))
let R = 6372.8 
return R * c} println(haversine(36.12, -86.67, 33.94, -118.40))
Translation of:
package require Tcl 8.5proc 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]]
distance=2887.3 km
'Sets decimal display to 32 places (0+.1^56)
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&
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
km OK 
Assemble with tasm /m /l; tlink /t
;.com files start here0100
;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
;push real onto FPU stack0107
D8 26 018Dr
;subtract real from top of stack (st(0) = st)010B
;(1.0-cos(st)) / 2.0010E
D9 06 0199r
;repeat for longitudes0112
D8 26 0195r
;st(1)=L st=Lons0119
D9 06 018Dr
;replace st with its cosine011F
D9 06 0191r
;st=cos(Lat2); st(1)=cos(Lat1); st(2)=L st(3)=Lons0125
;st=cos(Lat2)*cos(Lat1); st(1)=L st(2)=Lons0127
;st=cos(Lat2)*cos(Lat1)*L st(1)=Lons0129
;st=cos(Lat2)*cos(Lat1)*Lats + Lons012B
;replace st with its square root
;asin(x) = atan(x/sqrt(1-x^2))012D
;duplicate tos012F
;get 1.00133
;1 - x^20135
;sqrt(1-x^2)0137
;take atan(st(1)/st)0139
D8 0E 019Dr
;*2.0*Radius 
;Display value in FPU's top of stack (st)
;places before
; and after decimal point
;&=& allows scaler to be redefined, unlike equ
;repeat block &after& times
;scaler now = 10^after 013D
dword ptr scaler;use stack for convenient memory location0140
67| DA 0C 24
dword ptr [esp] ;st:= st*scaler0144
67| DB 1C 24
dword ptr [esp] ;round st to nearest integer0148
; and put it into eax 014A
66| BB 0000000A
;set up for idiv instruction0150
cx, before+after;set up loop counter0153
;co i.e: edx:= 00155
;eax:= edx:eax/ remainder in edx0158
;save least significant digit on stack0159
;cx--; loop back if not zero 015B
cl, before+after;(ch=0)015D
;used to suppress leading zeros015F
;get digit0160
;turn off suppression if not a zero0162
cl, after+1
;is digit immediately to left of decimal point?0165
;skip if not0167
;turn off leading zero suppression0168
;if leading zero then ' ' else add 0016A
bl, bl016C
al, ' '0170
;display character in al register0172
cl, after+1
;is digit immediately to left of decimal point?0175
;skip if not0177
;display decimal point0179
;loop until all digits displayed017D
;return to OS 017E
Haversine:
;return (1.0-Cos(Ang)) / 2.0 in st017E
D8 36 0189r
ret 0189
;36.12*pi/1800191
;33.94*pi/1800195
;-86.67*pi/1800199
;-118.40*pi/180019D
Radius2 dd
;6372.8 average radius of Earth (km) times 2
;(TASM isn't smart enough to do floating point constant calculations)
start 
include c:\cxpl\
\intrinsic 'code' declarations func real Haversine(Ang);real Areturn (1.0-Cos(Ang)) / 2.0; func real Dist(Lat1, Lat2, Lon1, Lon2); \Great circle distancereal 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.0.0;
\degrees to radiansRlOut(0, Dist(36.12*D2R, 33.94*D2R, -86.67*D2R, -118.40*D2R ));
Translation of:
haversine(36.12, -86.67, 33.94, -118.40).println(); fcn haversine(Lat1, Long1, Lat2, Long2){
const R = 6372.8;
Diff_Lat  := (Lat2
- Lat1) .toRad();
Diff_Long := (Long2 - Long1).toRad();
 := Lat1.toRad();
 := Lat2.toRad();
 := (Diff_Lat/2) .sin().pow(2) +
(Diff_Long/2).sin().pow(2) *
NLat.cos() * NLong.cos();
 := 2.0 * A.sqrt().asin();