function [d, dist] = calculateDistanceOnSphere(P_xyz, indx, r)

if nargin < 3
    r = 1;
end

n1 = P_xyz(indx(:, 1), :);
n2 = P_xyz(indx(:, 2), :);

% Great-circle distance with Spherical law of cosines:
dist = r * atan2(vecnorm(cross(n1, n2, 2), 2, 2), dot(n1, n2, 2));
d    = sum(dist);