地球上两点的经纬度计算他们距离的公式.doc

1、假设地球是一个标准球体，半径为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)+sinb*siny 因此AB两点的球面距离为 R*arccoscosb*cosy*cos(a-x)+sinb*siny 注：1.x,y,a,b都是角度，最后结果中给出的

2、arccos因为弧度形式。 2.所谓的“东经为正，西经为负，北纬为正，南纬为负”是为了计算的方便。 比如某点为西京145°，南纬36°，那么计算时可用(-145°,-36°) 3.AB对球心所张角的球法实际上是求<OA>和<OB>两向量的夹角K。 用公式<OA>*<OB>=|OA|*|OB|*cosK 可以得到其中地球平均半径为6371.004 km  假设地球是个标准的球体：半径可以查出来，假设是R:如图：要算出A到B的球面距离，先要求出A跟B的夹角，即角AOB，求角AOB可以先求AOB的最大边

3、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/6

4、0=1.8553258km=1855.3m 而每一分又有60秒,每一秒就代表1855.3m/60=30.92m 任意两点距离计算公式为 d111.12cos1/sinAsinB十cosAcosBcos(BA) 其中A点经度，纬度分别为A和A，B点的经度、纬度分别为B和B，d为距离。 至于比例尺计算就不废话了/这是主函数     double   CChartCtrl:CalcltDstns(float   fStarPtx,   float   fStarPty,   float   fEndPtx,  

5、 float   fEndPty)         /已知起始点坐标(fStartPtx,   fStartPty)及到达点坐标(fEndPtx,fEndPty)     /计算航程dbDstns.     /起始点,到达点坐标:经纬度     /航程:海里(1852米)     /Created   by   zhl     /2002.7.3     /precision:0.0001

6、60; 海里     /check   param     double   dbDir=CalcltDirct(fStarPtx,fStarPty,fEndPtx,fEndPty);         double   delta_fy=fEndPtx-fStarPtx;     double   delta_lnmg=fEndPty-fStarPty;         int   mk=(int)fEnd

7、Ptx*(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

8、  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-

9、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:Cal

10、cltDirct(float   fStarPtx,   float   fStarPty,   float   fEndPtx,   float   fEndPty)         /已知起始点坐标(fStartPtx,   fStartPty)及到达点坐标(fEndPtx,fEndPty)     /计算航向fDirect.     /起始点,到达点坐标:经纬度     /航向:角度    

11、 /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|f

12、EndPtx<-90.0f     |fEndPty>180.0f|fEndPty<-180.0f)         AfxMessageBox(strErr);     return   -1;         double   delta_fy=fEndPtx-fStarPtx;     double   delta_lnmg=fEndPty-fStarPty;     /经度差应小于

13、180度     if(delta_lnmg   <   -180.0)     delta_lnmg   +=   360.0;     if(delta_lnmg   >   180.0)     delta_lnmg   -=   360.0;     /delta_lnmg   >   0.0   从   西   ->   东

14、     delta_lnmg   <   0.0   从   东   ->   西     BOOL   bGoEast=FALSE,bGoNorth=FALSE;     if(delta_lnmg   >=   0.0)     bGoEast=TRUE;     else     bGoEast=FALSE;     /delta_fy &#

15、160; >   0.0   从   南   ->   北     delta_fy   <   0.0   从   北   ->   南     if(delta_fy>=0.0)     bGoNorth=TRUE;     else     bGoNorth=FALSE;         if(delta_fy=

16、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);/纬度渐

17、长率         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)*1