以下代码用python2运行

#!/usr/bin/python
#coding=utf-8
from __future__ import print_function
import sys
reload(sys)
sys.setdefaultencoding("utf-8")
import math

def haversine_distance(lon1, lat1, lon2, lat2):
    # 将纬度和经度从度转换为弧度
    lat1_rad = math.radians(lat1)
    lon1_rad = math.radians(lon1)
    lat2_rad = math.radians(lat2)
    lon2_rad = math.radians(lon2)

    # 计算差值
    dlon = lon2_rad - lon1_rad
    dlat = lat2_rad - lat1_rad

    # 应用 Haversine 公式
    a = math.sin(dlat / 2)**2 + math.cos(lat1_rad) * math.cos(lat2_rad) * math.sin(dlon / 2)**2
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
    distance = 6371 * c * 1000  # 地球平均半径约为 6371 千米

    return distance
    
def vincenty(lon1, lat1, lon2, lat2):
    # 定义椭球体参数
    a = 6378137.0  # 长半轴（单位：米）
    f = 1 / 298.257223563  # 扁率
    b = (1 - f) * a  # 短半轴（单位：米）

    # 将经纬度从度数转换为弧度
    lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])

    # 差值
    L = lon2 - lon1
    U1 = math.atan((1 - f) * math.tan(lat1))
    U2 = math.atan((1 - f) * math.tan(lat2))
    sinU1, cosU1 = math.sin(U1), math.cos(U1)
    sinU2, cosU2 = math.sin(U2), math.cos(U2)

    # 初始设定
    lamb = L
    iter_limit = 100  # 最大迭代次数

    for _ in range(iter_limit):
        sin_lamb = math.sin(lamb)
        cos_lamb = math.cos(lamb)
        sin_sigma = math.sqrt((cosU2 * sin_lamb) ** 2 + (cosU1 * sinU2 - sinU1 * cosU2 * cos_lamb) ** 2)
        if sin_sigma == 0:
            return 0  # 重合点
        cos_sigma = sinU1 * sinU2 + cosU1 * cosU2 * cos_lamb
        sigma = math.atan2(sin_sigma, cos_sigma)
        sin_alpha = cosU1 * cosU2 * sin_lamb / sin_sigma
        cos2_alpha = 1 - sin_alpha ** 2
        cos2_sigma_m = cos_sigma - 2 * sinU1 * sinU2 / cos2_alpha
        if math.isnan(cos2_sigma_m):
            cos2_sigma_m = 0  # 赤道情况
        C = f / 16 * cos2_alpha * (4 + f * (4 - 3 * cos2_alpha))
        lamb_prev = lamb
        lamb = L + (1 - C) * f * sin_alpha * (
            sigma + C * sin_sigma * (
                cos2_sigma_m + C * cos_sigma * (
                    -1 + 2 * cos2_sigma_m ** 2)))

        if abs(lamb - lamb_prev) < 1e-12:
            break
    else:
        raise ValueError("Vincenty formula failed to converge")

    u2 = cos2_alpha * (a ** 2 - b ** 2) / (b ** 2)
    A = 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)))
    B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)))
    delta_sigma = B * sin_sigma * (
        cos2_sigma_m + B / 4 * (
            cos_sigma * (-1 + 2 * cos2_sigma_m ** 2) -
            B / 6 * cos2_sigma_m * (-3 + 4 * sin_sigma ** 2) *
            (-3 + 4 * cos2_sigma_m ** 2)))

    s = b * A * (sigma - delta_sigma)

    return s  # 距离，单位：米

# 示例：计算北京和上海之间的距离
beijing_lon, beijing_lat = 116.407396, 39.904200
shanghai_lon, shanghai_lat = 121.473701, 31.230416

distance = vincenty(beijing_lon, beijing_lat, shanghai_lon, shanghai_lat)
print("北京和上海之间的距离是 {:.2f} 米".format(distance))

haversine算法与vincenty的差距很小，对于500多米的近距离，误差在1米多。

将haversine用sql语句表达结果为：

select 6371*1000 * acos(sin(radians(a.latitude)) * sin(radians(b.latitude)) + 
           cos(radians(a.latitude)) * cos(radians( b.latitude)) * 
           cos(radians(b.longitude) - radians(a.longitude))) AS distance