--- ray/src/common/fvect.c	2012/09/06 00:07:43	2.15
+++ ray/src/common/fvect.c	2014/12/08 23:51:12	2.21
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: fvect.c,v 2.15 2012/09/06 00:07:43 greg Exp $";
+static const char	RCSid[] = "$Id: fvect.c,v 2.21 2014/12/08 23:51:12 greg Exp $";
 #endif
 /*
  *  fvect.c - routines for floating-point vector calculations
@@ -7,11 +7,32 @@ static const char	RCSid[] = "$Id: fvect.c,v 2.15 2012/
 
 #include "copyright.h"
 
+#define _USE_MATH_DEFINES
 #include  <math.h>
 #include  "fvect.h"
+#include  "random.h"
 
+double
+Acos(double x)			/* insurance for touchy math library */
+{
+	if (x <= -1.+FTINY*FTINY)
+		return(M_PI);
+	if (x >= 1.-FTINY*FTINY)
+		return(.0);
+	return(acos(x));
+}
 
 double
+Asin(double x)			/* insurance for touchy math library */
+{
+	if (x <= -1.+FTINY*FTINY)
+		return(-M_PI/2.);
+	if (x >= 1.-FTINY*FTINY)
+		return(M_PI/2);
+	return(asin(x));
+}
+
+double
 fdot(				/* return the dot product of two vectors */
 const FVECT v1,
 const FVECT v2
@@ -29,9 +50,7 @@ const FVECT p2
 {
 	FVECT  delta;
 
-	delta[0] = p2[0] - p1[0];
-	delta[1] = p2[1] - p1[1];
-	delta[2] = p2[2] - p1[2];
+	VSUB(delta, p2, p1);
 
 	return(DOT(delta, delta));
 }
@@ -87,9 +106,13 @@ const FVECT v1,
 const FVECT v2
 )
 {
-	vres[0] = v1[1]*v2[2] - v1[2]*v2[1];
-	vres[1] = v1[2]*v2[0] - v1[0]*v2[2];
-	vres[2] = v1[0]*v2[1] - v1[1]*v2[0];
+	if ((vres == v1) | (vres == v2)) {
+		FVECT	vtmp;
+		VCROSS(vtmp, v1, v2);
+		VCOPY(vres, vtmp);
+		return;
+	}
+	VCROSS(vres, v1, v2);
 }
 
 
@@ -101,9 +124,7 @@ const FVECT v1,
 double f
 )
 {
-	vres[0] = v0[0] + f*v1[0];
-	vres[1] = v0[1] + f*v1[1];
-	vres[2] = v0[2] + f*v1[2];
+	VSUM(vres, v0, v1, f);
 }
 
 
@@ -135,6 +156,29 @@ FVECT  v
 
 
 int
+getperpendicular(		/* choose random perpedicular direction */
+FVECT vp,				/* returns normalized */
+const FVECT v				/* input vector must be normalized */
+)
+{
+	FVECT	v1;
+	int	i;
+					/* randomize other coordinates */
+	v1[0] = 0.5 - frandom();
+	v1[1] = 0.5 - frandom();
+	v1[2] = 0.5 - frandom();
+	for (i = 3; i--; )
+		if ((-0.6 < v[i]) & (v[i] < 0.6))
+			break;
+	if (i < 0)
+		return(0);
+	v1[i] = 1.0;
+	fcross(vp, v1, v);
+	return(normalize(vp) > 0.0);
+}
+
+
+int
 closestapproach(			/* closest approach of two rays */
 RREAL t[2],		/* returned distances along each ray */
 const FVECT rorg0,		/* first origin */
@@ -180,7 +224,7 @@ double theta		/* right-hand radians */
 	cost = cos(theta);
 	sint = sin(theta);
 	normprod = DOT(vorig, vnorm)*(1.-cost);
-	fcross(vperp, vnorm, vorig);
+	VCROSS(vperp, vnorm, vorig);
 	for (i = 0; i < 3; i++)
 		vres[i] = vorig[i]*cost + vnorm[i]*normprod + vperp[i]*sint;
 }
@@ -195,16 +239,16 @@ int meas		/* distance measure (radians, absolute, rela
 )
 {
 	FVECT	normtarg;
-	double	volen, dotprod, sint, cost;
+	double	volen, dotprod, sintr, cost;
 	int	i;
 
+	VCOPY(normtarg, vtarg);		/* in case vtarg==vres */
 	if (vres != vorig)
 		VCOPY(vres, vorig);
 	if (t == 0.0)
 		return(VLEN(vres));	/* no rotation requested */
 	if ((volen = normalize(vres)) == 0.0)
 		return(0.0);
-	VCOPY(normtarg, vtarg);
 	if (normalize(normtarg) == 0.0)
 		return(0.0);		/* target vector is zero */
 	dotprod = DOT(vres, normtarg);
@@ -222,10 +266,10 @@ int meas		/* distance measure (radians, absolute, rela
 	else if (meas == GEOD_REL)
 		t *= acos(dotprod);
 	cost = cos(t);
-	sint = sin(t);
+	sintr = sin(t) / sqrt(1. - dotprod*dotprod);
 	for (i = 0; i < 3; i++)
 		vres[i] = volen*( cost*vres[i] +
-				  sint*(normtarg[i] - dotprod*vres[i]) );
+				  sintr*(normtarg[i] - dotprod*vres[i]) );
 
 	return(volen);			/* return vector length */
 }