Haversine formula

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>

1. include <stdlib.h>
2. include <math.h>

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

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

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

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

return 0; }</lang>

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

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 rads = Math.PI / 180
 def R = 6372.8
 // In kilometers
 def dLat = (lat2 - lat1) * rads
 def dLon = (lon2 - lon1) * rads
 lat1 = lat1 * rads
 lat2 = lat2 * rads

 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

Works with: Free_Pascal

Library: Math

<lang pascal>Program HaversineDemo(output);

uses

 Math;

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

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

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

begin

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

end.</lang> Output:

Haversine distance: 2887.26 km.

Perl 6

<lang perl6>class EarthPoint {

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

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

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

       method haversine-dist(EarthPoint $p) {

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

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

               return $earth_radius * $c;
       }

}

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

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

Python

Translation of: Groovy

<lang python>>>> import math >>> def haversine(lat1, lon1, lat2, lon2):

 rads = math.pi / 180
 R = 6372.8
 # In kilometers
 dLat = (lat2 - lat1) * rads
 dLon = (lon2 - lon1) * rads
 lat1 = lat1 * rads
 lat2 = lat2 * rads

 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

Translation of: Groovy

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

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

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

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

}

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

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

distance=2887.3 km