根据2个经纬度点，计算这2个经纬度点之间的距离

原文 http://www.cnblogs.com/softfair/p/distance_of_two_latitude_and_longitude_points.html

根据 2 个经纬度点，计算这 2 个经纬度点之间的距离

球面上任意两点之间的距离计算公式可以参考维基百科上的下述文章。

值得一提的是，维基百科推荐使用Haversine公式，理由是Great-circle distance公式用到了大量余弦函数，而两点间距离很短时（比如地球表面上相距几百米的两点），余弦函数会得出0.999...的结果，会导致较大的舍入误差。而Haversine公式采用了正弦函数，即使距离很小，也能保持足够的有效数字。以前采用三角函数表计算时的确会有这个问题，但经过实际验证，采用计算机来计算时，两个公式的区别不大。稳妥起见，这里还是采用Haversine公式。

其中

R为地球半径，可取平均值 6371km；
φ ₁ , φ ₂ 表示两点的纬度；
Δλ 表示两点经度的差值。

根据 2 个经纬度坐标，距离计算函数

下面就是计算球面间两点(lat1, lon1) - (lat2, lon2)之间距离的函数。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace HarvenSin
{
    class Program
    {
  /// <summary>
  /// 根据经纬度，计算2个点之间的距离。
  /// </summary>
  /// <param name="args"></param>
  static void Main(string[] args)
  {
      //39.94607,116.32793  31.24063,121.42575
      Console.WriteLine(Distance(39.94607, 116.32793, 31.24063, 121.42575));
  }
  public static double HaverSin(double theta)
  {
      var v = Math.Sin(theta / 2);
      return v * v;
  }
  static double EARTH_RADIUS = 6371.0;//km 地球半径 平均值，千米
  /// <summary>
  /// 给定的经度1，纬度1；经度2，纬度2. 计算2个经纬度之间的距离。
  /// </summary>
  /// <param name="lat1">经度1</param>
  /// <param name="lon1">纬度1</param>
  /// <param name="lat2">经度2</param>
  /// <param name="lon2">纬度2</param>
  /// <returns>距离（公里、千米）</returns>
  public static double Distance(double lat1,double lon1, double lat2,double lon2)
  {
      //用haversine公式计算球面两点间的距离。
      //经纬度转换成弧度
      lat1 = ConvertDegreesToRadians(lat1);
      lon1 = ConvertDegreesToRadians(lon1);
      lat2 = ConvertDegreesToRadians(lat2);
      lon2 = ConvertDegreesToRadians(lon2);
      //差值
      var vLon = Math.Abs(lon1 - lon2);
      var vLat = Math.Abs(lat1 - lat2);
      //h is the great circle distance in radians, great circle就是一个球体上的切面，它的圆心即是球心的一个周长最大的圆。
      var h = HaverSin(vLat) + Math.Cos(lat1) * Math.Cos(lat2) * HaverSin(vLon);
      var distance = 2 * EARTH_RADIUS * Math.Asin(Math.Sqrt(h));
      return distance;
  }
  /// <summary>
  /// 将角度换算为弧度。
  /// </summary>
  /// <param name="degrees">角度</param>
  /// <returns>弧度</returns>
  public static double ConvertDegreesToRadians(double degrees)
  {
      return degrees * Math.PI / 180;
  }
  public static double ConvertRadiansToDegrees(double radian)
  {
      return radian * 180.0 / Math.PI;
  }
    }
}

公式来历：

VERSINE(F)=1-cos(F)

Haversine名字来历是Ha-VERSINE，即Half-Versine ，表示sin的一半的意思。

hav(A) = (1-cos(A))/2 = sin(A/2)* sin(A/2)

推倒过程：

如下一个半径为1 的圆，O是圆心，A、B是弦(chord)。角度AOB=theta。则角度AOC=theta/2。OC是垂直于AB的垂线(perpendicular)。AC长度是sin(theta/2)，AB长度是2*sin(theta/2)。

(图1)

如下地球图所示，假设半径R为1，O是球心，A (lat1,lon1) 和 B (lat2,lon2) 是我们感兴趣的2个点。2跟经度线 lon1，lon2相交于北极(north pole)N。EF所在的线是赤道(equator)。ACBD是平面上的等腰梯形的四个顶点(vertice)。AC和DB的弦（直线）在图上没有画出。CD的位置是：C (lat2,lon1) and D (lat1,lon2)。角度AOC是A点与C点的纬度差 dlat。角度EOF是经度E点和经度F点的差dlon。