1 原理公式

在地球上，计算两点之间的直线距离通常使用地理坐标系（例如WGS84）。计算两地直线距离的公式是根据经纬度之间的大圆距离（Great Circle Distance）来计算的。该公式基于球面三角学，常用的公式是 $H a v ers in e$ 公式。

形式一：
$d_{lon} = lon_2 - lon_1$
$d_{lat} = lat_2 - lat_1$
$sin²(d_{lat}/2) + cos(lat_1) * cos(lat_2) * sin²(d_{lon}/2)$
$atan^2(\sqrt{a}, \sqrt{(1-a)})$
$d = R * c$

形式二：
$d = R*acos(sin(lon_1)*sin(lon_2) + cos(lon_1)*cos(lon_2)*cos(lat_2-lat_1))$

其中：

$R$ 是地球的半径，约为6371千米。

在这里插入图片描述

请注意，这个公式假设地球是一个完美的球形。实际上，地球的形状更像一个椭球，因此使用更精确的地理信息系统（GIS）软件或库（如proj.4或GeographicLib）可能会得到更准确的结果。

此外，经纬度通常以度为单位，但上述公式中的角度应被视为弧度。如果经纬度是以度数形式给出的，需要将其转换为弧度。可以通过将度数乘以π/180并取整数部分得到弧度值。例如， $30°$ 转换为弧度为 $π /180 * 30$ 。

2 代码实现

2.1 JavaScript

function calculateDistance(lat1, lon1, lat2, lon2) {const R = 6371; // 地球半径，单位为千米const rad = (angle) => angle * Math.PI / 180; // 将角度转换为弧度const lat1Rad = rad(lat1);const lon1Rad = rad(lon1);const lat2Rad = rad(lat2);const lon2Rad = rad(lon2);const dLat = lat2Rad - lat1Rad;const dLon = lon2Rad - lon1Rad;const a = Math.sin(dLat / 2) * Math.sin(dLat / 2) +Math.cos(lat1Rad) * Math.cos(lat2Rad) *Math.sin(dLon / 2) * Math.sin(dLon / 2);const c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));const distance = R * c; // 返回单位为千米的距离return distance;
}// 示例用法
const lat1 = 39.9087; // 北京的经纬度
const lon1 = 116.4074;
const lat2 = 31.2304; // 上海的经纬度
const lon2 = 121.4737;const distance = calculateDistance(lat1, lon1, lat2, lon2);
console.log(distance); // 输出直线距离（单位：千米）

2.2 C++

#include <cmath>
#include <iostream>// 计算两个经纬度之间的距离（单位：千米）
double calculateDistance(double lat1, double lon1, double lat2, double lon2) {const double R = 6371; // 地球半径，单位为千米// 将经纬度转换为弧度double lat1Rad = std::atan(std::tan(lat1 * (M_PI / 180)) * std::cos(lon1 * (M_PI / 180)));double lon1Rad = lon1 * (M_PI / 180);double lat2Rad = std::atan(std::tan(lat2 * (M_PI / 180)) * std::cos(lon2 * (M_PI / 180)));double lon2Rad = lon2 * (M_PI / 180);// 计算两个经纬度之间的弧度差double dLat = lat2Rad - lat1Rad;double dLon = lon2Rad - lon1Rad;// 根据球面三角法公式计算距离double a = std::sin(dLat / 2) * std::sin(dLat / 2) +std::cos(lat1Rad) * std::cos(lat2Rad) *std::sin(dLon / 2) * std::sin(dLon / 2);double c = 2 * std::atan2(std::sqrt(a), std::sqrt(1 - a));// 返回距离return R * c;
}// 示例用法
int main() {double lat1 = 39.9087; // 北京的经纬度double lon1 = 116.4074;double lat2 = 31.2304; // 上海的经纬度double lon2 = 121.4737;double distance = calculateDistance(lat1, lon1, lat2, lon2);std::cout << "距离：" << distance << " 千米" << std::endl;return 0;
}

请注意，此代码使用了C++标准库中的数学函数和常量。此外，将经纬度转换为弧度的方法需要使用std::atan和std::tan函数。

2.3 Python

import mathdef calculate_distance(lat1, lon1, lat2, lon2):R = 6371  # 地球半径，单位为千米rad = lambda angle: angle * math.pi / 180  # 将角度转换为弧度lat1_rad = rad(lat1)lon1_rad = rad(lon1)lat2_rad = rad(lat2)lon2_rad = rad(lon2)dLat = lat2_rad - lat1_raddLon = lon2_rad - lon1_rada = math.sin(dLat / 2) ** 2 + math.cos(lat1_rad) * math.cos(lat2_rad) * math.sin(dLon / 2) ** 2c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))distance = R * c  # 返回单位为千米的距离return distance# 示例用法
lat1 = 39.9087  # 北京的经纬度
lon1 = 116.4074
lat2 = 31.2304  # 上海的经纬度
lon2 = 121.4737distance = calculate_distance(lat1, lon1, lat2, lon2)
print(distance)  # 输出直线距离（单位：千米）

请注意，由于Python和JavaScript之间的一些语法差异，需要使用math模块来进行数学计算。另外，由于Python中没有lambda函数的功能，因此需要使用普通的函数定义来代替。

2.4 MATLAB

function distance = calculateDistance(lat1, lon1, lat2, lon2)const R = 6371; % 地球半径，单位为千米lat1Rad = rad(lat1);lon1Rad = rad(lon1);lat2Rad = rad(lat2);lon2Rad = rad(lon2);dLat = lat2Rad - lat1Rad;dLon = lon2Rad - lon1Rad;a = sin(dLat / 2) .^ 2 + cos(lat1Rad) .* cos(lat2Rad) .* sin(dLon / 2) .^ 2;c = 2 * atan2(sqrt(a), sqrt(1 - a));distance = R * c; % 返回单位为千米的距离
end% 示例用法
lat1 = 39.9087; % 北京的经纬度
lon1 = 116.4074;
lat2 = 31.2304; % 上海的经纬度
lon2 = 121.4737;distance = calculateDistance(lat1, lon1, lat2, lon2);
disp(distance); % 输出直线距离（单位：千米）