Problem Statement

The haversine formula is an equation important in navigation, giving great-circle distances between two points on a sphere from their longitudes and latitudes. It is a special case of a more general formula in spherical trigonometry, the law of haversines, relating the sides and angles of spherical "triangles".

Implement a great-circle distance function, or use a library function, to show the great-circle distance between:

Most of the examples below adopted Kaimbridge's recommended value of 6372.8 km for the earth radius. However, the derivation of this ellipsoidal quadratic mean radius is wrong (the averaging over azimuth is biased). When applying these examples in real applications, it is better to use the mean earth radius, 6371 km. This value is recommended by the International Union of Geodesy and Geophysics and it minimizes the RMS relative error between the great circle and geodesic distance.

Let's start with the solution:

#include
"share/atspre_staload.hats"

staload "libc/SATS/math.sats"
staload _ = "libc/DATS/math.dats"
staload "libc/SATS/stdio.sats"
staload "libc/SATS/stdlib.sats"

#define R 6372.8
#define TO_RAD (3.1415926536 / 180)

typedef d = double

fun
dist
(
  th1: d, ph1: d, th2: d, ph2: d
) : d = let
  val ph1 = ph1 - ph2
  val ph1 = TO_RAD * ph1
  val th1 = TO_RAD * th1
  val th2 = TO_RAD * th2
  val dz = sin(th1) - sin(th2)
  val dx = cos(ph1) * cos(th1) - cos(th2)
  val dy = sin(ph1) * cos(th1)
in
  asin(sqrt(dx*dx + dy*dy + dz*dz)/2)*2*R
end // end of [dist]

implement
main0((*void*)) = let
  val d = dist(36.12, ~86.67, 33.94, ~118.4);
  /* Americans don't know kilometers */
in
  $extfcall(void, "printf", "dist: %.1f km (%.1f mi.)\n", d, d / 1.609344)
end // end of [main0]