The following code snippet is a modified version of the one posted by Rob Fowler.
/************************************************************************
** This function is due to Rob Fowler. Given dy and dx between 2 points
** A and B, we calculate a number in [0.0, 8.0) which is a monotonic
** function of the direction from A to B.
**
** (0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0) correspond to
** ( 0, 45, 90, 135, 180, 225, 270, 315, 360) degrees, measured
** counter-clockwise from the positive x axis.
*************************************************************************/
double FowlerAngle (double dy, double dx)
{
double adx, ady; /* Absolute Values of Dx and Dy */
int code; /* Angular Region Classification Code */
adx = (dx < 0) ? -dx : dx; /* Compute the absolute values. */
ady = (dy < 0) ? -dy : dy;
code = (adx < ady) ? 1 : 0;
if (dx < 0) code += 2;
if (dy < 0) code += 4;
switch (code)
{
case 0: return (dx==0) ? 0 : ady/adx; /* [ 0, 45] */
case 1: return (2.0 - (adx/ady)); /* ( 45, 90] */
case 3: return (2.0 + (adx/ady)); /* ( 90,135) */
case 2: return (4.0 - (ady/adx)); /* [135,180] */
case 6: return (4.0 + (ady/adx)); /* (180,225] */
case 7: return (6.0 - (adx/ady)); /* (225,270) */
case 5: return (6.0 + (adx/ady)); /* [270,315) */
case 4: return (8.0 - (ady/adx)); /* [315,360) */
}
}