在提供的代码中有一些错误，我在下面修复了它。

以上所有答案都假定地球是一个球体。然而，更精确的近似是扁球体。

a= 6378.137#equitorial radius in km
b= 6356.752#polar radius in km

def Distance(lat1, lons1, lat2, lons2):
    lat1=math.radians(lat1)
    lons1=math.radians(lons1)
    R1=(((((a**2)*math.cos(lat1))**2)+(((b**2)*math.sin(lat1))**2))/((a*math.cos(lat1))**2+(b*math.sin(lat1))**2))**0.5 #radius of earth at lat1
    x1=R1*math.cos(lat1)*math.cos(lons1)
    y1=R1*math.cos(lat1)*math.sin(lons1)
    z1=R1*math.sin(lat1)

    lat2=math.radians(lat2)
    lons2=math.radians(lons2)
    R2=(((((a**2)*math.cos(lat2))**2)+(((b**2)*math.sin(lat2))**2))/((a*math.cos(lat2))**2+(b*math.sin(lat2))**2))**0.5 #radius of earth at lat2
    x2=R2*math.cos(lat2)*math.cos(lons2)
    y2=R2*math.cos(lat2)*math.sin(lons2)
    z2=R2*math.sin(lat2)
    
    return ((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)**0.5

精确计算中长点之间距离所需的函数是复杂的，陷阱也很多。我不推荐哈弗辛或其他球形的解决方案，因为有很大的不准确性(地球不是一个完美的球体)。vincenty公式更好，但在某些情况下会抛出错误，即使编码正确。

与其自己编写函数，我建议使用geopy，它已经实现了非常精确的地理库来进行距离计算(论文来自作者)。

#pip install geopy
from geopy.distance import geodesic
NY = [40.71278,-74.00594]
Beijing = [39.90421,116.40739]
print("WGS84: ",geodesic(NY, Beijing).km) #WGS84 is Standard
print("Intl24: ",geodesic(NY, Beijing, ellipsoid='Intl 1924').km) #geopy includes different ellipsoids
print("Custom ellipsoid: ",geodesic(NY, Beijing, ellipsoid=(6377., 6356., 1 / 297.)).km) #custom ellipsoid

#supported ellipsoids:
#model             major (km)   minor (km)     flattening
#'WGS-84':        (6378.137,    6356.7523142,  1 / 298.257223563)
#'GRS-80':        (6378.137,    6356.7523141,  1 / 298.257222101)
#'Airy (1830)':   (6377.563396, 6356.256909,   1 / 299.3249646)
#'Intl 1924':     (6378.388,    6356.911946,   1 / 297.0)
#'Clarke (1880)': (6378.249145, 6356.51486955, 1 / 293.465)
#'GRS-67':        (6378.1600,   6356.774719,   1 / 298.25)

这个库的唯一缺点是它不支持向量化计算。 对于向量化计算，您可以使用新的gevectorslib。

#pip install geovectorslib
from geovectorslib import inverse
print(inverse(lats1,lons1,lats2,lons2)['s12'])

lat和lon是numpy数组。Geovectorslib是非常准确和非常快!我还没有找到改变椭球的方法。标准采用WGS84椭球，是大多数用途的最佳选择。

你可以使用CLLocationDistance中的构建来计算这个:

CLLocation *location1 = [[CLLocation alloc] initWithLatitude:latitude1 longitude:longitude1];
CLLocation *location2 = [[CLLocation alloc] initWithLatitude:latitude2 longitude:longitude2];
[self distanceInMetersFromLocation:location1 toLocation:location2]

- (int)distanceInMetersFromLocation:(CLLocation*)location1 toLocation:(CLLocation*)location2 {
    CLLocationDistance distanceInMeters = [location1 distanceFromLocation:location2];
    return distanceInMeters;
}

在你的例子中，如果你想要公里，只要除以1000。

在提供的代码中有一些错误，我在下面修复了它。

以上所有答案都假定地球是一个球体。然而，更精确的近似是扁球体。

a= 6378.137#equitorial radius in km
b= 6356.752#polar radius in km

def Distance(lat1, lons1, lat2, lons2):
    lat1=math.radians(lat1)
    lons1=math.radians(lons1)
    R1=(((((a**2)*math.cos(lat1))**2)+(((b**2)*math.sin(lat1))**2))/((a*math.cos(lat1))**2+(b*math.sin(lat1))**2))**0.5 #radius of earth at lat1
    x1=R1*math.cos(lat1)*math.cos(lons1)
    y1=R1*math.cos(lat1)*math.sin(lons1)
    z1=R1*math.sin(lat1)

    lat2=math.radians(lat2)
    lons2=math.radians(lons2)
    R2=(((((a**2)*math.cos(lat2))**2)+(((b**2)*math.sin(lat2))**2))/((a*math.cos(lat2))**2+(b*math.sin(lat2))**2))**0.5 #radius of earth at lat2
    x2=R2*math.cos(lat2)*math.cos(lons2)
    y2=R2*math.cos(lat2)*math.sin(lons2)
    z2=R2*math.sin(lat2)
    
    return ((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)**0.5

我在这里发布了我的工作示例。

在MySQL中列出表中指定点(我们使用一个随机点- lat:45.20327, long:23.7806)之间距离小于50 KM的所有点(表中字段为coord_lat和coord_long):

列出所有距离<50，单位:公里(地球半径6371公里):

SELECT denumire, (6371 * acos( cos( radians(45.20327) ) * cos( radians( coord_lat ) ) * cos( radians( 23.7806 ) - radians(coord_long) ) + sin( radians(45.20327) ) * sin( radians(coord_lat) ) )) AS distanta 
FROM obiective 
WHERE coord_lat<>'' 
    AND coord_long<>'' 
HAVING distanta<50 
ORDER BY distanta desc

上面的例子是在MySQL 5.0.95和5.5.16 (Linux)中测试的。

下面是另一个转换为Ruby代码的代码:

include Math
#Note: from/to = [lat, long]

def get_distance_in_km(from, to)
  radians = lambda { |deg| deg * Math.PI / 180 }
  radius = 6371 # Radius of the earth in kilometer
  dLat = radians[to[0]-from[0]]
  dLon = radians[to[1]-from[1]]

  cosines_product = Math.sin(dLat/2) * Math.sin(dLat/2) + Math.cos(radians[from[0]]) * Math.cos(radians[to[1]]) * Math.sin(dLon/2) * Math.sin(dLon/2)

  c = 2 * Math.atan2(Math.sqrt(cosines_product), Math.sqrt(1-cosines_product)) 
  return radius * c # Distance in kilometer
end