下面是c#语言(用纬度和弧度表示):

double CalculateGreatCircleDistance(double lat1, double long1, double lat2, double long2, double radius)
{
    return radius * Math.Acos(
        Math.Sin(lat1) * Math.Sin(lat2)
        + Math.Cos(lat1) * Math.Cos(lat2) * Math.Cos(long2 - long1));
}

如果你的纬度和长度是用角度表示的，那么除以180/PI就可以转换成弧度。

如果你使用的是。net，不要重新启动轮子。看到System.Device.Location。在另一个答案的评论中赞扬fnx。

using System.Device.Location;

double lat1 = 45.421527862548828D;
double long1 = -75.697189331054688D;
double lat2 = 53.64135D;
double long2 = -113.59273D;

GeoCoordinate geo1 = new GeoCoordinate(lat1, long1);
GeoCoordinate geo2 = new GeoCoordinate(lat2, long2);

double distance = geo1.GetDistanceTo(geo2);

下面是答案中的Swift实现

func degreesToRadians(degrees: Double) -> Double {
    return degrees * Double.pi / 180
}

func distanceInKmBetweenEarthCoordinates(lat1: Double, lon1: Double, lat2: Double, lon2: Double) -> Double {

    let earthRadiusKm: Double = 6371

    let dLat = degreesToRadians(degrees: lat2 - lat1)
    let dLon = degreesToRadians(degrees: lon2 - lon1)

    let lat1 = degreesToRadians(degrees: lat1)
    let lat2 = degreesToRadians(degrees: lat2)

    let a = sin(dLat/2) * sin(dLat/2) +
    sin(dLon/2) * sin(dLon/2) * cos(lat1) * cos(lat2)
    let c = 2 * atan2(sqrt(a), sqrt(1 - a))
    return earthRadiusKm * c
}

我认为R中的一个算法版本仍然缺失:

gpsdistance<-function(lat1,lon1,lat2,lon2){

# internal function to change deg to rad

degreesToRadians<- function (degrees) {
return (degrees * pi / 180)
}

R<-6371e3  #radius of Earth in meters

phi1<-degreesToRadians(lat1) # latitude 1
phi2<-degreesToRadians(lat2) # latitude 2
lambda1<-degreesToRadians(lon1) # longitude 1
lambda2<-degreesToRadians(lon2) # longitude 2

delta_phi<-phi1-phi2 # latitude-distance
delta_lambda<-lambda1-lambda2 # longitude-distance

a<-sin(delta_phi/2)*sin(delta_phi/2)+
cos(phi1)*cos(phi2)*sin(delta_lambda/2)*
sin(delta_lambda/2)

cc<-2*atan2(sqrt(a),sqrt(1-a))

distance<- R * cc

return(distance)  # in meters
}

在Python中，你可以使用geopy库使用WGS84椭球来计算测地线距离:

from geopy.distance import geodesic
newport_ri = (41.49008, -71.312796)
cleveland_oh = (41.499498, -81.695391)
print(geodesic(newport_ri, cleveland_oh).km)

这是我在Elixir中的实现

defmodule Geo do
  @earth_radius_km 6371
  @earth_radius_sm 3958.748
  @earth_radius_nm 3440.065
  @feet_per_sm 5280

  @d2r :math.pi / 180

  def deg_to_rad(deg), do: deg * @d2r

  def great_circle_distance(p1, p2, :km), do: haversine(p1, p2) * @earth_radius_km
  def great_circle_distance(p1, p2, :sm), do: haversine(p1, p2) * @earth_radius_sm
  def great_circle_distance(p1, p2, :nm), do: haversine(p1, p2) * @earth_radius_nm
  def great_circle_distance(p1, p2, :m), do: great_circle_distance(p1, p2, :km) * 1000
  def great_circle_distance(p1, p2, :ft), do: great_circle_distance(p1, p2, :sm) * @feet_per_sm

  @doc """
  Calculate the [Haversine](https://en.wikipedia.org/wiki/Haversine_formula)
  distance between two coordinates. Result is in radians. This result can be
  multiplied by the sphere's radius in any unit to get the distance in that unit.
  For example, multiple the result of this function by the Earth's radius in
  kilometres and you get the distance between the two given points in kilometres.
  """
  def haversine({lat1, lon1}, {lat2, lon2}) do
    dlat = deg_to_rad(lat2 - lat1)
    dlon = deg_to_rad(lon2 - lon1)

    radlat1 = deg_to_rad(lat1)
    radlat2 = deg_to_rad(lat2)

    a = :math.pow(:math.sin(dlat / 2), 2) +
        :math.pow(:math.sin(dlon / 2), 2) *
        :math.cos(radlat1) * :math.cos(radlat2)

    2 * :math.atan2(:math.sqrt(a), :math.sqrt(1 - a))
  end
end