球面两点间的球面距离的计算(2)

球面两点间的球面距离的计算(2)

 

1．计算

上篇中提到了两点间球面距离的计算的理论基础，在这篇中，进行计算。

 

同样，假设点A的经度为东经a1度，纬度为北纬b1度；点B的经度为东经a2度，纬度为南纬b2度。

所以α =  b1, β = b2， 设∠ACE的值为θ，那么θ = |a1 – a2|。

可以计算出

|OC| = R * sin(α)      |OD| = R * sin(β)

|AC| = R * cos(α)      |BD| = R * cos(β)

   

    由于|CE| = |BD|，|CD| = |OD| + |OC|

    由余弦定理可以知道

    |AE|² = |AC|² + |CE|² – 2 * |AC| *|CE| * cos(θ)

ð        |AE|² = R² * cos² (α) + R² * cos² (β) - 2 * R² * cos(α)cos(β) * cos(θ)

 

|BE| = |CD|

由|AB|² = |AE|² + |BE|²

ð        |AB|² = R² * cos² (α) + R² * cos² (β) - 2 * R² * cos(α)cos(β) * cos(θ) + R² * sin² (α) + R² * sin² (β) + 2 * R² * sin(α) sin(β)

ð        |AB|² = 2R² (1 + sin(α) sin(β) - cos² (α) cos² (β) cos(θ))

 

在三角形AOB中，亦可用余弦定理得到

ð        |OA|² + |OB|²– 2 |OA| * |OB| * cos(AOB) = |OB|²

ð        cos(AOB) = cos(α) * cos(β) * cos(θ) - sin(α) sin(β)

 

这样可以求出AOB的值(单位为弧度)，最后球面距离为R * (AOB)。

完毕，呵呵。

 

注意到，

（1）当A、B同为东经或同为西经时，θ = |a1 – a2|，那么当 A、B不是同为东经或同为西经时，θ = |a1 + a2|，若θ > 180度，θ = 360 -θ（假设单位是角度，而不是弧度）。

（2）当A、B同为南半球或同为北半球时，|CD| = | |OC| - |OD| |，当A、B在不同的南北半球时，才是|CD| = |OC| + |OD|。

 

 

2．编程实现

（1）函数接口定义

int CalGlobeDistance(const ST_GlobePnt a_stPntA,

const ST_GlobePnt a_stPntB,

const double a_dRadius,

double &a_dDistance)

 

返回值:

    0：执行OK，返回值正确；

1：a_stPntA 或a_stPntB 的东经、西经取值错误；

2：a_stPntA 或a_stPntB的南纬、北纬取值错误；

3：a_stPntA 或a_stPntB的经度范围越界；

4：a_stPntA 或a_stPntB的纬度范围越界；

10:半径取值错误(不能小于0)。

 

参数：

a_stPntA

    给定的第一个点的坐标，关于结构体ST_GlobePnt的请参照以下的ST_GlobePnt定义；

 

a_stPntB

    给定的第二个点的坐标，关于结构体ST_GlobePnt的请参照以下的ST_GlobePnt定义；

 

a_dRadius

给定的球体的半径，取值不能小于0；

 

a_dDistance

最后求得的球面距离，如果在计算过程中出现错误，则其值为负数。

 

结构体ST_GlobePnt

 

typedef struct

{

    byte btEorW;      // 取值为1,东经; 0,西经;

    byte btNorS;      // 取值为1,北纬; 0,南纬;   

    double dLong;     // 经度，取值范围[0.0，180.0)

    double dLat;      // 纬度，取值范围[0.0，90.0]

}ST_GlobePnt;

 

 

（2）程序代码

 

常量定义

typedef struct

{

char btEorW;        // 取值为1,东经; 0,西经;

char btNorS;        // 取值为1,北纬; 0,南纬;   

double dLong;       // 经度，取值范围[0.0，180.0)

double dLat;        // 纬度，取值范围[0.0，90.0]

}ST_GlobePnt;

 

const char EAST = 1;

const char WEST = 0;

const char NORTH = 1;

const char SOUTH = 0;

 

const int EXC_OK = 0;        // 0：执行OK，返回值正确；

const int EorW_ERROR = 1;    // 1：东经、西经取值错误；

const int NorS_ERROR = 2;    // 2：南纬、北纬取值错误；

const int LONG_ERROR = 3;    // 3：经度范围越界；

const int LAT_ERROR = 4;     // 4: 纬度范围越界；

const int RADIUS_ERROR = 10;  // 10:半径取值错误(不能小于0)。

 

const double ACCURACY = 0.00001;    // 精确度

 

