--- ray/src/common/tcos.c	2003/02/22 02:07:22	3.2
+++ ray/src/common/tcos.c	2013/02/07 20:19:52	3.7
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: tcos.c,v 3.2 2003/02/22 02:07:22 greg Exp $";
+static const char	RCSid[] = "$Id: tcos.c,v 3.7 2013/02/07 20:19:52 greg Exp $";
 #endif
 /*
  * Table-based cosine approximation.
@@ -9,82 +9,22 @@ static const char	RCSid[] = "$Id: tcos.c,v 3.2 2003/02
  *
  * No interpolation in this version.
  *
- * External symbols declared in standard.h
+ * External symbols declared in rtmath.h
  */
 
-/* ====================================================================
- * The Radiance Software License, Version 1.0
- *
- * Copyright (c) 1990 - 2002 The Regents of the University of California,
- * through Lawrence Berkeley National Laboratory.   All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- *
- * 1. Redistributions of source code must retain the above copyright
- *         notice, this list of conditions and the following disclaimer.
- *
- * 2. Redistributions in binary form must reproduce the above copyright
- *       notice, this list of conditions and the following disclaimer in
- *       the documentation and/or other materials provided with the
- *       distribution.
- *
- * 3. The end-user documentation included with the redistribution,
- *           if any, must include the following acknowledgment:
- *             "This product includes Radiance software
- *                 (http://radsite.lbl.gov/)
- *                 developed by the Lawrence Berkeley National Laboratory
- *               (http://www.lbl.gov/)."
- *       Alternately, this acknowledgment may appear in the software itself,
- *       if and wherever such third-party acknowledgments normally appear.
- *
- * 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
- *       and "The Regents of the University of California" must
- *       not be used to endorse or promote products derived from this
- *       software without prior written permission. For written
- *       permission, please contact radiance@radsite.lbl.gov.
- *
- * 5. Products derived from this software may not be called "Radiance",
- *       nor may "Radiance" appear in their name, without prior written
- *       permission of Lawrence Berkeley National Laboratory.
- *
- * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
- * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
- * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
- * DISCLAIMED.   IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
- * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
- * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
- * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
- * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
- * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- * ====================================================================
- *
- * This software consists of voluntary contributions made by many
- * individuals on behalf of Lawrence Berkeley National Laboratory.   For more
- * information on Lawrence Berkeley National Laboratory, please see
- * <http://www.lbl.gov/>.
- */
+#include "copyright.h"
 
 #include <math.h>
 
+#include "rtmath.h"
+
 #ifndef NCOSENTRY
-#define NCOSENTRY	256
+#define NCOSENTRY	1024
 #endif
 
-#ifdef M_PI
-#define PI		((double)M_PI)
-#else
-#define PI		3.14159265358979323846
-#endif
 
-
 double
-tcos(x)				/* approximate cosine */
-register double	x;
+tcos(double x)				/* approximate cosine */
 {
 	static double	costab[NCOSENTRY+1];
 	register int	i;
@@ -96,8 +36,8 @@ register double	x;
 	if (x < 0.)
 		x = -x;
 	i = (NCOSENTRY*2./PI) * x  +  0.5;
-	if (i >= 4*NCOSENTRY)
-		i %= 4*NCOSENTRY;
+	while (i >= 4*NCOSENTRY)
+		i -= 4*NCOSENTRY;
 	switch (i / NCOSENTRY) {
 	case 0:
 		return(costab[i]);
@@ -109,4 +49,24 @@ register double	x;
 		return(costab[(4*NCOSENTRY)-i]);
 	}
 	return(0.);		/* should never be reached */
+}
+
+
+/* Fast arctangent approximation due to Rajan et al. 2006 */
+double
+aatan2(double y, double x)
+{
+	double	ratio, aratio, val;
+
+	if (x == 0)
+		return (y > 0) ? PI/2. : 3./2.*PI;
+
+	aratio = (ratio = y/x) >= 0 ? ratio : -ratio;
+
+	if (aratio > 1.01)
+		return PI/2. - aatan2(x, y);
+
+	val = PI/4.*ratio - ratio*(aratio - 1.)*(0.2447 + 0.0663*aratio);
+
+	return val + PI*(x < 0);
 }