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 haversine distance function, or use a library function, to show the haversine 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).
C
<lang 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; }</lang>
D
<lang d>import std.stdio, std.math;
real haversineDistance(in real dth1, in real dph1,
in real dth2, in real dph2) pure nothrow {
enum real R = 6371; enum real TO_RAD = PI / 180;
alias immutable(real) imr; imr ph1d = dph1 - dph2; imr ph1 = ph1d * TO_RAD; imr th1 = dth1 * TO_RAD; imr th2 = dth2 * TO_RAD;
imr dz = sin(th1) - sin(th2); imr dx = cos(ph1) * cos(th1) - cos(th2); imr dy = sin(ph1) * cos(th1); return asin(sqrt(dx ^^ 2 + dy ^^ 2 + dz ^^ 2) / 2) * 2 * R;
}
void main() {
immutable d = haversineDistance(36.12, -86.67, 33.94, -118.4);
writefln("Haversine distance: %.1f km", d);
}</lang> Output:
Haversine distance: 2886.4 km
Forth
<lang forth>: s>f s>d d>f ;
- deg>rad 174532925199433e-16 f* ;
- difference f- deg>rad 2 s>f f/ fsin fdup f* ;
- haversine ( lat1 lon1 lat2 lon2 -- haversine)
frot difference ( lat1 lat2 dLon^2) frot frot fover fover ( dLon^2 lat1 lat2 lat1 lat2) fswap difference ( dLon^2 lat1 lat2 dLat^2) fswap deg>rad fcos ( dLon^2 lat1 dLat^2 lat2) frot deg>rad fcos f* ( dLon^2 dLat2 lat1*lat2) frot f* f+ ( lat1*lat2*dLon^2+dLat^2) fsqrt fasin 127456 s>f f* 10 s>f f/ ( haversine)
36.12e -86.67e 33.94e -118.40e haversine cr f.</lang> Output:
2887.25995060711
Go
<lang go>package main
import (
"fmt" "math"
)
func haversine(θ float64) float64 {
return .5 * (1 - math.Cos(θ))
}
type pos struct {
φ float64 // latitude, radians ψ float64 // longitude, radians
}
func degPos(lat, lon float64) pos {
return pos{lat * math.Pi / 180, lon * math.Pi / 180}
}
const rEarth = 6372.8 // km
func hsDist(p1, p2 pos) float64 {
return 2 * rEarth * math.Asin(math.Sqrt(haversine(p2.φ-p1.φ)+
math.Cos(p1.φ)*math.Cos(p2.φ)*haversine(p2.ψ-p1.ψ)))
}
func main() {
fmt.Println(hsDist(degPos(36.12, -86.67), degPos(33.94, -118.40)))
}</lang> Output:
2887.2599506071097
Groovy
<lang Groovy>def haversine(lat1, lon1, lat2, lon2) {
def R = 6372.8 // In kilometers def dLat = Math.toRadians(lat2 - lat1) def dLon = Math.toRadians(lon2 - lon1) lat1 = Math.toRadians(lat1) lat2 = Math.toRadians(lat2)
def a = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.sin(dLon / 2) * Math.sin(dLon / 2) * Math.cos(lat1) * Math.cos(lat2) def c = 2 * Math.asin(Math.sqrt(a)) R * c
}
haversine(36.12, -86.67, 33.94, -118.40)
> 2887.25995060711</lang>
Pascal
<lang pascal>Program HaversineDemo(output);
uses
Math;
function haversineDist(th1, ph1, th2, ph2: double): double;
const diameter = 2 * 6372.8; var dx, dy, dz: double; begin ph1 := degtorad(ph1 - ph2); th1 := degtorad(th1); th2 := degtorad(th2); dz := sin(th1) - sin(th2); dx := cos(ph1) * cos(th1) - cos(th2); dy := sin(ph1) * cos(th1); haversineDist := arcsin(sqrt(dx**2 + dy**2 + dz**2) / 2) * diameter; end;
begin
writeln ('Haversine distance: ', haversineDist(36.12, -86.67, 33.94, -118.4):7:2, ' km.');
end.</lang> Output:
Haversine distance: 2887.26 km.
Perl 6
<lang perl6>class EarthPoint {
has $.lat; # latitude
has $.lon; # longitude
has $earth_radius = 6371; # mean earth radius
has $radian_ratio = pi / 180;
# accessors for radians
method latR { $.lat * $radian_ratio }
method lonR { $.lon * $radian_ratio }
method haversine-dist(EarthPoint $p) {
my EarthPoint $arc .= new(
lat => $!lat - $p.lat,
lon => $!lon - $p.lon );
my $a = sin($arc.latR/2) ** 2 + sin($arc.lonR/2) ** 2
* cos($.latR) * cos($p.latR);
my $c = 2 * asin( sqrt($a) );
return $earth_radius * $c;
}
}
my EarthPoint $BNA .= new(lat => 36.12, lon => -86.67); my EarthPoint $LAX .= new(lat => 33.94, lon => -118.4);
say $BNA.haversine-dist($LAX); # 2886.44444099822</lang>
Python
<lang python>>>> import math >>> def haversine(lat1, lon1, lat2, lon2):
R = 6372.8 # In kilometers dLat = math.radians(lat2 - lat1) dLon = math.radians(lon2 - lon1) lat1 = math.radians(lat1) lat2 = math.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
>>> haversine(36.12, -86.67, 33.94, -118.40) 2887.2599506071106 >>> </lang>
Tcl
<lang tcl>package require Tcl 8.5 proc haversineFormula {lat1 lon1 lat2 lon2} {
set rads [expr atan2(0,-1)/180] set R 6372.8 ;# In kilometers
set dLat [expr {($lat2-$lat1) * $rads}]
set dLon [expr {($lon2-$lon1) * $rads}]
set lat1 [expr {$lat1 * $rads}]
set lat2 [expr {$lat2 * $rads}]
set a [expr {sin($dLat/2)**2 + sin($dLon/2)**2*cos($lat1)*cos($lat2)}]
set c [expr {2*asin(sqrt($a))}]
return [expr {$R * $c}]
}
- Don't bother with too much inappropriate accuracy!
puts [format "distance=%.1f km" [haversineFormula 36.12 -86.67 33.94 -118.40]]</lang> Output:
distance=2887.3 km