--- ray/src/common/fvect.c	1989/02/02 10:34:32	1.1
+++ ray/src/common/fvect.c	2003/02/25 02:47:21	2.7
@@ -1,20 +1,16 @@
-/* Copyright (c) 1986 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: fvect.c,v 2.7 2003/02/25 02:47:21 greg Exp $";
 #endif
-
 /*
- *  fvect.c - routines for float vector calculations
- *
- *     8/14/85
+ *  fvect.c - routines for floating-point vector calculations
  */
 
+#include "copyright.h"
+
+#include  <math.h>
 #include  "fvect.h"
 
-#define  FTINY		1e-7
 
-
 double
 fdot(v1, v2)			/* return the dot product of two vectors */
 register FVECT  v1, v2;
@@ -27,11 +23,12 @@ double
 dist2(p1, p2)			/* return square of distance between points */
 register FVECT  p1, p2;
 {
-	static FVECT  delta;
+	FVECT  delta;
 
 	delta[0] = p2[0] - p1[0];
 	delta[1] = p2[1] - p1[1];
 	delta[2] = p2[2] - p1[2];
+
 	return(DOT(delta, delta));
 }
 
@@ -41,13 +38,13 @@ dist2line(p, ep1, ep2)		/* return square of distance t
 FVECT  p;		/* the point */
 FVECT  ep1, ep2;	/* points on the line */
 {
-	static double  d, d1, d2;
+	register double  d, d1, d2;
 
 	d = dist2(ep1, ep2);
 	d1 = dist2(ep1, p);
-	d2 = dist2(ep2, p);
+	d2 = d + d1 - dist2(ep2, p);
 
-	return(d1 - (d+d1-d2)*(d+d1-d2)/d/4);
+	return(d1 - 0.25*d2*d2/d);
 }
 
 
@@ -56,7 +53,7 @@ dist2lseg(p, ep1, ep2)		/* return square of distance t
 FVECT  p;		/* the point */
 FVECT  ep1, ep2;	/* the end points */
 {
-	static double  d, d1, d2;
+	register double  d, d1, d2;
 
 	d = dist2(ep1, ep2);
 	d1 = dist2(ep1, p);
@@ -69,11 +66,13 @@ FVECT  ep1, ep2;	/* the end points */
 		if (d1 - d2 > d)
 			return(d2);
 	}
+	d2 = d + d1 - d2;
 
-	return(d1 - (d+d1-d2)*(d+d1-d2)/d/4);	/* distance to line */
+	return(d1 - 0.25*d2*d2/d);	/* distance to line */
 }
 
 
+void
 fcross(vres, v1, v2)		/* vres = v1 X v2 */
 register FVECT  vres, v1, v2;
 {
@@ -83,24 +82,59 @@ register FVECT  vres, v1, v2;
 }
 
 
+void
+fvsum(vres, v0, v1, f)		/* vres = v0 + f*v1 */
+register FVECT  vres, v0, v1;
+register double  f;
+{
+	vres[0] = v0[0] + f*v1[0];
+	vres[1] = v0[1] + f*v1[1];
+	vres[2] = v0[2] + f*v1[2];
+}
+
+
 double
 normalize(v)			/* normalize a vector, return old magnitude */
 register FVECT  v;
 {
-	static double  len;
+	register double  len, d;
 	
-	len = DOT(v, v);
+	d = DOT(v, v);
 	
-	if (len <= FTINY*FTINY)
+	if (d <= 0.0)
 		return(0.0);
 	
-	if (len >= (1.0-FTINY)*(1.0-FTINY) &&
-			len <= (1.0+FTINY)*(1.0+FTINY))
-		return(1.0);
+	if (d <= 1.0+FTINY && d >= 1.0-FTINY)
+		len = 0.5 + 0.5*d;	/* first order approximation */
+	else
+		len = sqrt(d);
 
-	len = sqrt(len);
-	v[0] /= len;
-	v[1] /= len;
-	v[2] /= len;
+	v[0] *= d = 1.0/len;
+	v[1] *= d;
+	v[2] *= d;
+
 	return(len);
+}
+
+
+void
+spinvector(vres, vorig, vnorm, theta)	/* rotate vector around normal */
+FVECT  vres, vorig, vnorm;
+double  theta;
+{
+	double  sint, cost, normprod;
+	FVECT  vperp;
+	register int  i;
+	
+	if (theta == 0.0) {
+		if (vres != vorig)
+			VCOPY(vres, vorig);
+		return;
+	}
+	cost = cos(theta);
+	sint = sin(theta);
+	normprod = DOT(vorig, vnorm)*(1.-cost);
+	fcross(vperp, vnorm, vorig);
+	for (i = 0; i < 3; i++)
+		vres[i] = vorig[i]*cost + vnorm[i]*normprod + vperp[i]*sint;
 }