[Neftali [Germán.Estévez]]

En su momento utilicé esta con bastates buenos resultados.
No es de mi autoría, pero no recuerdo de dónde salió para poner el link. La he tenido que buscar en el código del proyecto.

Código Delphi [-]
// distance in meters onto the WGS84 geoid between two Coordinates
// some data are exageretely precise (way beyond GPS's capacity) so I
// added the AlmostNul switch, as arccos(1) returns an error when two points
// are too close (and therefore distance is 0)
function WGS84Distance(plon1, plat1, plon2, plat2: Double;
  Almostnul: boolean = False): Double;
var
  d, s, c: Double;
  lat1: Double;
  lon1: Double;
  lat2: Double;
  lon2: Double;
  almost1: Double;
Const
  R = 6378137; // Medium earth radius in meter
begin
  // d = R x arcos [ sin(lat1) x sin(lat2) + cos(lat1) x cos(lat2) x cos(lon2-lon1) ]
  lat1 := DegToRad(plat1);
  lon1 := DegToRad(plon1);
  lat2 := DegToRad(plat2);
  lon2 := DegToRad(plon2);
  if (lon1 = lon2) AND (lat1 = lat2) then
    d := 0
  else begin
    s := sin(lat1) * sin(lat2);
    c := cos(lat1) * cos(lat2) * cos(lon2 - lon1);

    if Almostnul then
      almost1 := round((s + c) * 100000000) / 100000000
    else
      almost1 := (s + c);

    if almost1 <> 1 then
      d := R * arccos(s + c)
    else
      d := 0;
  end;
  Result := d;
end;