谁知道根据地球上两点的经纬度计算他们距离的公式或者算法？

Question

谁知道根据地球上两点的经纬度计算他们距离的公式或者算法？

Answer 1

//这是主函数
double CChartCtrl::CalcltDstns(float fStarPtx, float fStarPty, float fEndPtx, float fEndPty)
{
//已知起始点坐标(fStartPtx, fStartPty)及到达点坐标(fEndPtx,fEndPty)
//计算航程dbDstns.
//起始点,到达点坐标:经纬度
//航程:海里(1852米)
//Created by zhl
//2002.7.3
//precision:0.0001 海里
//check param
double dbDir=CalcltDirct(fStarPtx,fStarPty,fEndPtx,fEndPty); double delta_fy=fEndPtx-fStarPtx;
double delta_lnmg=fEndPty-fStarPty; int mk=(int)fEndPtx*(int)fStarPtx;
double fy_m,dbDstns;//
if(mk>=0)
{//不跨赤道航行
fy_m=(fStarPtx+fEndPtx)/2;
}else
{//跨赤道航行
fy_m=fabs(fStarPtx)>fabs(fEndPtx)?fStarPtx/2:fEndPtx/2;
}
double delta_l=(1852.2-9.3*cos(fy_m*M_PI/180+fy_m*M_PI/180))
*delta_fy*60/1852; if((dbDir>80&&dbDir<100)||(dbDir>260&&dbDir<280))
{//东西向
dbDstns=1.00181*delta_lnmg*60*cos(fy_m*M_PI/180)*
sqrt(1-e2*sin(fy_m*M_PI/180)*sin(fy_m*M_PI/180))/
sin(dbDir*M_PI/180);
}else
{//南北向
dbDstns=delta_l/cos(dbDir*M_PI/180);
}
return dbDstns;
}

Answer 2

//这是计算两点间航向的函数
double CChartCtrl::CalcltDirct(float fStarPtx, float fStarPty, float fEndPtx, float fEndPty)
{
//已知起始点坐标(fStartPtx, fStartPty)及到达点坐标(fEndPtx,fEndPty)
//计算航向fDirect.
//起始点,到达点坐标:经纬度
//航向:角度
//Created by zhl
//2002.7.2
//check param
CString strErr;
strErr.LoadString(IDS_CHK_15002);
if(fStarPtx>90.0f||fStarPtx<-90.0f||fStarPty>180.0f
||fStarPty<-180.0f||fEndPtx>90.0f||fEndPtx<-90.0f
||fEndPty>180.0f||fEndPty<-180.0f)
{
AfxMessageBox(strErr);
return -1;
}
double delta_fy=fEndPtx-fStarPtx;
double delta_lnmg=fEndPty-fStarPty;
//经度差应小于180度
if(delta_lnmg < -180.0)
delta_lnmg += 360.0;
if(delta_lnmg > 180.0)
delta_lnmg -= 360.0;
//delta_lnmg > 0.0 从 西 ---> 东  delta_lnmg < 0.0 从 东 ---> 西
BOOL bGoEast=FALSE,bGoNorth=FALSE;
if(delta_lnmg >= 0.0)
bGoEast=TRUE;
else
bGoEast=FALSE;
//delta_fy > 0.0 从 南 ---> 北  delta_fy < 0.0 从 北 ---> 南
if(delta_fy>=0.0)
bGoNorth=TRUE;
else
bGoNorth=FALSE; if(delta_fy==0)
{
if(delta_lnmg==0)return 0;
return bGoEast?90:270;
}
double d1=7915.7045*(e/2*log10((1-e*sin(fStarPtx*M_PI/180))
/(1+e*sin(fStarPtx*M_PI/180)))
+log10(tan((45+fStarPtx/2)*M_PI/180.0)));//纬度渐长率 double d2=7915.7045*(e/2*log10((1-e*sin(fEndPtx*M_PI/180))/
(1+e*sin(fEndPtx*M_PI/180)))
+log10(tan((45+fEndPtx/2)*M_PI/180.0)));//纬度渐长率 double delta_d=d2-d1;////纬度渐长率差(分) double dbDir=atan(delta_lnmg*60/delta_d)*180/M_PI;
if(!bGoEast&&bGoNorth)dbDir=360+dbDir;
if(!bGoEast&&!bGoNorth)dbDir=180+dbDir;
if(bGoEast&&!bGoNorth)dbDir=180+dbDir;
return dbDir;
}