Garfield++ v1r0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Numerics.cc
Go to the documentation of this file.
1#include <iostream>
2#include <cmath>
3
4#include "Numerics.hh"
5
6namespace Garfield {
7
8namespace Numerics {
9
10void Dfact(const int n, std::vector<std::vector<double> >& a,
11 std::vector<int>& ir, int& ifail, double& det, int& jfail) {
12
13 const double g1 = 1.e-19;
14 const double g2 = 1.e-19;
15
16 double tf, p, q, t, s11, s12;
17 int k;
18
19 ifail = jfail = 0;
20
21 int nxch = 0;
22 det = 1.;
23
24 for (int j = 1; j <= n; ++j) {
25 k = j;
26 p = fabs(a[j - 1][j - 1]);
27 if (j == n) {
28 if (p <= 0.) {
29 det = 0.;
30 ifail = -1;
31 jfail = 0;
32 return;
33 }
34 det *= a[j - 1][j - 1];
35 a[j - 1][j - 1] = 1. / a[j - 1][j - 1];
36 t = fabs(det);
37 if (t < g1) {
38 det = 0.;
39 if (jfail == 0) jfail = -1;
40 } else if (t > g2) {
41 det = 1.;
42 if (jfail == 0) jfail = +1;
43 }
44 continue;
45 }
46 for (int i = j + 1; i <= n; ++i) {
47 q = fabs(a[i - 1][j - 1]);
48 if (q <= p) continue;
49 k = i;
50 p = q;
51 }
52 if (k != j) {
53 for (int l = 1; l <= n; ++l) {
54 tf = a[j - 1][l - 1];
55 a[j - 1][l - 1] = a[k - 1][l - 1];
56 a[k - 1][l - 1] = tf;
57 }
58 ++nxch;
59 ir[nxch - 1] = j * 4096 + k;
60 } else if (p <= 0.) {
61 det = 0.;
62 ifail = -1;
63 jfail = 0;
64 return;
65 }
66 det *= a[j - 1][j - 1];
67 a[j - 1][j - 1] = 1. / a[j - 1][j - 1];
68 t = fabs(det);
69 if (t < g1) {
70 det = 0.;
71 if (jfail == 0) jfail = -1;
72 } else if (t > g2) {
73 det = 1.;
74 if (jfail == 0) jfail = +1;
75 }
76 for (k = j + 1; k <= n; ++k) {
77 s11 = -a[j - 1][k - 1];
78 s12 = -a[k - 1][j];
79 if (j == 1) {
80 a[j - 1][k - 1] = -s11 * a[j - 1][j - 1];
81 a[k - 1][j] = -(s12 + a[j - 1][j] * a[k - 1][j - 1]);
82 continue;
83 }
84 for (int i = 1; i <= j - 1; ++i) {
85 s11 += a[i - 1][k - 1] * a[j - 1][i - 1];
86 s12 += a[i - 1][j] * a[k - 1][i - 1];
87 }
88 a[j - 1][k - 1] = -s11 * a[j - 1][j - 1];
89 a[k - 1][j] = -a[j - 1][j] * a[k - 1][j - 1] - s12;
90 }
91 }
92
93 if (nxch % 2 != 0) det = -det;
94 if (jfail != 0) det = 0.;
95 ir[n - 1] = nxch;
96}
97
98void Dfeqn(const int n, std::vector<std::vector<double> >& a,
99 std::vector<int>& ir, std::vector<double>& b) {
100
101 double te;
102 double s21, s22;
103
104 if (n <= 0) return;
105
106 int nxch = ir[n - 1];
107 if (nxch != 0) {
108 for (int m = 1; m <= nxch; ++m) {
109 int ij = ir[m - 1];
110 int i = ij / 4096;
111 int j = ij % 4096;
112 te = b[i - 1];
113 b[i - 1] = b[j - 1];
114 b[j - 1] = te;
115 }
116 }
117
118 b[0] *= a[0][0];
119 if (n == 1) return;
120
121 for (int i = 2; i <= n; ++i) {
122 s21 = -b[i - 1];
123 for (int j = 1; j <= i - 1; ++j) {
124 s21 += a[i - 1][j - 1] * b[j - 1];
125 }
126 b[i - 1] = -a[i - 1][i - 1] * s21;
127 }
128
129 for (int i = 1; i <= n - 1; ++i) {
130 s22 = -b[n - i - 1];
131 for (int j = 1; j <= i; ++j) {
132 s22 += a[n - i - 1][n - j] * b[n - j];
133 }
134 b[n - i - 1] = -s22;
135 }
136}
137
138void Dfinv(const int n, std::vector<std::vector<double> >& a,
139 std::vector<int>& ir) {
140
141 double ti;
142 double s31, s32, s33, s34;
143
144 if (n <= 1) return;
145 a[1][0] = -a[1][1] * a[0][0] * a[1][0];
146 a[0][1] = -a[0][1];
147 if (n > 2) {
148 for (int i = 3; i <= n; ++i) {
149 for (int j = 1; j <= i - 2; ++j) {
150 s31 = 0.;
151 s32 = a[j - 1][i - 1];
152 for (int k = j; k <= i - 2; ++k) {
153 s31 += a[k - 1][j - 1] * a[i - 1][k - 1];
154 s32 += a[j - 1][k] * a[k][i - 1];
155 }
156 a[i - 1][j - 1] =
157 -a[i - 1][i - 1] * (s31 + a[i - 2][j - 1] * a[i - 1][i - 2]);
158 a[j - 1][i - 1] = -s32;
159 }
160 a[i - 1][i - 2] = -a[i - 1][i - 1] * a[i - 2][i - 2] * a[i - 1][i - 2];
161 a[i - 2][i - 1] = -a[i - 2][i - 1];
162 }
163 }
164
165 for (int i = 1; i <= n - 1; ++i) {
166 for (int j = 1; j <= i; ++j) {
167 s33 = a[i - 1][j - 1];
168 for (int k = 1; k <= n - i; ++k) {
169 s33 += a[i + k - 1][j - 1] * a[i - 1][i + k - 1];
170 }
171 a[i - 1][j - 1] = s33;
172 }
173 for (int j = 1; j <= n - i; ++j) {
174 s34 = 0.;
175 for (int k = j; k <= n - i; ++k) {
176 s34 += a[i + k - 1][i + j - 1] * a[i - 1][i + k - 1];
177 }
178 a[i - 1][i + j - 1] = s34;
179 }
180 }
181
182 int nxch = ir[n - 1];
183 if (nxch == 0) return;
184
185 for (int m = 1; m <= nxch; ++m) {
186 int k = nxch - m + 1;
187 int ij = ir[k - 1];
188 int i = ij / 4096;
189 int j = ij % 4096;
190 for (k = 1; k <= n; ++k) {
191 ti = a[k - 1][i - 1];
192 a[k - 1][i - 1] = a[k - 1][j - 1];
193 a[k - 1][j - 1] = ti;
194 }
195 }
196}
197
198// ******************************************************************
199//
200// REPLACES B BY THE SOLUTION X OF A*X=B, AND REPLACES A BY ITS IN-
201// VERSE.
202//
203// n ORDER OF THE SQUARE MATRIX IN ARRAY A.
204// A (DOUBLE PRECISION) TWO-DIMENSIONAL ARRAY CONTAINING
205// AN n BY n MATRIX.
206//
207// IFAIL OUTPUT PARAMETER. IFAIL= 0 ... NORMAL EXIT.
208// IFAIL=-1 ... SINGULAR MATRIX.
209//
210// B (DOUBLE PRECISION) ONE-DIMENSIONAL ARRAY
211//
212// CALLS ... DFACT, DFINV.
213//
214// ******************************************************************
215
216void Deqinv(const int n, std::vector<std::vector<double> >& a, int& ifail,
217 std::vector<double>& b) {
218
219 double t1, t2, t3;
220 double det, temp, s;
221 double b1, b2, c11, c12, c13, c21, c22, c23, c31, c32, c33;
222
223 std::vector<int> ir;
224 ir.clear();
225 ir.resize(n);
226 for (int i = 0; i < n; ++i) ir[i] = 0;
227
228 // TEST FOR PARAMETER ERRORS.
229 if (n < 1) {
230 ifail = +1;
231 return;
232 }
233
234 ifail = 0;
235 int jfail = 0;
236
237 if (n > 3) {
238 // n > 3 CASES.
239 // FACTORIZE MATRIX, INVERT AND SOLVE SYSTEM.
240 Dfact(n, a, ir, ifail, det, jfail);
241 if (ifail != 0) return;
242 Dfeqn(n, a, ir, b);
243 Dfinv(n, a, ir);
244 } else if (n == 3) {
245 // n = 3 CASE.
246 // COMPUTE COFACTORS.
247 c11 = a[1][1] * a[2][2] - a[1][2] * a[2][1];
248 c12 = a[1][2] * a[2][0] - a[1][0] * a[2][2];
249 c13 = a[1][0] * a[2][1] - a[1][1] * a[2][0];
250 c21 = a[2][1] * a[0][2] - a[2][2] * a[0][1];
251 c22 = a[2][2] * a[0][0] - a[2][0] * a[0][2];
252 c23 = a[2][0] * a[0][1] - a[2][1] * a[0][0];
253 c31 = a[0][1] * a[1][2] - a[0][2] * a[1][1];
254 c32 = a[0][2] * a[1][0] - a[0][0] * a[1][2];
255 c33 = a[0][0] * a[1][1] - a[0][1] * a[1][0];
256 t1 = fabs(a[0][0]);
257 t2 = fabs(a[1][0]);
258 t3 = fabs(a[2][0]);
259
260 // SET temp = PIVOT AND det = PIVOT * det.
261 if (t1 >= t2) {
262 if (t3 >= t1) {
263 // PIVOT IS A31
264 temp = a[2][0];
265 det = c23 * c12 - c22 * c13;
266 } else {
267 // PIVOT IS A11
268 temp = a[0][0];
269 det = c22 * c33 - c23 * c32;
270 }
271 } else {
272 if (t3 >= t2) {
273 // PIVOT IS A31
274 temp = a[2][0];
275 det = c23 * c12 - c22 * c13;
276 } else {
277 // PIVOT IS A21
278 temp = a[1][0];
279 det = c13 * c32 - c12 * c33;
280 }
281 }
282
283 // SET ELEMENTS OF INVERSE IN A.
284 if (det == 0.) {
285 ifail = -1;
286 return;
287 }
288 s = temp / det;
289 a[0][0] = s * c11;
290 a[0][1] = s * c21;
291 a[0][2] = s * c31;
292 a[1][0] = s * c12;
293 a[1][1] = s * c22;
294 a[1][2] = s * c32;
295 a[2][0] = s * c13;
296 a[2][1] = s * c23;
297 a[2][2] = s * c33;
298
299 // REPLACE B BY AINV*B.
300 b1 = b[0];
301 b2 = b[1];
302 b[0] = a[0][0] * b1 + a[0][1] * b2 + a[0][2] * b[2];
303 b[1] = a[1][0] * b1 + a[1][1] * b2 + a[1][2] * b[2];
304 b[2] = a[2][0] * b1 + a[2][1] * b2 + a[2][2] * b[2];
305 } else if (n == 2) {
306 // n = 2 CASE BY CRAMERS RULE.
307 det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
308 if (det == 0.) {
309 ifail = -1;
310 return;
311 }
312 s = 1. / det;
313 c11 = s * a[1][1];
314 a[0][1] = -s * a[0][1];
315 a[1][0] = -s * a[1][0];
316 a[1][1] = s * a[0][0];
317 a[0][0] = c11;
318
319 b1 = b[0];
320 b[0] = c11 * b1 + a[0][1] * b[1];
321 b[1] = a[1][0] * b1 + a[1][1] * b[1];
322 } else {
323 // n = 1 CASE.
324 if (a[0][0] == 0.) {
325 ifail = -1;
326 return;
327 }
328 a[0][0] = 1. / a[0][0];
329 b[0] = a[0][0] * b[0];
330 }
331}
332
333void Cfact(const int n, std::vector<std::vector<std::complex<double> > >& a,
334 std::vector<int>& ir, int& ifail, std::complex<double>& det,
335 int& jfail) {
336
337 const double g1 = 1.e-19;
338 const double g2 = 1.e-19;
339
340 std::complex<double> tf;
341 double p, q, t;
342 std::complex<double> s11, s12;
343 int k;
344
345 ifail = jfail = 0;
346
347 int nxch = 0;
348 det = std::complex<double>(1., 0.);
349
350 for (int j = 1; j <= n; ++j) {
351 k = j;
352 p = std::max(fabs(real(a[j - 1][j - 1])), fabs(imag(a[j - 1][j - 1])));
353 if (j == n) {
354 if (p <= 0.) {
355 det = std::complex<double>(0., 0.);
356 ifail = -1;
357 jfail = 0;
358 return;
359 }
360 det *= a[j - 1][j - 1];
361 a[j - 1][j - 1] = std::complex<double>(1., 0.) / a[j - 1][j - 1];
362 t = std::max(fabs(real(det)), fabs(imag(det)));
363 if (t < g1) {
364 det = std::complex<double>(0., 0.);
365 if (jfail == 0) jfail = -1;
366 } else if (t > g2) {
367 det = std::complex<double>(1., 0.);
368 if (jfail == 0) jfail = +1;
369 }
370 continue;
371 }
372 for (int i = j + 1; i <= n; ++i) {
373 q = std::max(fabs(real(a[i - 1][j - 1])), fabs(imag(a[i - 1][j - 1])));
374 if (q <= p) continue;
375 k = i;
376 p = q;
377 }
378 if (k != j) {
379 for (int l = 1; l <= n; ++l) {
380 tf = a[j - 1][l - 1];
381 a[j - 1][l - 1] = a[k - 1][l - 1];
382 a[k - 1][l - 1] = tf;
383 }
384 ++nxch;
385 ir[nxch - 1] = j * 4096 + k;
386 } else if (p <= 0.) {
387 det = std::complex<double>(0., 0.);
388 ifail = -1;
389 jfail = 0;
390 return;
391 }
392 det *= a[j - 1][j - 1];
393 a[j - 1][j - 1] = 1. / a[j - 1][j - 1];
394 t = std::max(fabs(real(det)), fabs(imag(det)));
395 if (t < g1) {
396 det = std::complex<double>(0., 0.);
397 if (jfail == 0) jfail = -1;
398 } else if (t > g2) {
399 det = std::complex<double>(1., 0.);
400 if (jfail == 0) jfail = +1;
401 }
402 for (k = j + 1; k <= n; ++k) {
403 s11 = -a[j - 1][k - 1];
404 s12 = -a[k - 1][j];
405 if (j == 1) {
406 a[j - 1][k - 1] = -s11 * a[j - 1][j - 1];
407 a[k - 1][j] = -(s12 + a[j - 1][j] * a[k - 1][j - 1]);
408 continue;
409 }
410 for (int i = 1; i <= j - 1; ++i) {
411 s11 += a[i - 1][k - 1] * a[j - 1][i - 1];
412 s12 += a[i - 1][j] * a[k - 1][i - 1];
413 }
414 a[j - 1][k - 1] = -s11 * a[j - 1][j - 1];
415 a[k - 1][j] = -a[j - 1][j] * a[k - 1][j - 1] - s12;
416 }
417 }
418
419 if (nxch % 2 != 0) det = -det;
420 if (jfail != 0) det = std::complex<double>(0., 0.);
421 ir[n - 1] = nxch;
422}
423
424void Cfinv(const int n, std::vector<std::vector<std::complex<double> > >& a,
425 std::vector<int>& ir) {
426
427 std::complex<double> ti;
428 std::complex<double> s31, s32, s33, s34;
429
430 if (n <= 1) return;
431 a[1][0] = -a[1][1] * a[0][0] * a[1][0];
432 a[0][1] = -a[0][1];
433 if (n > 2) {
434 for (int i = 3; i <= n; ++i) {
435 for (int j = 1; j <= i - 2; ++j) {
436 s31 = std::complex<double>(0., 0.);
437 s32 = a[j - 1][i - 1];
438 for (int k = j; k <= i - 2; ++k) {
439 s31 += a[k - 1][j - 1] * a[i - 1][k - 1];
440 s32 += a[j - 1][k] * a[k][i - 1];
441 }
442 a[i - 1][j - 1] =
443 -a[i - 1][i - 1] * (s31 + a[i - 2][j - 1] * a[i - 1][i - 2]);
444 a[j - 1][i - 1] = -s32;
445 }
446 a[i - 1][i - 2] = -a[i - 1][i - 1] * a[i - 2][i - 2] * a[i - 1][i - 2];
447 a[i - 2][i - 1] = -a[i - 2][i - 1];
448 }
449 }
450
451 for (int i = 1; i <= n - 1; ++i) {
452 for (int j = 1; j <= i; ++j) {
453 s33 = a[i - 1][j - 1];
454 for (int k = 1; k <= n - i; ++k) {
455 s33 += a[i + k - 1][j - 1] * a[i - 1][i + k - 1];
456 }
457 a[i - 1][j - 1] = s33;
458 }
459 for (int j = 1; j <= n - i; ++j) {
460 s34 = std::complex<double>(0., 0.);
461 for (int k = j; k <= n - i; ++k) {
462 s34 += a[i + k - 1][i + j - 1] * a[i - 1][i + k - 1];
463 }
464 a[i - 1][i + j - 1] = s34;
465 }
466 }
467
468 int nxch = ir[n - 1];
469 if (nxch == 0) return;
470
471 for (int m = 1; m <= nxch; ++m) {
472 int k = nxch - m + 1;
473 int ij = ir[k - 1];
474 int i = ij / 4096;
475 int j = ij % 4096;
476 for (k = 1; k <= n; ++k) {
477 ti = a[k - 1][i - 1];
478 a[k - 1][i - 1] = a[k - 1][j - 1];
479 a[k - 1][j - 1] = ti;
480 }
481 }
482}
483
484// ******************************************************************
485//
486// REPLACES A BY ITS INVERSE.
487//
488// (PARAMETERS AS FOR CEQINV.)
489//
490// CALLS ... CFACT, CFINV.
491//
492// ******************************************************************
493
494void Cinv(const int n, std::vector<std::vector<std::complex<double> > >& a,
495 int& ifail) {
496
497 double t1, t2, t3;
498 std::complex<double> det, temp, s;
499 std::complex<double> c11, c12, c13, c21, c22, c23, c31, c32, c33;
500
501 std::vector<int> ir;
502 ir.clear();
503 ir.resize(n);
504 for (int i = 0; i < n; ++i) ir[i] = 0;
505
506 // TEST FOR PARAMETER ERRORS.
507 if (n < 1) {
508 ifail = 1;
509 return;
510 }
511
512 ifail = 0;
513 int jfail = 0;
514
515 if (n > 3) {
516 // n > 3 CASES.
517 // FACTORIZE MATRIX AND INVERT.
518 Cfact(n, a, ir, ifail, det, jfail);
519 if (ifail != 0) return;
520 Cfinv(n, a, ir);
521 } else if (n == 3) {
522 // n = 3 CASE.
523 // COMPUTE COFACTORS.
524 c11 = a[1][1] * a[2][2] - a[1][2] * a[2][1];
525 c12 = a[1][2] * a[2][0] - a[1][0] * a[2][2];
526 c13 = a[1][0] * a[2][1] - a[1][1] * a[2][0];
527 c21 = a[2][1] * a[0][2] - a[2][2] * a[0][1];
528 c22 = a[2][2] * a[0][0] - a[2][0] * a[0][2];
529 c23 = a[2][0] * a[0][1] - a[2][1] * a[0][0];
530 c31 = a[0][1] * a[1][2] - a[0][2] * a[1][1];
531 c32 = a[0][2] * a[1][0] - a[0][0] * a[1][2];
532 c33 = a[0][0] * a[1][1] - a[0][1] * a[1][0];
533 t1 = fabs(real(a[0][0])) + fabs(imag(a[0][0]));
534 t2 = fabs(real(a[1][0])) + fabs(imag(a[1][0]));
535 t3 = fabs(real(a[2][0])) + fabs(imag(a[2][0]));
536
537 // SET temp = PIVOT AND det = PIVOT * det.
538 if (t1 >= t2) {
539 if (t3 >= t1) {
540 // PIVOT IS A31
541 temp = a[2][0];
542 det = c23 * c12 - c22 * c13;
543 } else {
544 // PIVOT IS A11
545 temp = a[0][0];
546 det = c22 * c33 - c23 * c32;
547 }
548 } else {
549 if (t3 >= t2) {
550 // PIVOT IS A31
551 temp = a[2][0];
552 det = c23 * c12 - c22 * c13;
553 } else {
554 // PIVOT IS A21
555 temp = a[1][0];
556 det = c13 * c32 - c12 * c33;
557 }
558 }
559 // SET ELEMENTS OF INVERSE IN A.
560 if (real(det) == 0. && imag(det) == 0.) {
561 ifail = -1;
562 return;
563 }
564 s = temp / det;
565 a[0][0] = s * c11;
566 a[0][1] = s * c21;
567 a[0][2] = s * c31;
568 a[1][0] = s * c12;
569 a[1][1] = s * c22;
570 a[1][2] = s * c32;
571 a[2][0] = s * c13;
572 a[2][1] = s * c23;
573 a[2][2] = s * c33;
574 } else if (n == 2) {
575 // n=2 CASE BY CRAMERS RULE.
576 det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
577 if (real(det) == 0. && imag(det) == 0.) {
578 ifail = -1;
579 return;
580 }
581 s = std::complex<double>(1., 0.) / det;
582 c11 = s * a[1][1];
583 a[0][1] = -s * a[0][1];
584 a[1][0] = -s * a[1][0];
585 a[1][1] = s * a[0][0];
586 a[0][0] = c11;
587 } else {
588 // n = 1 CASE.
589 if (real(a[0][0]) == 0. && imag(a[0][0]) == 0.) {
590 ifail = -1;
591 return;
592 }
593 a[0][0] = std::complex<double>(1., 0.) / a[0][0];
594 }
595}
596
597// Numerical integration using 15-point Gauss-Kronrod algorithm
598// Origin: QUADPACK
599double GaussKronrod15(double (*f)(const double), const double a,
600 const double b) {
601
602 // Abscissae of the 15-point Kronrod rule
603 // xGK[1], xGK[3], ... abscissae of the 7-point Gauss rule
604 // xGK[0], xGK[2], ... abscissae which are optimally added
605 // to the 7-point Gauss rule
606 const double xGK[8] = {9.914553711208126e-01, 9.491079123427585e-01,
607 8.648644233597691e-01, 7.415311855993944e-01,
608 5.860872354676911e-01, 4.058451513773972e-01,
609 2.077849550078985e-01, 0.0e+00};
610 // Weights of the 15-point Kronrod rule
611 const double wGK[8] = {2.293532201052922e-02, 6.309209262997855e-02,
612 1.047900103222502e-01, 1.406532597155259e-01,
613 1.690047266392679e-01, 1.903505780647854e-01,
614 2.044329400752989e-01, 2.094821410847278e-01};
615 // Weights of the 7-point Gauss rule
616 const double wG[4] = {1.294849661688697e-01, 2.797053914892767e-01,
617 3.818300505051189e-01, 4.179591836734694e-01};
618
619 // Mid-point of the interval
620 const double center = 0.5 * (a + b);
621 // Half-length of the interval
622 const double halfLength = 0.5 * (b - a);
623
624 double fC = f(center);
625 // Result of the 7-point Gauss formula
626 double resG = fC * wG[3];
627 // Result of the 15-point Kronrod formula
628 double resK = fC * wGK[7];
629
630 for (int j = 0; j < 3; ++j) {
631 const int i = j * 2 + 1;
632 // Abscissa
633 const double x = halfLength * xGK[i];
634 // Function value
635 const double fSum = f(center - x) + f(center + x);
636 resG += wG[j] * fSum;
637 resK += wGK[i] * fSum;
638 }
639
640 for (int j = 0; j < 4; ++j) {
641 const int i = j * 2;
642 const double x = halfLength * xGK[i];
643 const double fSum = f(center - x) + f(center + x);
644 resK += wGK[i] * fSum;
645 }
646
647 return resK * halfLength;
648}
649
650double Divdif(const std::vector<double>& f, const std::vector<double>& a,
651 int nn, double x, int mm) {
652
653 // C++ version of DIVDIF (CERN program library E105) which performs
654 // tabular interpolation using symmetrically placed argument points.
655
656 double t[20], d[20];
657
658 const int mmax = 10;
659
660 // Check the arguments.
661 if (nn < 2) {
662 std::cerr << "Divdif:\n";
663 std::cerr << " Array length < 2.\n";
664 return 0.;
665 }
666 if (mm < 1) {
667 std::cerr << "Divdif:\n";
668 std::cerr << " Interpolation order < 1.\n";
669 return 0.;
670 }
671
672 // Deal with the case that X is located at A(1) or A(N).
673 if (fabs(x - a[0]) < 1.e-6 * (fabs(a[0]) + fabs(a[nn - 1]))) {
674 return f[0];
675 }
676 if (fabs(x - a[nn - 1]) < 1.e-6 * (fabs(a[0]) + fabs(a[nn - 1]))) {
677 return f[nn - 1];
678 }
679
680 // Find subscript IX of X in array A.
681 int n = nn;
682 int m;
683 if (mm <= mmax && mm <= n - 1) {
684 m = mm;
685 } else {
686 if (mmax <= n - 1) {
687 m = mmax;
688 } else {
689 m = n - 1;
690 }
691 }
692 int mplus = m + 1;
693 int ix = 0;
694 int iy = n + 1;
695 int mid;
696 if (a[0] > a[n - 1]) {
697 // Search decreasing arguments.
698 do {
699 mid = (ix + iy) / 2;
700 if (x > a[mid - 1]) {
701 iy = mid;
702 } else {
703 ix = mid;
704 }
705 } while (iy - ix > 1);
706 } else {
707 // Search increasing arguments.
708 do {
709 mid = (ix + iy) / 2;
710 if (x < a[mid - 1]) {
711 iy = mid;
712 } else {
713 ix = mid;
714 }
715 } while (iy - ix > 1);
716 }
717 // Copy reordered interpolation points into (T[I],D[I]), setting
718 // EXTRA to True if M+2 points to be used.
719 int npts = m + 2 - (m % 2);
720 int ip = 0;
721 int l = 0;
722 int isub;
723 do {
724 isub = ix + l;
725 if ((1 > isub) || (isub > n)) {
726 // Skip point.
727 npts = mplus;
728 } else {
729 // Insert point.
730 ip++;
731 t[ip - 1] = a[isub - 1];
732 d[ip - 1] = f[isub - 1];
733 }
734 if (ip < npts) {
735 l = -l;
736 if (l >= 0) {
737 l++;
738 }
739 }
740 } while (ip < npts);
741
742 bool extra = npts != mplus;
743 // Replace d by the leading diagonal of a divided-difference table,
744 // supplemented by an extra line if EXTRA is True.
745 for (int l = 1; l <= m; l++) {
746 if (extra) {
747 isub = mplus - l;
748 d[m + 1] = (d[m + 1] - d[m - 1]) / (t[m + 1] - t[isub - 1]);
749 }
750 int i = mplus;
751 for (int j = l; j <= m; j++) {
752 isub = i - l;
753 d[i - 1] = (d[i - 1] - d[i - 1 - 1]) / (t[i - 1] - t[isub - 1]);
754 i--;
755 }
756 }
757 // Evaluate the Newton interpolation formula at X, averaging two values
758 // of last difference if EXTRA is True.
759 double sum = d[mplus - 1];
760 if (extra) {
761 sum = 0.5 * (sum + d[m + 1]);
762 }
763 int j = m;
764 for (int l = 1; l <= m; l++) {
765 sum = d[j - 1] + (x - t[j - 1]) * sum;
766 j--;
767 }
768 return sum;
769}
770
771bool Boxin3(std::vector<std::vector<std::vector<double> > >& value,
772 std::vector<double>& xAxis, std::vector<double>& yAxis,
773 std::vector<double>& zAxis, int nx, int ny, int nz, double xx,
774 double yy, double zz, double& f, int iOrder) {
775
776 // std::cout << nx << ", " << ny << ", " << nz << "\n";
777 //-----------------------------------------------------------------------
778 // BOXIN3 - interpolation of order 1 and 2 in an irregular rectangular
779 // 3-dimensional grid.
780 // (Last changed on 13/ 2/00.)
781 //-----------------------------------------------------------------------
782
783 int iX0 = 0, iX1 = 0;
784 int iY0 = 0, iY1 = 0;
785 int iZ0 = 0, iZ1 = 0;
786 double fX[4], fY[4], fZ[4];
787
788 // Ensure we are in the grid.
789 const double x = std::min(std::max(xx, std::min(xAxis[0], xAxis[nx - 1])),
790 std::max(xAxis[0], xAxis[nx - 1]));
791 const double y = std::min(std::max(yy, std::min(yAxis[0], yAxis[ny - 1])),
792 std::max(yAxis[0], yAxis[ny - 1]));
793 const double z = std::min(std::max(zz, std::min(zAxis[0], zAxis[nz - 1])),
794 std::max(zAxis[0], zAxis[nz - 1]));
795
796 // Make sure we have enough points.
797 if (iOrder < 0 || iOrder > 2 || nx < 1 || ny < 1 || nz < 1) {
798 std::cerr << "Boxin3:\n";
799 std::cerr << " Incorrect order or number of points.\n";
800 std::cerr << " No interpolation.\n";
801 f = 0.;
802 return false;
803 }
804
805 if (iOrder == 0 || nx == 1) {
806 // Zeroth order interpolation in x.
807 // Find the nearest node.
808 double dist = fabs(x - xAxis[0]);
809 int iNode = 0;
810 for (int i = 1; i < nx; i++) {
811 if (fabs(x - xAxis[i]) < dist) {
812 dist = fabs(x - xAxis[i]);
813 iNode = i;
814 }
815 }
816 // Set the summing range.
817 iX0 = iNode;
818 iX1 = iNode;
819 // Establish the shape functions.
820 fX[0] = 0.;
821 fX[1] = 0.;
822 fX[2] = 0.;
823 fX[3] = 0.;
824 } else if (iOrder == 1 || nx == 2) {
825 // First order interpolation in x.
826 // Find the grid segment containing this point.
827 int iGrid = 0;
828 for (int i = 1; i < nx; i++) {
829 if ((xAxis[i - 1] - x) * (x - xAxis[i]) >= 0.) {
830 iGrid = i;
831 }
832 }
833 // Ensure there won't be divisions by zero.
834 if (xAxis[iGrid] == xAxis[iGrid - 1]) {
835 std::cerr << "Boxin3:\n";
836 std::cerr << " Incorrect grid; no interpolation.\n";
837 f = 0.;
838 return false;
839 }
840 // Compute local coordinates.
841 const double xLocal =
842 (x - xAxis[iGrid - 1]) / (xAxis[iGrid] - xAxis[iGrid - 1]);
843 // Set the summing range.
844 iX0 = iGrid - 1;
845 iX1 = iGrid;
846 // Set the shape functions.
847 fX[0] = 1. - xLocal;
848 fX[1] = xLocal;
849 fX[2] = 0.;
850 fX[3] = 0.;
851 } else if (iOrder == 2) {
852 // Second order interpolation in x.
853 // Find the grid segment containing this point.
854 int iGrid = 0;
855 for (int i = 1; i < nx; i++) {
856 if ((xAxis[i - 1] - x) * (x - xAxis[i]) >= 0.) {
857 iGrid = i;
858 }
859 }
860 // Compute the local coordinate for this grid segment.
861 const double xLocal =
862 (x - xAxis[iGrid - 1]) / (xAxis[iGrid] - xAxis[iGrid - 1]);
863 // Set the summing range and shape functions.
864 if (iGrid == 1) {
865 iX0 = iGrid - 1;
866 iX1 = iGrid + 1;
867 if (xAxis[iX0] == xAxis[iX0 + 1] || xAxis[iX0] == xAxis[iX0 + 2] ||
868 xAxis[iX0 + 1] == xAxis[iX0 + 2]) {
869 std::cerr << "Boxin3:\n";
870 std::cerr << " One or more grid points in x coincide.\n";
871 std::cerr << " No interpolation.\n";
872 f = 0.;
873 return false;
874 }
875 fX[0] = (x - xAxis[iX0 + 1]) * (x - xAxis[iX0 + 2]) /
876 ((xAxis[iX0] - xAxis[iX0 + 1]) * (xAxis[iX0] - xAxis[iX0 + 2]));
877 fX[1] =
878 (x - xAxis[iX0]) * (x - xAxis[iX0 + 2]) /
879 ((xAxis[iX0 + 1] - xAxis[iX0]) * (xAxis[iX0 + 1] - xAxis[iX0 + 2]));
880 fX[2] =
881 (x - xAxis[iX0]) * (x - xAxis[iX0 + 1]) /
882 ((xAxis[iX0 + 2] - xAxis[iX0]) * (xAxis[iX0 + 2] - xAxis[iX0 + 1]));
883 } else if (iGrid == nx - 1) {
884 iX0 = iGrid - 2;
885 iX1 = iGrid;
886 if (xAxis[iX0] == xAxis[iX0 + 1] || xAxis[iX0] == xAxis[iX0 + 2] ||
887 xAxis[iX0 + 1] == xAxis[iX0 + 2]) {
888 std::cerr << "Boxin3:\n";
889 std::cerr << " One or more grid points in x coincide.\n";
890 std::cerr << " No interpolation.\n";
891 f = 0.;
892 return false;
893 }
894 fX[0] = (x - xAxis[iX0 + 1]) * (x - xAxis[iX0 + 2]) /
895 ((xAxis[iX0] - xAxis[iX0 + 1]) * (xAxis[iX0] - xAxis[iX0 + 2]));
896 fX[1] =
897 (x - xAxis[iX0]) * (x - xAxis[iX0 + 2]) /
898 ((xAxis[iX0 + 1] - xAxis[iX0]) * (xAxis[iX0 + 1] - xAxis[iX0 + 2]));
899 fX[2] =
900 (x - xAxis[iX0]) * (x - xAxis[iX0 + 1]) /
901 ((xAxis[iX0 + 2] - xAxis[iX0]) * (xAxis[iX0 + 2] - xAxis[iX0 + 1]));
902 } else {
903 iX0 = iGrid - 2;
904 iX1 = iGrid + 1;
905 if (xAxis[iX0] == xAxis[iX0 + 1] || xAxis[iX0] == xAxis[iX0 + 2] ||
906 xAxis[iX0] == xAxis[iX0 + 3] || xAxis[iX0 + 1] == xAxis[iX0 + 2] ||
907 xAxis[iX0 + 1] == xAxis[iX0 + 3] ||
908 xAxis[iX0 + 2] == xAxis[iX0 + 3]) {
909 std::cerr << "Boxin3:\n";
910 std::cerr << " One or more grid points in x coincide.\n";
911 std::cerr << " No interpolation.\n";
912 f = 0.;
913 return false;
914 }
915 fX[0] = (x - xAxis[iX0 + 1]) * (x - xAxis[iX0 + 2]) /
916 ((xAxis[iX0] - xAxis[iX0 + 1]) * (xAxis[iX0] - xAxis[iX0 + 2]));
917 fX[1] =
918 (x - xAxis[iX0]) * (x - xAxis[iX0 + 2]) /
919 ((xAxis[iX0 + 1] - xAxis[iX0]) * (xAxis[iX0 + 1] - xAxis[iX0 + 2]));
920 fX[2] =
921 (x - xAxis[iX0]) * (x - xAxis[iX0 + 1]) /
922 ((xAxis[iX0 + 2] - xAxis[iX0]) * (xAxis[iX0 + 2] - xAxis[iX0 + 1]));
923 fX[0] *= (1. - xLocal);
924 fX[1] = fX[1] * (1. - xLocal) + xLocal * (x - xAxis[iX0 + 2]) *
925 (x - xAxis[iX0 + 3]) /
926 ((xAxis[iX0 + 1] - xAxis[iX0 + 2]) *
927 (xAxis[iX0 + 1] - xAxis[iX0 + 3]));
928 fX[2] = fX[2] * (1. - xLocal) + xLocal * (x - xAxis[iX0 + 1]) *
929 (x - xAxis[iX0 + 3]) /
930 ((xAxis[iX0 + 2] - xAxis[iX0 + 1]) *
931 (xAxis[iX0 + 2] - xAxis[iX0 + 3]));
932 fX[3] = xLocal * (x - xAxis[iX0 + 1]) * (x - xAxis[iX0 + 2]) /
933 ((xAxis[iX0 + 3] - xAxis[iX0 + 1]) *
934 (xAxis[iX0 + 3] - xAxis[iX0 + 2]));
935 }
936 }
937
938 if (iOrder == 0 || ny == 1) {
939 // Zeroth order interpolation in y.
940 // Find the nearest node.
941 double dist = fabs(y - yAxis[0]);
942 int iNode = 0;
943 for (int i = 1; i < ny; i++) {
944 if (fabs(y - yAxis[i]) < dist) {
945 dist = fabs(y - yAxis[i]);
946 iNode = i;
947 }
948 }
949 // Set the summing range.
950 iY0 = iNode;
951 iY1 = iNode;
952 // Establish the shape functions.
953 fY[0] = 1.;
954 fY[1] = 0.;
955 fY[2] = 0.;
956 } else if (iOrder == 1 || ny == 2) {
957 // First order interpolation in y.
958 // Find the grid segment containing this point.
959 int iGrid = 0;
960 for (int i = 1; i < ny; i++) {
961 if ((yAxis[i - 1] - y) * (y - yAxis[i]) >= 0.) {
962 iGrid = i;
963 }
964 }
965 // Ensure there won't be divisions by zero.
966 if (yAxis[iGrid] == yAxis[iGrid - 1]) {
967 std::cerr << "Boxin3:\n";
968 std::cerr << " Incorrect grid; no interpolation.\n";
969 f = 0.;
970 return false;
971 }
972 // Compute local coordinates.
973 const double yLocal =
974 (y - yAxis[iGrid - 1]) / (yAxis[iGrid] - yAxis[iGrid - 1]);
975 // Set the summing range.
976 iY0 = iGrid - 1;
977 iY1 = iGrid;
978 // Set the shape functions.
979 fY[0] = 1. - yLocal;
980 fY[1] = yLocal;
981 fY[2] = 0.;
982 } else if (iOrder == 2) {
983 // Second order interpolation in y.
984 // Find the grid segment containing this point.
985 int iGrid = 0;
986 for (int i = 1; i < ny; i++) {
987 if ((yAxis[i - 1] - y) * (y - yAxis[i]) >= 0.) {
988 iGrid = i;
989 }
990 }
991 // Compute the local coordinate for this grid segment.
992 const double yLocal =
993 (y - yAxis[iGrid - 1]) / (yAxis[iGrid] - yAxis[iGrid - 1]);
994 // Set the summing range and shape functions.
995 // These assignments are shared by all of the following conditions,
996 // so it's easier to take them out.
997 fY[0] = (y - yAxis[iY0 + 1]) * (y - yAxis[iY0 + 2]) /
998 ((yAxis[iY0] - yAxis[iY0 + 1]) * (yAxis[iY0] - yAxis[iY0 + 2]));
999 fY[1] = (y - yAxis[iY0]) * (y - yAxis[iY0 + 2]) /
1000 ((yAxis[iY0 + 1] - yAxis[iY0]) * (yAxis[iY0 + 1] - yAxis[iY0 + 2]));
1001 fY[2] = (y - yAxis[iY0]) * (y - yAxis[iY0 + 1]) /
1002 ((yAxis[iY0 + 2] - yAxis[iY0]) * (yAxis[iY0 + 2] - yAxis[iY0 + 1]));
1003
1004 if (iGrid == 1) {
1005 iY0 = iGrid - 1;
1006 iY1 = iGrid + 1;
1007 if (yAxis[iY0] == yAxis[iY0 + 1] || yAxis[iY0] == yAxis[iY0 + 2] ||
1008 yAxis[iY0 + 1] == yAxis[iY0 + 2]) {
1009 std::cerr << "Boxin3:\n";
1010 std::cerr << " One or more grid points in y coincide.\n";
1011 std::cerr << " No interpolation.\n";
1012 f = 0.;
1013 return false;
1014 }
1015 fY[0] = (y - yAxis[iY0 + 1]) * (y - yAxis[iY0 + 2]) /
1016 ((yAxis[iY0] - yAxis[iY0 + 1]) * (yAxis[iY0] - yAxis[iY0 + 2]));
1017 fY[1] =
1018 (y - yAxis[iY0]) * (y - yAxis[iY0 + 2]) /
1019 ((yAxis[iY0 + 1] - yAxis[iY0]) * (yAxis[iY0 + 1] - yAxis[iY0 + 2]));
1020 fY[2] =
1021 (y - yAxis[iY0]) * (y - yAxis[iY0 + 1]) /
1022 ((yAxis[iY0 + 2] - yAxis[iY0]) * (yAxis[iY0 + 2] - yAxis[iY0 + 1]));
1023 } else if (iGrid == ny - 1) {
1024 iY0 = iGrid - 2;
1025 iY1 = iGrid;
1026 if (yAxis[iY0] == yAxis[iY0 + 1] || yAxis[iY0] == yAxis[iY0 + 2] ||
1027 yAxis[iY0 + 1] == yAxis[iY0 + 2]) {
1028 std::cerr << "Boxin3:\n";
1029 std::cerr << " One or more grid points in y coincide.\n";
1030 std::cerr << " No interpolation.\n";
1031 f = 0.;
1032 return false;
1033 }
1034 fY[0] = (y - yAxis[iY0 + 1]) * (y - yAxis[iY0 + 2]) /
1035 ((yAxis[iY0] - yAxis[iY0 + 1]) * (yAxis[iY0] - yAxis[iY0 + 2]));
1036 fY[1] =
1037 (y - yAxis[iY0]) * (y - yAxis[iY0 + 2]) /
1038 ((yAxis[iY0 + 1] - yAxis[iY0]) * (yAxis[iY0 + 1] - yAxis[iY0 + 2]));
1039 fY[2] =
1040 (y - yAxis[iY0]) * (y - yAxis[iY0 + 1]) /
1041 ((yAxis[iY0 + 2] - yAxis[iY0]) * (yAxis[iY0 + 2] - yAxis[iY0 + 1]));
1042 } else {
1043 iY0 = iGrid - 2;
1044 iY1 = iGrid + 1;
1045 if (yAxis[iY0] == yAxis[iY0 + 1] || yAxis[iY0] == yAxis[iY0 + 2] ||
1046 yAxis[iY0] == yAxis[iY0 + 3] || yAxis[iY0 + 1] == yAxis[iY0 + 2] ||
1047 yAxis[iY0 + 1] == yAxis[iY0 + 3] ||
1048 yAxis[iY0 + 2] == yAxis[iY0 + 3]) {
1049 std::cerr << "Boxin3:\n";
1050 std::cerr << " One or more grid points in y coincide.\n";
1051 std::cerr << " No interpolation.\n";
1052 f = 0.;
1053 return false;
1054 }
1055 fY[0] = (y - yAxis[iY0 + 1]) * (y - yAxis[iY0 + 2]) /
1056 ((yAxis[iY0] - yAxis[iY0 + 1]) * (yAxis[iY0] - yAxis[iY0 + 2]));
1057 fY[1] =
1058 (y - yAxis[iY0]) * (y - yAxis[iY0 + 2]) /
1059 ((yAxis[iY0 + 1] - yAxis[iY0]) * (yAxis[iY0 + 1] - yAxis[iY0 + 2]));
1060 fY[2] =
1061 (y - yAxis[iY0]) * (y - yAxis[iY0 + 1]) /
1062 ((yAxis[iY0 + 2] - yAxis[iY0]) * (yAxis[iY0 + 2] - yAxis[iY0 + 1]));
1063
1064 fY[0] *= (1. - yLocal);
1065 fY[1] = fY[1] * (1. - yLocal) + yLocal * (y - yAxis[iY0 + 2]) *
1066 (y - yAxis[iY0 + 3]) /
1067 ((yAxis[iY0 + 1] - yAxis[iY0 + 2]) *
1068 (yAxis[iY0 + 1] - yAxis[iY0 + 3]));
1069 fY[2] = fY[2] * (1. - yLocal) + yLocal * (y - yAxis[iY0 + 1]) *
1070 (y - yAxis[iY0 + 3]) /
1071 ((yAxis[iY0 + 2] - yAxis[iY0 + 1]) *
1072 (yAxis[iY0 + 2] - yAxis[iY0 + 3]));
1073 fY[3] = yLocal * (y - yAxis[iY0 + 1]) * (y - yAxis[iY0 + 2]) /
1074 ((yAxis[iY0 + 3] - yAxis[iY0 + 1]) *
1075 (yAxis[iY0 + 3] - yAxis[iY0 + 2]));
1076 }
1077 }
1078
1079 if (iOrder == 0 || nz == 1) {
1080 // Zeroth order interpolation in z.
1081 // Find the nearest node.
1082 double dist = fabs(z - zAxis[0]);
1083 int iNode = 0;
1084 for (int i = 1; i < nz; i++) {
1085 if (fabs(z - zAxis[i]) < dist) {
1086 dist = fabs(z - zAxis[i]);
1087 iNode = i;
1088 }
1089 }
1090 // Set the summing range.
1091 iZ0 = iNode;
1092 iZ1 = iNode;
1093 // Establish the shape functions.
1094 fZ[0] = 1.;
1095 fZ[1] = 0.;
1096 fZ[2] = 0.;
1097 } else if (iOrder == 1 || nz == 2) {
1098 // First order interpolation in z.
1099 // Find the grid segment containing this point.
1100 int iGrid = 0;
1101 for (int i = 1; i < nz; i++) {
1102 if ((zAxis[i - 1] - z) * (z - zAxis[i]) >= 0.) {
1103 iGrid = i;
1104 }
1105 }
1106 // Ensure there won't be divisions by zero.
1107 if (zAxis[iGrid] == zAxis[iGrid - 1]) {
1108 std::cerr << "Boxin3:\n";
1109 std::cerr << " Incorrect grid; no interpolation.\n";
1110 f = 0.;
1111 return false;
1112 }
1113 // Compute local coordinates.
1114 const double zLocal =
1115 (z - zAxis[iGrid - 1]) / (zAxis[iGrid] - zAxis[iGrid - 1]);
1116 // Set the summing range.
1117 iZ0 = iGrid - 1;
1118 iZ1 = iGrid;
1119 // Set the shape functions.
1120 fZ[0] = 1. - zLocal;
1121 fZ[1] = zLocal;
1122 fZ[2] = 0.;
1123 } else if (iOrder == 2) {
1124 // Second order interpolation in z.
1125 // Find the grid segment containing this point.
1126 int iGrid = 0;
1127 for (int i = 1; i < nz; i++) {
1128 if ((zAxis[i - 1] - z) * (z - zAxis[i]) >= 0.) {
1129 iGrid = i;
1130 }
1131 }
1132 // Compute the local coordinate for this grid segment.
1133 const double zLocal =
1134 (z - zAxis[iGrid - 1]) / (zAxis[iGrid] - zAxis[iGrid - 1]);
1135 // Set the summing range and shape functions.
1136 // These assignments are shared by all of the following conditions,
1137 // so it's easier to take them out.
1138 fZ[0] = (z - zAxis[iZ0 + 1]) * (z - zAxis[iZ0 + 2]) /
1139 ((zAxis[iZ0] - zAxis[iZ0 + 1]) * (zAxis[iZ0] - zAxis[iZ0 + 2]));
1140 fZ[1] = (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 2]) /
1141 ((zAxis[iZ0 + 1] - zAxis[iZ0]) * (zAxis[iZ0 + 1] - zAxis[iZ0 + 2]));
1142 fZ[2] = (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 1]) /
1143 ((zAxis[iZ0 + 2] - zAxis[iZ0]) * (zAxis[iZ0 + 2] - zAxis[iZ0 + 1]));
1144
1145 if (iGrid == 1) {
1146 iZ0 = iGrid - 1;
1147 iZ1 = iGrid + 1;
1148 if (zAxis[iZ0] == zAxis[iZ0 + 1] || zAxis[iZ0] == zAxis[iZ0 + 2] ||
1149 zAxis[iZ0 + 1] == zAxis[iZ0 + 2]) {
1150 std::cerr << "Boxin3:\n";
1151 std::cerr << " One or more grid points in z coincide.\n";
1152 std::cerr << " No interpolation.\n";
1153 f = 0.;
1154 return false;
1155 }
1156 fZ[0] = (z - zAxis[iZ0 + 1]) * (z - zAxis[iZ0 + 2]) /
1157 ((zAxis[iZ0] - zAxis[iZ0 + 1]) * (zAxis[iZ0] - zAxis[iZ0 + 2]));
1158 fZ[1] =
1159 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 2]) /
1160 ((zAxis[iZ0 + 1] - zAxis[iZ0]) * (zAxis[iZ0 + 1] - zAxis[iZ0 + 2]));
1161 fZ[2] =
1162 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 1]) /
1163 ((zAxis[iZ0 + 2] - zAxis[iZ0]) * (zAxis[iZ0 + 2] - zAxis[iZ0 + 1]));
1164 } else if (iGrid == nz - 1) {
1165 iZ0 = iGrid - 2;
1166 iZ1 = iGrid;
1167 if (zAxis[iZ0] == zAxis[iZ0 + 1] || zAxis[iZ0] == zAxis[iZ0 + 2] ||
1168 zAxis[iZ0 + 1] == zAxis[iZ0 + 2]) {
1169 std::cerr << "Boxin3:\n";
1170 std::cerr << " One or more grid points in z coincide.\n";
1171 std::cerr << " No interpolation.\n";
1172 f = 0.;
1173 return false;
1174 }
1175 fZ[0] = (z - zAxis[iZ0 + 1]) * (z - zAxis[iZ0 + 2]) /
1176 ((zAxis[iZ0] - zAxis[iZ0 + 1]) * (zAxis[iZ0] - zAxis[iZ0 + 2]));
1177 fZ[1] =
1178 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 2]) /
1179 ((zAxis[iZ0 + 1] - zAxis[iZ0]) * (zAxis[iZ0 + 1] - zAxis[iZ0 + 2]));
1180 fZ[2] =
1181 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 1]) /
1182 ((zAxis[iZ0 + 2] - zAxis[iZ0]) * (zAxis[iZ0 + 2] - zAxis[iZ0 + 1]));
1183 } else {
1184 iZ0 = iGrid - 2;
1185 iZ1 = iGrid + 1;
1186
1187 if (zAxis[iZ0] == zAxis[iZ0 + 1] || zAxis[iZ0] == zAxis[iZ0 + 2] ||
1188 zAxis[iZ0] == zAxis[iZ0 + 3] || zAxis[iZ0 + 1] == zAxis[iZ0 + 2] ||
1189 zAxis[iZ0 + 1] == zAxis[iZ0 + 3] ||
1190 zAxis[iZ0 + 2] == zAxis[iZ0 + 3]) {
1191 std::cerr << "Boxin3:\n";
1192 std::cerr << " One or more grid points in z coincide.\n";
1193 std::cerr << " No interpolation.\n";
1194 f = 0.;
1195 return false;
1196 }
1197
1198 fZ[0] = (z - zAxis[iZ0 + 1]) * (z - zAxis[iZ0 + 2]) /
1199 ((zAxis[iZ0] - zAxis[iZ0 + 1]) * (zAxis[iZ0] - zAxis[iZ0 + 2]));
1200 fZ[1] =
1201 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 2]) /
1202 ((zAxis[iZ0 + 1] - zAxis[iZ0]) * (zAxis[iZ0 + 1] - zAxis[iZ0 + 2]));
1203 fZ[2] =
1204 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 1]) /
1205 ((zAxis[iZ0 + 2] - zAxis[iZ0]) * (zAxis[iZ0 + 2] - zAxis[iZ0 + 1]));
1206
1207 fZ[0] *= (1. - zLocal);
1208 fZ[1] = fZ[1] * (1. - zLocal) + zLocal * (z - zAxis[iZ0 + 2]) *
1209 (z - zAxis[iZ0 + 3]) /
1210 ((zAxis[iZ0 + 1] - zAxis[iZ0 + 2]) *
1211 (zAxis[iZ0 + 1] - zAxis[iZ0 + 3]));
1212 fZ[2] = fZ[2] * (1. - zLocal) + zLocal * (z - zAxis[iZ0 + 1]) *
1213 (z - zAxis[iZ0 + 3]) /
1214 ((zAxis[iZ0 + 2] - zAxis[iZ0 + 1]) *
1215 (zAxis[iZ0 + 2] - zAxis[iZ0 + 3]));
1216 fZ[3] = zLocal * (z - zAxis[iZ0 + 1]) * (z - zAxis[iZ0 + 2]) /
1217 ((zAxis[iZ0 + 3] - zAxis[iZ0 + 1]) *
1218 (zAxis[iZ0 + 3] - zAxis[iZ0 + 2]));
1219 }
1220 }
1221
1222 f = 0.;
1223 for (int i = iX0; i <= iX1; ++i) {
1224 for (int j = iY0; j <= iY1; ++j) {
1225 for (int k = iZ0; k <= iZ1; ++k) {
1226 // std::cout << "i = " << i << ", j = " << j << ", k = " << k << "\n";
1227 // std::cout << "value: " << value[i][j][k] << "\n";
1228 // std::cout << "fX = " << fX[i - iX0] << ", fY = " << fY[j - iY0] << ",
1229 // fZ = " << fZ[k - iZ0] << "\n";
1230 f += value[i][j][k] * fX[i - iX0] * fY[j - iY0] * fZ[k - iZ0];
1231 }
1232 }
1233 }
1234 // std::cout << f << std::endl;
1235 return true;
1236}
1237}
1238}
fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:616
Garfield::Numerics::Dfinv
void Dfinv(const int n, std::vector< std::vector< double > > &a, std::vector< int > &ir)
Definition: Numerics.cc:138
Garfield::Numerics::Dfact
void Dfact(const int n, std::vector< std::vector< double > > &a, std::vector< int > &ir, int &ifail, double &det, int &jfail)
Definition: Numerics.cc:10
Garfield::Numerics::Cfact
void Cfact(const int n, std::vector< std::vector< std::complex< double > > > &a, std::vector< int > &ir, int &ifail, std::complex< double > &det, int &jfail)
Definition: Numerics.cc:333
Garfield::Numerics::Boxin3
bool Boxin3(std::vector< std::vector< std::vector< double > > > &value, std::vector< double > &xAxis, std::vector< double > &yAxis, std::vector< double > &zAxis, int nx, int ny, int nz, double xx, double yy, double zz, double &f, int iOrder)
Definition: Numerics.cc:771
Garfield::Numerics::GaussKronrod15
double GaussKronrod15(double(*f)(const double), const double a, const double b)
Definition: Numerics.cc:599
Garfield::Numerics::Cinv
void Cinv(const int n, std::vector< std::vector< std::complex< double > > > &a, int &ifail)
Definition: Numerics.cc:494
Garfield::Numerics::Divdif
double Divdif(const std::vector< double > &f, const std::vector< double > &a, int nn, double x, int mm)
Definition: Numerics.cc:650
Garfield::Numerics::Deqinv
void Deqinv(const int n, std::vector< std::vector< double > > &a, int &ifail, std::vector< double > &b)
Definition: Numerics.cc:216
Garfield::Numerics::Cfinv
void Cfinv(const int n, std::vector< std::vector< std::complex< double > > > &a, std::vector< int > &ir)
Definition: Numerics.cc:424
Garfield::Numerics::Dfeqn
void Dfeqn(const int n, std::vector< std::vector< double > > &a, std::vector< int > &ir, std::vector< double > &b)
Definition: Numerics.cc:98