/* atan - arctangent functions */

/*
functions used to compute the arctangent of two numbers (the sine term
and the cosine term.

This stolen from the c floating point library at yale cs dept vax.

Based on Hart & Cheney (#5077 , 19.56D)
*/

static	double	sq2p1	= 2.414213562373095048802e0;
static	double	sq2m1	= 0.414213562373095048802e0;
static	double	pio2	= 1.570796326794896619231e0;
static	double	pio4	= 0.785398163397448309615e0;
static	double	p4	= 0.161536412982230228262e2;
static	double	p3	= 0.268425481955039737941e3;
static	double	p2	= 0.115302935154048501154e4;
static	double	p1	= 0.178040631643319697105e4;
static	double	p0	= 0.896785974036638619599e3;
static	double	q4	= 0.589569705084446222279e2;
static	double	q3	= 0.536265374031215315104e3;
static	double	q2	= 0.166678381488163371845e4;
static	double	q1	= 0.207933497444540981287e4;
static	double	q0	= 0.896785974036638619624e3;

static	double	xatan(arg)

double	arg;
{
	double	argsq;
	double	value;

	argsq = arg*arg;
	value = ((((p4*argsq+p3)*argsq+p2)*argsq+p1)*argsq+p0);
	value = value/(((((argsq+q4)*argsq+q3)*argsq+q2)*argsq+q1)*argsq+q0);
	return(value*arg);
}

/* satan reduces range to single quadrant and calls xatan to do work */

static	double	satan(arg)

double	arg;
{
	if (arg < sq2m1)
	    return(xatan(arg));
	else
	    if (arg > sq2p1)
		return(pio2-xatan(1.0/arg));
	    else
		return(pio4+xatan((arg-1.0)/(arg+1.0)));
}
/* atan2 figures out the quadrant and returns correct atan */

double	atan2(arg2,arg1)

double	arg2,arg1;
{
	if ((arg1+arg2) == arg1)
	{
	    if (arg1 == 0.0) return(0.0);
	    if (arg1 > 0.0)
		return(pio2);
	    else
		return(-pio2);
	}
	else
	{
	    if (arg2 < 0.0)
	    {
		if (arg1 >= 0.0)
		    return(pio2+pio2-satan(-arg1/arg2));
		else
		    return(-pio2-pio2+satan(arg1/arg2));
	    }
	    else
	    {
		if (arg1 > 0.0)
		    return(satan(arg1/arg2));
		else
		    return(-satan(-arg1/arg2));
}	}   }