|
发表于 2005-9-15 20:04:07
|
显示全部楼层
引用 远风 在 2005-9-15 14:21 时的帖子:
我用这个公司写了程序算地球表面两点间距离,但是发现答案不正确啊,望各位高手指点。
请参考:球面两点间的球面距离的计算(2) http://www.zahui.com/html/9/41685.htm
里面有计算公式的推导和编程实现。
转载:
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|2 = |AC|2 + |CE|2 – 2 * |AC| *|CE| * cos(θ)
→ |AE|2 = R2 * cos2 (α) + R2 * cos2 (β) - 2 * R2 * cos(α)cos(β) * cos(θ)
|BE| = |CD|
由|AB|2 = |AE|2 + |BE|2
→ |AB|2 = R2 * cos2 (α) + R2 * cos2 (β) - 2 * R2 * cos(α)cos(β) * cos(θ) + R2 * sin2 (α) + R2 * sin2 (β) + 2 * R2 * sin(α) sin(β)
→ |AB|2 = 2R2 (1 + sin(α) sin(β) - cos2 (α) cos2 (β) cos(θ))
在三角形AOB中,亦可用余弦定理得到
→ |OA|2 + |OB|2– 2 |OA| * |OB| * cos(AOB)
→ 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; |
|