function d = sphdist (b1, a1, b2, a2) ## b1, b2: latitudes; a1, a2: longitudes; all in radians. d = acos (sin (b1) .* sin (b2) + cos (b1) .* cos (b2) .* cos (a1 .- a2)); endfunction