The haversine formula

from math import radians, cos, sin, asin, sqrt

def haversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees)
    """
    # convert decimal degrees to radians
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])

    # haversine formula
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a))
    r = 6371 # Radius of earth in kilometers. Use 3956 for miles
    return c * r

# The haversine formula is for calculating the great circle distance between two
# latitude/longitude coordinate pairs

haversine = function(lon1, lat1, lon2, lat2){
    # convert decimal degrees to radians
    lon1 = lon1 * pi / 180
    lon2 = lon2 * pi / 180
    lat1 = lat1 * pi / 180
    lat2 = lat2 * pi / 180
    # haversine formula
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * atan2(sqrt(a), sqrt(1-a))
    km = 6367 * c
    km
}

haversine.wrapper = function(lang.data,
      longitude="longitude",
      latitude="latitude"){
  # A version of the haversine function, which accepts a data.frame
  # as argument rather than a series of coordinate pairs
  # Usage example:
  # geodist <- log10(combn(nrow(langdata), 2, haversine.wrapper(langdata)))
  function(indices){
    haversine(lang.data[[longitude]][indices[1]],
      lang.data[[latitude]][indices[1]],
      lang.data[[longitude]][indices[2]],
      lang.data[[latitude]][indices[2]])
  }
}