/*
   FUNCTION:  double atan2 (y, x)

   PURPOSE:
	To compute the double precision arctangent of y/x in radians.

   INPUT PARAMETERS:
	double y - The numerator (rise) of y/x.
	double x - The denomenator (run) of y/x.

   OUTPUT PARAMETERS:
	Returned - Arctan(y/x) in radians.

   DESCRIPTION:
	This routine computes the arctangent via polynomial evaluation. See
	the Arctangent algorithm 5076 from Hart's "Computer Approximations".

   WRITTEN BY:  Scott D. Roth		DATE:  March 27, 1981
*/


#define WHOLE_PI	3.1415926535897932d
#define HALF_PI		1.5707963267948966d
#define FORTH_PI	0.78539816339744831d
#define TAN_8TH_PI	0.41421356d
#define TAN_3_8THS_PI	2.4142136d

static double q[4] = {		/* Coefficients of polynomial q */
	.44541340059290680E2d,
	.92324801072300975E2d,
	.62835930511032377E2d,
	.15503977551421988E2d};
static double p[5] = {		/* Coefficients of polynomial p */
	.44541340059290680E2d,
	.77477687719204209E2d,
	.40969264832102256E2d,
	.66605790170092627E1d,
	.15897402884823070E0d};

double atan2 (y, x)
    double y, x;
	{
	double yx, adjust, atan_yx, s, p_val, q_val;
	register int y_pos, x_pos;			/* Booleans */

	y_pos = (y >= 0);		/* Remember signs for later */
	if (!(x_pos = (x > 0)))
	    if (x == 0)			/* Special 0 denominator case */
		return(y_pos ? HALF_PI : -HALF_PI);

	yx = y / x;
	if (y_pos ^ x_pos)		/* Make yx positive */
	    yx = -yx;
	adjust = 0;
	if (yx >= TAN_8TH_PI)		/* Adjust value to be within */
	    if (yx >= TAN_3_8THS_PI)	/* tan(PI/8) and tan(3*PI/8) */
		{
		yx = -1.0d / yx;
		adjust = HALF_PI;
		}
	    else
		{
		yx = (yx - 1.0d) / (yx + 1.0d);
		adjust = FORTH_PI;
		}		
	s = yx * yx;
					/* Evaluate P polynomial */
	p_val = ((((p[4] * s + p[3]) * s + p[2]) * s + p[1]) * s + p[0]) * yx;
					/* Evaluate Q polynomial */
	q_val = (((s + q[3]) * s + q[2]) * s + q[1]) * s + q[0];

					/* Atan = P(yx**2) / Q(yx**2) */
	atan_yx = p_val / q_val + adjust;

	if (!x_pos)
	    atan_yx = WHOLE_PI - atan_yx;
	return (y_pos ? atan_yx : -atan_yx);
	}