/************************************************************************
 *									*
 *   Replacement atan2() function for PML library. Original function    *
 *   from Fred Fish does not conform to Standard C.  It has arguments   *
 *   in a reverse order.  It also does not complain if both arguments   *
 *   are 0.0                                                            *
 *									*
 *   Michal Jaegermann - December 1990
 *									*
 ************************************************************************
 */
/************************************************************************
 *									*
 *				N O T I C E				*
 *									*
 *			Copyright Abandoned, 1987, Fred Fish		*
 *									*
 *	This previously copyrighted work has been placed into the	*
 *	public domain by the author (Fred Fish) and may be freely used	*
 *	for any purpose, private or commercial.  I would appreciate	*
 *	it, as a courtesy, if this notice is left in all copies and	*
 *	derivative works.  Thank you, and enjoy...			*
 *									*
 *	The author makes no warranty of any kind with respect to this	*
 *	product and explicitly disclaims any implied warranties of	*
 *	merchantability or fitness for any particular purpose.		*
 *									*
 ************************************************************************
 */


/*
 *  FUNCTION
 *
 *	atan2   double precision arc tangent of two arguments
 *
 *  KEY WORDS
 *
 *	atan2
 *	machine independent routines
 *	trigonometric functions
 *	math libraries
 *
 *  DESCRIPTION
 *
 *	Returns double precision arc tangent of two
 *	double precision floating point arguments ( atan(Y/X) ).
 *
 *  USAGE
 *
 *	double atan2(y, x)
 *	double y;
 *	double x;
 *
 *  (This order of arguments really makes sense if you will notice
 *   that this function is mainly used during a conversion from
 *   rectangular to polar coordinates)
 *
 *  REFERENCES
 *
 *	Fortran 77 user's guide, Digital Equipment Corp. pp B-4.
 *
 *  RESTRICTIONS
 *
 *	For precision information refer to documentation of the
 *	other floating point library routines called.
 *	
 *  PROGRAMMER
 *
 *	Fred Fish
 *	Modified by Michal Jaegermann
 *
 *  INTERNALS
 *
 *	Computes atan2(y, x) from:
 *
 *		1.	If x = 0 and y != 0 then
 *			atan2(y, x) = PI/2 * (sign(y))
 *                      otherwise error signaled and
 *                      atan2(y ,x) is 0.
 *
 *		2.	If x > 0 then
 *			atan2(y, x) = atan(y/x)
 *
 *		3.	If x < 0 then atan2(y, x) =
 *			PI*(sign(y)) + atan(y/x)
 *
 */

#include <stdio.h>
#include <pmluser.h>
#include "pml.h"

static char funcname[] = "atan2";

double atan2 (y, x)
double y;
double x;
{
    double result;
    extern double sign();
    extern double atan();
    struct exception xcpt;

    ENTER ("atan2");
    DEBUG4 ("atan2in", "x = %le y = %le", x, y);
    if (x == 0.0) {
	if (y == 0.0) {
	    result = 0.0;
	    xcpt.type = DOMAIN;
	    xcpt.name = funcname;
	    xcpt.arg1 = x;
	    xcpt.retval = result;
	    if (!matherr(&xcpt)) {
		fprintf (stderr, "%s: DOMAIN error\n", funcname);
		errno = EDOM;
	    }
	} else {
	    result = sign (HALFPI, y);
	}
    } else {
	result = atan (y/x);
	if (x < 0.0) {
	    result += sign (PI, y);
	}
    }
    DEBUG3 ("atan2out", "result %le", result);
    LEAVE ();
    return (result);
}
