你可以在f#的fssnip中找到这个实现(有一些很好的解释)

以下是重要的部分:

let GreatCircleDistance<[&ltMeasure>] 'u> (R : float<'u>) (p1 : Location) (p2 : Location) =
    let degToRad (x : float&ltdeg>) = System.Math.PI * x / 180.0&ltdeg/rad>

    let sq x = x * x
    // take the sin of the half and square the result
    let sinSqHf (a : float&ltrad>) = (System.Math.Sin >> sq) (a / 2.0&ltrad>)
    let cos (a : float&ltdeg>) = System.Math.Cos (degToRad a / 1.0&ltrad>)

    let dLat = (p2.Latitude - p1.Latitude) |> degToRad
    let dLon = (p2.Longitude - p1.Longitude) |> degToRad

    let a = sinSqHf dLat + cos p1.Latitude * cos p2.Latitude * sinSqHf dLon
    let c = 2.0 * System.Math.Atan2(System.Math.Sqrt(a), System.Math.Sqrt(1.0-a))

    R * c

对于java

public static double degreesToRadians(double degrees) {
    return degrees * Math.PI / 180;
}

public static double distanceInKmBetweenEarthCoordinates(Location location1, Location location2) {
    double earthRadiusKm = 6371;

    double dLat = degreesToRadians(location2.getLatitude()-location1.getLatitude());
    double dLon = degreesToRadians(location2.getLongitude()-location1.getLongitude());

    double lat1 = degreesToRadians(location1.getLatitude());
    double lat2 = degreesToRadians(location2.getLatitude());

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

如果你需要更准确的数据，可以看看这个。

Vincenty's formulae are two related iterative methods used in geodesy to calculate the distance between two points on the surface of a spheroid, developed by Thaddeus Vincenty (1975a) They are based on the assumption that the figure of the Earth is an oblate spheroid, and hence are more accurate than methods such as great-circle distance which assume a spherical Earth. The first (direct) method computes the location of a point which is a given distance and azimuth (direction) from another point. The second (inverse) method computes the geographical distance and azimuth between two given points. They have been widely used in geodesy because they are accurate to within 0.5 mm (0.020″) on the Earth ellipsoid.

打印稿版本

export const degreeToRadian = (degree: number) => {
  return degree * Math.PI / 180;
}

export const distanceBetweenEarthCoordinatesInKm = (lat1: number, lon1: number, lat2: number, lon2: number) => {
    const earthRadiusInKm = 6371;

    const dLat = degreeToRadian(lat2 - lat1);
    const dLon = degreeToRadian(lon2 - lon1);

    lat1 = degreeToRadian(lat1);
    lat2 = degreeToRadian(lat2);

    const a = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.sin(dLon / 2) * Math.sin(dLon / 2) * Math.cos(lat1) * Math.cos(lat2);
    const c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));

    return earthRadiusInKm * c;
}

    private double deg2rad(double deg)
    {
        return (deg * Math.PI / 180.0);
    }

    private double rad2deg(double rad)
    {
        return (rad / Math.PI * 180.0);
    }

    private double GetDistance(double lat1, double lon1, double lat2, double lon2)
    {
        //code for Distance in Kilo Meter
        double theta = lon1 - lon2;
        double dist = Math.Sin(deg2rad(lat1)) * Math.Sin(deg2rad(lat2)) + Math.Cos(deg2rad(lat1)) * Math.Cos(deg2rad(lat2)) * Math.Cos(deg2rad(theta));
        dist = Math.Abs(Math.Round(rad2deg(Math.Acos(dist)) * 60 * 1.1515 * 1.609344 * 1000, 0));
        return (dist);
    }

    private double GetDirection(double lat1, double lon1, double lat2, double lon2)
    {
        //code for Direction in Degrees
        double dlat = deg2rad(lat1) - deg2rad(lat2);
        double dlon = deg2rad(lon1) - deg2rad(lon2);
        double y = Math.Sin(dlon) * Math.Cos(lat2);
        double x = Math.Cos(deg2rad(lat1)) * Math.Sin(deg2rad(lat2)) - Math.Sin(deg2rad(lat1)) * Math.Cos(deg2rad(lat2)) * Math.Cos(dlon);
        double direct = Math.Round(rad2deg(Math.Atan2(y, x)), 0);
        if (direct < 0)
            direct = direct + 360;
        return (direct);
    }

    private double GetSpeed(double lat1, double lon1, double lat2, double lon2, DateTime CurTime, DateTime PrevTime)
    {
        //code for speed in Kilo Meter/Hour
        TimeSpan TimeDifference = CurTime.Subtract(PrevTime);
        double TimeDifferenceInSeconds = Math.Round(TimeDifference.TotalSeconds, 0);
        double theta = lon1 - lon2;
        double dist = Math.Sin(deg2rad(lat1)) * Math.Sin(deg2rad(lat2)) + Math.Cos(deg2rad(lat1)) * Math.Cos(deg2rad(lat2)) * Math.Cos(deg2rad(theta));
        dist = rad2deg(Math.Acos(dist)) * 60 * 1.1515 * 1.609344;
        double Speed = Math.Abs(Math.Round((dist / Math.Abs(TimeDifferenceInSeconds)) * 60 * 60, 0));
        return (Speed);
    }

    private double GetDuration(DateTime CurTime, DateTime PrevTime)
    {
        //code for speed in Kilo Meter/Hour
        TimeSpan TimeDifference = CurTime.Subtract(PrevTime);
        double TimeDifferenceInSeconds = Math.Abs(Math.Round(TimeDifference.TotalSeconds, 0));
        return (TimeDifferenceInSeconds);
    }

这是我在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