great circle - Distance between lat/long points

I am attempting to calculate the distance between two latitude/longitude points. I have a piece of code that mostly works that I yanked from this post but I do not really understand how it works.

Here is the code:

// POINT 1
$thisLat = deg2rad(44.638);
$thisLong = deg2rad(-63.587);

// POINT 2
$otherLat = deg2rad(44.644);
$otherLong = deg2rad(-63.911);

$MeanRadius = 6378 - 21 * sin($lat1);


$xa = (Cos($thisLat)) * (Cos($thisLong));
$ya = (Cos($thisLat)) * (Sin($thisLong));
$za = (Sin($thisLat));

$xb = (Cos($otherLat)) * (Cos($otherLong));
$yb = (Cos($otherLat)) * (Sin($otherLong));
$zb = (Sin($otherLat));

$distance = $MeanRadius * Acos($xa * $xb + $ya * $yb + $za * $zb);


echo $distance;
?>

I have a couple questions:

1. what are xa, ya, za? I understand that they are points on a 3D cartesian plane but where are they relative to? The center of the earth?


2. How does this cos($xa * $xb + $ya * $yb + $za * $zb) calculate the distance between the points? I know that in 2D I would do this:

Pythagorean Theorem 

distance^2 = b^2 + a^2
distance = sqr((y2-y1)^2 + (x2 - x1)^2)

3. How accurate will this be? There was some discussion about that on the other page. But I specifically want to use the distance to tell if users are within the something like 10m, 20m or 50m of each other. Will I be able to do this with good accuracy?


4. What should I use for $MeanRadius? Is that a reasonable value? I think that that value assumes that the earth is a elipse.

Answer

This is terrible code for general-purpose use because it can give erroneous results or even fail altogether for short distances. Use the Haversine Formula instead.

(The formula on which your code is based converts two points on the sphere (not an ellipse) into their 3D Cartesian coordinates (xa,ya,za) and (xb,yb,zb) on the unit sphere and forms their dot product, which will equal the cosine of the angle between them. The ACos function returns that angle, which when scaled by the earth's radius will estimate the distance. The problem is that the cosine of a small angle, say of size 'e' in radians, differs from 1 by an amount close to e^2/2. This disappears into the floating point error cloud when e is smaller than the square root of twice the floating point precision. If you're computing in single precision this means values of e less than 0.001 -- about one kilometer -- will be confused with zero! In double precision the cutoff is around e = 10^-8, but by the time e = 10^-4 or so (about 10 meters) you potentially can lose so much precision that you need to worry, unless the implementation of the ACos function is unusually accurate (e.g., has some high-precision internal computations built in)).