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)).
Yes, you will get these kinds of errors with a global Mercator projection: it is accurate at the equator and the distortion increases exponentially with latitude away from the equator. The distance distortion is exactly 2 (100%) at 60 degrees latitude. At your test latitudes (64.14 degrees) I compute a distortion of 2.294, exactly agreeing with the ratio 904/394 = 2.294. (Earlier I computed 2.301 but that was based on a sphere, not the WGS84 ellipsoid. The difference (of 0.3%) gives us a sense of the accuracy you might gain from using an ellipsoid-based projection versus the sphere-based Haversine formula.)
There is no such thing as a global projection that yields highly accurate distances everywhere. That's one reason the UTM zone system is used!
One solution is to use spherical geometry for all your calculations, but you have rejected that (which is reasonable if you're going to be doing complex operations, but the decision might be worth revisiting).
Another solution is to adapt the projection to the points being compared. For example, you could safely use a transverse Mercator (as in the UTM system) with a meridian lying near the center of the region of interest. Moving the meridian is a simple thing to do: just subtract the meridian's longitude from all the longitudes and use a single TM projection centered at the Prime Meridian (with a scale factor of 1, rather than the 0.9996 of the UTM system). For your work this will tend to be more accurate than using UTM itself. It will give correct angles (TM is conformal) and will be remarkably accurate for points separated by only a few tens of kilometers: expect better than six-digit accuracy. In fact, I would be inclined to attribute any small differences between these adapted-TM distances and the Haversine distances to the difference between the ellipsoid (used for the TM projection) and the sphere (used by Haversine), rather than to distortion in the projection.
Best Answer
As mentioned by @mkkenedy in a comment on the question: once you've converted your latitude/longitude coordinates to Mercator coordinates, they are on a Euclidean plane, where you can use the Pythagorean theorem.
Specifically, if your Mercator coordinates are
(x1, y1)and(x2, y2), the distance is: