两经纬度之间的距离计算

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

 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 }