Calculate distance between two points on a globe(http://www.codecodex.com/wiki/Calculate_distance_between_two_points_on_a_globe)

c++

/// @brief The usual PI/180 constant
static const double DEG_TO_RAD = 0.017453292519943295769236907684886;
/// @brief Earth's quatratic mean radius for WGS-84
static const double EARTH_RADIUS_IN_METERS = 6372797.560856;

/** @brief Computes the arc, in radian, between two WGS-84 positions.
  *
  * The result is equal to <code>Distance(from,to)/EARTH_RADIUS_IN_METERS</code>
  *    <code>= 2*asin(sqrt(h(d/EARTH_RADIUS_IN_METERS )))</code>
  *
  * where:<ul>
  *    <li>d is the distance in meters between 'from' and 'to' positions.</li>
  *    <li>h is the haversine function: <code>h(x)=sin²(x/2)</code></li>
  * </ul>
  *
  * The haversine formula gives:
  *    <code>h(d/R) = h(from.lat-to.lat)+h(from.lon-to.lon)+cos(from.lat)*cos(to.lat)</code>
  *
  * @sa http://en.wikipedia.org/wiki/Law_of_haversines
  */
double ArcInRadians(const Position& from, const Position& to) {
    double latitudeArc  = (from.lat - to.lat) * DEG_TO_RAD;
    double longitudeArc = (from.lon - to.lon) * DEG_TO_RAD;
    double latitudeH = sin(latitudeArc * 0.5);
    latitudeH *= latitudeH;
    double lontitudeH = sin(longitudeArc * 0.5);
    lontitudeH *= lontitudeH;
    double tmp = cos(from.lat*DEG_TO_RAD) * cos(to.lat*DEG_TO_RAD);
    return 2.0 * asin(sqrt(latitudeH + tmp*lontitudeH));
}

/** @brief Computes the distance, in meters, between two WGS-84 positions.
  *
  * The result is equal to <code>EARTH_RADIUS_IN_METERS*ArcInRadians(from,to)</code>
  *
  * @sa ArcInRadians
  */
double DistanceInMeters(const Position& from, const Position& to) {
    return EARTH_RADIUS_IN_METERS*ArcInRadians(from, to);
}

DistanceCalculator.java

import com.google.android.maps.GeoPoint;

public class DistanceCalculator {

   private double Radius;

   // R = earth's radius (mean radius = 6,371km)
   // Constructor
   DistanceCalculator(double R) {
      Radius = R;
   }

   public double CalculationByDistance(GeoPoint StartP, GeoPoint EndP) {
      double lat1 = StartP.getLatitudeE6()/1E6;
      double lat2 = EndP.getLatitudeE6()/1E6;
      double lon1 = StartP.getLongitudeE6()/1E6;
      double lon2 = EndP.getLongitudeE6()/1E6;
      double dLat = Math.toRadians(lat2-lat1);
      double dLon = Math.toRadians(lon2-lon1);
      double a = Math.sin(dLat/2) * Math.sin(dLat/2) +
         Math.cos(Math.toRadians(lat1)) * Math.cos(Math.toRadians(lat2)) *
         Math.sin(dLon/2) * Math.sin(dLon/2);
      double c = 2 * Math.asin(Math.sqrt(a));
      return Radius * c;
   }
}

haversine.rb

# haversine.rb
#
# haversine formula to compute the great circle distance between two points given their latitude and longitudes
#
# Copyright (C) 2008, 360VL, Inc
# Copyright (C) 2008, Landon Cox
#
# http://www.esawdust.com (Landon Cox)
# contact:
# http://www.esawdust.com/blog/businesscard/businesscard.html
#
# LICENSE: GNU Affero GPL v3
# The ruby implementation of the Haversine formula is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License version 3 as published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
# implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Affero General Public
# License version 3 for more details.  http://www.gnu.org/licenses/
#
# Landon Cox - 9/25/08
#
# Notes:
#
# translated into Ruby based on information contained in:
#   http://mathforum.org/library/drmath/view/51879.html  Doctors Rick and Peterson - 4/20/99
#   http://www.movable-type.co.uk/scripts/latlong.html
#   http://en.wikipedia.org/wiki/Haversine_formula
#
# This formula can compute accurate distances between two points given latitude and longitude, even for
# short distances.

# PI = 3.1415926535
RAD_PER_DEG = 0.017453293  #  PI/180

# the great circle distance d will be in whatever units R is in

Rmiles = 3956           # radius of the great circle in miles
Rkm = 6371              # radius in kilometers...some algorithms use 6367
Rfeet = Rmiles * 5282   # radius in feet
Rmeters = Rkm * 1000    # radius in meters

@distances = Hash.new   # this is global because if computing lots of track point distances, it didn't make
                        # sense to new a Hash each time over potentially 100's of thousands of points

=begin rdoc
  given two lat/lon points, compute the distance between the two points using the haversine formula
  the result will be a Hash of distances which are key'd by 'mi','km','ft', and 'm'
=end

def haversine_distance( lat1, lon1, lat2, lon2 )

  dlon = lon2 - lon1
  dlat = lat2 - lat1

  dlon_rad = dlon * RAD_PER_DEG
  dlat_rad = dlat * RAD_PER_DEG

  lat1_rad = lat1 * RAD_PER_DEG
  lon1_rad = lon1 * RAD_PER_DEG

  lat2_rad = lat2 * RAD_PER_DEG
  lon2_rad = lon2 * RAD_PER_DEG

  # puts "dlon: #{dlon}, dlon_rad: #{dlon_rad}, dlat: #{dlat}, dlat_rad: #{dlat_rad}"

  a = Math.sin(dlat_rad/2)**2 + Math.cos(lat1_rad) * Math.cos(lat2_rad) * Math.sin(dlon_rad/2)**2
  c = 2 * Math.asin( Math.sqrt(a))

  dMi = Rmiles * c          # delta between the two points in miles
  dKm = Rkm * c             # delta in kilometers
  dFeet = Rfeet * c         # delta in feet
  dMeters = Rmeters * c     # delta in meters

  @distances["mi"] = dMi
  @distances["km"] = dKm
  @distances["ft"] = dFeet
  @distances["m"] = dMeters
end

def test_haversine

  lon1 = -104.88544
  lat1 = 39.06546

  lon2 = -104.80
  lat2 = lat1

  haversine_distance( lat1, lon1, lat2, lon2 )

  puts "the distance from  #{lat1}, #{lon1} to #{lat2}, #{lon2} is:"
  puts "#{@distances['mi']} mi"
  puts "#{@distances['km']} km"
  puts "#{@distances['ft']} ft"
  puts "#{@distances['m']} m"

  if @distances['km'].to_s =~ /7\.376*/
    puts "Test: Success"
  else
    puts "Test: Failed"
  end

end

test_haversine

javascript

var R = 6371; // km
var dLat = (lat2-lat1)*Math.PI/180;
var dLon = (lon2-lon1)*Math.PI/180;
var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
        Math.cos(lat1*Math.PI/180) * Math.cos(lat2*Math.PI/180) *
        Math.sin(dLon/2) * Math.sin(dLon/2);
var c = 2 * Math.asin(Math.sqrt(a));
var d = R * c;

Objective C

// adapted from C++ code on this page
+ (CGFloat)directMetersFromCoordinate:(CLLocationCoordinate2D)from toCoordinate:(CLLocationCoordinate2D)to {

    static const double DEG_TO_RAD = 0.017453292519943295769236907684886;
    static const double EARTH_RADIUS_IN_METERS = 6372797.560856;

    double latitudeArc  = (from.latitude - to.latitude) * DEG_TO_RAD;
    double longitudeArc = (from.longitude - to.longitude) * DEG_TO_RAD;
    double latitudeH = sin(latitudeArc * 0.5);
    latitudeH *= latitudeH;
    double lontitudeH = sin(longitudeArc * 0.5);
    lontitudeH *= lontitudeH;
    double tmp = cos(from.latitude*DEG_TO_RAD) * cos(to.latitude*DEG_TO_RAD);
    return EARTH_RADIUS_IN_METERS * 2.0 * asin(sqrt(latitudeH + tmp*lontitudeH));
}

PostgreSQL PL

create or replace function distance_in_km (numeric,numeric,numeric,numeric) returns numeric as
$body$
use strict;
use Math::Trig qw(great_circle_distance deg2rad);

my $lat1 = shift(@_);
my $lon1 = shift(@_);
my $lat2 = shift(@_);
my $lon2 = shift(@_);

# Notice the 90 - latitude: phi zero is at the North Pole.
sub NESW { deg2rad($_[0]), deg2rad(90 - $_[1]) }
my @L = NESW( $lat1, $lon1 );
my @T = NESW( $lat2, $lon2 );
my $km = great_circle_distance(@L, @T, 6378);

return $km;
$body$
language plperlu strict security definer;

SELECT *, (
    6371 * acos(cos(radians(search_lat)) * cos(radians(lat) ) *
    cos(radians(lng) - radians(search_lng)) + sin(radians(search_lat)) * sin(radians(lat)))
) AS distance
FROM table
WHERE lat != search_lat AND lng != search_lng AND distance < 25
ORDER BY distance
FETCH 10 ONLY

Python

#coding:UTF-8
"""
  Python implementation of Haversine formula
  Copyright (C) <2009>  Bartek Górny <[email protected]>

  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program.  If not, see <http://www.gnu.org/licenses/>.
"""

import math

def recalculate_coordinate(val,  _as=None):
  """
    Accepts a coordinate as a tuple (degree, minutes, seconds)
    You can give only one of them (e.g. only minutes as a floating point number) and it will be duly
    recalculated into degrees, minutes and seconds.
    Return value can be specified as 'deg', 'min' or 'sec'; default return value is a proper coordinate tuple.
  """
  deg,  min,  sec = val
  # pass outstanding values from right to left
  min = (min or 0) + int(sec) / 60
  sec = sec % 60
  deg = (deg or 0) + int(min) / 60
  min = min % 60
  # pass decimal part from left to right
  dfrac,  dint = math.modf(deg)
  min = min + dfrac * 60
  deg = dint
  mfrac,  mint = math.modf(min)
  sec = sec + mfrac * 60
  min = mint
  if _as:
    sec = sec + min * 60 + deg * 3600
    if _as == 'sec': return sec
    if _as == 'min': return sec / 60
    if _as == 'deg': return sec / 3600
  return deg,  min,  sec


def points2distance(start,  end):
  """
    Calculate distance (in kilometers) between two points given as (long, latt) pairs
    based on Haversine formula (http://en.wikipedia.org/wiki/Haversine_formula).
    Implementation inspired by JavaScript implementation from http://www.movable-type.co.uk/scripts/latlong.html
    Accepts coordinates as tuples (deg, min, sec), but coordinates can be given in any form - e.g.
    can specify only minutes:
    (0, 3133.9333, 0)
    is interpreted as
    (52.0, 13.0, 55.998000000008687)
    which, not accidentally, is the lattitude of Warsaw, Poland.
  """
  start_long = math.radians(recalculate_coordinate(start[0],  'deg'))
  start_latt = math.radians(recalculate_coordinate(start[1],  'deg'))
  end_long = math.radians(recalculate_coordinate(end[0],  'deg'))
  end_latt = math.radians(recalculate_coordinate(end[1],  'deg'))
  d_latt = end_latt - start_latt
  d_long = end_long - start_long
  a = math.sin(d_latt/2)**2 + math.cos(start_latt) * math.cos(end_latt) * math.sin(d_long/2)**2
  c = 2 * math.asin(math.sqrt(a))
  return 6371 * c


if __name__ == '__main__':
 warsaw = ((21,  0,  30),  (52, 13, 56))
 cracow = ((19, 56, 18),  (50, 3, 41))
 print points2distance(warsaw,  cracow)