ARTS 2.5.10 (git: 2f1c442c)
interpolation_lagrange.cc
Go to the documentation of this file.
1#include "interpolation_lagrange.h"
2
3namespace Interpolation {
4
5Array<Lagrange> LagrangeVector(
6 const ConstVectorView& xs, const ConstVectorView& xi, const Index polyorder,
7 const Numeric extrapol, const bool do_derivs, const GridType type, const std::pair<Numeric, Numeric> cycle) {
8 if (xs.size() == 1) {
9 check_lagrange_interpolation(xi, polyorder, xs[0], extrapol, type, cycle);
10 } else {
11 check_lagrange_interpolation(xi, polyorder, {min(xs), max(xs)}, extrapol, type, cycle);
12 }
13 Array<Lagrange> out;
14 out.reserve(xs.size());
15 bool has_one = false;
16 for (auto x : xs) {
17 if (has_one) {
18 out.emplace_back(out.back().pos, x, xi, polyorder, do_derivs, type, cycle);
19 } else {
20 out.emplace_back(start_pos_finder(x, xi), x, xi, polyorder, do_derivs, type, cycle);
21 has_one = true;
22 }
23 }
24 return out;
25}
26
30
34
35Vector interpweights(const Lagrange& dim0) { return dim0.lx; }
36
37Grid<Vector, 1> interpweights(const Array<Lagrange>& dim0) {
38 Grid<Vector, 1> out(dim0.size());
39 for (std::size_t i = 0; i < dim0.size(); i++) out(i) = interpweights(dim0[i]);
40 return out;
41}
42
46
47Vector dinterpweights(const Lagrange& dim0) { return dim0.dlx; }
48
49Grid<Vector, 1> dinterpweights(const Array<Lagrange>& dim0) {
50 Grid<Vector, 1> out(dim0.size());
51 for (std::size_t i = 0; i < dim0.size(); i++)
52 out(i) = dinterpweights(dim0[i]);
53 return out;
54}
55
59
60Numeric interp(const ConstVectorView& yi, const ConstVectorView& iw,
61 const Lagrange& dim0) {
62 Numeric out(0.0);
63 const Index I = yi.size();
64 for (Index i = 0; i < dim0.size(); i++)
65 out += iw[i] * yi[cycler(i + dim0.pos, I)];
66 return out;
67}
68
72
73void reinterp(VectorView out, const ConstVectorView& iy,
74 const Grid<Vector, 1>& iw, const Array<Lagrange>& dim0) {
75 for (std::size_t i = 0; i < dim0.size(); i++)
76 out[i] = interp(iy, iw(i), dim0[i]);
77}
78
79Vector reinterp(const ConstVectorView& iy, const Grid<Vector, 1>& iw,
80 const Array<Lagrange>& dim0) {
81 Vector out(dim0.size());
82 for (std::size_t i = 0; i < dim0.size(); i++)
83 out[i] = interp(iy, iw(i), dim0[i]);
84 return out;
85}
86
90
94
95Matrix interpweights(const Lagrange& dim0, const Lagrange& dim1) {
96 Matrix out(dim0.size(), dim1.size());
97 for (Index i = 0; i < dim0.size(); i++)
98 for (Index j = 0; j < dim1.size(); j++) out(i, j) = dim0.lx[i] * dim1.lx[j];
99 return out;
100}
101
102Grid<Matrix, 2> interpweights(const Array<Lagrange>& dim0,
103 const Array<Lagrange>& dim1) {
104 Grid<Matrix, 2> out(dim0.size(), dim1.size());
105 for (std::size_t i = 0; i < dim0.size(); i++)
106 for (std::size_t j = 0; j < dim1.size(); j++)
107 out(i, j) = interpweights(dim0[i], dim1[j]);
108 return out;
109}
110
114
115Matrix dinterpweights(const Lagrange& dim0, const Lagrange& dim1, Index dim) {
116 Matrix out(dim0.size(), dim1.size());
117 for (Index i = 0; i < dim0.size(); i++)
118 for (Index j = 0; j < dim1.size(); j++)
119 out(i, j) = (dim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
120 (dim == 1 ? dim1.dlx[j] : dim1.lx[j]);
121 return out;
122}
123
124Grid<Matrix, 2> dinterpweights(const Array<Lagrange>& dim0,
125 const Array<Lagrange>& dim1, Index dim) {
126 Grid<Matrix, 2> out(dim0.size(), dim1.size());
127 for (std::size_t i = 0; i < dim0.size(); i++)
128 for (std::size_t j = 0; j < dim1.size(); j++)
129 out(i, j) = dinterpweights(dim0[i], dim1[j], dim);
130 return out;
131}
132
136
137Numeric interp(const ConstMatrixView& yi, const ConstMatrixView& iw,
138 const Lagrange& dim0, const Lagrange& dim1) {
139 Numeric out(0.0);
140 const Index I = yi.nrows();
141 const Index J = yi.ncols();
142 for (Index i = 0; i < dim0.size(); i++)
143 for (Index j = 0; j < dim1.size(); j++)
144 out += iw(i, j) * yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J));
145 return out;
146}
147
151
152void reinterp(MatrixView out, const ConstMatrixView& iy,
153 const Grid<Matrix, 2>& iw, const Array<Lagrange>& dim0,
154 const Array<Lagrange>& dim1) {
155 for (std::size_t i = 0; i < dim0.size(); i++)
156 for (std::size_t j = 0; j < dim1.size(); j++)
157 out(i, j) = interp(iy, iw(i, j), dim0[i], dim1[j]);
158}
159
160Matrix reinterp(const ConstMatrixView& iy, const Grid<Matrix, 2>& iw,
161 const Array<Lagrange>& dim0,
162 const Array<Lagrange>& dim1) {
163 Matrix out(dim0.size(), dim1.size());
164 for (std::size_t i = 0; i < dim0.size(); i++)
165 for (std::size_t j = 0; j < dim1.size(); j++)
166 out(i, j) = interp(iy, iw(i, j), dim0[i], dim1[j]);
167 return out;
168}
169
173
177
178Tensor3 interpweights(const Lagrange& dim0, const Lagrange& dim1,
179 const Lagrange& dim2) {
180 Tensor3 out(dim0.size(), dim1.size(), dim2.size());
181 for (Index i = 0; i < dim0.size(); i++)
182 for (Index j = 0; j < dim1.size(); j++)
183 for (Index k = 0; k < dim2.size(); k++)
184 out(i, j, k) = dim0.lx[i] * dim1.lx[j] * dim2.lx[k];
185 return out;
186}
187
188Grid<Tensor3, 3> interpweights(const Array<Lagrange>& dim0,
189 const Array<Lagrange>& dim1,
190 const Array<Lagrange>& dim2) {
191 Grid<Tensor3, 3> out(dim0.size(), dim1.size(), dim2.size());
192 for (std::size_t i = 0; i < dim0.size(); i++)
193 for (std::size_t j = 0; j < dim1.size(); j++)
194 for (std::size_t k = 0; k < dim2.size(); k++)
195 out(i, j, k) = interpweights(dim0[i], dim1[j], dim2[k]);
196 return out;
197}
198
202
203Tensor3 dinterpweights(const Lagrange& dim0, const Lagrange& dim1,
204 const Lagrange& dim2, Index dim) {
205 Tensor3 out(dim0.size(), dim1.size(), dim2.size());
206 for (Index i = 0; i < dim0.size(); i++)
207 for (Index j = 0; j < dim1.size(); j++)
208 for (Index k = 0; k < dim2.size(); k++)
209 out(i, j, k) = (dim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
210 (dim == 1 ? dim1.dlx[j] : dim1.lx[j]) *
211 (dim == 2 ? dim2.dlx[k] : dim2.lx[k]);
212 return out;
213}
214
215Grid<Tensor3, 3> dinterpweights(const Array<Lagrange>& dim0,
216 const Array<Lagrange>& dim1,
217 const Array<Lagrange>& dim2, Index dim) {
218 Grid<Tensor3, 3> out(dim0.size(), dim1.size(), dim2.size());
219 for (std::size_t i = 0; i < dim0.size(); i++)
220 for (std::size_t j = 0; j < dim1.size(); j++)
221 for (std::size_t k = 0; k < dim2.size(); k++)
222 out(i, j, k) = dinterpweights(dim0[i], dim1[j], dim2[k], dim);
223 return out;
224}
225
229
230Numeric interp(const ConstTensor3View& yi, const ConstTensor3View& iw,
231 const Lagrange& dim0, const Lagrange& dim1,
232 const Lagrange& dim2) {
233 Numeric out(0.0);
234 const Index I = yi.npages();
235 const Index J = yi.nrows();
236 const Index K = yi.ncols();
237 for (Index i = 0; i < dim0.size(); i++)
238 for (Index j = 0; j < dim1.size(); j++)
239 for (Index k = 0; k < dim2.size(); k++)
240 out += iw(i, j, k) * yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J), cycler(k + dim2.pos, K));
241 return out;
242}
243
247
248void reinterp(Tensor3View out, const ConstTensor3View& iy,
249 const Grid<Tensor3, 3>& iw, const Array<Lagrange>& dim0,
250 const Array<Lagrange>& dim1,
251 const Array<Lagrange>& dim2) {
252 for (std::size_t i = 0; i < dim0.size(); i++)
253 for (std::size_t j = 0; j < dim1.size(); j++)
254 for (std::size_t k = 0; k < dim2.size(); k++)
255 out(i, j, k) = interp(iy, iw(i, j, k), dim0[i], dim1[j], dim2[k]);
256}
257
258Tensor3 reinterp(const ConstTensor3View& iy, const Grid<Tensor3, 3>& iw,
259 const Array<Lagrange>& dim0,
260 const Array<Lagrange>& dim1,
261 const Array<Lagrange>& dim2) {
262 Tensor3 out(dim0.size(), dim1.size(), dim2.size());
263 for (std::size_t i = 0; i < dim0.size(); i++)
264 for (std::size_t j = 0; j < dim1.size(); j++)
265 for (std::size_t k = 0; k < dim2.size(); k++)
266 out(i, j, k) = interp(iy, iw(i, j, k), dim0[i], dim1[j], dim2[k]);
267 return out;
268}
269
273
277
278Tensor4 interpweights(const Lagrange& dim0, const Lagrange& dim1,
279 const Lagrange& dim2, const Lagrange& dim3) {
280 Tensor4 out(dim0.size(), dim1.size(), dim2.size(), dim3.size());
281 for (Index i = 0; i < dim0.size(); i++)
282 for (Index j = 0; j < dim1.size(); j++)
283 for (Index k = 0; k < dim2.size(); k++)
284 for (Index l = 0; l < dim3.size(); l++)
285 out(i, j, k, l) = dim0.lx[i] * dim1.lx[j] * dim2.lx[k] * dim3.lx[l];
286 return out;
287}
288
289Grid<Tensor4, 4> interpweights(const Array<Lagrange>& dim0,
290 const Array<Lagrange>& dim1,
291 const Array<Lagrange>& dim2,
292 const Array<Lagrange>& dim3) {
293 Grid<Tensor4, 4> out(dim0.size(), dim1.size(), dim2.size(), dim3.size());
294 for (std::size_t i = 0; i < dim0.size(); i++)
295 for (std::size_t j = 0; j < dim1.size(); j++)
296 for (std::size_t k = 0; k < dim2.size(); k++)
297 for (std::size_t l = 0; l < dim3.size(); l++)
298 out(i, j, k, l) = interpweights(dim0[i], dim1[j], dim2[k], dim3[l]);
299 return out;
300}
301
305
306Tensor4 dinterpweights(const Lagrange& dim0, const Lagrange& dim1,
307 const Lagrange& dim2, const Lagrange& dim3, Index dim) {
308 Tensor4 out(dim0.size(), dim1.size(), dim2.size(), dim3.size());
309 for (Index i = 0; i < dim0.size(); i++)
310 for (Index j = 0; j < dim1.size(); j++)
311 for (Index k = 0; k < dim2.size(); k++)
312 for (Index l = 0; l < dim3.size(); l++)
313 out(i, j, k, l) = (dim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
314 (dim == 1 ? dim1.dlx[j] : dim1.lx[j]) *
315 (dim == 2 ? dim2.dlx[k] : dim2.lx[k]) *
316 (dim == 3 ? dim3.dlx[l] : dim3.lx[l]);
317 return out;
318}
319
320Grid<Tensor4, 4> dinterpweights(const Array<Lagrange>& dim0,
321 const Array<Lagrange>& dim1,
322 const Array<Lagrange>& dim2,
323 const Array<Lagrange>& dim3, Index dim) {
324 Grid<Tensor4, 4> out(dim0.size(), dim1.size(), dim2.size(), dim3.size());
325 for (std::size_t i = 0; i < dim0.size(); i++)
326 for (std::size_t j = 0; j < dim1.size(); j++)
327 for (std::size_t k = 0; k < dim2.size(); k++)
328 for (std::size_t l = 0; l < dim3.size(); l++)
329 out(i, j, k, l) =
330 dinterpweights(dim0[i], dim1[j], dim2[k], dim3[l], dim);
331 return out;
332}
333
337
338Numeric interp(const ConstTensor4View& yi, const ConstTensor4View& iw,
339 const Lagrange& dim0, const Lagrange& dim1, const Lagrange& dim2,
340 const Lagrange& dim3) {
341 Numeric out(0.0);
342 const Index I = yi.nbooks();
343 const Index J = yi.npages();
344 const Index K = yi.nrows();
345 const Index L = yi.ncols();
346 for (Index i = 0; i < dim0.size(); i++)
347 for (Index j = 0; j < dim1.size(); j++)
348 for (Index k = 0; k < dim2.size(); k++)
349 for (Index l = 0; l < dim3.size(); l++)
350 out += iw(i, j, k, l) *
351 yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J),
352 cycler(k + dim2.pos, K), cycler(l + dim3.pos, L));
353 return out;
354}
355
359
360void reinterp(Tensor4View out, const ConstTensor4View& iy,
361 const Grid<Tensor4, 4>& iw, const Array<Lagrange>& dim0,
362 const Array<Lagrange>& dim1,
363 const Array<Lagrange>& dim2,
364 const Array<Lagrange>& dim3) {
365 for (std::size_t i = 0; i < dim0.size(); i++)
366 for (std::size_t j = 0; j < dim1.size(); j++)
367 for (std::size_t k = 0; k < dim2.size(); k++)
368 for (std::size_t l = 0; l < dim3.size(); l++)
369 out(i, j, k, l) =
370 interp(iy, iw(i, j, k, l), dim0[i], dim1[j], dim2[k], dim3[l]);
371}
372
373Tensor4 reinterp(const ConstTensor4View& iy, const Grid<Tensor4, 4>& iw,
374 const Array<Lagrange>& dim0,
375 const Array<Lagrange>& dim1,
376 const Array<Lagrange>& dim2,
377 const Array<Lagrange>& dim3) {
378 Tensor4 out(dim0.size(), dim1.size(), dim2.size(), dim3.size());
379 for (std::size_t i = 0; i < dim0.size(); i++)
380 for (std::size_t j = 0; j < dim1.size(); j++)
381 for (std::size_t k = 0; k < dim2.size(); k++)
382 for (std::size_t l = 0; l < dim3.size(); l++)
383 out(i, j, k, l) =
384 interp(iy, iw(i, j, k, l), dim0[i], dim1[j], dim2[k], dim3[l]);
385 return out;
386}
387
391
395
396Tensor5 interpweights(const Lagrange& dim0, const Lagrange& dim1,
397 const Lagrange& dim2, const Lagrange& dim3,
398 const Lagrange& dim4) {
399 Tensor5 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size());
400 for (Index i = 0; i < dim0.size(); i++)
401 for (Index j = 0; j < dim1.size(); j++)
402 for (Index k = 0; k < dim2.size(); k++)
403 for (Index l = 0; l < dim3.size(); l++)
404 for (Index m = 0; m < dim4.size(); m++)
405 out(i, j, k, l, m) =
406 dim0.lx[i] * dim1.lx[j] * dim2.lx[k] * dim3.lx[l] * dim4.lx[m];
407 return out;
408}
409
410Grid<Tensor5, 5> interpweights(const Array<Lagrange>& dim0,
411 const Array<Lagrange>& dim1,
412 const Array<Lagrange>& dim2,
413 const Array<Lagrange>& dim3,
414 const Array<Lagrange>& dim4) {
415 Grid<Tensor5, 5> out(dim0.size(), dim1.size(), dim2.size(), dim3.size(),
416 dim4.size());
417 for (std::size_t i = 0; i < dim0.size(); i++)
418 for (std::size_t j = 0; j < dim1.size(); j++)
419 for (std::size_t k = 0; k < dim2.size(); k++)
420 for (std::size_t l = 0; l < dim3.size(); l++)
421 for (std::size_t m = 0; m < dim4.size(); m++)
422 out(i, j, k, l, m) =
423 interpweights(dim0[i], dim1[j], dim2[k], dim3[l], dim4[m]);
424 return out;
425}
426
430
431Tensor5 dinterpweights(const Lagrange& dim0, const Lagrange& dim1,
432 const Lagrange& dim2, const Lagrange& dim3,
433 const Lagrange& dim4, Index dim) {
434 Tensor5 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size());
435 for (Index i = 0; i < dim0.size(); i++)
436 for (Index j = 0; j < dim1.size(); j++)
437 for (Index k = 0; k < dim2.size(); k++)
438 for (Index l = 0; l < dim3.size(); l++)
439 for (Index m = 0; m < dim4.size(); m++)
440 out(i, j, k, l, m) = (dim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
441 (dim == 1 ? dim1.dlx[j] : dim1.lx[j]) *
442 (dim == 2 ? dim2.dlx[k] : dim2.lx[k]) *
443 (dim == 3 ? dim3.dlx[l] : dim3.lx[l]) *
444 (dim == 4 ? dim4.dlx[m] : dim4.lx[m]);
445 return out;
446}
447
448Grid<Tensor5, 5> dinterpweights(const Array<Lagrange>& dim0,
449 const Array<Lagrange>& dim1,
450 const Array<Lagrange>& dim2,
451 const Array<Lagrange>& dim3,
452 const Array<Lagrange>& dim4, Index dim) {
453 Grid<Tensor5, 5> out(dim0.size(), dim1.size(), dim2.size(), dim3.size(),
454 dim4.size());
455 for (std::size_t i = 0; i < dim0.size(); i++)
456 for (std::size_t j = 0; j < dim1.size(); j++)
457 for (std::size_t k = 0; k < dim2.size(); k++)
458 for (std::size_t l = 0; l < dim3.size(); l++)
459 for (std::size_t m = 0; m < dim4.size(); m++)
460 out(i, j, k, l, m) = dinterpweights(dim0[i], dim1[j], dim2[k],
461 dim3[l], dim4[m], dim);
462 return out;
463}
464
468
469Numeric interp(const ConstTensor5View& yi, const ConstTensor5View& iw,
470 const Lagrange& dim0, const Lagrange& dim1, const Lagrange& dim2,
471 const Lagrange& dim3, const Lagrange& dim4) {
472 Numeric out(0.0);
473 const Index I = yi.nshelves();
474 const Index J = yi.nbooks();
475 const Index K = yi.npages();
476 const Index L = yi.nrows();
477 const Index M = yi.ncols();
478 for (Index i = 0; i < dim0.size(); i++)
479 for (Index j = 0; j < dim1.size(); j++)
480 for (Index k = 0; k < dim2.size(); k++)
481 for (Index l = 0; l < dim3.size(); l++)
482 for (Index m = 0; m < dim4.size(); m++)
483 out += iw(i, j, k, l, m) *
484 yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J),
485 cycler(k + dim2.pos, K), cycler(l + dim3.pos, L),
486 cycler(m + dim4.pos, M));
487 return out;
488}
489
493
494void reinterp(Tensor5View out, const ConstTensor5View& iy,
495 const Grid<Tensor5, 5>& iw, const Array<Lagrange>& dim0,
496 const Array<Lagrange>& dim1,
497 const Array<Lagrange>& dim2,
498 const Array<Lagrange>& dim3,
499 const Array<Lagrange>& dim4) {
500 for (std::size_t i = 0; i < dim0.size(); i++)
501 for (std::size_t j = 0; j < dim1.size(); j++)
502 for (std::size_t k = 0; k < dim2.size(); k++)
503 for (std::size_t l = 0; l < dim3.size(); l++)
504 for (std::size_t m = 0; m < dim4.size(); m++)
505 out(i, j, k, l, m) = interp(iy, iw(i, j, k, l, m), dim0[i], dim1[j],
506 dim2[k], dim3[l], dim4[m]);
507}
508
509Tensor5 reinterp(const ConstTensor5View& iy, const Grid<Tensor5, 5>& iw,
510 const Array<Lagrange>& dim0,
511 const Array<Lagrange>& dim1,
512 const Array<Lagrange>& dim2,
513 const Array<Lagrange>& dim3,
514 const Array<Lagrange>& dim4) {
515 Tensor5 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size());
516 for (std::size_t i = 0; i < dim0.size(); i++)
517 for (std::size_t j = 0; j < dim1.size(); j++)
518 for (std::size_t k = 0; k < dim2.size(); k++)
519 for (std::size_t l = 0; l < dim3.size(); l++)
520 for (std::size_t m = 0; m < dim4.size(); m++)
521 out(i, j, k, l, m) = interp(iy, iw(i, j, k, l, m), dim0[i], dim1[j],
522 dim2[k], dim3[l], dim4[m]);
523 return out;
524}
525
529
533
534Tensor6 interpweights(const Lagrange& dim0, const Lagrange& dim1,
535 const Lagrange& dim2, const Lagrange& dim3,
536 const Lagrange& dim4, const Lagrange& dim5) {
537 Tensor6 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
538 dim5.size());
539 for (Index i = 0; i < dim0.size(); i++)
540 for (Index j = 0; j < dim1.size(); j++)
541 for (Index k = 0; k < dim2.size(); k++)
542 for (Index l = 0; l < dim3.size(); l++)
543 for (Index m = 0; m < dim4.size(); m++)
544 for (Index n = 0; n < dim5.size(); n++)
545 out(i, j, k, l, m, n) = dim0.lx[i] * dim1.lx[j] * dim2.lx[k] *
546 dim3.lx[l] * dim4.lx[m] * dim5.lx[n];
547 return out;
548}
549
550Grid<Tensor6, 6> interpweights(const Array<Lagrange>& dim0,
551 const Array<Lagrange>& dim1,
552 const Array<Lagrange>& dim2,
553 const Array<Lagrange>& dim3,
554 const Array<Lagrange>& dim4,
555 const Array<Lagrange>& dim5) {
556 Grid<Tensor6, 6> out(dim0.size(), dim1.size(), dim2.size(), dim3.size(),
557 dim4.size(), dim5.size());
558 for (std::size_t i = 0; i < dim0.size(); i++)
559 for (std::size_t j = 0; j < dim1.size(); j++)
560 for (std::size_t k = 0; k < dim2.size(); k++)
561 for (std::size_t l = 0; l < dim3.size(); l++)
562 for (std::size_t m = 0; m < dim4.size(); m++)
563 for (std::size_t n = 0; n < dim5.size(); n++)
564 out(i, j, k, l, m, n) = interpweights(dim0[i], dim1[j], dim2[k],
565 dim3[l], dim4[m], dim5[n]);
566 return out;
567}
568
572
573Tensor6 dinterpweights(const Lagrange& dim0, const Lagrange& dim1,
574 const Lagrange& dim2, const Lagrange& dim3,
575 const Lagrange& dim4, const Lagrange& dim5, Index dim) {
576 Tensor6 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
577 dim5.size());
578 for (Index i = 0; i < dim0.size(); i++)
579 for (Index j = 0; j < dim1.size(); j++)
580 for (Index k = 0; k < dim2.size(); k++)
581 for (Index l = 0; l < dim3.size(); l++)
582 for (Index m = 0; m < dim4.size(); m++)
583 for (Index n = 0; n < dim5.size(); n++)
584 out(i, j, k, l, m, n) = (dim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
585 (dim == 1 ? dim1.dlx[j] : dim1.lx[j]) *
586 (dim == 2 ? dim2.dlx[k] : dim2.lx[k]) *
587 (dim == 3 ? dim3.dlx[l] : dim3.lx[l]) *
588 (dim == 4 ? dim4.dlx[m] : dim4.lx[m]) *
589 (dim == 5 ? dim5.dlx[n] : dim5.lx[n]);
590 return out;
591}
592
593Grid<Tensor6, 6> dinterpweights(const Array<Lagrange>& dim0,
594 const Array<Lagrange>& dim1,
595 const Array<Lagrange>& dim2,
596 const Array<Lagrange>& dim3,
597 const Array<Lagrange>& dim4,
598 const Array<Lagrange>& dim5, Index dim) {
599 Grid<Tensor6, 6> out(dim0.size(), dim1.size(), dim2.size(), dim3.size(),
600 dim4.size(), dim5.size());
601 for (std::size_t i = 0; i < dim0.size(); i++)
602 for (std::size_t j = 0; j < dim1.size(); j++)
603 for (std::size_t k = 0; k < dim2.size(); k++)
604 for (std::size_t l = 0; l < dim3.size(); l++)
605 for (std::size_t m = 0; m < dim4.size(); m++)
606 for (std::size_t n = 0; n < dim5.size(); n++)
607 out(i, j, k, l, m, n) = dinterpweights(
608 dim0[i], dim1[j], dim2[k], dim3[l], dim4[m], dim5[n], dim);
609 return out;
610}
611
615
616Numeric interp(const ConstTensor6View& yi, const ConstTensor6View& iw,
617 const Lagrange& dim0, const Lagrange& dim1, const Lagrange& dim2,
618 const Lagrange& dim3, const Lagrange& dim4,
619 const Lagrange& dim5) {
620 Numeric out(0.0);
621 const Index I = yi.nvitrines();
622 const Index J = yi.nshelves();
623 const Index K = yi.nbooks();
624 const Index L = yi.npages();
625 const Index M = yi.nrows();
626 const Index N = yi.ncols();
627 for (Index i = 0; i < dim0.size(); i++)
628 for (Index j = 0; j < dim1.size(); j++)
629 for (Index k = 0; k < dim2.size(); k++)
630 for (Index l = 0; l < dim3.size(); l++)
631 for (Index m = 0; m < dim4.size(); m++)
632 for (Index n = 0; n < dim5.size(); n++)
633 out += iw(i, j, k, l, m, n) *
634 yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J),
635 cycler(k + dim2.pos, K), cycler(l + dim3.pos, L),
636 cycler(m + dim4.pos, M), cycler(n + dim5.pos, N));
637 return out;
638}
639
643
644void reinterp(Tensor6View out, const ConstTensor6View& iy,
645 const Grid<Tensor6, 6>& iw, const Array<Lagrange>& dim0,
646 const Array<Lagrange>& dim1,
647 const Array<Lagrange>& dim2,
648 const Array<Lagrange>& dim3,
649 const Array<Lagrange>& dim4,
650 const Array<Lagrange>& dim5) {
651 for (std::size_t i = 0; i < dim0.size(); i++)
652 for (std::size_t j = 0; j < dim1.size(); j++)
653 for (std::size_t k = 0; k < dim2.size(); k++)
654 for (std::size_t l = 0; l < dim3.size(); l++)
655 for (std::size_t m = 0; m < dim4.size(); m++)
656 for (std::size_t n = 0; n < dim5.size(); n++)
657 out(i, j, k, l, m, n) =
658 interp(iy, iw(i, j, k, l, m, n), dim0[i], dim1[j], dim2[k],
659 dim3[l], dim4[m], dim5[n]);
660}
661
662Tensor6 reinterp(const ConstTensor6View& iy, const Grid<Tensor6, 6>& iw,
663 const Array<Lagrange>& dim0,
664 const Array<Lagrange>& dim1,
665 const Array<Lagrange>& dim2,
666 const Array<Lagrange>& dim3,
667 const Array<Lagrange>& dim4,
668 const Array<Lagrange>& dim5) {
669 Tensor6 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
670 dim5.size());
671 for (std::size_t i = 0; i < dim0.size(); i++)
672 for (std::size_t j = 0; j < dim1.size(); j++)
673 for (std::size_t k = 0; k < dim2.size(); k++)
674 for (std::size_t l = 0; l < dim3.size(); l++)
675 for (std::size_t m = 0; m < dim4.size(); m++)
676 for (std::size_t n = 0; n < dim5.size(); n++)
677 out(i, j, k, l, m, n) =
678 interp(iy, iw(i, j, k, l, m, n), dim0[i], dim1[j], dim2[k],
679 dim3[l], dim4[m], dim5[n]);
680 return out;
681}
682
686
690
691Tensor7 interpweights(const Lagrange& dim0, const Lagrange& dim1,
692 const Lagrange& dim2, const Lagrange& dim3,
693 const Lagrange& dim4, const Lagrange& dim5,
694 const Lagrange& dim6) {
695 Tensor7 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
696 dim5.size(), dim6.size());
697 for (Index i = 0; i < dim0.size(); i++)
698 for (Index j = 0; j < dim1.size(); j++)
699 for (Index k = 0; k < dim2.size(); k++)
700 for (Index l = 0; l < dim3.size(); l++)
701 for (Index m = 0; m < dim4.size(); m++)
702 for (Index n = 0; n < dim5.size(); n++)
703 for (Index o = 0; o < dim6.size(); o++)
704 out(i, j, k, l, m, n, o) = dim0.lx[i] * dim1.lx[j] *
705 dim2.lx[k] * dim3.lx[l] *
706 dim4.lx[m] * dim5.lx[n] * dim6.lx[o];
707 return out;
708}
709
710Grid<Tensor7, 7> interpweights(const Array<Lagrange>& dim0,
711 const Array<Lagrange>& dim1,
712 const Array<Lagrange>& dim2,
713 const Array<Lagrange>& dim3,
714 const Array<Lagrange>& dim4,
715 const Array<Lagrange>& dim5,
716 const Array<Lagrange>& dim6) {
717 Grid<Tensor7, 7> out(dim0.size(), dim1.size(), dim2.size(), dim3.size(),
718 dim4.size(), dim5.size(), dim6.size());
719 for (std::size_t i = 0; i < dim0.size(); i++)
720 for (std::size_t j = 0; j < dim1.size(); j++)
721 for (std::size_t k = 0; k < dim2.size(); k++)
722 for (std::size_t l = 0; l < dim3.size(); l++)
723 for (std::size_t m = 0; m < dim4.size(); m++)
724 for (std::size_t n = 0; n < dim5.size(); n++)
725 for (std::size_t o = 0; o < dim6.size(); o++)
726 out(i, j, k, l, m, n, o) =
727 interpweights(dim0[i], dim1[j], dim2[k], dim3[l], dim4[m],
728 dim5[n], dim6[o]);
729 return out;
730}
731
735
736Tensor7 dinterpweights(const Lagrange& dim0, const Lagrange& dim1,
737 const Lagrange& dim2, const Lagrange& dim3,
738 const Lagrange& dim4, const Lagrange& dim5,
739 const Lagrange& dim6, Index dim) {
740 Tensor7 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
741 dim5.size(), dim6.size());
742 for (Index i = 0; i < dim0.size(); i++)
743 for (Index j = 0; j < dim1.size(); j++)
744 for (Index k = 0; k < dim2.size(); k++)
745 for (Index l = 0; l < dim3.size(); l++)
746 for (Index m = 0; m < dim4.size(); m++)
747 for (Index n = 0; n < dim5.size(); n++)
748 for (Index o = 0; o < dim6.size(); o++)
749 out(i, j, k, l, m, n, o) =
750 (dim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
751 (dim == 1 ? dim1.dlx[j] : dim1.lx[j]) *
752 (dim == 2 ? dim2.dlx[k] : dim2.lx[k]) *
753 (dim == 3 ? dim3.dlx[l] : dim3.lx[l]) *
754 (dim == 4 ? dim4.dlx[m] : dim4.lx[m]) *
755 (dim == 5 ? dim5.dlx[n] : dim5.lx[n]) *
756 (dim == 6 ? dim6.dlx[o] : dim6.lx[o]);
757 return out;
758}
759
760Grid<Tensor7, 7> dinterpweights(const Array<Lagrange>& dim0,
761 const Array<Lagrange>& dim1,
762 const Array<Lagrange>& dim2,
763 const Array<Lagrange>& dim3,
764 const Array<Lagrange>& dim4,
765 const Array<Lagrange>& dim5,
766 const Array<Lagrange>& dim6, Index dim) {
767 Grid<Tensor7, 7> out(dim0.size(), dim1.size(), dim2.size(), dim3.size(),
768 dim4.size(), dim5.size(), dim6.size());
769 for (std::size_t i = 0; i < dim0.size(); i++)
770 for (std::size_t j = 0; j < dim1.size(); j++)
771 for (std::size_t k = 0; k < dim2.size(); k++)
772 for (std::size_t l = 0; l < dim3.size(); l++)
773 for (std::size_t m = 0; m < dim4.size(); m++)
774 for (std::size_t n = 0; n < dim5.size(); n++)
775 for (std::size_t o = 0; o < dim6.size(); o++)
776 out(i, j, k, l, m, n, o) =
777 dinterpweights(dim0[i], dim1[j], dim2[k], dim3[l], dim4[m],
778 dim5[n], dim6[o], dim);
779 return out;
780}
781
785
786Numeric interp(const ConstTensor7View& yi, const ConstTensor7View& iw,
787 const Lagrange& dim0, const Lagrange& dim1, const Lagrange& dim2,
788 const Lagrange& dim3, const Lagrange& dim4, const Lagrange& dim5,
789 const Lagrange& dim6) {
790 Numeric out(0.0);
791 const Index I = yi.nlibraries();
792 const Index J = yi.nvitrines();
793 const Index K = yi.nshelves();
794 const Index L = yi.nbooks();
795 const Index M = yi.npages();
796 const Index N = yi.nrows();
797 const Index O = yi.ncols();
798 for (Index i = 0; i < dim0.size(); i++)
799 for (Index j = 0; j < dim1.size(); j++)
800 for (Index k = 0; k < dim2.size(); k++)
801 for (Index l = 0; l < dim3.size(); l++)
802 for (Index m = 0; m < dim4.size(); m++)
803 for (Index n = 0; n < dim5.size(); n++)
804 for (Index o = 0; o < dim6.size(); o++)
805 out += iw(i, j, k, l, m, n, o) *
806 yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J),
807 cycler(k + dim2.pos, K), cycler(l + dim3.pos, L),
808 cycler(m + dim4.pos, M), cycler(n + dim5.pos, N),
809 cycler(o + dim6.pos, O));
810 return out;
811}
812
816
817void reinterp(Tensor7View out, const ConstTensor7View& iy,
818 const Grid<Tensor7, 7>& iw, const Array<Lagrange>& dim0,
819 const Array<Lagrange>& dim1,
820 const Array<Lagrange>& dim2,
821 const Array<Lagrange>& dim3,
822 const Array<Lagrange>& dim4,
823 const Array<Lagrange>& dim5,
824 const Array<Lagrange>& dim6) {
825 for (std::size_t i = 0; i < dim0.size(); i++)
826 for (std::size_t j = 0; j < dim1.size(); j++)
827 for (std::size_t k = 0; k < dim2.size(); k++)
828 for (std::size_t l = 0; l < dim3.size(); l++)
829 for (std::size_t m = 0; m < dim4.size(); m++)
830 for (std::size_t n = 0; n < dim5.size(); n++)
831 for (std::size_t o = 0; o < dim6.size(); o++)
832 out(i, j, k, l, m, n, o) =
833 interp(iy, iw(i, j, k, l, m, n, o), dim0[i], dim1[j],
834 dim2[k], dim3[l], dim4[m], dim5[n], dim6[o]);
835}
836
837Tensor7 reinterp(const ConstTensor7View& iy, const Grid<Tensor7, 7>& iw,
838 const Array<Lagrange>& dim0,
839 const Array<Lagrange>& dim1,
840 const Array<Lagrange>& dim2,
841 const Array<Lagrange>& dim3,
842 const Array<Lagrange>& dim4,
843 const Array<Lagrange>& dim5,
844 const Array<Lagrange>& dim6) {
845 Tensor7 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
846 dim5.size(), dim6.size());
847 for (std::size_t i = 0; i < dim0.size(); i++)
848 for (std::size_t j = 0; j < dim1.size(); j++)
849 for (std::size_t k = 0; k < dim2.size(); k++)
850 for (std::size_t l = 0; l < dim3.size(); l++)
851 for (std::size_t m = 0; m < dim4.size(); m++)
852 for (std::size_t n = 0; n < dim5.size(); n++)
853 for (std::size_t o = 0; o < dim6.size(); o++)
854 out(i, j, k, l, m, n, o) =
855 interp(iy, iw(i, j, k, l, m, n, o), dim0[i], dim1[j],
856 dim2[k], dim3[l], dim4[m], dim5[n], dim6[o]);
857 return out;
858}
859} // namespace Interpolation
max
base max(const Array< base > &x)
Max function.
Definition: array.h:145
min
base min(const Array< base > &x)
Min function.
Definition: array.h:161
Array
This can be used to make arrays out of anything.
Definition: array.h:48
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:1065
ConstMatrixView::nrows
Index nrows() const noexcept
Definition: matpackI.h:1079
ConstMatrixView::ncols
Index ncols() const noexcept
Definition: matpackI.h:1080
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:133
ConstTensor3View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:145
ConstTensor3View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:148
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:151
ConstTensor4View
A constant view of a Tensor4.
Definition: matpackIV.h:133
ConstTensor4View::ncols
Index ncols() const noexcept
Definition: matpackIV.h:146
ConstTensor4View::nrows
Index nrows() const noexcept
Definition: matpackIV.h:145
ConstTensor4View::nbooks
Index nbooks() const noexcept
Definition: matpackIV.h:143
ConstTensor4View::npages
Index npages() const noexcept
Definition: matpackIV.h:144
ConstTensor5View
A constant view of a Tensor5.
Definition: matpackV.h:144
ConstTensor5View::nrows
Index nrows() const noexcept
Definition: matpackV.h:157
ConstTensor5View::ncols
Index ncols() const noexcept
Definition: matpackV.h:158
ConstTensor5View::npages
Index npages() const noexcept
Definition: matpackV.h:156
ConstTensor5View::nbooks
Index nbooks() const noexcept
Definition: matpackV.h:155
ConstTensor5View::nshelves
Index nshelves() const noexcept
Definition: matpackV.h:154
ConstTensor6View
A constant view of a Tensor6.
Definition: matpackVI.h:150
ConstTensor6View::nbooks
Index nbooks() const noexcept
Definition: matpackVI.h:160
ConstTensor6View::nvitrines
Index nvitrines() const noexcept
Definition: matpackVI.h:158
ConstTensor6View::ncols
Index ncols() const noexcept
Definition: matpackVI.h:163
ConstTensor6View::npages
Index npages() const noexcept
Definition: matpackVI.h:161
ConstTensor6View::nshelves
Index nshelves() const noexcept
Definition: matpackVI.h:159
ConstTensor6View::nrows
Index nrows() const noexcept
Definition: matpackVI.h:162
ConstTensor7View
A constant view of a Tensor7.
Definition: matpackVII.h:148
ConstTensor7View::ncols
Index ncols() const noexcept
Definition: matpackVII.h:164
ConstTensor7View::npages
Index npages() const noexcept
Definition: matpackVII.h:162
ConstTensor7View::nrows
Index nrows() const noexcept
Definition: matpackVII.h:163
ConstTensor7View::nlibraries
Index nlibraries() const noexcept
Definition: matpackVII.h:158
ConstTensor7View::nvitrines
Index nvitrines() const noexcept
Definition: matpackVII.h:159
ConstTensor7View::nshelves
Index nshelves() const noexcept
Definition: matpackVII.h:160
ConstTensor7View::nbooks
Index nbooks() const noexcept
Definition: matpackVII.h:161
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:521
ConstVectorView::size
Index size() const noexcept
Definition: matpackI.h:550
Grid
Row-major grid creation.
Definition: grids.h:52
MatrixView
The MatrixView class.
Definition: matpackI.h:1188
Matrix
The Matrix class.
Definition: matpackI.h:1285
Tensor3View
The Tensor3View class.
Definition: matpackIII.h:251
Tensor3
The Tensor3 class.
Definition: matpackIII.h:352
Tensor4View
The Tensor4View class.
Definition: matpackIV.h:296
Tensor4
The Tensor4 class.
Definition: matpackIV.h:435
Tensor5View
The Tensor5View class.
Definition: matpackV.h:348
Tensor5
The Tensor5 class.
Definition: matpackV.h:524
Tensor6View
The Tensor6View class.
Definition: matpackVI.h:635
Tensor6
The Tensor6 class.
Definition: matpackVI.h:1105
Tensor7View
The Tensor7View class.
Definition: matpackVII.h:1308
Tensor7
The Tensor7 class.
Definition: matpackVII.h:2407
VectorView
The VectorView class.
Definition: matpackI.h:674
Vector
The Vector class.
Definition: matpackI.h:910
GridType
GridType
Definition: gridded_fields.h:118
interpolation_lagrange.h
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
Interpolation::cycler
constexpr Index cycler(const Index n, const Index N) noexcept
Definition: interpolation_lagrange.h:26
Interpolation::interp
Numeric interp(const ConstVectorView &yi, const ConstVectorView &iw, const Lagrange &dim0)
Definition: interpolation_lagrange.cc:60
Interpolation::l
constexpr Numeric l(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Index j, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
Definition: interpolation_lagrange.h:214
Interpolation::check_lagrange_interpolation
void check_lagrange_interpolation(const SortedVectorType &xi, const Index polyorder=1, const std::pair< Numeric, Numeric > x={std::numeric_limits< Numeric >::infinity(), -std::numeric_limits< Numeric >::infinity()}, const Numeric extrapol=0.5, const GridType type=GridType::Standard, const std::pair< Numeric, Numeric > cycle={-180, 180})
Definition: interpolation_lagrange.h:403
Interpolation::reinterp
void reinterp(VectorView out, const ConstVectorView &iy, const Grid< Vector, 1 > &iw, const Array< Lagrange > &dim0)
Definition: interpolation_lagrange.cc:73
Interpolation::dinterpweights
Vector dinterpweights(const Lagrange &dim0)
Definition: interpolation_lagrange.cc:47
Interpolation::start_pos_finder
constexpr Index start_pos_finder(const Numeric x, const SortedVectorType &xvec) noexcept
Definition: interpolation_lagrange.h:85
Interpolation::interpweights
Vector interpweights(const Lagrange &dim0)
Definition: interpolation_lagrange.cc:35
Interpolation::LagrangeVector
Array< Lagrange > LagrangeVector(const ConstVectorView &xs, const ConstVectorView &xi, const Index polyorder, const Numeric extrapol, const bool do_derivs, const GridType type, const std::pair< Numeric, Numeric > cycle)
Definition: interpolation_lagrange.cc:5
N
#define N
Definition: rng.cc:164
M
#define M
Definition: rng.cc:165
Interpolation::Lagrange
A Lagrange interpolation computer.
Definition: interpolation_lagrange.h:449
Interpolation::Lagrange::size
Index size() const noexcept
Definition: interpolation_lagrange.h:460