--- ray/src/common/fvect.c 2012/09/06 00:07:43 2.15 +++ ray/src/common/fvect.c 2013/06/29 21:03:44 2.19 @@ -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.19 2013/06/29 21:03:44 greg Exp $"; #endif /* * fvect.c - routines for floating-point vector calculations @@ -7,11 +7,31 @@ static const char RCSid[] = "$Id: fvect.c,v 2.15 2012/ #include "copyright.h" +#define _USE_MATH_DEFINES #include #include "fvect.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 +49,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 +105,7 @@ 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]; + VCROSS(vres, v1, v2); } @@ -101,9 +117,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); } @@ -180,7 +194,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 +209,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 +236,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 */ }