假设地球是一个标准球体，半径为R,并且假设东经为正，西经为负，北纬为正，南纬为负，
则A(x,y)的坐标可表示为（R*cosy*cosx, R*cosy*sinx,R*siny）
B(a,b)可表示为(R*cosb*cosa ,R*cosb*sina,R*sinb)
于是，AB对于球心所张的角的余弦大小为
cosb*cosy*(cosa*cosx+sina*sinx)+sinb*siny=cosb*cosy*cos(a-x)+s inb*siny
因此AB两点的球面距离为
R*{arccos[cosb*cosy*cos(a-x)+sinb*siny]}
注：1.x,y,a,b都是角度，最后结果中给出的arccos因为弧度形式。
2.所谓的“东经为正，西经为负，北纬为正，南纬为负”是为了计算的方便。
比如某点为西京145°，南纬36°，那么计算时可用(-145°,-36°)
3.AB对球心所张角的球法实际上是求
用公式
可以得到
其中地球平均半径为6371.004 km
假设地球是个标准的球体：半径可以查出来，假设是R:
如图：
要算出A到B的球面距离，先要求出A跟B的夹角，即角AOB，
求角AOB可以先求AOB的最大边AB的长度。在根据余弦定律可以求夹角。
AB在三角形AQB中，AQ的长度可以根据AB的纬度之差计算。
BQ在三角形BPQ中，BP和PQ可求，角BPQ可以根据两者的经度求出，这样BQ的长度也可以求出来，
所以AB的长度是可以求出来的。因为三角形ABQ是直角三角形，已经得到两个边
知道了角AOB后，AB的弧长是可以求的。
这样推出其公式就不难了
关于用经纬度计算距离：
地球赤道上环绕地球一周走一圈共40075.04公里,而@一圈分成360°,而每1°(度)有60,每一度一秒在赤道上的长度计算如下：
40075.04km/360°=111.31955km
111.31955km/60=1.8553258km=1855.3m
而每一分又有60秒,每一秒就代表1855.3m/60=30.92m
任意两点距离计算公式为
d＝111.12cos{1/[sinΦAsinΦB十cosΦAcosΦBcos(λB—λA)]}
其中A点经度，纬度分别为λA和ΦA，B点的经度、纬度分别为λB和ΦB，d为距离。至于比例尺计算就不废话了
//这是主函数
double CChartCtrl::CalcltDstns(float fStarPtx, float fStarPty, float fEndPt x, 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;
}
//这是计算两点间航向的函数
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;
}
M_PI就是pi的值