%% Haversine 公式
function distance = haversine(lat1, lon1, lat2, lon2)
    % 将角度转换为弧度
    R = 6371; % 地球半径，单位为千米
    lat1 = deg2rad(lat1);
    lon1 = deg2rad(lon1);
    lat2 = deg2rad(lat2);
    lon2 = deg2rad(lon2);

    % Haversine公式
    dlat = lat2 - lat1;
    dlon = lon2 - lon1;
    a = sin(dlat/2)^2 + cos(lat1) * cos(lat2) * sin(dlon/2)^2;
    c = 2 * atan2(sqrt(a), sqrt(1-a));

    % 距离
    distance = R * c;
end