两经纬度之间的距离计算

Posted 2020-06-14

以下是计算两经纬度之间距离的代码，分为：头文件、源代码和测试代码三部分。
具体如下：

 1 // LatLonDistanceDlg.h : 头文件  
 2 //  
 3   
 4 #pragma once  
 5   
 6 typedef struct  
 7 {  
 8     double dLongitude;  
 9     double dLatitude;  
10 }MyLatLong_T,*pMyLatLong_T;  
11   
12 // CLatLonDistanceDlg 对话框  
13 class CLatLonDistanceDlg : public CDialog  
14 {  
15     ......  
16   
17     void DistanceOfTwoPositions(MyLatLong_T ll_PosStart, MyLatLong_T ll_PosEnd, double &dis);  
18 // 实现  
19 protected:  
20     ......  
21 };  
22   
23 // LatLonDistanceDlg.cpp : 实现文件  
24 //  
25   
26 // 常量定义  
27 const double gConstPI = 3.1415926536;  
28 // WGS-84 长轴半径  
29 #define EARTH_LONG_RADIUS_WGS84 6378137.0       // 单位: 米  
30 // WGS-84 地球扁率(earth flattening)  
31 #define EARTH_FLATTENING_WGS84  298.257223563  
32 // meters per sea mile  
33 #define METERS_PER_SEA_MILE     1852.0  
34   
35 // 输入参数: MyLatLong_T 是个结构，包括经度和纬度  
36 // 输出参数: dDistance 就是返回的大圆距离  
37 void CLatLonDistanceDlg::DistanceOfTwoPositions(MyLatLong_T StartLatLong, MyLatLong_T EndLatLong, double &dDistance)  
38 {  
39     double DInRadians;  
40     double dTmpVal = 0.0;  
41     double dFi = 0.0;  
42     double dFi2 = 0.0;  
43     double dDrda = EndLatLong.dLongitude - StartLatLong.dLongitude;  
44     double dTmp = 0.0;  
45     double dTmp2 = 0.0;  
46   
47     dDrda = dDrda * gConstPI / 180.0;   // in radians  
48     dFi = StartLatLong.dLatitude;  
49     dFi2 = EndLatLong.dLatitude;  
50   
51     dFi = dFi * gConstPI / 180.0;       // in radians  
52     dFi2 = dFi2 * gConstPI / 180.0;     // in radians  
53   
54     dTmpVal = sin(dFi) * sin(dFi2) + cos(dFi) * cos(dFi2) * cos(dDrda);  
55     if(fabs(dTmpVal) > 1.0)  
56     {  
57         AfxMessageBox(L"Invalidate value of arccos!");      // Use Unicode Character Set  
58         return ;  
59     }  
60     DInRadians = acos(dTmpVal);         // in radians  
61     dTmp = (sin(dFi) + sin(dFi2));  
62     dTmp2 = (sin(dFi) - sin(dFi2));  
63     dTmpVal = ((3 * sin(DInRadians) - DInRadians) * dTmp * dTmp ) / (1 + cos(DInRadians));  
64     dTmpVal = dTmpVal - ((3 * sin(DInRadians) + DInRadians) * dTmp2 * dTmp2) / (1 - cos(DInRadians));  
65   
66     dDistance = EARTH_LONG_RADIUS_WGS84 * DInRadians + (EARTH_LONG_RADIUS_WGS84 / (4 * EARTH_FLATTENING_WGS84)) * dTmpVal;//in meters  
67 }  
68   
69 // 测试代码  
70 // TODO: 在此添加额外的初始化代码  
71 {  
72     // 深圳百合酒店 纬度:22.601369,经度114.115145  
73     // 深圳百合星城 纬度:22.601334,经度114.115807  
74     // 深圳布吉农批 纬度:22.583596,经度114.112227  
75     MyLatLong_T LatLon;  
76     MyLatLong_T LatLon2;  
77     MyLatLong_T LatLon3;  
78     double dDis = 0.0;  
79   
80     LatLon.dLatitude = 22.601369;  
81     LatLon.dLongitude = 114.115145;  
82     LatLon2.dLatitude = 22.601334;  
83     LatLon2.dLongitude = 114.115807;  
84     LatLon3.dLatitude = 22.583596;  
85     LatLon3.dLongitude = 114.112227;  
86   
87     DistanceOfTwoPositions(LatLon,LatLon2,dDis);  
88     TRACE("%f\r\n",dDis);   // 68.177865  
89     DistanceOfTwoPositions(LatLon,LatLon3,dDis);  
90     TRACE("%f\r\n",dDis);   // 1990.891295  
91 }