postgresql计算球面距离的一个应用

开发有一个需求，需要计算两个不同经纬度之间的球面距离。写成函数如下：

create FUNCTION getdistance
(  
    i_latbegin REAL,i_lngbegin REAL,i_latend REAL,i_lngend REAL  
) 
 
RETURNS FLOAT  
AS  

$BODY$

--author  :kenyon
--created :2012-06-12
--purpose :calc the distance between diff lat and lng

declare

v_distance real;
v_earth_radius real;
v_radlatbegin real;
v_radlatend real;
v_radlatdiff real;
v_radlngdiff real;

BEGIN  

    --地球半径
    v_earth_radius:=6378.137;
     
    v_radlatbegin := i_latbegin * pi()/180.0;
    v_radlatend   := i_latend * pi()/180.0;
    v_radlatdiff  := v_radlatbegin - v_radlatend;
    v_radlngdiff  := i_lngbegin * pi()/180.0 - i_lngend * pi()/180.0;

    v_distance := 2 * asin(sqrt(power(sin(v_radlatdiff / 2),2) + cos(v_radlatbegin) * cos(v_radlatend) * power(sin(v_radlngdiff/2),2)));
    v_distance := v_distance * v_earth_radius*1000;
      
    RETURN v_distance; 
END;
$BODY$
  
LANGUAGE 'plpgsql' VOLATILE;

目前只有这个应用用到计算空间距离，故暂时没去安装和应用postgis了。上述lng和lat分别是经度(longitude)和纬度(latitude)的简写。

参考的函数原型：

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

EARTH_RADIUS=6371           # 地球平均半径，6371km

def hav(theta):
    s = sin(theta / 2)
    return s * s

def get_distance_hav(lat0,lng0,lat1,lng1):
    "用haversine公式计算球面两点间的距离。"
    # 经纬度转换成弧度
    lat0 = radians(lat0)
    lat1 = radians(lat1)
    lng0 = radians(lng0)
    lng1 = radians(lng1)

    dlng = fabs(lng0 - lng1)
    dlat = fabs(lat0 - lat1)
    h = hav(dlat) + cos(lat0) * cos(lat1) * hav(dlng)
    distance = 2 * EARTH_RADIUS * asin(sqrt(h))

    return distance