地球上兩點的經緯度計算他們距離的公式.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的長度。在根據(jù)余弦定律可以求夾角。AB在三角形AQB中，AQ的長度可以根據(jù)AB的緯度之差計算。BQ在三角形BPQ中，BP和PQ可求，角BPQ可以根據(jù)兩者的經度求出，這樣BQ的長度也可以求出來，所以AB的長度是可以求出來的。因為三角形ABQ是直角三角形，已經得到兩個邊知道了角AOB后，AB的弧長是可以求的。這樣推出其公式就不難了 關于用經緯度計算距離： 地球赤道上環(huán)繞地球一周走一圈共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為距離。 至于比例尺計算就不廢話了/這是主函數(shù)     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;        /這是計算兩點間航向的函數(shù)     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