7 |
|
|
8 |
|
#include "copyright.h" |
9 |
|
|
10 |
+ |
#define _USE_MATH_DEFINES |
11 |
|
#include <math.h> |
12 |
|
#include "fvect.h" |
13 |
+ |
#include "random.h" |
14 |
|
|
15 |
+ |
double |
16 |
+ |
Acos(double x) /* insurance for touchy math library */ |
17 |
+ |
{ |
18 |
+ |
if (x <= -1.+FTINY*FTINY) |
19 |
+ |
return(M_PI); |
20 |
+ |
if (x >= 1.-FTINY*FTINY) |
21 |
+ |
return(.0); |
22 |
+ |
return(acos(x)); |
23 |
+ |
} |
24 |
|
|
25 |
|
double |
26 |
+ |
Asin(double x) /* insurance for touchy math library */ |
27 |
+ |
{ |
28 |
+ |
if (x <= -1.+FTINY*FTINY) |
29 |
+ |
return(-M_PI/2.); |
30 |
+ |
if (x >= 1.-FTINY*FTINY) |
31 |
+ |
return(M_PI/2); |
32 |
+ |
return(asin(x)); |
33 |
+ |
} |
34 |
+ |
|
35 |
+ |
double |
36 |
|
fdot( /* return the dot product of two vectors */ |
37 |
|
const FVECT v1, |
38 |
|
const FVECT v2 |
50 |
|
{ |
51 |
|
FVECT delta; |
52 |
|
|
53 |
< |
delta[0] = p2[0] - p1[0]; |
33 |
< |
delta[1] = p2[1] - p1[1]; |
34 |
< |
delta[2] = p2[2] - p1[2]; |
53 |
> |
VSUB(delta, p2, p1); |
54 |
|
|
55 |
|
return(DOT(delta, delta)); |
56 |
|
} |
106 |
|
const FVECT v2 |
107 |
|
) |
108 |
|
{ |
109 |
< |
vres[0] = v1[1]*v2[2] - v1[2]*v2[1]; |
110 |
< |
vres[1] = v1[2]*v2[0] - v1[0]*v2[2]; |
111 |
< |
vres[2] = v1[0]*v2[1] - v1[1]*v2[0]; |
109 |
> |
if ((vres == v1) | (vres == v2)) { |
110 |
> |
FVECT vtmp; |
111 |
> |
VCROSS(vtmp, v1, v2); |
112 |
> |
VCOPY(vres, vtmp); |
113 |
> |
return; |
114 |
> |
} |
115 |
> |
VCROSS(vres, v1, v2); |
116 |
|
} |
117 |
|
|
118 |
|
|
124 |
|
double f |
125 |
|
) |
126 |
|
{ |
127 |
< |
vres[0] = v0[0] + f*v1[0]; |
105 |
< |
vres[1] = v0[1] + f*v1[1]; |
106 |
< |
vres[2] = v0[2] + f*v1[2]; |
127 |
> |
VSUM(vres, v0, v1, f); |
128 |
|
} |
129 |
|
|
130 |
|
|
156 |
|
|
157 |
|
|
158 |
|
int |
159 |
+ |
getperpendicular( /* choose perpedicular direction */ |
160 |
+ |
FVECT vp, /* returns normalized */ |
161 |
+ |
const FVECT v, /* input vector must be normalized */ |
162 |
+ |
int randomize /* randomize orientation */ |
163 |
+ |
) |
164 |
+ |
{ |
165 |
+ |
int ord[3]; |
166 |
+ |
FVECT v1; |
167 |
+ |
int i; |
168 |
+ |
|
169 |
+ |
if (randomize) { /* randomize coordinates? */ |
170 |
+ |
v1[0] = 0.5 - frandom(); |
171 |
+ |
v1[1] = 0.5 - frandom(); |
172 |
+ |
v1[2] = 0.5 - frandom(); |
173 |
+ |
switch ((int)(frandom()*5.99999)) { |
174 |
+ |
case 0: ord[0] = 0; ord[1] = 1; ord[2] = 2; break; |
175 |
+ |
case 1: ord[0] = 0; ord[1] = 2; ord[2] = 1; break; |
176 |
+ |
case 2: ord[0] = 1; ord[1] = 0; ord[2] = 2; break; |
177 |
+ |
case 3: ord[0] = 1; ord[1] = 2; ord[2] = 0; break; |
178 |
+ |
case 4: ord[0] = 2; ord[1] = 0; ord[2] = 1; break; |
179 |
+ |
case 5: ord[0] = 2; ord[1] = 1; ord[2] = 0; break; |
180 |
+ |
} |
181 |
+ |
} else { |
182 |
+ |
v1[0] = v1[1] = v1[2] = 0.0; |
183 |
+ |
ord[0] = 0; ord[1] = 1; ord[2] = 2; |
184 |
+ |
} |
185 |
+ |
|
186 |
+ |
for (i = 3; i--; ) |
187 |
+ |
if ((-0.6 < v[ord[i]]) & (v[ord[i]] < 0.6)) |
188 |
+ |
break; |
189 |
+ |
if (i < 0) |
190 |
+ |
return(0); |
191 |
+ |
|
192 |
+ |
v1[ord[i]] = 1.0; |
193 |
+ |
fcross(vp, v1, v); |
194 |
+ |
|
195 |
+ |
return(normalize(vp) > 0.0); |
196 |
+ |
} |
197 |
+ |
|
198 |
+ |
|
199 |
+ |
int |
200 |
|
closestapproach( /* closest approach of two rays */ |
201 |
|
RREAL t[2], /* returned distances along each ray */ |
202 |
|
const FVECT rorg0, /* first origin */ |
242 |
|
cost = cos(theta); |
243 |
|
sint = sin(theta); |
244 |
|
normprod = DOT(vorig, vnorm)*(1.-cost); |
245 |
< |
fcross(vperp, vnorm, vorig); |
245 |
> |
VCROSS(vperp, vnorm, vorig); |
246 |
|
for (i = 0; i < 3; i++) |
247 |
|
vres[i] = vorig[i]*cost + vnorm[i]*normprod + vperp[i]*sint; |
248 |
|
} |
257 |
|
) |
258 |
|
{ |
259 |
|
FVECT normtarg; |
260 |
< |
double volen, dotprod, sint, cost; |
260 |
> |
double volen, dotprod, sintr, cost; |
261 |
|
int i; |
262 |
|
|
263 |
|
VCOPY(normtarg, vtarg); /* in case vtarg==vres */ |
284 |
|
else if (meas == GEOD_REL) |
285 |
|
t *= acos(dotprod); |
286 |
|
cost = cos(t); |
287 |
< |
sint = sin(t); |
287 |
> |
sintr = sin(t) / sqrt(1. - dotprod*dotprod); |
288 |
|
for (i = 0; i < 3; i++) |
289 |
|
vres[i] = volen*( cost*vres[i] + |
290 |
< |
sint*(normtarg[i] - dotprod*vres[i]) ); |
290 |
> |
sintr*(normtarg[i] - dotprod*vres[i]) ); |
291 |
|
|
292 |
|
return(volen); /* return vector length */ |
293 |
|
} |