It does work that way. This function places the two points in a polar coordinates system (radius and angle). The resulting system is in fact centered in the first point, and the function computes the angle coordinate of the second point (with the Z coordinates ignored upon the translation).
In 3D, the problem is obviously expanded to three dimensions - more specifically a spherical coordinate system (radius, elevation angle, and azimuth angle). Because of the two distinct angles, separate functions would be needed to compute their values, with the radius remaining the classical euclidean distance between the two (the Azimuth Angle is the same as the 2D angle).
The
Wiki article on the subject is very neat, and offers much more information, including conversion formulas (just in case you need something like this for your 3D solution).