const double PI = 3.14159265359;

 

函数实现

辅助函数int CheckGPnt(const ST_GlobePnt &a_stPnt)

    功能判断点是否符合要求，实现如下：

 

int CheckGPnt(const ST_GlobePnt &a_stPnt)

{

    int iRet = EXC_OK;

   

    // check if east or west

    if(a_stPnt.btEorW != EAST && a_stPnt.btEorW != WEST)

    {

        iRet = EorW_ERROR;

        return iRet;

    }

   

    // check if nort or south

    if(a_stPnt.btNorS != NORTH && a_stPnt.btNorS != SOUTH)

    {

        iRet = EorW_ERROR;

        return iRet;

    }

   

    // check the range of Long.

    if(a_stPnt.dLong - 180.0 > ACCURACY || a_stPnt.dLong <= - ACCURACY)

    {

        iRet = LONG_ERROR;

        return iRet;       

    }

   

    // check the range of Lat

    if(a_stPnt.dLat - 90.0 > ACCURACY || a_stPnt.dLat <= -ACCURACY)

    {

        iRet = LAT_ERROR;

        return iRet;

    }   

   

    return iRet;

   

}

 

 

接口函数如下：

int CalGlobeDistance(const ST_GlobePnt &a_stPntA,

                    const ST_GlobePnt &a_stPntB,

                    const double &a_dRadius,

                    double &a_dDistance)

{

    int iRet = EXC_OK;   

    a_dDistance = -1.0;

 

     double dTheta;    // the angle of the two Longs

     double dAlpha;    // the angle of a_stPntA's Lat

     double dBeta;     // the angle of a_stPntB's Lat  

     double dResult;   // the angle of the two lines

                       //     which one is through the given point and the center

     double dCosResult; // = cos(dResult);

 

   // Check if the given points is OK --Begin-------   

    iRet = CheckGPnt(a_stPntA);

    if(iRet != EXC_OK)

    {

        return iRet;

    }

   

    iRet = CheckGPnt(a_stPntB);

    if(iRet != EXC_OK)

    {

        return iRet;

    }

   // Check if the given points is OK --End--------- 

 

   // check the range of Radius

    if(a_dRadius < -ACCURACY)

    {

        iRet = RADIUS_ERROR;

        return iRet;

    }

   

    dAlpha = a_stPntA.dLat;

    dBeta = a_stPntB.dLat;

   

   

    // Calculate the the angle of the two Longs --- Begin --------

    if(a_stPntA.btEorW == a_stPntB.btEorW)

    {

        dTheta = fabs(a_stPntA.dLong - a_stPntB.dLong);   

    }else // !=

    {

        dTheta = a_stPntA.dLong + a_stPntB.dLong;

    }

       

    if(dTheta - 180.00 > -ACCURACY)

    {

        dTheta = 360.0 - dTheta;

    }

    // Calculate the the angle of the two Longs --- End -------- 

    // Change DEGREE to RADIAN     --- Begin ------

    dAlpha = dAlpha / 180.0 * PI;

dBeta = dBeta / 180.0 * PI;

dTheta = dTheta / 180.0 * PI;

    // Change DEGREE to RADIAN     --- End ------ 

 

    // Calculate the angle of the two lines  --- Begin ---------

    if(a_stPntA.btNorS != a_stPntB.btNorS)

    {

        dCosResult = cos(dAlpha) * cos(dBeta) * cos(dTheta) - sin(dAlpha) * sin(dBeta);

    }else // !=

    {

        dCosResult = cos(dAlpha) * cos(dBeta) * cos(dTheta) + sin(dAlpha) * sin(dBeta);

    }

   

    dResult = acos(dCosResult);//  dResult ranges [0,Pi],needn't change

    // Calculate the angle of the two lines  --- End ---------

   

    // Get the distance around the globe

    a_dDistance = a_dRadius * dResult;           

   

    return iRet;

}

 

 

（3）补充说明

利用程序计算，可能存在一定的误差，减少误差的方法可能用到计算方法中提到的算法，呼呼，先告一段落，暂不考虑这个吧。

3．关于本文

    本文的版本为 Version 1.2 ，于 2007 年 5 月 8 日，在 1.1 的基础上稍作修改。 Version 1.1 ，完稿于 2005 年 1 月 12 日。

    你可以到我的个人主页http://zijinshi.cn(网站服务器在国外，可能需要代理服务器才能访问)中找到本文原文和相关源代码，当然，你也可以从我的blog：http://blog.csdn.net/zijinshi 上找到这篇文章。

    我的Email:frank@zijinshi.cn, 欢迎不吝赐教。☺

(完)