Garfield++ v2r0
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 do {
723 const int isub = ix + l;
724 if ((1 > isub) || (isub > n)) {
725 // Skip point.
726 npts = mplus;
727 } else {
728 // Insert point.
729 ip++;
730 t[ip - 1] = a[isub - 1];
731 d[ip - 1] = f[isub - 1];
732 }
733 if (ip < npts) {
734 l = -l;
735 if (l >= 0) {
736 l++;
737 }
738 }
739 } while (ip < npts);
740
741 bool extra = npts != mplus;
742 // Replace d by the leading diagonal of a divided-difference table,
743 // supplemented by an extra line if EXTRA is True.
744 for (l = 1; l <= m; l++) {
745 if (extra) {
746 const int isub = mplus - l;
747 d[m + 1] = (d[m + 1] - d[m - 1]) / (t[m + 1] - t[isub - 1]);
748 }
749 int i = mplus;
750 for (int j = l; j <= m; j++) {
751 const int isub = i - l;
752 d[i - 1] = (d[i - 1] - d[i - 1 - 1]) / (t[i - 1] - t[isub - 1]);
753 i--;
754 }
755 }
756 // Evaluate the Newton interpolation formula at X, averaging two values
757 // of last difference if EXTRA is True.
758 double sum = d[mplus - 1];
759 if (extra) {
760 sum = 0.5 * (sum + d[m + 1]);
761 }
762 int j = m;
763 for (l = 1; l <= m; l++) {
764 sum = d[j - 1] + (x - t[j - 1]) * sum;
765 j--;
766 }
767 return sum;
768}
769
770bool Boxin3(const std::vector<std::vector<std::vector<double> > >& value,
771 const std::vector<double>& xAxis,
772 const std::vector<double>& yAxis,
773 const std::vector<double>& zAxis,
774 const int nx, const int ny, const int nz,
775 const double xx, const double yy, const double zz,
776 double& f, const int iOrder) {
777
778 // std::cout << nx << ", " << ny << ", " << nz << "\n";
779 //-----------------------------------------------------------------------
780 // BOXIN3 - interpolation of order 1 and 2 in an irregular rectangular
781 // 3-dimensional grid.
782 // (Last changed on 13/ 2/00.)
783 //-----------------------------------------------------------------------
784
785 int iX0 = 0, iX1 = 0;
786 int iY0 = 0, iY1 = 0;
787 int iZ0 = 0, iZ1 = 0;
788 double fX[4], fY[4], fZ[4];
789
790 // Ensure we are in the grid.
791 const double x = std::min(std::max(xx, std::min(xAxis[0], xAxis[nx - 1])),
792 std::max(xAxis[0], xAxis[nx - 1]));
793 const double y = std::min(std::max(yy, std::min(yAxis[0], yAxis[ny - 1])),
794 std::max(yAxis[0], yAxis[ny - 1]));
795 const double z = std::min(std::max(zz, std::min(zAxis[0], zAxis[nz - 1])),
796 std::max(zAxis[0], zAxis[nz - 1]));
797
798 // Make sure we have enough points.
799 if (iOrder < 0 || iOrder > 2 || nx < 1 || ny < 1 || nz < 1) {
800 std::cerr << "Boxin3:\n";
801 std::cerr << " Incorrect order or number of points.\n";
802 std::cerr << " No interpolation.\n";
803 f = 0.;
804 return false;
805 }
806
807 if (iOrder == 0 || nx == 1) {
808 // Zeroth order interpolation in x.
809 // Find the nearest node.
810 double dist = fabs(x - xAxis[0]);
811 int iNode = 0;
812 for (int i = 1; i < nx; i++) {
813 if (fabs(x - xAxis[i]) < dist) {
814 dist = fabs(x - xAxis[i]);
815 iNode = i;
816 }
817 }
818 // Set the summing range.
819 iX0 = iNode;
820 iX1 = iNode;
821 // Establish the shape functions.
822 fX[0] = 1.;
823 fX[1] = 0.;
824 fX[2] = 0.;
825 fX[3] = 0.;
826 } else if (iOrder == 1 || nx == 2) {
827 // First order interpolation in x.
828 // Find the grid segment containing this point.
829 int iGrid = 0;
830 for (int i = 1; i < nx; i++) {
831 if ((xAxis[i - 1] - x) * (x - xAxis[i]) >= 0.) {
832 iGrid = i;
833 }
834 }
835 // Ensure there won't be divisions by zero.
836 if (xAxis[iGrid] == xAxis[iGrid - 1]) {
837 std::cerr << "Boxin3:\n";
838 std::cerr << " Incorrect grid; no interpolation.\n";
839 f = 0.;
840 return false;
841 }
842 // Compute local coordinates.
843 const double xLocal =
844 (x - xAxis[iGrid - 1]) / (xAxis[iGrid] - xAxis[iGrid - 1]);
845 // Set the summing range.
846 iX0 = iGrid - 1;
847 iX1 = iGrid;
848 // Set the shape functions.
849 fX[0] = 1. - xLocal;
850 fX[1] = xLocal;
851 fX[2] = 0.;
852 fX[3] = 0.;
853 } else if (iOrder == 2) {
854 // Second order interpolation in x.
855 // Find the grid segment containing this point.
856 int iGrid = 0;
857 for (int i = 1; i < nx; i++) {
858 if ((xAxis[i - 1] - x) * (x - xAxis[i]) >= 0.) {
859 iGrid = i;
860 }
861 }
862 // Compute the local coordinate for this grid segment.
863 const double xLocal =
864 (x - xAxis[iGrid - 1]) / (xAxis[iGrid] - xAxis[iGrid - 1]);
865 // Set the summing range and shape functions.
866 if (iGrid == 1) {
867 iX0 = iGrid - 1;
868 iX1 = iGrid + 1;
869 if (xAxis[iX0] == xAxis[iX0 + 1] || xAxis[iX0] == xAxis[iX0 + 2] ||
870 xAxis[iX0 + 1] == xAxis[iX0 + 2]) {
871 std::cerr << "Boxin3:\n";
872 std::cerr << " One or more grid points in x coincide.\n";
873 std::cerr << " No interpolation.\n";
874 f = 0.;
875 return false;
876 }
877 fX[0] = (x - xAxis[iX0 + 1]) * (x - xAxis[iX0 + 2]) /
878 ((xAxis[iX0] - xAxis[iX0 + 1]) * (xAxis[iX0] - xAxis[iX0 + 2]));
879 fX[1] =
880 (x - xAxis[iX0]) * (x - xAxis[iX0 + 2]) /
881 ((xAxis[iX0 + 1] - xAxis[iX0]) * (xAxis[iX0 + 1] - xAxis[iX0 + 2]));
882 fX[2] =
883 (x - xAxis[iX0]) * (x - xAxis[iX0 + 1]) /
884 ((xAxis[iX0 + 2] - xAxis[iX0]) * (xAxis[iX0 + 2] - xAxis[iX0 + 1]));
885 } else if (iGrid == nx - 1) {
886 iX0 = iGrid - 2;
887 iX1 = iGrid;
888 if (xAxis[iX0] == xAxis[iX0 + 1] || xAxis[iX0] == xAxis[iX0 + 2] ||
889 xAxis[iX0 + 1] == xAxis[iX0 + 2]) {
890 std::cerr << "Boxin3:\n";
891 std::cerr << " One or more grid points in x coincide.\n";
892 std::cerr << " No interpolation.\n";
893 f = 0.;
894 return false;
895 }
896 fX[0] = (x - xAxis[iX0 + 1]) * (x - xAxis[iX0 + 2]) /
897 ((xAxis[iX0] - xAxis[iX0 + 1]) * (xAxis[iX0] - xAxis[iX0 + 2]));
898 fX[1] =
899 (x - xAxis[iX0]) * (x - xAxis[iX0 + 2]) /
900 ((xAxis[iX0 + 1] - xAxis[iX0]) * (xAxis[iX0 + 1] - xAxis[iX0 + 2]));
901 fX[2] =
902 (x - xAxis[iX0]) * (x - xAxis[iX0 + 1]) /
903 ((xAxis[iX0 + 2] - xAxis[iX0]) * (xAxis[iX0 + 2] - xAxis[iX0 + 1]));
904 } else {
905 iX0 = iGrid - 2;
906 iX1 = iGrid + 1;
907 if (xAxis[iX0] == xAxis[iX0 + 1] || xAxis[iX0] == xAxis[iX0 + 2] ||
908 xAxis[iX0] == xAxis[iX0 + 3] || xAxis[iX0 + 1] == xAxis[iX0 + 2] ||
909 xAxis[iX0 + 1] == xAxis[iX0 + 3] ||
910 xAxis[iX0 + 2] == xAxis[iX0 + 3]) {
911 std::cerr << "Boxin3:\n";
912 std::cerr << " One or more grid points in x coincide.\n";
913 std::cerr << " No interpolation.\n";
914 f = 0.;
915 return false;
916 }
917 fX[0] = (x - xAxis[iX0 + 1]) * (x - xAxis[iX0 + 2]) /
918 ((xAxis[iX0] - xAxis[iX0 + 1]) * (xAxis[iX0] - xAxis[iX0 + 2]));
919 fX[1] =
920 (x - xAxis[iX0]) * (x - xAxis[iX0 + 2]) /
921 ((xAxis[iX0 + 1] - xAxis[iX0]) * (xAxis[iX0 + 1] - xAxis[iX0 + 2]));
922 fX[2] =
923 (x - xAxis[iX0]) * (x - xAxis[iX0 + 1]) /
924 ((xAxis[iX0 + 2] - xAxis[iX0]) * (xAxis[iX0 + 2] - xAxis[iX0 + 1]));
925 fX[0] *= (1. - xLocal);
926 fX[1] = fX[1] * (1. - xLocal) + xLocal * (x - xAxis[iX0 + 2]) *
927 (x - xAxis[iX0 + 3]) /
928 ((xAxis[iX0 + 1] - xAxis[iX0 + 2]) *
929 (xAxis[iX0 + 1] - xAxis[iX0 + 3]));
930 fX[2] = fX[2] * (1. - xLocal) + xLocal * (x - xAxis[iX0 + 1]) *
931 (x - xAxis[iX0 + 3]) /
932 ((xAxis[iX0 + 2] - xAxis[iX0 + 1]) *
933 (xAxis[iX0 + 2] - xAxis[iX0 + 3]));
934 fX[3] = xLocal * (x - xAxis[iX0 + 1]) * (x - xAxis[iX0 + 2]) /
935 ((xAxis[iX0 + 3] - xAxis[iX0 + 1]) *
936 (xAxis[iX0 + 3] - xAxis[iX0 + 2]));
937 }
938 }
939
940 if (iOrder == 0 || ny == 1) {
941 // Zeroth order interpolation in y.
942 // Find the nearest node.
943 double dist = fabs(y - yAxis[0]);
944 int iNode = 0;
945 for (int i = 1; i < ny; i++) {
946 if (fabs(y - yAxis[i]) < dist) {
947 dist = fabs(y - yAxis[i]);
948 iNode = i;
949 }
950 }
951 // Set the summing range.
952 iY0 = iNode;
953 iY1 = iNode;
954 // Establish the shape functions.
955 fY[0] = 1.;
956 fY[1] = 0.;
957 fY[2] = 0.;
958 } else if (iOrder == 1 || ny == 2) {
959 // First order interpolation in y.
960 // Find the grid segment containing this point.
961 int iGrid = 0;
962 for (int i = 1; i < ny; i++) {
963 if ((yAxis[i - 1] - y) * (y - yAxis[i]) >= 0.) {
964 iGrid = i;
965 }
966 }
967 // Ensure there won't be divisions by zero.
968 if (yAxis[iGrid] == yAxis[iGrid - 1]) {
969 std::cerr << "Boxin3:\n";
970 std::cerr << " Incorrect grid; no interpolation.\n";
971 f = 0.;
972 return false;
973 }
974 // Compute local coordinates.
975 const double yLocal =
976 (y - yAxis[iGrid - 1]) / (yAxis[iGrid] - yAxis[iGrid - 1]);
977 // Set the summing range.
978 iY0 = iGrid - 1;
979 iY1 = iGrid;
980 // Set the shape functions.
981 fY[0] = 1. - yLocal;
982 fY[1] = yLocal;
983 fY[2] = 0.;
984 } else if (iOrder == 2) {
985 // Second order interpolation in y.
986 // Find the grid segment containing this point.
987 int iGrid = 0;
988 for (int i = 1; i < ny; i++) {
989 if ((yAxis[i - 1] - y) * (y - yAxis[i]) >= 0.) {
990 iGrid = i;
991 }
992 }
993 // Compute the local coordinate for this grid segment.
994 const double yLocal =
995 (y - yAxis[iGrid - 1]) / (yAxis[iGrid] - yAxis[iGrid - 1]);
996 // Set the summing range and shape functions.
997 // These assignments are shared by all of the following conditions,
998 // so it's easier to take them out.
999 fY[0] = (y - yAxis[iY0 + 1]) * (y - yAxis[iY0 + 2]) /
1000 ((yAxis[iY0] - yAxis[iY0 + 1]) * (yAxis[iY0] - yAxis[iY0 + 2]));
1001 fY[1] = (y - yAxis[iY0]) * (y - yAxis[iY0 + 2]) /
1002 ((yAxis[iY0 + 1] - yAxis[iY0]) * (yAxis[iY0 + 1] - yAxis[iY0 + 2]));
1003 fY[2] = (y - yAxis[iY0]) * (y - yAxis[iY0 + 1]) /
1004 ((yAxis[iY0 + 2] - yAxis[iY0]) * (yAxis[iY0 + 2] - yAxis[iY0 + 1]));
1005
1006 if (iGrid == 1) {
1007 iY0 = iGrid - 1;
1008 iY1 = iGrid + 1;
1009 if (yAxis[iY0] == yAxis[iY0 + 1] || yAxis[iY0] == yAxis[iY0 + 2] ||
1010 yAxis[iY0 + 1] == yAxis[iY0 + 2]) {
1011 std::cerr << "Boxin3:\n";
1012 std::cerr << " One or more grid points in y coincide.\n";
1013 std::cerr << " No interpolation.\n";
1014 f = 0.;
1015 return false;
1016 }
1017 fY[0] = (y - yAxis[iY0 + 1]) * (y - yAxis[iY0 + 2]) /
1018 ((yAxis[iY0] - yAxis[iY0 + 1]) * (yAxis[iY0] - yAxis[iY0 + 2]));
1019 fY[1] =
1020 (y - yAxis[iY0]) * (y - yAxis[iY0 + 2]) /
1021 ((yAxis[iY0 + 1] - yAxis[iY0]) * (yAxis[iY0 + 1] - yAxis[iY0 + 2]));
1022 fY[2] =
1023 (y - yAxis[iY0]) * (y - yAxis[iY0 + 1]) /
1024 ((yAxis[iY0 + 2] - yAxis[iY0]) * (yAxis[iY0 + 2] - yAxis[iY0 + 1]));
1025 } else if (iGrid == ny - 1) {
1026 iY0 = iGrid - 2;
1027 iY1 = iGrid;
1028 if (yAxis[iY0] == yAxis[iY0 + 1] || yAxis[iY0] == yAxis[iY0 + 2] ||
1029 yAxis[iY0 + 1] == yAxis[iY0 + 2]) {
1030 std::cerr << "Boxin3:\n";
1031 std::cerr << " One or more grid points in y coincide.\n";
1032 std::cerr << " No interpolation.\n";
1033 f = 0.;
1034 return false;
1035 }
1036 fY[0] = (y - yAxis[iY0 + 1]) * (y - yAxis[iY0 + 2]) /
1037 ((yAxis[iY0] - yAxis[iY0 + 1]) * (yAxis[iY0] - yAxis[iY0 + 2]));
1038 fY[1] =
1039 (y - yAxis[iY0]) * (y - yAxis[iY0 + 2]) /
1040 ((yAxis[iY0 + 1] - yAxis[iY0]) * (yAxis[iY0 + 1] - yAxis[iY0 + 2]));
1041 fY[2] =
1042 (y - yAxis[iY0]) * (y - yAxis[iY0 + 1]) /
1043 ((yAxis[iY0 + 2] - yAxis[iY0]) * (yAxis[iY0 + 2] - yAxis[iY0 + 1]));
1044 } else {
1045 iY0 = iGrid - 2;
1046 iY1 = iGrid + 1;
1047 if (yAxis[iY0] == yAxis[iY0 + 1] || yAxis[iY0] == yAxis[iY0 + 2] ||
1048 yAxis[iY0] == yAxis[iY0 + 3] || yAxis[iY0 + 1] == yAxis[iY0 + 2] ||
1049 yAxis[iY0 + 1] == yAxis[iY0 + 3] ||
1050 yAxis[iY0 + 2] == yAxis[iY0 + 3]) {
1051 std::cerr << "Boxin3:\n";
1052 std::cerr << " One or more grid points in y coincide.\n";
1053 std::cerr << " No interpolation.\n";
1054 f = 0.;
1055 return false;
1056 }
1057 fY[0] = (y - yAxis[iY0 + 1]) * (y - yAxis[iY0 + 2]) /
1058 ((yAxis[iY0] - yAxis[iY0 + 1]) * (yAxis[iY0] - yAxis[iY0 + 2]));
1059 fY[1] =
1060 (y - yAxis[iY0]) * (y - yAxis[iY0 + 2]) /
1061 ((yAxis[iY0 + 1] - yAxis[iY0]) * (yAxis[iY0 + 1] - yAxis[iY0 + 2]));
1062 fY[2] =
1063 (y - yAxis[iY0]) * (y - yAxis[iY0 + 1]) /
1064 ((yAxis[iY0 + 2] - yAxis[iY0]) * (yAxis[iY0 + 2] - yAxis[iY0 + 1]));
1065
1066 fY[0] *= (1. - yLocal);
1067 fY[1] = fY[1] * (1. - yLocal) + yLocal * (y - yAxis[iY0 + 2]) *
1068 (y - yAxis[iY0 + 3]) /
1069 ((yAxis[iY0 + 1] - yAxis[iY0 + 2]) *
1070 (yAxis[iY0 + 1] - yAxis[iY0 + 3]));
1071 fY[2] = fY[2] * (1. - yLocal) + yLocal * (y - yAxis[iY0 + 1]) *
1072 (y - yAxis[iY0 + 3]) /
1073 ((yAxis[iY0 + 2] - yAxis[iY0 + 1]) *
1074 (yAxis[iY0 + 2] - yAxis[iY0 + 3]));
1075 fY[3] = yLocal * (y - yAxis[iY0 + 1]) * (y - yAxis[iY0 + 2]) /
1076 ((yAxis[iY0 + 3] - yAxis[iY0 + 1]) *
1077 (yAxis[iY0 + 3] - yAxis[iY0 + 2]));
1078 }
1079 }
1080
1081 if (iOrder == 0 || nz == 1) {
1082 // Zeroth order interpolation in z.
1083 // Find the nearest node.
1084 double dist = fabs(z - zAxis[0]);
1085 int iNode = 0;
1086 for (int i = 1; i < nz; i++) {
1087 if (fabs(z - zAxis[i]) < dist) {
1088 dist = fabs(z - zAxis[i]);
1089 iNode = i;
1090 }
1091 }
1092 // Set the summing range.
1093 iZ0 = iNode;
1094 iZ1 = iNode;
1095 // Establish the shape functions.
1096 fZ[0] = 1.;
1097 fZ[1] = 0.;
1098 fZ[2] = 0.;
1099 } else if (iOrder == 1 || nz == 2) {
1100 // First order interpolation in z.
1101 // Find the grid segment containing this point.
1102 int iGrid = 0;
1103 for (int i = 1; i < nz; i++) {
1104 if ((zAxis[i - 1] - z) * (z - zAxis[i]) >= 0.) {
1105 iGrid = i;
1106 }
1107 }
1108 // Ensure there won't be divisions by zero.
1109 if (zAxis[iGrid] == zAxis[iGrid - 1]) {
1110 std::cerr << "Boxin3:\n";
1111 std::cerr << " Incorrect grid; no interpolation.\n";
1112 f = 0.;
1113 return false;
1114 }
1115 // Compute local coordinates.
1116 const double zLocal =
1117 (z - zAxis[iGrid - 1]) / (zAxis[iGrid] - zAxis[iGrid - 1]);
1118 // Set the summing range.
1119 iZ0 = iGrid - 1;
1120 iZ1 = iGrid;
1121 // Set the shape functions.
1122 fZ[0] = 1. - zLocal;
1123 fZ[1] = zLocal;
1124 fZ[2] = 0.;
1125 } else if (iOrder == 2) {
1126 // Second order interpolation in z.
1127 // Find the grid segment containing this point.
1128 int iGrid = 0;
1129 for (int i = 1; i < nz; i++) {
1130 if ((zAxis[i - 1] - z) * (z - zAxis[i]) >= 0.) {
1131 iGrid = i;
1132 }
1133 }
1134 // Compute the local coordinate for this grid segment.
1135 const double zLocal =
1136 (z - zAxis[iGrid - 1]) / (zAxis[iGrid] - zAxis[iGrid - 1]);
1137 // Set the summing range and shape functions.
1138 // These assignments are shared by all of the following conditions,
1139 // so it's easier to take them out.
1140 fZ[0] = (z - zAxis[iZ0 + 1]) * (z - zAxis[iZ0 + 2]) /
1141 ((zAxis[iZ0] - zAxis[iZ0 + 1]) * (zAxis[iZ0] - zAxis[iZ0 + 2]));
1142 fZ[1] = (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 2]) /
1143 ((zAxis[iZ0 + 1] - zAxis[iZ0]) * (zAxis[iZ0 + 1] - zAxis[iZ0 + 2]));
1144 fZ[2] = (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 1]) /
1145 ((zAxis[iZ0 + 2] - zAxis[iZ0]) * (zAxis[iZ0 + 2] - zAxis[iZ0 + 1]));
1146
1147 if (iGrid == 1) {
1148 iZ0 = iGrid - 1;
1149 iZ1 = iGrid + 1;
1150 if (zAxis[iZ0] == zAxis[iZ0 + 1] || zAxis[iZ0] == zAxis[iZ0 + 2] ||
1151 zAxis[iZ0 + 1] == zAxis[iZ0 + 2]) {
1152 std::cerr << "Boxin3:\n";
1153 std::cerr << " One or more grid points in z coincide.\n";
1154 std::cerr << " No interpolation.\n";
1155 f = 0.;
1156 return false;
1157 }
1158 fZ[0] = (z - zAxis[iZ0 + 1]) * (z - zAxis[iZ0 + 2]) /
1159 ((zAxis[iZ0] - zAxis[iZ0 + 1]) * (zAxis[iZ0] - zAxis[iZ0 + 2]));
1160 fZ[1] =
1161 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 2]) /
1162 ((zAxis[iZ0 + 1] - zAxis[iZ0]) * (zAxis[iZ0 + 1] - zAxis[iZ0 + 2]));
1163 fZ[2] =
1164 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 1]) /
1165 ((zAxis[iZ0 + 2] - zAxis[iZ0]) * (zAxis[iZ0 + 2] - zAxis[iZ0 + 1]));
1166 } else if (iGrid == nz - 1) {
1167 iZ0 = iGrid - 2;
1168 iZ1 = iGrid;
1169 if (zAxis[iZ0] == zAxis[iZ0 + 1] || zAxis[iZ0] == zAxis[iZ0 + 2] ||
1170 zAxis[iZ0 + 1] == zAxis[iZ0 + 2]) {
1171 std::cerr << "Boxin3:\n";
1172 std::cerr << " One or more grid points in z coincide.\n";
1173 std::cerr << " No interpolation.\n";
1174 f = 0.;
1175 return false;
1176 }
1177 fZ[0] = (z - zAxis[iZ0 + 1]) * (z - zAxis[iZ0 + 2]) /
1178 ((zAxis[iZ0] - zAxis[iZ0 + 1]) * (zAxis[iZ0] - zAxis[iZ0 + 2]));
1179 fZ[1] =
1180 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 2]) /
1181 ((zAxis[iZ0 + 1] - zAxis[iZ0]) * (zAxis[iZ0 + 1] - zAxis[iZ0 + 2]));
1182 fZ[2] =
1183 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 1]) /
1184 ((zAxis[iZ0 + 2] - zAxis[iZ0]) * (zAxis[iZ0 + 2] - zAxis[iZ0 + 1]));
1185 } else {
1186 iZ0 = iGrid - 2;
1187 iZ1 = iGrid + 1;
1188
1189 if (zAxis[iZ0] == zAxis[iZ0 + 1] || zAxis[iZ0] == zAxis[iZ0 + 2] ||
1190 zAxis[iZ0] == zAxis[iZ0 + 3] || zAxis[iZ0 + 1] == zAxis[iZ0 + 2] ||
1191 zAxis[iZ0 + 1] == zAxis[iZ0 + 3] ||
1192 zAxis[iZ0 + 2] == zAxis[iZ0 + 3]) {
1193 std::cerr << "Boxin3:\n";
1194 std::cerr << " One or more grid points in z coincide.\n";
1195 std::cerr << " No interpolation.\n";
1196 f = 0.;
1197 return false;
1198 }
1199
1200 fZ[0] = (z - zAxis[iZ0 + 1]) * (z - zAxis[iZ0 + 2]) /
1201 ((zAxis[iZ0] - zAxis[iZ0 + 1]) * (zAxis[iZ0] - zAxis[iZ0 + 2]));
1202 fZ[1] =
1203 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 2]) /
1204 ((zAxis[iZ0 + 1] - zAxis[iZ0]) * (zAxis[iZ0 + 1] - zAxis[iZ0 + 2]));
1205 fZ[2] =
1206 (z - zAxis[iZ0]) * (z - zAxis[iZ0 + 1]) /
1207 ((zAxis[iZ0 + 2] - zAxis[iZ0]) * (zAxis[iZ0 + 2] - zAxis[iZ0 + 1]));
1208
1209 fZ[0] *= (1. - zLocal);
1210 fZ[1] = fZ[1] * (1. - zLocal) + zLocal * (z - zAxis[iZ0 + 2]) *
1211 (z - zAxis[iZ0 + 3]) /
1212 ((zAxis[iZ0 + 1] - zAxis[iZ0 + 2]) *
1213 (zAxis[iZ0 + 1] - zAxis[iZ0 + 3]));
1214 fZ[2] = fZ[2] * (1. - zLocal) + zLocal * (z - zAxis[iZ0 + 1]) *
1215 (z - zAxis[iZ0 + 3]) /
1216 ((zAxis[iZ0 + 2] - zAxis[iZ0 + 1]) *
1217 (zAxis[iZ0 + 2] - zAxis[iZ0 + 3]));
1218 fZ[3] = zLocal * (z - zAxis[iZ0 + 1]) * (z - zAxis[iZ0 + 2]) /
1219 ((zAxis[iZ0 + 3] - zAxis[iZ0 + 1]) *
1220 (zAxis[iZ0 + 3] - zAxis[iZ0 + 2]));
1221 }
1222 }
1223
1224 f = 0.;
1225 for (int i = iX0; i <= iX1; ++i) {
1226 for (int j = iY0; j <= iY1; ++j) {
1227 for (int k = iZ0; k <= iZ1; ++k) {
1228 // std::cout << "i = " << i << ", j = " << j << ", k = " << k << "\n";
1229 // std::cout << "value: " << value[i][j][k] << "\n";
1230 // std::cout << "fX = " << fX[i - iX0] << ", fY = " << fY[j - iY0] << ",
1231 // fZ = " << fZ[k - iZ0] << "\n";
1232 f += value[i][j][k] * fX[i - iX0] * fY[j - iY0] * fZ[k - iZ0];
1233 }
1234 }
1235 }
1236 // std::cout << f << std::endl;
1237 return true;
1238}
1239}
1240}
Garfield::Numerics::Boxin3
bool Boxin3(const std::vector< std::vector< std::vector< double > > > &value, const std::vector< double > &xAxis, const std::vector< double > &yAxis, const std::vector< double > &zAxis, const int nx, const int ny, const int nz, const double xx, const double yy, const double zz, double &f, const int iOrder)
Definition: Numerics.cc:770
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::GaussKronrod15
double GaussKronrod15(double(*f)(const double), const double a, const double b)
Numerical integration using 15-point Gauss-Kronrod algorithm.
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