1 |
|
|
/* |
2 |
|
|
Teem: Tools to process and visualize scientific data and images . |
3 |
|
|
Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago |
4 |
|
|
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann |
5 |
|
|
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah |
6 |
|
|
|
7 |
|
|
This library is free software; you can redistribute it and/or |
8 |
|
|
modify it under the terms of the GNU Lesser General Public License |
9 |
|
|
(LGPL) as published by the Free Software Foundation; either |
10 |
|
|
version 2.1 of the License, or (at your option) any later version. |
11 |
|
|
The terms of redistributing and/or modifying this software also |
12 |
|
|
include exceptions to the LGPL that facilitate static linking. |
13 |
|
|
|
14 |
|
|
This library is distributed in the hope that it will be useful, |
15 |
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of |
16 |
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
17 |
|
|
Lesser General Public License for more details. |
18 |
|
|
|
19 |
|
|
You should have received a copy of the GNU Lesser General Public License |
20 |
|
|
along with this library; if not, write to Free Software Foundation, Inc., |
21 |
|
|
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
22 |
|
|
*/ |
23 |
|
|
|
24 |
|
|
#include "echo.h" |
25 |
|
|
#include "privateEcho.h" |
26 |
|
|
|
27 |
|
|
echoPos_t |
28 |
|
|
_echo_SuperquadX_v(echoPos_t x, echoPos_t y, echoPos_t z, |
29 |
|
|
echoPos_t A, echoPos_t B) { |
30 |
|
|
echoPos_t xxb, yya, zza; |
31 |
|
|
|
32 |
|
|
xxb = pow(x*x, 1/B); yya = pow(y*y, 1/A); zza = pow(z*z, 1/A); |
33 |
|
|
return pow(yya + zza, A/B) + xxb - 1; |
34 |
|
|
} |
35 |
|
|
|
36 |
|
|
echoPos_t |
37 |
|
|
_echo_SuperquadY_v(echoPos_t x, echoPos_t y, echoPos_t z, |
38 |
|
|
echoPos_t A, echoPos_t B) { |
39 |
|
|
echoPos_t xxa, yyb, zza; |
40 |
|
|
|
41 |
|
|
xxa = pow(x*x, 1/A); yyb = pow(y*y, 1/B); zza = pow(z*z, 1/A); |
42 |
|
|
return pow(xxa + zza, A/B) + yyb - 1; |
43 |
|
|
} |
44 |
|
|
|
45 |
|
|
echoPos_t |
46 |
|
|
_echo_SuperquadZ_v(echoPos_t x, echoPos_t y, echoPos_t z, |
47 |
|
|
echoPos_t A, echoPos_t B) { |
48 |
|
|
echoPos_t xxa, yya, zzb; |
49 |
|
|
|
50 |
|
|
xxa = pow(x*x, 1/A); yya = pow(y*y, 1/A); zzb = pow(z*z, 1/B); |
51 |
|
|
return pow(xxa + yya, A/B) + zzb - 1; |
52 |
|
|
} |
53 |
|
|
|
54 |
|
|
/* -------------------------------------------------------- */ |
55 |
|
|
|
56 |
|
|
echoPos_t |
57 |
|
|
_echo_SuperquadX_vg(echoPos_t grad[3], |
58 |
|
|
echoPos_t x, echoPos_t y, echoPos_t z, |
59 |
|
|
echoPos_t A, echoPos_t B) { |
60 |
|
|
echoPos_t R, xxb, yya, zza; |
61 |
|
|
|
62 |
|
|
xxb = pow(x*x, 1/B); yya = pow(y*y, 1/A); zza = pow(z*z, 1/A); |
63 |
|
|
R = pow(yya + zza, (A/B)-1); |
64 |
|
|
ELL_3V_SET(grad, 2*xxb/(B*x), 2*R*yya/(B*y), 2*R*zza/(B*z)); |
65 |
|
|
return pow(yya + zza, A/B) + xxb - 1; |
66 |
|
|
} |
67 |
|
|
|
68 |
|
|
echoPos_t |
69 |
|
|
_echo_SuperquadY_vg(echoPos_t grad[3], |
70 |
|
|
echoPos_t x, echoPos_t y, echoPos_t z, |
71 |
|
|
echoPos_t A, echoPos_t B) { |
72 |
|
|
echoPos_t R, xxa, yyb, zza; |
73 |
|
|
|
74 |
|
|
xxa = pow(x*x, 1/A); yyb = pow(y*y, 1/B); zza = pow(z*z, 1/A); |
75 |
|
|
R = pow(xxa + zza, (A/B)-1); |
76 |
|
|
ELL_3V_SET(grad, 2*R*xxa/(B*x), 2*yyb/(B*y), 2*R*zza/(B*z)); |
77 |
|
|
return pow(xxa + zza, A/B) + yyb - 1; |
78 |
|
|
} |
79 |
|
|
|
80 |
|
|
echoPos_t |
81 |
|
|
_echo_SuperquadZ_vg(echoPos_t grad[3], |
82 |
|
|
echoPos_t x, echoPos_t y, echoPos_t z, |
83 |
|
|
echoPos_t A, echoPos_t B) { |
84 |
|
|
echoPos_t R, xxa, yya, zzb; |
85 |
|
|
|
86 |
|
|
xxa = pow(x*x, 1/A); yya = pow(y*y, 1/A); zzb = pow(z*z, 1/B); |
87 |
|
|
R = pow(xxa + yya, (A/B)-1); |
88 |
|
|
ELL_3V_SET(grad, 2*R*xxa/(B*x), 2*R*yya/(B*y), 2*zzb/(B*z)); |
89 |
|
|
return pow(xxa + yya, A/B) + zzb - 1; |
90 |
|
|
} |
91 |
|
|
|
92 |
|
|
/* -------------------------------------------------------- */ |
93 |
|
|
|
94 |
|
|
echoPos_t |
95 |
|
|
_echo_SuperquadX_lvg(echoPos_t grad[3], |
96 |
|
|
echoPos_t x, echoPos_t y, echoPos_t z, |
97 |
|
|
echoPos_t A, echoPos_t B) { |
98 |
|
|
echoPos_t R, xxb, yya, zza, larg; |
99 |
|
|
|
100 |
|
|
echoPos_t ret; |
101 |
|
|
|
102 |
|
|
xxb = pow(x*x, 1/B); yya = pow(y*y, 1/A); zza = pow(z*z, 1/A); |
103 |
|
|
R = pow(yya + zza, 1-(A/B))*xxb; |
104 |
|
|
ELL_3V_SET(grad, |
105 |
|
|
2/(B*x*(1 + pow(yya + zza, A/B)/xxb)), |
106 |
|
|
2*yya/(B*y*(yya + zza + R)), |
107 |
|
|
2*zza/(B*z*(yya + zza + R))); |
108 |
|
|
larg = pow(yya + zza, A/B) + xxb; |
109 |
|
|
ret= larg > 0 ? log(larg) : ECHO_POS_MIN; |
110 |
|
|
/* |
111 |
|
|
if (!( AIR_EXISTS(grad[0]) && AIR_EXISTS(grad[1]) && AIR_EXISTS(grad[2]) )) { |
112 |
|
|
fprintf(stderr, "_echo_SuperquadX_lvg: PANIC\n"); |
113 |
|
|
fprintf(stderr, "x = %g, y = %g, z = %g, A = %g, B = %g\n", |
114 |
|
|
x, y, z, A, B); |
115 |
|
|
fprintf(stderr, "pow(%g * %g = %g, 1/%g = %g) = %g\n", |
116 |
|
|
x, x, x*x, B, 1/B, pow(x*x, 1/B)); |
117 |
|
|
fprintf(stderr, "xxb = %g, yya = %g, zza = %g\n", |
118 |
|
|
xxb, yya, zza); |
119 |
|
|
fprintf(stderr, "R: pow(%g + %g = %g, 1-(%g/%g = %g) = %g) = %g*%g = %g\n", |
120 |
|
|
yya, zza, yya + zza, |
121 |
|
|
A, B, A/B, 1-(A/B), |
122 |
|
|
pow(yya + zza, 1-(A/B)), xxb, |
123 |
|
|
pow(yya + zza, 1-(A/B))*xxb); |
124 |
|
|
fprintf(stderr, "grad[0]: 2/(%g * %g * (1 + pow(%g + %g = %g, %g/%g = %g)/%g = %g)) = %g\n", |
125 |
|
|
B, x, yya, zza, yya+zza, A, B, A/B, xxb, |
126 |
|
|
pow(yya + zza, A/B)/xxb, grad[0]); |
127 |
|
|
fprintf(stderr, "grad[1]: 2*%g/(%g*%g*(%g + %g + %g = %g) = %g) = %g\n", |
128 |
|
|
yya, B, y, yya, zza, R, yya + zza + R, |
129 |
|
|
B*y*(yya + zza + R), |
130 |
|
|
2*yya/(B*y*(yya + zza + R))); |
131 |
|
|
|
132 |
|
|
fprintf(stderr, "log(pow(%g + %g = %g, %g) = %g + %g) = %g\n", |
133 |
|
|
yya, zza, yya+zza, A/B, pow(yya + zza, A/B), xxb, ret); |
134 |
|
|
fprintf(stderr, "grad = %g %g %g\n", grad[0], grad[1], grad[2]); |
135 |
|
|
fprintf(stderr, "\n----------\n\n"); |
136 |
|
|
} |
137 |
|
|
*/ |
138 |
|
|
return ret; |
139 |
|
|
} |
140 |
|
|
|
141 |
|
|
echoPos_t |
142 |
|
|
_echo_SuperquadY_lvg(echoPos_t grad[3], |
143 |
|
|
echoPos_t x, echoPos_t y, echoPos_t z, |
144 |
|
|
echoPos_t A, echoPos_t B) { |
145 |
|
|
echoPos_t R, xxa, yyb, zza, larg; |
146 |
|
|
|
147 |
|
|
xxa = pow(x*x, 1/A); yyb = pow(y*y, 1/B); zza = pow(z*z, 1/A); |
148 |
|
|
R = pow(xxa + zza, 1-(A/B))*yyb; |
149 |
|
|
ELL_3V_SET(grad, |
150 |
|
|
2*xxa/(B*x*(xxa + zza + R)), |
151 |
|
|
2/(B*y*(1 + pow(xxa + zza, A/B)/yyb)), |
152 |
|
|
2*zza/(B*z*(xxa + zza + R))); |
153 |
|
|
larg = pow(xxa + zza, A/B) + yyb; |
154 |
|
|
return larg > 0 ? log(larg) : ECHO_POS_MIN; |
155 |
|
|
} |
156 |
|
|
|
157 |
|
|
echoPos_t |
158 |
|
|
_echo_SuperquadZ_lvg(echoPos_t grad[3], |
159 |
|
|
echoPos_t x, echoPos_t y, echoPos_t z, |
160 |
|
|
echoPos_t A, echoPos_t B) { |
161 |
|
|
echoPos_t R, xxa, yya, zzb, larg; |
162 |
|
|
|
163 |
|
|
echoPos_t ret; |
164 |
|
|
|
165 |
|
|
xxa = pow(x*x, 1/A); yya = pow(y*y, 1/A); zzb = pow(z*z, 1/B); |
166 |
|
|
R = pow(xxa + yya, 1-(A/B))*zzb; |
167 |
|
|
ELL_3V_SET(grad, |
168 |
|
|
2*xxa/(B*x*(xxa + yya + R)), |
169 |
|
|
2*yya/(B*y*(xxa + yya + R)), |
170 |
|
|
2/(B*z*(1 + pow(xxa + yya, A/B)/zzb))); |
171 |
|
|
larg = pow(xxa + yya, A/B) + zzb; |
172 |
|
|
|
173 |
|
|
ret = larg > 0 ? log(larg) : ECHO_POS_MIN; |
174 |
|
|
/* |
175 |
|
|
if (!AIR_EXISTS(ret)) { |
176 |
|
|
fprintf(stderr, "_echo_SuperquadZ_lvg: PANIC\n"); |
177 |
|
|
fprintf(stderr, "x = %g, y = %g, z = %g, A = %g, B = %g\n", |
178 |
|
|
x, y, z, A, B); |
179 |
|
|
fprintf(stderr, "pow(%g*%g = %g, %g) = %g\n", |
180 |
|
|
x, x, x*x, 1/A, xxa); |
181 |
|
|
fprintf(stderr, "pow(%g*%g = %g, %g) = %g\n", |
182 |
|
|
y, y, y*y, 1/A, yya); |
183 |
|
|
fprintf(stderr, "log(pow(%g, %g) = %g + %g) = %g\n", |
184 |
|
|
xxa + yya, A/B, pow(xxa + yya, A/B), zzb, ret); |
185 |
|
|
exit(0); |
186 |
|
|
} |
187 |
|
|
*/ |
188 |
|
|
return ret; |
189 |
|
|
} |
190 |
|
|
|
191 |
|
|
/* -------------------------------------------------------- */ |
192 |
|
|
|
193 |
|
|
int |
194 |
|
|
_echoRayIntx_Superquad(RAYINTX_ARGS(Superquad)) { |
195 |
|
|
char me[]="_echoRayIntx_Superquad"; |
196 |
|
|
echoPos_t TT=0, Tmin, Tmax, t0, t1, t2, t3, v1, v2, diff, tol, |
197 |
|
|
saveTmin, Vmin, Vmax, VV=0, dV, dVmin, dVmax, tmp, |
198 |
|
|
(*v)(echoPos_t, echoPos_t, echoPos_t, |
199 |
|
|
echoPos_t, echoPos_t), |
200 |
|
|
(*vg)(echoPos_t[3], |
201 |
|
|
echoPos_t, echoPos_t, echoPos_t, |
202 |
|
|
echoPos_t, echoPos_t), |
203 |
|
|
(*lvg)(echoPos_t[3], |
204 |
|
|
echoPos_t, echoPos_t, echoPos_t, |
205 |
|
|
echoPos_t, echoPos_t), |
206 |
|
|
from[3], grad[3], pos[3]; /* these two used only by macros */ |
207 |
|
|
int iter; |
208 |
|
|
|
209 |
|
|
if (!_echoRayIntx_CubeSolid(&Tmin, &Tmax, |
210 |
|
|
-1-2*ECHO_EPSILON, 1+2*ECHO_EPSILON, |
211 |
|
|
-1-2*ECHO_EPSILON, 1+2*ECHO_EPSILON, |
212 |
|
|
-1-2*ECHO_EPSILON, 1+2*ECHO_EPSILON, ray)) { |
213 |
|
|
return AIR_FALSE; |
214 |
|
|
} |
215 |
|
|
switch(obj->axis) { |
216 |
|
|
case 0: |
217 |
|
|
v = _echo_SuperquadX_v; |
218 |
|
|
vg = _echo_SuperquadX_vg; |
219 |
|
|
lvg = _echo_SuperquadX_lvg; |
220 |
|
|
break; |
221 |
|
|
case 1: |
222 |
|
|
v = _echo_SuperquadY_v; |
223 |
|
|
vg = _echo_SuperquadY_vg; |
224 |
|
|
lvg = _echo_SuperquadY_lvg; |
225 |
|
|
break; |
226 |
|
|
case 2: default: |
227 |
|
|
v = _echo_SuperquadZ_v; |
228 |
|
|
vg = _echo_SuperquadZ_vg; |
229 |
|
|
lvg = _echo_SuperquadZ_lvg; |
230 |
|
|
break; |
231 |
|
|
} |
232 |
|
|
if (tstate->verbose) { |
233 |
|
|
fprintf(stderr, "%s%s: Tmin, Tmax = %g, %g, ax = %d\n", |
234 |
|
|
_echoDot(tstate->depth), me, Tmin, Tmax, obj->axis); |
235 |
|
|
} |
236 |
|
|
|
237 |
|
|
#define VAL(TT) \ |
238 |
|
|
(ELL_3V_SCALE_ADD2(pos, 1, from, (TT), ray->dir), \ |
239 |
|
|
v(pos[0], pos[1], pos[2], obj->A, obj->B)) |
240 |
|
|
|
241 |
|
|
#define VALGRAD(VV, DV, TT) \ |
242 |
|
|
ELL_3V_SCALE_ADD2(pos, 1, from, (TT), ray->dir); \ |
243 |
|
|
(VV) = vg(grad, pos[0], pos[1], pos[2], obj->A, obj->B); \ |
244 |
|
|
(DV) = ELL_3V_DOT(grad, ray->dir) |
245 |
|
|
|
246 |
|
|
#define LVALGRAD(VV, DV, TT) \ |
247 |
|
|
ELL_3V_SCALE_ADD2(pos, 1, from, (TT), ray->dir); \ |
248 |
|
|
(VV) = lvg(grad, pos[0], pos[1], pos[2], obj->A, obj->B); \ |
249 |
|
|
(DV) = ELL_3V_DOT(grad, ray->dir) |
250 |
|
|
|
251 |
|
|
#define RR 0.61803399 |
252 |
|
|
#define CC (1.0-RR) |
253 |
|
|
#define SHIFT3(a,b,c,d) (a)=(b); (b)=(c); (c)=(d) |
254 |
|
|
#define SHIFT2(a,b,c) (a)=(b); (b)=(c) |
255 |
|
|
|
256 |
|
|
/* testing */ |
257 |
|
|
ELL_3V_SCALE_ADD2(from, 1, ray->from, Tmin, ray->dir); |
258 |
|
|
saveTmin = Tmin; |
259 |
|
|
Tmin = 0; |
260 |
|
|
Tmax -= saveTmin; |
261 |
|
|
/* |
262 |
|
|
ELL_3V_COPY(from, ray->from); |
263 |
|
|
saveTmin = 0; |
264 |
|
|
*/ |
265 |
|
|
|
266 |
|
|
/* evaluate value and derivatives at Tmin and Tmax */ |
267 |
|
|
VALGRAD(Vmin, dVmin, Tmin); |
268 |
|
|
VALGRAD(Vmax, dVmax, Tmax); |
269 |
|
|
|
270 |
|
|
/* if the segment start and end are both positive or both negative, |
271 |
|
|
and if the derivatives also don't change sign, there's no root. |
272 |
|
|
Also, due to convexity, if values at start and end are negative, |
273 |
|
|
then there is no root */ |
274 |
|
|
if ( (Vmin*Vmax >= 0 && dVmin*dVmax >= 0) |
275 |
|
|
|| (Vmin <= 0 && Vmax <= 0) ) { |
276 |
|
|
return AIR_FALSE; |
277 |
|
|
} |
278 |
|
|
if (tstate->verbose) { |
279 |
|
|
fprintf(stderr, "%s%s: dVmin = %g, dVmax = %g, Vmin = %g, Vmax = %g\n", |
280 |
|
|
_echoDot(tstate->depth), me, dVmin, dVmax, Vmin, Vmax); |
281 |
|
|
} |
282 |
|
|
|
283 |
|
|
/* either the value changed sign, or the derivative changed sign, or |
284 |
|
|
both. If, as is common, the derivatives changed sign, but the |
285 |
|
|
values didn't (which means they are both positive, due to a test |
286 |
|
|
above), we need to limit the interval by minimizing the value |
287 |
|
|
until either we see a negative value, or until the minimization |
288 |
|
|
converged. Based on Golden Section Search, NR pg.401 */ |
289 |
|
|
if (dVmin*dVmax < 0 && Vmin*Vmax >= 0) { |
290 |
|
|
t0 = Tmin; |
291 |
|
|
t1 = RR*Tmin + CC*Tmax; |
292 |
|
|
t2 = CC*Tmin + RR*Tmax; |
293 |
|
|
t3 = Tmax; |
294 |
|
|
v1 = VAL(t1); |
295 |
|
|
v2 = VAL(t2); |
296 |
|
|
if (tstate->verbose) { |
297 |
|
|
fprintf(stderr, "%s%s: \n" |
298 |
|
|
" t0 = % 31.15f\n" |
299 |
|
|
" t1 = % 31.15f -> v1 = % 31.15f\n" |
300 |
|
|
" t2 = % 31.15f -> v2 = % 31.15f\n" |
301 |
|
|
" t3 = % 31.15f\n", |
302 |
|
|
_echoDot(tstate->depth), me, |
303 |
|
|
t0, t1, v1, t2, v2, t3); |
304 |
|
|
} |
305 |
|
|
tol = sqrt(ECHO_POS_EPS); |
306 |
|
|
while ( (t3-t0 > tol*(t1+t2)) /* still haven't converged */ |
307 |
|
|
&& (v1 > 0 && v2 > 0) ) { /* v1 and v2 both positive */ |
308 |
|
|
diff = v2 - v1; |
309 |
|
|
if (v1 < v2) { |
310 |
|
|
SHIFT3(t3, t2, t1, CC*t0 + RR*t2); |
311 |
|
|
SHIFT2(v2, v1, VAL(t1)); |
312 |
|
|
} else { |
313 |
|
|
SHIFT3(t0, t1, t2, RR*t1 + CC*t3); |
314 |
|
|
SHIFT2(v1, v2, VAL(t2)); |
315 |
|
|
} |
316 |
|
|
if (tstate->verbose) { |
317 |
|
|
fprintf(stderr, "%s%s: %s ---> \n" |
318 |
|
|
" t0 = % 31.15f\n" |
319 |
|
|
" t1 = % 31.15f -> v1 = % 31.15f\n" |
320 |
|
|
" t2 = % 31.15f -> v2 = % 31.15f\n" |
321 |
|
|
" t3 = % 31.15f\n", |
322 |
|
|
_echoDot(tstate->depth), me, |
323 |
|
|
diff > 0 ? "v1 < v2" : "v1 > v2", |
324 |
|
|
t0, t1, v1, t2, v2, t3); |
325 |
|
|
} |
326 |
|
|
} |
327 |
|
|
if (v1 > 0 && v2 > 0) { |
328 |
|
|
/* the minimization stopped, yet both v1 and v2 are still positive, |
329 |
|
|
so there's no way we can have a root */ |
330 |
|
|
if (tstate->verbose) { |
331 |
|
|
fprintf(stderr, "%s%s: minimization found no root\n", |
332 |
|
|
_echoDot(tstate->depth), me); |
333 |
|
|
} |
334 |
|
|
return AIR_FALSE; |
335 |
|
|
} |
336 |
|
|
/* else either v1 or v2 <= 0, so there is a root (actually two). |
337 |
|
|
By construction, f(t0) is positive, so we can now bracket the |
338 |
|
|
root between t0 and t1 or t2, whichever one is associated with |
339 |
|
|
a smaller value */ |
340 |
|
|
Tmin = t0; |
341 |
|
|
/* HEY: shouldn't I just be using whichever one is closer? */ |
342 |
|
|
Tmax = v1 < v2 ? t1 : t2; |
343 |
|
|
Vmin = VAL(Tmin); |
344 |
|
|
Vmax = VAL(Tmax); |
345 |
|
|
} |
346 |
|
|
|
347 |
|
|
/* the value isn't necessarily monotonic between Tmin and Tmax, but |
348 |
|
|
we know that there is only one root. Find it with newton-raphson, |
349 |
|
|
using the log of function, both for values and for derivatives */ |
350 |
|
|
iter = 0; |
351 |
|
|
TT = (Tmin + Tmax)/2; |
352 |
|
|
LVALGRAD(VV, dV, TT); |
353 |
|
|
|
354 |
|
|
while (iter < parm->sqNRI && AIR_ABS(VV) > parm->sqTol |
355 |
|
|
&& AIR_EXISTS(VV) && AIR_EXISTS(dV)) { |
356 |
|
|
if (tstate->verbose) { |
357 |
|
|
fprintf(stderr, "%s%s: iter = %d: TT = %g, VV = %g, dV = %g\n", |
358 |
|
|
_echoDot(tstate->depth), me, iter, TT, VV, dV); |
359 |
|
|
} |
360 |
|
|
TT -= VV/dV; |
361 |
|
|
if (!AIR_IN_OP(Tmin, TT, Tmax)) { |
362 |
|
|
/* newton-raphson sent us flying out of bounds; find a tighter |
363 |
|
|
[Tmin,Tmax] bracket with bisection and try again */ |
364 |
|
|
TT = (Tmin + Tmax)/2; |
365 |
|
|
VV = VAL(TT); |
366 |
|
|
if (Vmin*VV < 0) { |
367 |
|
|
Tmax = TT; |
368 |
|
|
Vmax = VV; |
369 |
|
|
} else { |
370 |
|
|
Tmin = TT; |
371 |
|
|
Vmin = VV; |
372 |
|
|
} |
373 |
|
|
TT = (Tmin + Tmax)/2; |
374 |
|
|
} |
375 |
|
|
LVALGRAD(VV, dV, TT); |
376 |
|
|
iter++; |
377 |
|
|
} |
378 |
|
|
|
379 |
|
|
if (!( AIR_EXISTS(VV) && AIR_EXISTS(dV) )) { |
380 |
|
|
/* we bailed out of the loop above because |
381 |
|
|
we got screwed by numerical errors |
382 |
|
|
--> pretend that there was no intersection, |
383 |
|
|
and HEY this will have to be debugged later */ |
384 |
|
|
if (tstate->verbose) { |
385 |
|
|
fprintf(stderr, "%s%s: SORRY, numerical problems!\n", |
386 |
|
|
_echoDot(tstate->depth), me); |
387 |
|
|
} |
388 |
|
|
return AIR_FALSE; |
389 |
|
|
} |
390 |
|
|
|
391 |
|
|
/* else we succeeded in finding the intersection */ |
392 |
|
|
intx->t = TT + saveTmin; |
393 |
|
|
VALGRAD(VV, dV, TT); /* puts gradient into grad */ |
394 |
|
|
ELL_3V_NORM(intx->norm, grad, tmp); |
395 |
|
|
intx->obj = OBJECT(obj); |
396 |
|
|
/* set in intx: |
397 |
|
|
yes: t, norm |
398 |
|
|
no: u, v |
399 |
|
|
*/ |
400 |
|
|
return AIR_TRUE; |
401 |
|
|
} |