Calculation of Distance from GPS Coordinates in DIAdem

Issue Details

I have an application, in which my hardware measures GPS-coordinates of two points in less than 10 centimeter distance, namely their corresponding latitude and longitude values. My problem is, that when I am using the Spherical Law of Cosines to calculate the actual distance, I get always exactly zero as a result, although the points are several centimeters apart.

The formulas look like the following in DIAdem:

Zeta=acos(sin(Latitude1)*sin(Latitude2)+cos(Latitude1)*cos(Latitude2)*cos(Longitude1-Longitude2))
CalculatedDistance=R*Zeta

Here R=6371 km, the radius of the Earth.

I use the Double datatype to store the raw latitude and longitude values, and also to store the resulting distance value.

 

Solution

The problem comes from the finite resolution of the Double datatype. According to the IEEE 754 standard the so called significand precision is 53 bits, so the maximum relative rounding error is 2^-53. It means, that considering decimal  numbers, there will be maximum 53*log(2)~16 significant digits, i.e. nonzero digits stored by the fraction part of the variable.
Let us consider the Spherical Law of Cosines in DIAdem now.
If the physical distance is less than 6.371 cm-s, Zeta will be in reality less than 10^-8 radians,
so in cos(Zeta)~1-Zeta^2/2 only the 17th significant digit is different from 9, it will be 5. After rounding we get for the argument of the acos-function exactly 1, so Zeta will be calculated to be exactly zero radian, and the calculated distance will be 0 cm.

To avoid this rounding error for distances under 1 km use the  so called Equirectangular approximation:

Zeta=sqrt((Latitude1-Latitude2)^2+(Longitude1-Longitude2)^2*(cos((Latitude1+Latitude2)/2))^2)
CalculatedDistance=R*Zeta

If the actual distance is under 1 km, this approximation causes negligible relative errors in the calculated distance, well under 1 ppm.

To get geometrically exact results use the Haversine formula.