ARTS 2.5.4 (git: 4c0d3b4d)
interpolation_lagrange.h
Go to the documentation of this file.
1#ifndef interpolation_lagrange_h
2#define interpolation_lagrange_h
3
4#include "constants.h"
5#include "debug.h"
6#include "enums.h"
7#include "grids.h"
8#include "matpackVII.h"
9#include "nonstd.h"
10#include <algorithm>
11#include <array>
12#include <functional>
13#include <limits>
14#include <memory>
15#include <numeric>
16#include <type_traits>
17#include <vector>
18
19namespace Interpolation {
26constexpr Index cycler(const Index n, const Index N) noexcept { return n >= N ? n - N : n; }
27
39constexpr Numeric cyclic_clamp(const Numeric x,
40 const std::pair<Numeric, Numeric> xlim) noexcept {
41 if (x < xlim.first)
42 return cyclic_clamp(x + xlim.second - xlim.first, xlim);
43 if (x >= xlim.second)
44 return cyclic_clamp(x - xlim.second + xlim.first, xlim);
45 return x;
46}
47
54constexpr Numeric min_cyclic(const Numeric x,
55 const std::pair<Numeric, Numeric> xlim) noexcept {
56 const Numeric dx = xlim.second - xlim.first;
57 const bool lo = nonstd::abs(x) < nonstd::abs(x - dx);
58 const bool hi = nonstd::abs(x) < nonstd::abs(x + dx);
59 const bool me = nonstd::abs(x + dx) < nonstd::abs(x - dx);
60 return (lo and hi) ? x : me ? x + dx : x - dx;
61}
62
70constexpr bool full_cycle(const Numeric xo,
71 const Numeric xn,
72 const std::pair<Numeric, Numeric> xlim) noexcept {
73 return (xo == xlim.first and xn == xlim.second) or
74 (xo == xlim.second and xn == xlim.first);
75}
76
84template <class SortedVectorType>
85constexpr Index start_pos_finder(const Numeric x, const SortedVectorType& xvec) noexcept {
86 if (const Index n = xvec.size(); n > 1) {
87 const Numeric x0 = xvec[ 0];
88 const Numeric x1 = xvec[n - 1];
89 const Numeric frac = (x - x0) / (x1 - x0);
90 const auto start_pos = Index(frac * (Numeric)(n - 2));
91 return start_pos > 0 ? (start_pos < n ? start_pos : n - 1) : 0;
92 }
93 return 0;
94}
95
97
113constexpr Index IMAX(const Index a, const Index b) noexcept { return a > b ? a : b; }
114
116
132constexpr Index IMIN(const Index a, const Index b) noexcept { return a < b ? a : b; }
133
148template <class SortedVectorType>
149constexpr Index pos_finder(const Index pos0, const Numeric x, const SortedVectorType& xi,
150 const Index polyorder, const bool cyclic,
151 const bool ascending,
152 const std::pair<Numeric, Numeric> cycle = {
153 -180, 180}) noexcept {
154 if (cyclic and (x < cycle.first or x > cycle.second))
155 // We are below or above the cycle so we must redo calculations after clamping x to the cycle
156 return pos_finder(pos0, cyclic_clamp(x, cycle), xi, polyorder, cyclic, ascending, cycle);
157 const Index N = xi.size()-1;
158 Index p0=pos0;
159
160 // Loops to find the first position with a neighbor larger or smaller
161 if (ascending) {
162 while (p0 < N and xi[p0] < x) ++p0;
163 while (p0 > 0 and xi[p0] > x) --p0;
164 } else {
165 while (p0 < N and xi[p0] > x) ++p0;
166 while (p0 > 0 and xi[p0] < x) --p0;
167 }
168
169 // Adjustment for higher and lower polynominal orders so that x is in the middle
170 if (polyorder) {
171 if (cyclic) // Max N since we can overstep bounds
172 return IMIN(IMAX(p0 - polyorder / 2, 0), N);
173 // Max N-polyorder since we cannot overstep bounds
174 return IMIN(IMAX(p0 - polyorder / 2, 0), N-polyorder);
175 }
176 // In the nearest neighbor case, we mus choose the closest neighbor
177 if (p0 < N and (nonstd::abs(xi[p0] - x) >= nonstd::abs(xi[p0 + 1] - x)))
178 return p0 + 1;
179 return p0;
180}
181
193 Standard, /* 1-to-1 interpolation grid */
194 Cyclic, /* Cyclic interpolation grid */
195 Log, /* Natural logarithm interpolation grid */
196 Log10, /* 10-base logarithm interpolation grid */
197 Log2, /* 2-base logarithm interpolation grid */
198 SinDeg, /* Cosine in degrees interpolation grid, grid only defined [-90, 90] */
199 SinRad, /* Cosine in radians interpolation grid, grid only defined [-PI/2, PI/2] */
200 CosDeg, /* Cosine in degrees interpolation grid, grid only defined [0, 180] */
201 CosRad /* Cosine in radians interpolation grid, grid only defined [0, PI] */
202 );
203
213template <GridType type, class SortedVectorType>
214constexpr Numeric l(const Index p0, const Index n, const Numeric x,
215 const SortedVectorType& xi, const Index j,
216 [[maybe_unused]] const std::pair<Numeric, Numeric> cycle = {
217 -180, 180}) noexcept {
218 Numeric val = 1.0;
219 for (Index m = 0; m < n; m++) {
220 if (m not_eq j) {
221 if constexpr (type == GridType::Log) {
222 // Natural log weights
223 val *= (std::log(x) - std::log(xi[m + p0])) /
224 (std::log(xi[j + p0]) - std::log(xi[m + p0]));
225 } else if constexpr (type == GridType::Log10) {
226 // Common log weights
227 val *= (std::log10(x) - std::log10(xi[m + p0])) /
228 (std::log10(xi[j + p0]) - std::log10(xi[m + p0]));
229 } else if constexpr (type == GridType::Log2) {
230 // Binary log weights
231 val *= (std::log2(x) - std::log2(xi[m + p0])) /
232 (std::log2(xi[j + p0]) - std::log2(xi[m + p0]));
233 } else if constexpr (type == GridType::SinDeg) {
234 // Sine in degrees weights
235 using Conversion::sind;
236 val *= (sind(x) - sind(xi[m + p0])) / (sind(xi[j + p0]) - sind(xi[m + p0]));
237 } else if constexpr (type == GridType::SinRad) {
238 // Sine in radians weights
239 using std::sin;
240 val *= (sin(x) - sin(xi[m + p0])) / (sin(xi[j + p0]) - sin(xi[m + p0]));
241 } else if constexpr (type == GridType::CosDeg) {
242 // Cosine in degrees weights (nb. order changed)
243 using Conversion::cosd;
244 val *= (cosd(xi[m + p0]) - cosd(x)) / (cosd(xi[m + p0]) - cosd(xi[j + p0]));
245 } else if constexpr (type == GridType::CosRad) {
246 // Cosine in radians weights (nb. order changed)
247 using std::cos;
248 val *= (cos(xi[m + p0]) - cos(x)) / (cos(xi[m + p0]) - cos(xi[j + p0]));
249 } else if constexpr (type == GridType::Standard) {
250 // Linear weights, simple and straightforward
251 val *= (x - xi[m + p0]) / (xi[j + p0] - xi[m + p0]);
252 } else if constexpr (type == GridType::Cyclic) {
253 // Cyclic weights
254 // We have to ensure that all weights are cyclic (e.g., 355 degrees < -6 degrees)
255 const Index N = Index(xi.size());
256 const Index m_pos = cycler(m + p0, N);
257 const Index j_pos = cycler(j + p0, N);
258 const Numeric x_val = cyclic_clamp(x, cycle);
259 if (full_cycle(xi[0], xi[N-1], cycle)) {
260 // We ignore the last point in full cycles
261 if (j_pos == N - 1)
262 return 0;
263 if (m_pos not_eq N - 1)
264 val *= min_cyclic(x_val - xi[m_pos], cycle) / min_cyclic(xi[j_pos] - xi[m_pos], cycle);
265 } else {
266 val *= min_cyclic(x_val - xi[m_pos], cycle) / min_cyclic(xi[j_pos] - xi[m_pos], cycle);
267 }
268 }
269 }
270 }
271 return val;
272}
273
288template <GridType type, class SortedVectorType,
289 class LagrangeVectorType>
290constexpr double dl_dval(
291 const Index p0, const Index n, const Numeric x, const SortedVectorType& xi,
292 [[maybe_unused]] const LagrangeVectorType& li, const Index j, const Index i,
293 [[maybe_unused]] const std::pair<Numeric, Numeric> cycle) noexcept {
294 if constexpr (type == GridType::Standard) {
295 // Linear weights, simple and straightforward
296 if (x not_eq xi[i + p0]) {
297 // A simple case when x is not on the grid
298 return li[j] / (x - xi[i + p0]);
299 }
300 // We have to resort to the full recalculations
301 Numeric val = 1.0 / (xi[j + p0] - xi[i + p0]);
302 for (Index m = 0; m < n; m++) {
303 if (m not_eq j and m not_eq i) {
304 val *= (x - xi[m + p0]) / (xi[j + p0] - xi[m + p0]);
305 }
306 }
307 return val;
308 } else if constexpr (type == GridType::Cyclic) {
309 // Cyclic weights
310 // We have to ensure that all weights are cyclic (e.g., 355 degrees < -6 degrees)
311 const decltype(i + p0) N = xi.size();
312 const Index i_pos = cycler(i + p0, N);
313 const Index j_pos = cycler(j + p0, N);
314 const Numeric x_val = cyclic_clamp(x, cycle);
315 if (full_cycle(xi[0], xi[N-1], cycle)) {
316 // We ignore the last point in full cycles
317 if (i_pos == N - 1 or j_pos == N - 1)
318 return 0;
319 if (x_val not_eq xi[i_pos])
320 // A simple case when x is not on the grid
321 return li[j] / min_cyclic(x_val - xi[i_pos], cycle);
322
323 // We have to resort to the full recalculations
324 Numeric val = 1.0 / min_cyclic(xi[j_pos] - xi[i_pos], cycle);
325 for (Index m = 0; m < n; m++) {
326 if (m not_eq j and m not_eq i) {
327 const Index m_pos = cycler(m + p0, N);
328 if (m_pos not_eq N - 1) {
329 val *= min_cyclic(x_val - xi[m_pos], cycle) / min_cyclic(xi[j_pos] - xi[m_pos], cycle);
330 }
331 }
332 }
333 return val;
334 }
335 if (x_val not_eq xi[i_pos])
336 // A simple case when x is not on the grid
337 return li[j] / min_cyclic(x_val - xi[i_pos], cycle);
338 // We have to resort to the full recalculations
339 Numeric val = 1.0 / min_cyclic(xi[j_pos] - xi[i_pos], cycle);
340 for (Index m = 0; m < n; m++) {
341 if (m not_eq j and m not_eq i) {
342 const Index m_pos = cycler(m + p0, N);
343 val *= min_cyclic(x_val - xi[m_pos], cycle) / min_cyclic(xi[j_pos] - xi[m_pos], cycle);
344 }
345 }
346 return val;
347 } else /*if any other case */ {
348 // All other cases have to use full calculations because we need a linear derivative
349 Numeric val = 1.0 / (xi[j + p0] - xi[i + p0]);
350 for (Index m = 0; m < n; m++) {
351 if (m not_eq j and m not_eq i) {
352 val *= (x - xi[m + p0]) / (xi[j + p0] - xi[m + p0]);
353 }
354 }
355 return val;
356 }
357}
358
369template <GridType type, class SortedVectorType,
370 class LagrangeVectorType>
371constexpr Numeric dl(const Index p0, const Index n, const Numeric x,
372 const SortedVectorType& xi, const LagrangeVectorType& li,
373 const Index j,
374 const std::pair<Numeric, Numeric> cycle = {-180,
375 180}) noexcept {
376 Numeric dval = 0.0;
377 for (Index i = 0; i < n; i++) {
378 if (i not_eq j) {
379 dval += dl_dval<type>(p0, n, x, xi, li, j, i, cycle);
380 }
381 }
382 return dval;
383}
384
386template <class SortedVectorType>
387constexpr bool is_ascending(const SortedVectorType& xi) noexcept {
388 if (xi.nelem() > 1)
389 return xi[0] < xi[1];
390 return false;
391}
392
402template <class SortedVectorType>
403void check_lagrange_interpolation([[maybe_unused]] const SortedVectorType& xi,
404 [[maybe_unused]] const Index polyorder = 1,
405 [[maybe_unused]] const std::pair<Numeric, Numeric> x={std::numeric_limits<Numeric>::infinity(), -std::numeric_limits<Numeric>::infinity()},
406 [[maybe_unused]] const Numeric extrapol = 0.5,
407 [[maybe_unused]] const GridType type = GridType::Standard,
408 [[maybe_unused]] const std::pair<Numeric, Numeric> cycle = {-180, 180}) {
409 const Index n = Index(xi.size());
410
411 ARTS_USER_ERROR_IF (polyorder >= n,"Interpolation setup has failed!\n"
412 "\tRequesting greater interpolation order than possible with given input grid")
413 ARTS_USER_ERROR_IF(type == GridType::Cyclic and cycle.first >= cycle.second,
414 "Interpolation setup has failed!\n" "\tBad cycle, must be [first, second)")
415 if (polyorder and extrapol > 0 and GridType::Cyclic not_eq type) {
416 const bool ascending = is_ascending(xi);
417 const Numeric xmin = ascending ? xi[0 ] - extrapol * nonstd::abs(xi[1 ] - xi[0 ]) :
418 xi[n-1] - extrapol * nonstd::abs(xi[n-2] - xi[n-1]) ;
419 const Numeric xmax = ascending ? xi[n-1] + extrapol * nonstd::abs(xi[n-2] - xi[n-1]) :
420 xi[0 ] + extrapol * nonstd::abs(xi[1 ] - xi[0 ]) ;
421 ARTS_USER_ERROR_IF(x.first < xmin or x.second > xmax,
422 "Interpolation setup has failed!\n"
423 "\tThe new grid has limits: ", x.first, ' ', x.second, '\n',
424 "\tThe old grid has limits: ", xmin, ' ', xmax)
425 }
426}
427
428
438template <class SortedVectorType>
439void check_lagrange_interpolation([[maybe_unused]] const SortedVectorType& xi,
440 [[maybe_unused]] const Index polyorder,
441 [[maybe_unused]] const Numeric x,
442 [[maybe_unused]] const Numeric extrapol = 0.5,
443 [[maybe_unused]] const GridType type = GridType::Standard,
444 [[maybe_unused]] const std::pair<Numeric, Numeric> cycle = {-180, 180}) {
445 check_lagrange_interpolation(xi, polyorder, {x, x}, extrapol, type, cycle);
446}
447
449struct Lagrange {
452
455
458
459 /* Number of weights */
460 [[nodiscard]] Index size() const noexcept { return lx.size(); }
461
463 Lagrange(Lagrange&& l) noexcept : pos(l.pos), lx(std::move(l.lx)), dlx(std::move(l.dlx)) {}
464
466 Lagrange(const Lagrange& l) = default;
467
470 pos = l.pos;
471 lx = std::move(l.lx);
472 dlx = std::move(l.dlx);
473 return *this;
474 }
475
477 Lagrange& operator=(const Lagrange& l) = default;
478
489 template <class SortedVectorType>
490 Lagrange(const Index pos0, const Numeric x, const SortedVectorType& xi,
491 const Index polyorder = 1,
492 const bool do_derivs = false,
493 const GridType type = GridType::Standard,
494 const std::pair<Numeric, Numeric> cycle = {-180, 180}) noexcept {
495 const Index n = xi.size();
496 if (n > 1) {
497 const Index p = polyorder + 1;
498 const bool ascending = is_ascending(xi);
499
500 pos = pos_finder(pos0, x, xi, polyorder, type == GridType::Cyclic, ascending, cycle);
501 lx = lx_finder(pos, p, x, xi, type, cycle);
502 if (do_derivs) dlx = dlx_finder(pos, p, x, xi, lx, type, cycle);
503 } else {
504 pos = 0;
505 lx = Array<Numeric>(1, 1);
506 dlx = Array<Numeric>(1, 0);
507 }
508 }
509
511 Lagrange() noexcept : lx(1, 1), dlx(1, 0) {}
512
514 friend std::ostream& operator<<(std::ostream& os, const Lagrange& l) {
515 os << "pos: " << l.pos << '\n' << "lx: ";
516 for (auto x : l.lx) os << ' ' << x;
517 os << '\n' << "dlx: ";
518 for (auto x : l.dlx) os << ' ' << x;
519 return os << '\n';
520 }
521
522 private:
532 template <class SortedVectorType>
534 const Index p0, const Index n, const Numeric x,
535 const SortedVectorType& xi, const GridType type,
536 const std::pair<Numeric, Numeric> cycle) noexcept {
537 Array<Numeric> out(n);
538 switch (type) {
539 case GridType::Standard:
540 for (Index j = 0; j < n; j++)
541 out[j] = l<GridType::Standard>(p0, n, x, xi, j);
542 break;
543 case GridType::Log:
544 for (Index j = 0; j < n; j++)
545 out[j] = l<GridType::Log>(p0, n, x, xi, j);
546 break;
547 case GridType::Log10:
548 for (Index j = 0; j < n; j++)
549 out[j] = l<GridType::Log10>(p0, n, x, xi, j);
550 break;
551 case GridType::Log2:
552 for (Index j = 0; j < n; j++)
553 out[j] = l<GridType::Log2>(p0, n, x, xi, j);
554 break;
555 case GridType::Cyclic:
556 for (Index j = 0; j < n; j++)
557 out[j] = l<GridType::Cyclic>(p0, n, x, xi, j, cycle);
558 break;
559 case GridType::SinDeg:
560 for (Index j = 0; j < n; j++)
561 out[j] = l<GridType::SinDeg>(p0, n, x, xi, j);
562 break;
563 case GridType::SinRad:
564 for (Index j = 0; j < n; j++)
565 out[j] = l<GridType::SinRad>(p0, n, x, xi, j);
566 break;
567 case GridType::CosDeg:
568 for (Index j = 0; j < n; j++)
569 out[j] = l<GridType::CosDeg>(p0, n, x, xi, j);
570 break;
571 case GridType::CosRad:
572 for (Index j = 0; j < n; j++)
573 out[j] = l<GridType::CosRad>(p0, n, x, xi, j);
574 break;
575 case GridType::FINAL: { /* pass */
576 }
577 }
578 return out;
579 }
580
591 template <class SortedVectorType>
593 const Index p0, const Index n, const Numeric x,
594 const SortedVectorType& xi, const Array<Numeric>& li,
595 const GridType type,
596 const std::pair<Numeric, Numeric> cycle) noexcept {
597 Array<Numeric> out(n);
598 switch (type) {
599 case GridType::Standard:
600 for (Index j = 0; j < n; j++)
601 out[j] = dl<GridType::Standard>(p0, n, x, xi, li, j);
602 break;
603 case GridType::Log:
604 for (Index j = 0; j < n; j++)
605 out[j] = dl<GridType::Log>(p0, n, x, xi, li, j);
606 break;
607 case GridType::Log10:
608 for (Index j = 0; j < n; j++)
609 out[j] = dl<GridType::Log10>(p0, n, x, xi, li, j);
610 break;
611 case GridType::Log2:
612 for (Index j = 0; j < n; j++)
613 out[j] = dl<GridType::Log2>(p0, n, x, xi, li, j);
614 break;
615 case GridType::Cyclic:
616 for (Index j = 0; j < n; j++)
617 out[j] = dl<GridType::Cyclic>(p0, n, x, xi, li, j, cycle);
618 break;
619 case GridType::SinDeg:
620 for (Index j = 0; j < n; j++)
621 out[j] = dl<GridType::SinDeg>(p0, n, x, xi, li, j);
622 break;
623 case GridType::SinRad:
624 for (Index j = 0; j < n; j++)
625 out[j] = dl<GridType::SinRad>(p0, n, x, xi, li, j);
626 break;
627 case GridType::CosDeg:
628 for (Index j = 0; j < n; j++)
629 out[j] = dl<GridType::CosDeg>(p0, n, x, xi, li, j);
630 break;
631 case GridType::CosRad:
632 for (Index j = 0; j < n; j++)
633 out[j] = dl<GridType::CosRad>(p0, n, x, xi, li, j);
634 break;
635 case GridType::FINAL: { /* pass */
636 }
637 }
638 return out;
639 }
640}; // Lagrange
641
643template <std::size_t PolyOrder>
647
649 std::array<Numeric, PolyOrder + 1> lx;
650
652 std::array<Numeric, PolyOrder + 1> dlx;
653
654 /* Number of weights */
655 static constexpr Index size() noexcept { return PolyOrder + 1; }
656
658 constexpr FixedLagrange(FixedLagrange&& l) noexcept : pos(l.pos), lx(std::move(l.lx)), dlx(std::move(l.dlx)) {}
659
661 constexpr FixedLagrange& operator=(FixedLagrange&& l) noexcept {
662 pos = l.pos;
663 lx = std::move(l.lx);
664 dlx = std::move(l.dlx);
665 return *this;
666 }
667
677 template <class SortedVectorType>
678 constexpr FixedLagrange(const Index p0, const Numeric x,
679 const SortedVectorType& xi,
680 const bool do_derivs = false,
681 const GridType type = GridType::Standard,
682 const std::pair<Numeric, Numeric> cycle = {-180, 180}) noexcept
683 : pos(pos_finder(p0, x, xi, PolyOrder,
684 type == GridType::Cyclic,
685 xi.size() > 1 ? xi[0] < xi[1] : false, cycle)),
686 lx(lx_finder(pos, x, xi, type, cycle)),
687 dlx(do_derivs ? dlx_finder(pos, x, xi, lx, type, cycle)
688 : std::array<Numeric, PolyOrder + 1>{}) {}
689
691 friend std::ostream& operator<<(std::ostream& os, const FixedLagrange& l) {
692 os << "pos: " << l.pos << '\n' << "lx: ";
693 for (auto x : l.lx) os << ' ' << x;
694 os << '\n' << "dlx: ";
695 for (auto x : l.dlx) os << ' ' << x;
696 return os << '\n';
697 }
698
699 private:
708 template <class SortedVectorType>
709 static constexpr std::array<Numeric, PolyOrder + 1> lx_finder(
710 const Index p0, const Numeric x, const SortedVectorType& xi,
711 const GridType type,
712 const std::pair<Numeric, Numeric> cycle) noexcept {
713 std::array<Numeric, PolyOrder + 1> out{};
714 constexpr Index n = PolyOrder + 1;
715 switch (type) {
716 case GridType::Standard:
717 for (Index j = 0; j < n; j++)
718 out[j] = l<GridType::Standard>(p0, n, x, xi, j);
719 break;
720 case GridType::Log:
721 for (Index j = 0; j < n; j++)
722 out[j] = l<GridType::Log>(p0, n, x, xi, j);
723 break;
724 case GridType::Log10:
725 for (Index j = 0; j < n; j++)
726 out[j] = l<GridType::Log10>(p0, n, x, xi, j);
727 break;
728 case GridType::Log2:
729 for (Index j = 0; j < n; j++)
730 out[j] = l<GridType::Log2>(p0, n, x, xi, j);
731 break;
732 case GridType::Cyclic:
733 for (Index j = 0; j < n; j++)
734 out[j] = l<GridType::Cyclic>(p0, n, x, xi, j, cycle);
735 break;
736 case GridType::SinDeg:
737 for (Index j = 0; j < n; j++)
738 out[j] = l<GridType::SinDeg>(p0, n, x, xi, j);
739 break;
740 case GridType::SinRad:
741 for (Index j = 0; j < n; j++)
742 out[j] = l<GridType::SinRad>(p0, n, x, xi, j);
743 break;
744 case GridType::CosDeg:
745 for (Index j = 0; j < n; j++)
746 out[j] = l<GridType::CosDeg>(p0, n, x, xi, j);
747 break;
748 case GridType::CosRad:
749 for (Index j = 0; j < n; j++)
750 out[j] = l<GridType::CosRad>(p0, n, x, xi, j);
751 break;
752 case GridType::FINAL: { /* pass */
753 }
754 }
755 return out;
756 }
757
767 template <class SortedVectorType>
768 static constexpr std::array<Numeric, PolyOrder + 1> dlx_finder(
769 const Index p0, const Numeric x, const SortedVectorType& xi,
770 const std::array<Numeric, PolyOrder + 1>& li, const GridType type,
771 const std::pair<Numeric, Numeric> cycle) noexcept {
772 std::array<Numeric, PolyOrder + 1> out{};
773 constexpr Index n = PolyOrder + 1;
774 switch (type) {
775 case GridType::Standard:
776 for (Index j = 0; j < n; j++)
777 out[j] = dl<GridType::Standard>(p0, n, x, xi, li, j);
778 break;
779 case GridType::Log:
780 for (Index j = 0; j < n; j++)
781 out[j] = dl<GridType::Log>(p0, n, x, xi, li, j);
782 break;
783 case GridType::Log10:
784 for (Index j = 0; j < n; j++)
785 out[j] = dl<GridType::Log10>(p0, n, x, xi, li, j);
786 break;
787 case GridType::Log2:
788 for (Index j = 0; j < n; j++)
789 out[j] = dl<GridType::Log2>(p0, n, x, xi, li, j);
790 break;
791 case GridType::SinDeg:
792 for (Index j = 0; j < n; j++)
793 out[j] = dl<GridType::SinDeg>(p0, n, x, xi, li, j);
794 break;
795 case GridType::SinRad:
796 for (Index j = 0; j < n; j++)
797 out[j] = dl<GridType::SinRad>(p0, n, x, xi, li, j);
798 break;
799 case GridType::CosDeg:
800 for (Index j = 0; j < n; j++)
801 out[j] = dl<GridType::CosDeg>(p0, n, x, xi, li, j);
802 break;
803 case GridType::CosRad:
804 for (Index j = 0; j < n; j++)
805 out[j] = dl<GridType::CosRad>(p0, n, x, xi, li, j);
806 break;
807 case GridType::Cyclic:
808 for (Index j = 0; j < n; j++)
809 out[j] = dl<GridType::Cyclic>(p0, n, x, xi, li, j, cycle);
810 break;
811 case GridType::FINAL: { /* pass */
812 }
813 }
814 return out;
815 }
816}; // FixedLagrange
817
821
836 const ConstVectorView& x, const ConstVectorView& xi, const Index polyorder,
837 const Numeric extrapol=0.5, const bool do_derivs=false, const GridType type=GridType::Standard, const std::pair<Numeric, Numeric> cycle={-180, 180});
838
850template <std::size_t PolyOrder, class UnsortedVectorType,
851 class SortedVectorType>
853 const UnsortedVectorType& xs, const SortedVectorType& xi,
854 const bool do_derivs=false, const GridType type=GridType::Standard, const std::pair<Numeric, Numeric> cycle={-180, 180}) {
856 out.reserve(xs.size());
857 bool has_one = false;
858 for (auto x : xs) {
859 if (has_one) {
860 out.emplace_back(out.back().pos, x, xi, do_derivs, type, cycle);
861 } else {
862 out.emplace_back(start_pos_finder(x, xi), x, xi, do_derivs, type, cycle);
863 has_one = true;
864 }
865 }
866 return out;
867}
868
872
876
882Vector interpweights(const Lagrange& dim0);
883
890
896template <std::size_t PolyOrder>
897constexpr const std::array<Numeric, PolyOrder + 1>& interpweights(
898 const FixedLagrange<PolyOrder>& dim0) {
899 return dim0.lx;
900}
901
907template <std::size_t PolyOrder>
909 const Array<FixedLagrange<PolyOrder>>& dim0) {
910 Grid<std::array<Numeric, PolyOrder + 1>, 1> out(dim0.size());
911 for (std::size_t i = 0; i < dim0.size(); i++) out(i) = interpweights(dim0[i]);
912 return out;
913}
914
918
924Vector dinterpweights(const Lagrange& dim0);
925
932
938template <std::size_t PolyOrder>
939constexpr std::array<Numeric, PolyOrder + 1> dinterpweights(
940 const FixedLagrange<PolyOrder>& dim0) {
941 return dim0.dlx;
942}
943
949template <std::size_t PolyOrder>
951 const Array<FixedLagrange<PolyOrder>>& dim0) {
952 Grid<std::array<Numeric, PolyOrder + 1>, 1> out(dim0.size());
953 for (std::size_t i = 0; i < dim0.size(); i++)
954 out(i) = dinterpweights(dim0[i]);
955 return out;
956}
957
961
970Numeric interp(const ConstVectorView& yi, const ConstVectorView& iw,
971 const Lagrange& dim0);
972
981template <std::size_t PolyOrder, class VectorType>
982constexpr Numeric interp(const VectorType& yi,
983 const std::array<Numeric, PolyOrder + 1>& iw,
984 const FixedLagrange<PolyOrder>& dim0) {
985 Numeric out(0.0);
986 const Index I = yi.size();
987 for (Index i = 0; i < dim0.size(); i++)
988 out += iw[i] * yi[cycler(i + dim0.pos, I)];
989 return out;
990}
991
995
1005void reinterp(VectorView out, const ConstVectorView& iy,
1006 const Grid<Vector, 1>& iw, const Array<Lagrange>& dim0);
1007
1016Vector reinterp(const ConstVectorView& iy, const Grid<Vector, 1>& iw,
1017 const Array<Lagrange>& dim0);
1018
1027template <std::size_t PolyOrder>
1029 const Grid<std::array<Numeric, PolyOrder + 1>, 1>& iw,
1030 const Array<FixedLagrange<PolyOrder>>& dim0) {
1031 Vector out(dim0.size());
1032 for (std::size_t i = 0; i < dim0.size(); i++)
1033 out[i] = interp(iy, iw(i), dim0[i]);
1034 return out;
1035}
1036
1040
1044
1051Matrix interpweights(const Lagrange& dim0, const Lagrange& dim1);
1052
1060 const Array<Lagrange>& dim1);
1067template <std::size_t PolyOrder0, std::size_t PolyOrder1>
1069 const FixedLagrange<PolyOrder0>& dim0,
1070 const FixedLagrange<PolyOrder1>& dim1) {
1072 for (Index i = 0; i < dim0.size(); i++)
1073 for (Index j = 0; j < dim1.size(); j++) out(i, j) = dim0.lx[i] * dim1.lx[j];
1074 return out;
1075}
1076
1083template <std::size_t PolyOrder0, std::size_t PolyOrder1>
1085 const Array<FixedLagrange<PolyOrder0>>& dim0,
1086 const Array<FixedLagrange<PolyOrder1>>& dim1) {
1088 dim1.size());
1089 for (std::size_t i = 0; i < dim0.size(); i++)
1090 for (std::size_t j = 0; j < dim1.size(); j++)
1091 out(i, j) = interpweights(dim0[i], dim1[j]);
1092 return out;
1093}
1094
1098
1106Matrix dinterpweights(const Lagrange& dim0, const Lagrange& dim1, Index dim);
1107
1116 const Array<Lagrange>& dim1, Index dim);
1117
1125template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1126 std::size_t DerivativeDim>
1128 const FixedLagrange<PolyOrder0>& dim0,
1129 const FixedLagrange<PolyOrder1>& dim1) {
1130 static_assert(DerivativeDim < 2);
1131
1133 for (Index i = 0; i < dim0.size(); i++)
1134 for (Index j = 0; j < dim1.size(); j++)
1135 out(i, j) = (DerivativeDim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
1136 (DerivativeDim == 1 ? dim1.dlx[j] : dim1.lx[j]);
1137 return out;
1138}
1139
1147template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1148 std::size_t DerivativeDim>
1150 const Array<FixedLagrange<PolyOrder0>>& dim0,
1151 const Array<FixedLagrange<PolyOrder1>>& dim1) {
1152 static_assert(DerivativeDim < 2);
1153
1155 dim1.size());
1156 for (std::size_t i = 0; i < dim0.size(); i++)
1157 for (std::size_t j = 0; j < dim1.size(); j++)
1158 out(i, j) = dinterpweights<PolyOrder0, PolyOrder1, DerivativeDim>(
1159 dim0[i], dim1[j]);
1160 return out;
1161}
1162
1166
1176Numeric interp(const ConstMatrixView& yi, const ConstMatrixView& iw,
1177 const Lagrange& dim0, const Lagrange& dim1);
1178
1188template <std::size_t PolyOrder0, std::size_t PolyOrder1, class MatrixType>
1189constexpr Numeric interp(
1190 const MatrixType& yi,
1192 const FixedLagrange<PolyOrder0>& dim0,
1193 const FixedLagrange<PolyOrder1>& dim1) {
1194 Numeric out(0.0);
1195 const Index I = yi.nrows();
1196 const Index J = yi.ncols();
1197 for (Index i = 0; i < dim0.size(); i++)
1198 for (Index j = 0; j < dim1.size(); j++)
1199 out += iw(i, j) * yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J));
1200 return out;
1201}
1202
1206
1217void reinterp(MatrixView out, const ConstMatrixView& iy,
1218 const Grid<Matrix, 2>& iw, const Array<Lagrange>& dim0,
1219 const Array<Lagrange>& dim1);
1220
1230Matrix reinterp(const ConstMatrixView& iy, const Grid<Matrix, 2>& iw,
1231 const Array<Lagrange>& dim0,
1232 const Array<Lagrange>& dim1);
1233
1243template <std::size_t PolyOrder0, std::size_t PolyOrder1>
1245 const ConstMatrixView& iy,
1247 const Array<FixedLagrange<PolyOrder0>>& dim0,
1248 const Array<FixedLagrange<PolyOrder1>>& dim1) {
1249 Matrix out(dim0.size(), dim1.size());
1250 for (std::size_t i = 0; i < dim0.size(); i++)
1251 for (std::size_t j = 0; j < dim1.size(); j++)
1252 out(i, j) = interp(iy, iw(i, j), dim0[i], dim1[j]);
1253 return out;
1254}
1255
1259
1263
1271Tensor3 interpweights(const Lagrange& dim0, const Lagrange& dim1,
1272 const Lagrange& dim2);
1273
1282 const Array<Lagrange>& dim1,
1283 const Array<Lagrange>& dim2);
1284
1292template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1293 std::size_t PolyOrder2>
1296 const FixedLagrange<PolyOrder1>& dim1,
1297 const FixedLagrange<PolyOrder2>& dim2) {
1299 for (Index i = 0; i < dim0.size(); i++)
1300 for (Index j = 0; j < dim1.size(); j++)
1301 for (Index k = 0; k < dim2.size(); k++)
1302 out(i, j, k) = dim0.lx[i] * dim1.lx[j] * dim2.lx[k];
1303 return out;
1304}
1305
1313template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1314 std::size_t PolyOrder2>
1317 const Array<FixedLagrange<PolyOrder1>>& dim1,
1318 const Array<FixedLagrange<PolyOrder2>>& dim2) {
1320 out(dim0.size(), dim1.size(), dim2.size());
1321 for (std::size_t i = 0; i < dim0.size(); i++)
1322 for (std::size_t j = 0; j < dim1.size(); j++)
1323 for (std::size_t k = 0; k < dim2.size(); k++)
1324 out(i, j, k) = interpweights(dim0[i], dim1[j], dim2[k]);
1325 return out;
1326}
1327
1331
1340Tensor3 dinterpweights(const Lagrange& dim0, const Lagrange& dim1,
1341 const Lagrange& dim2, Index dim);
1342
1352 const Array<Lagrange>& dim1,
1353 const Array<Lagrange>& dim2, Index dim);
1354
1363template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1364 std::size_t PolyOrder2, std::size_t DerivativeDim>
1367 const FixedLagrange<PolyOrder1>& dim1,
1368 const FixedLagrange<PolyOrder2>& dim2) {
1369 static_assert(DerivativeDim < 3);
1370
1372 for (Index i = 0; i < dim0.size(); i++)
1373 for (Index j = 0; j < dim1.size(); j++)
1374 for (Index k = 0; k < dim2.size(); k++)
1375 out(i, j, k) = (DerivativeDim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
1376 (DerivativeDim == 1 ? dim1.dlx[j] : dim1.lx[j]) *
1377 (DerivativeDim == 2 ? dim2.dlx[k] : dim2.lx[k]);
1378 return out;
1379}
1380
1389template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1390 std::size_t PolyOrder2, std::size_t DerivativeDim>
1393 const Array<FixedLagrange<PolyOrder1>>& dim1,
1394 const Array<FixedLagrange<PolyOrder2>>& dim2) {
1396 out(dim0.size(), dim0.size(), dim2.size());
1397 for (std::size_t i = 0; i < dim0.size(); i++)
1398 for (std::size_t j = 0; j < dim1.size(); j++)
1399 for (std::size_t k = 0; k < dim2.size(); k++)
1400 out(i, j, k) =
1401 dinterpweights<PolyOrder0, PolyOrder1, PolyOrder2, DerivativeDim>(
1402 dim0[i], dim1[j], dim2[k]);
1403 return out;
1404}
1405
1409
1420Numeric interp(const ConstTensor3View& yi, const ConstTensor3View& iw,
1421 const Lagrange& dim0, const Lagrange& dim1,
1422 const Lagrange& dim2);
1423
1434template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1435 std::size_t PolyOrder2, class Tensor3Type>
1436constexpr Numeric interp(const Tensor3Type& yi,
1437 const FixedGrid<Numeric, PolyOrder0 + 1,
1438 PolyOrder1 + 1, PolyOrder2 + 1>& iw,
1439 const FixedLagrange<PolyOrder0>& dim0,
1440 const FixedLagrange<PolyOrder1>& dim1,
1441 const FixedLagrange<PolyOrder2>& dim2) {
1442 Numeric out(0.0);
1443 const Index I = yi.npages();
1444 const Index J = yi.nrows();
1445 const Index K = yi.ncols();
1446 for (Index i = 0; i < dim0.size(); i++)
1447 for (Index j = 0; j < dim1.size(); j++)
1448 for (Index k = 0; k < dim2.size(); k++)
1449 out +=
1450 iw(i, j, k) * yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J),
1451 cycler(k + dim2.pos, K));
1452 return out;
1453}
1454
1458
1470void reinterp(Tensor3View out, const ConstTensor3View& iy,
1471 const Grid<Tensor3, 3>& iw, const Array<Lagrange>& dim0,
1472 const Array<Lagrange>& dim1,
1473 const Array<Lagrange>& dim2);
1474
1486 const Array<Lagrange>& dim0,
1487 const Array<Lagrange>& dim1,
1488 const Array<Lagrange>& dim2);
1489
1500template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1501 std::size_t PolyOrder2>
1503 const ConstTensor3View& iy,
1504 const Grid<
1506 iw,
1507 const Array<FixedLagrange<PolyOrder0>>& dim0,
1508 const Array<FixedLagrange<PolyOrder1>>& dim1,
1509 const Array<FixedLagrange<PolyOrder2>>& dim2) {
1510 Tensor3 out(dim0.size(), dim1.size(), dim2.size());
1511 for (std::size_t i = 0; i < dim0.size(); i++)
1512 for (std::size_t j = 0; j < dim1.size(); j++)
1513 for (std::size_t k = 0; k < dim2.size(); k++)
1514 out(i, j, k) = interp(iy, iw(i, j, k), dim0[i], dim1[j], dim2[k]);
1515 return out;
1516}
1517
1521
1525
1534Tensor4 interpweights(const Lagrange& dim0, const Lagrange& dim1,
1535 const Lagrange& dim2, const Lagrange& dim3);
1536
1546 const Array<Lagrange>& dim1,
1547 const Array<Lagrange>& dim2,
1548 const Array<Lagrange>& dim3);
1549
1558template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1559 std::size_t PolyOrder2, std::size_t PolyOrder3>
1560constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1561 PolyOrder3 + 1>
1563 const FixedLagrange<PolyOrder1>& dim1,
1564 const FixedLagrange<PolyOrder2>& dim2,
1565 const FixedLagrange<PolyOrder3>& dim3) {
1566 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1567 PolyOrder3 + 1>
1568 out;
1569 for (Index i = 0; i < dim0.size(); i++)
1570 for (Index j = 0; j < dim1.size(); j++)
1571 for (Index k = 0; k < dim2.size(); k++)
1572 for (Index l = 0; l < dim3.size(); l++)
1573 out(i, j, k, l) = dim0.lx[i] * dim1.lx[j] * dim2.lx[k] * dim3.lx[l];
1574 return out;
1575}
1576
1585template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1586 std::size_t PolyOrder2, std::size_t PolyOrder3>
1587Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1588 PolyOrder3 + 1>,
1589 4>
1591 const Array<FixedLagrange<PolyOrder1>>& dim1,
1592 const Array<FixedLagrange<PolyOrder2>>& dim2,
1593 const Array<FixedLagrange<PolyOrder3>>& dim3) {
1594 Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1595 PolyOrder3 + 1>,
1596 4>
1597 out(dim0.size(), dim1.size(), dim2.size(), dim3.size());
1598 for (std::size_t i = 0; i < dim0.size(); i++)
1599 for (std::size_t j = 0; j < dim1.size(); j++)
1600 for (std::size_t k = 0; k < dim2.size(); k++)
1601 for (std::size_t l = 0; l < dim3.size(); l++)
1602 out(i, j, k, l) = interpweights(dim0[i], dim1[j], dim2[k], dim3[l]);
1603 return out;
1604}
1605
1609
1619Tensor4 dinterpweights(const Lagrange& dim0, const Lagrange& dim1,
1620 const Lagrange& dim2, const Lagrange& dim3, Index dim);
1621
1632 const Array<Lagrange>& dim1,
1633 const Array<Lagrange>& dim2,
1634 const Array<Lagrange>& dim3, Index dim);
1635
1645template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1646 std::size_t PolyOrder2, std::size_t PolyOrder3,
1647 std::size_t DerivativeDim>
1648constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1649 PolyOrder3 + 1>
1651 const FixedLagrange<PolyOrder1>& dim1,
1652 const FixedLagrange<PolyOrder2>& dim2,
1653 const FixedLagrange<PolyOrder3>& dim3) {
1654 static_assert(DerivativeDim < 4);
1655
1656 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1657 PolyOrder3 + 1>
1658 out;
1659 for (Index i = 0; i < dim0.size(); i++)
1660 for (Index j = 0; j < dim1.size(); j++)
1661 for (Index k = 0; k < dim2.size(); k++)
1662 for (Index l = 0; l < dim3.size(); l++)
1663 out(i, j, k, l) = (DerivativeDim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
1664 (DerivativeDim == 1 ? dim1.dlx[j] : dim1.lx[j]) *
1665 (DerivativeDim == 2 ? dim2.dlx[k] : dim2.lx[k]) *
1666 (DerivativeDim == 3 ? dim3.dlx[l] : dim3.lx[l]);
1667 return out;
1668}
1669
1679template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1680 std::size_t PolyOrder2, std::size_t PolyOrder3,
1681 std::size_t DerivativeDim>
1682Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1683 PolyOrder3 + 1>,
1684 4>
1686 const Array<FixedLagrange<PolyOrder1>>& dim1,
1687 const Array<FixedLagrange<PolyOrder2>>& dim2,
1688 const Array<FixedLagrange<PolyOrder3>>& dim3) {
1689 Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1690 PolyOrder3 + 1>,
1691 4>
1692 out(dim0.size(), dim1.size(), dim2.size(), dim3.size());
1693 for (std::size_t i = 0; i < dim0.size(); i++)
1694 for (std::size_t j = 0; j < dim1.size(); j++)
1695 for (std::size_t k = 0; k < dim2.size(); k++)
1696 for (std::size_t l = 0; l < dim3.size(); l++)
1697 out(i, j, k, l) =
1698 dinterpweights<PolyOrder0, PolyOrder1, PolyOrder2, PolyOrder3,
1699 DerivativeDim>(dim0[i], dim1[j], dim2[k], dim3[l]);
1700 return out;
1701}
1702
1706
1718Numeric interp(const ConstTensor4View& yi, const ConstTensor4View& iw,
1719 const Lagrange& dim0, const Lagrange& dim1, const Lagrange& dim2,
1720 const Lagrange& dim3);
1721
1733template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1734 std::size_t PolyOrder2, std::size_t PolyOrder3, class Tensor4Type>
1735constexpr Numeric interp(
1736 const Tensor4Type& yi,
1737 const FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1738 PolyOrder3 + 1>& iw,
1739 const FixedLagrange<PolyOrder0>& dim0,
1740 const FixedLagrange<PolyOrder1>& dim1,
1741 const FixedLagrange<PolyOrder2>& dim2,
1742 const FixedLagrange<PolyOrder3>& dim3) {
1743 Numeric out(0.0);
1744 const Index I = yi.nbooks();
1745 const Index J = yi.npages();
1746 const Index K = yi.nrows();
1747 const Index L = yi.ncols();
1748 for (Index i = 0; i < dim0.size(); i++)
1749 for (Index j = 0; j < dim1.size(); j++)
1750 for (Index k = 0; k < dim2.size(); k++)
1751 for (Index l = 0; l < dim3.size(); l++)
1752 out += iw(i, j, k, l) *
1753 yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J),
1754 cycler(k + dim2.pos, K), cycler(l + dim3.pos, L));
1755 return out;
1756}
1757
1761
1774void reinterp(Tensor4View out, const ConstTensor4View& iy,
1775 const Grid<Tensor4, 4>& iw, const Array<Lagrange>& dim0,
1776 const Array<Lagrange>& dim1,
1777 const Array<Lagrange>& dim2,
1778 const Array<Lagrange>& dim3);
1779
1792 const Array<Lagrange>& dim0,
1793 const Array<Lagrange>& dim1,
1794 const Array<Lagrange>& dim2,
1795 const Array<Lagrange>& dim3);
1796
1808template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1809 std::size_t PolyOrder2, std::size_t PolyOrder3>
1811 const Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1,
1812 PolyOrder2 + 1, PolyOrder3 + 1>,
1813 4>& iw,
1814 const Array<FixedLagrange<PolyOrder0>>& dim0,
1815 const Array<FixedLagrange<PolyOrder1>>& dim1,
1816 const Array<FixedLagrange<PolyOrder2>>& dim2,
1817 const Array<FixedLagrange<PolyOrder3>>& dim3) {
1818 Tensor4 out(dim0.size(), dim1.size(), dim2.size(), dim3.size());
1819 for (std::size_t i = 0; i < dim0.size(); i++)
1820 for (std::size_t j = 0; j < dim1.size(); j++)
1821 for (std::size_t k = 0; k < dim2.size(); k++)
1822 for (std::size_t l = 0; l < dim3.size(); l++)
1823 out(i, j, k, l) =
1824 interp(iy, iw(i, j, k, l), dim0[i], dim1[j], dim2[k], dim3[l]);
1825 return out;
1826}
1827
1831
1835
1845Tensor5 interpweights(const Lagrange& dim0, const Lagrange& dim1,
1846 const Lagrange& dim2, const Lagrange& dim3,
1847 const Lagrange& dim4);
1848
1859 const Array<Lagrange>& dim1,
1860 const Array<Lagrange>& dim2,
1861 const Array<Lagrange>& dim3,
1862 const Array<Lagrange>& dim4);
1863
1873template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1874 std::size_t PolyOrder2, std::size_t PolyOrder3,
1875 std::size_t PolyOrder4>
1876constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1877 PolyOrder3 + 1, PolyOrder4 + 1>
1879 const FixedLagrange<PolyOrder1>& dim1,
1880 const FixedLagrange<PolyOrder2>& dim2,
1881 const FixedLagrange<PolyOrder3>& dim3,
1882 const FixedLagrange<PolyOrder4>& dim4) {
1883 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1884 PolyOrder3 + 1, PolyOrder4 + 1>
1885 out;
1886 for (Index i = 0; i < dim0.size(); i++)
1887 for (Index j = 0; j < dim1.size(); j++)
1888 for (Index k = 0; k < dim2.size(); k++)
1889 for (Index l = 0; l < dim3.size(); l++)
1890 for (Index m = 0; m < dim4.size(); m++)
1891 out(i, j, k, l, m) =
1892 dim0.lx[i] * dim1.lx[j] * dim2.lx[k] * dim3.lx[l] * dim4.lx[m];
1893 return out;
1894}
1895
1905template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1906 std::size_t PolyOrder2, std::size_t PolyOrder3,
1907 std::size_t PolyOrder4>
1908Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1909 PolyOrder3 + 1, PolyOrder4 + 1>,
1910 5>
1912 const Array<FixedLagrange<PolyOrder1>>& dim1,
1913 const Array<FixedLagrange<PolyOrder2>>& dim2,
1914 const Array<FixedLagrange<PolyOrder3>>& dim3,
1915 const Array<FixedLagrange<PolyOrder4>>& dim4) {
1916 Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1917 PolyOrder3 + 1, PolyOrder4 + 1>,
1918 5>
1919 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size());
1920 for (std::size_t i = 0; i < dim0.size(); i++)
1921 for (std::size_t j = 0; j < dim1.size(); j++)
1922 for (std::size_t k = 0; k < dim2.size(); k++)
1923 for (std::size_t l = 0; l < dim3.size(); l++)
1924 for (std::size_t m = 0; m < dim4.size(); m++)
1925 out(i, j, k, l, m) =
1926 interpweights(dim0[i], dim1[j], dim2[k], dim3[l], dim4[m]);
1927 return out;
1928}
1929
1933
1944Tensor5 dinterpweights(const Lagrange& dim0, const Lagrange& dim1,
1945 const Lagrange& dim2, const Lagrange& dim3,
1946 const Lagrange& dim4, Index dim);
1947
1959 const Array<Lagrange>& dim1,
1960 const Array<Lagrange>& dim2,
1961 const Array<Lagrange>& dim3,
1962 const Array<Lagrange>& dim4, Index dim);
1963
1974template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1975 std::size_t PolyOrder2, std::size_t PolyOrder3,
1976 std::size_t PolyOrder4, std::size_t DerivativeDim>
1977constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1978 PolyOrder3 + 1, PolyOrder4 + 1>
1980 const FixedLagrange<PolyOrder1>& dim1,
1981 const FixedLagrange<PolyOrder2>& dim2,
1982 const FixedLagrange<PolyOrder3>& dim3,
1983 const FixedLagrange<PolyOrder4>& dim4) {
1984 static_assert(DerivativeDim < 5);
1985
1986 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1987 PolyOrder3 + 1, PolyOrder4 + 1>
1988 out;
1989 for (Index i = 0; i < dim0.size(); i++)
1990 for (Index j = 0; j < dim1.size(); j++)
1991 for (Index k = 0; k < dim2.size(); k++)
1992 for (Index l = 0; l < dim3.size(); l++)
1993 for (Index m = 0; m < dim4.size(); m++)
1994 out(i, j, k, l, m) =
1995 (DerivativeDim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
1996 (DerivativeDim == 1 ? dim1.dlx[j] : dim1.lx[j]) *
1997 (DerivativeDim == 2 ? dim2.dlx[k] : dim2.lx[k]) *
1998 (DerivativeDim == 3 ? dim3.dlx[l] : dim3.lx[l]) *
1999 (DerivativeDim == 4 ? dim4.dlx[m] : dim4.lx[m]);
2000 return out;
2001}
2002
2013template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2014 std::size_t PolyOrder2, std::size_t PolyOrder3,
2015 std::size_t PolyOrder4, std::size_t DerivativeDim>
2016Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2017 PolyOrder3 + 1, PolyOrder4 + 1>,
2018 5>
2020 const Array<FixedLagrange<PolyOrder1>>& dim1,
2021 const Array<FixedLagrange<PolyOrder2>>& dim2,
2022 const Array<FixedLagrange<PolyOrder3>>& dim3,
2023 const Array<FixedLagrange<PolyOrder4>>& dim4) {
2024 Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2025 PolyOrder3 + 1, PolyOrder4 + 1>,
2026 5>
2027 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size());
2028 for (std::size_t i = 0; i < dim0.size(); i++)
2029 for (std::size_t j = 0; j < dim1.size(); j++)
2030 for (std::size_t k = 0; k < dim2.size(); k++)
2031 for (std::size_t l = 0; l < dim3.size(); l++)
2032 for (std::size_t m = 0; m < dim4.size(); m++)
2033 out(i, j, k, l, m) =
2034 dinterpweights<PolyOrder0, PolyOrder1, PolyOrder2, PolyOrder3,
2035 PolyOrder4, DerivativeDim>(
2036 dim0[i], dim1[j], dim2[k], dim3[l], dim4[m]);
2037 return out;
2038}
2039
2043
2056Numeric interp(const ConstTensor5View& yi, const ConstTensor5View& iw,
2057 const Lagrange& dim0, const Lagrange& dim1, const Lagrange& dim2,
2058 const Lagrange& dim3, const Lagrange& dim4);
2059
2072template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2073 std::size_t PolyOrder2, std::size_t PolyOrder3,
2074 std::size_t PolyOrder4, class Tensor5Type>
2075constexpr Numeric interp(
2076 const Tensor5Type& yi,
2077 const FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2078 PolyOrder3 + 1, PolyOrder4 + 1>& iw,
2079 const FixedLagrange<PolyOrder0>& dim0,
2080 const FixedLagrange<PolyOrder1>& dim1,
2081 const FixedLagrange<PolyOrder2>& dim2,
2082 const FixedLagrange<PolyOrder3>& dim3,
2083 const FixedLagrange<PolyOrder4>& dim4) {
2084 Numeric out(0.0);
2085 const Index I = yi.nshelves();
2086 const Index J = yi.nbooks();
2087 const Index K = yi.npages();
2088 const Index L = yi.nrows();
2089 const Index M = yi.ncols();
2090 for (Index i = 0; i < dim0.size(); i++)
2091 for (Index j = 0; j < dim1.size(); j++)
2092 for (Index k = 0; k < dim2.size(); k++)
2093 for (Index l = 0; l < dim3.size(); l++)
2094 for (Index m = 0; m < dim4.size(); m++)
2095 out += iw(i, j, k, l, m) *
2096 yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J),
2097 cycler(k + dim2.pos, K), cycler(l + dim3.pos, L),
2098 cycler(m + dim4.pos, M));
2099 return out;
2100}
2101
2105
2119void reinterp(Tensor5View out, const ConstTensor5View& iy,
2120 const Grid<Tensor5, 5>& iw, const Array<Lagrange>& dim0,
2121 const Array<Lagrange>& dim1,
2122 const Array<Lagrange>& dim2,
2123 const Array<Lagrange>& dim3,
2124 const Array<Lagrange>& dim4);
2125
2139 const Array<Lagrange>& dim0,
2140 const Array<Lagrange>& dim1,
2141 const Array<Lagrange>& dim2,
2142 const Array<Lagrange>& dim3,
2143 const Array<Lagrange>& dim4);
2144
2157template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2158 std::size_t PolyOrder2, std::size_t PolyOrder3,
2159 std::size_t PolyOrder4>
2161 const ConstTensor5View& iy,
2162 const Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1,
2163 PolyOrder2 + 1, PolyOrder3 + 1, PolyOrder4 + 1>,
2164 5>& iw,
2165 const FixedLagrange<PolyOrder0>& dim0,
2166 const FixedLagrange<PolyOrder1>& dim1,
2167 const FixedLagrange<PolyOrder2>& dim2,
2168 const FixedLagrange<PolyOrder3>& dim3,
2169 const FixedLagrange<PolyOrder4>& dim4) {
2170 Tensor5 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size());
2171 for (std::size_t i = 0; i < dim0.size(); i++)
2172 for (std::size_t j = 0; j < dim1.size(); j++)
2173 for (std::size_t k = 0; k < dim2.size(); k++)
2174 for (std::size_t l = 0; l < dim3.size(); l++)
2175 for (std::size_t m = 0; m < dim4.size(); m++)
2176 out(i, j, k, l, m) = interp(iy, iw(i, j, k, l, m), dim0[i], dim1[j],
2177 dim2[k], dim3[l], dim4[m]);
2178 return out;
2179}
2180
2184
2188
2199Tensor6 interpweights(const Lagrange& dim0, const Lagrange& dim1,
2200 const Lagrange& dim2, const Lagrange& dim3,
2201 const Lagrange& dim4, const Lagrange& dim5);
2202
2214 const Array<Lagrange>& dim1,
2215 const Array<Lagrange>& dim2,
2216 const Array<Lagrange>& dim3,
2217 const Array<Lagrange>& dim4,
2218 const Array<Lagrange>& dim5);
2219
2230template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2231 std::size_t PolyOrder2, std::size_t PolyOrder3,
2232 std::size_t PolyOrder4, std::size_t PolyOrder5>
2233constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2234 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>
2236 const FixedLagrange<PolyOrder1>& dim1,
2237 const FixedLagrange<PolyOrder2>& dim2,
2238 const FixedLagrange<PolyOrder3>& dim3,
2239 const FixedLagrange<PolyOrder4>& dim4,
2240 const FixedLagrange<PolyOrder5>& dim5) {
2241 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2242 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>
2243 out;
2244 for (Index i = 0; i < dim0.size(); i++)
2245 for (Index j = 0; j < dim1.size(); j++)
2246 for (Index k = 0; k < dim2.size(); k++)
2247 for (Index l = 0; l < dim3.size(); l++)
2248 for (Index m = 0; m < dim4.size(); m++)
2249 for (Index n = 0; n < dim5.size(); n++)
2250 out(i, j, k, l, m, n) = dim0.lx[i] * dim1.lx[j] * dim2.lx[k] *
2251 dim3.lx[l] * dim4.lx[m] * dim5.lx[n];
2252 return out;
2253}
2254
2265template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2266 std::size_t PolyOrder2, std::size_t PolyOrder3,
2267 std::size_t PolyOrder4, std::size_t PolyOrder5>
2268Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2269 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>,
2270 6>
2272 const Array<FixedLagrange<PolyOrder1>>& dim1,
2273 const Array<FixedLagrange<PolyOrder2>>& dim2,
2274 const Array<FixedLagrange<PolyOrder3>>& dim3,
2275 const Array<FixedLagrange<PolyOrder4>>& dim4,
2276 const Array<FixedLagrange<PolyOrder5>>& dim5) {
2277 Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2278 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>,
2279 6>
2280 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
2281 dim5.size());
2282 for (std::size_t i = 0; i < dim0.size(); i++)
2283 for (std::size_t j = 0; j < dim1.size(); j++)
2284 for (std::size_t k = 0; k < dim2.size(); k++)
2285 for (std::size_t l = 0; l < dim3.size(); l++)
2286 for (std::size_t m = 0; m < dim4.size(); m++)
2287 for (std::size_t n = 0; n < dim5.size(); n++)
2288 out(i, j, k, l, m, n) = interpweights(dim0[i], dim1[j], dim2[k],
2289 dim3[l], dim4[m], dim5[n]);
2290 return out;
2291}
2292
2296
2308Tensor6 dinterpweights(const Lagrange& dim0, const Lagrange& dim1,
2309 const Lagrange& dim2, const Lagrange& dim3,
2310 const Lagrange& dim4, const Lagrange& dim5, Index dim);
2311
2324 const Array<Lagrange>& dim1,
2325 const Array<Lagrange>& dim2,
2326 const Array<Lagrange>& dim3,
2327 const Array<Lagrange>& dim4,
2328 const Array<Lagrange>& dim5, Index dim);
2329
2341template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2342 std::size_t PolyOrder2, std::size_t PolyOrder3,
2343 std::size_t PolyOrder4, std::size_t PolyOrder5,
2344 std::size_t DerivativeDim>
2345constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2346 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>
2348 const FixedLagrange<PolyOrder1>& dim1,
2349 const FixedLagrange<PolyOrder2>& dim2,
2350 const FixedLagrange<PolyOrder3>& dim3,
2351 const FixedLagrange<PolyOrder4>& dim4,
2352 const FixedLagrange<PolyOrder5>& dim5) {
2353 static_assert(DerivativeDim < 6);
2354
2355 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2356 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>
2357 out;
2358 for (Index i = 0; i < dim0.size(); i++)
2359 for (Index j = 0; j < dim1.size(); j++)
2360 for (Index k = 0; k < dim2.size(); k++)
2361 for (Index l = 0; l < dim3.size(); l++)
2362 for (Index m = 0; m < dim4.size(); m++)
2363 for (Index n = 0; n < dim5.size(); n++)
2364 out(i, j, k, l, m, n) =
2365 (DerivativeDim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
2366 (DerivativeDim == 1 ? dim1.dlx[j] : dim1.lx[j]) *
2367 (DerivativeDim == 2 ? dim2.dlx[k] : dim2.lx[k]) *
2368 (DerivativeDim == 3 ? dim3.dlx[l] : dim3.lx[l]) *
2369 (DerivativeDim == 4 ? dim4.dlx[m] : dim4.lx[m]) *
2370 (DerivativeDim == 5 ? dim5.dlx[n] : dim5.lx[n]);
2371 return out;
2372}
2373
2385template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2386 std::size_t PolyOrder2, std::size_t PolyOrder3,
2387 std::size_t PolyOrder4, std::size_t PolyOrder5,
2388 std::size_t DerivativeDim>
2389Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2390 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>,
2391 6>
2393 const Array<FixedLagrange<PolyOrder1>>& dim1,
2394 const Array<FixedLagrange<PolyOrder2>>& dim2,
2395 const Array<FixedLagrange<PolyOrder3>>& dim3,
2396 const Array<FixedLagrange<PolyOrder4>>& dim4,
2397 const Array<FixedLagrange<PolyOrder5>>& dim5) {
2398 Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2399 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>,
2400 6>
2401 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
2402 dim5.size());
2403 for (std::size_t i = 0; i < dim0.size(); i++)
2404 for (std::size_t j = 0; j < dim1.size(); j++)
2405 for (std::size_t k = 0; k < dim2.size(); k++)
2406 for (std::size_t l = 0; l < dim3.size(); l++)
2407 for (std::size_t m = 0; m < dim4.size(); m++)
2408 for (std::size_t n = 0; n < dim5.size(); n++)
2409 out(i, j, k, l, m, n) =
2410 dinterpweights<PolyOrder0, PolyOrder1, PolyOrder2, PolyOrder3,
2411 PolyOrder4, PolyOrder5, DerivativeDim>(
2412 dim0[i], dim1[j], dim2[k], dim3[l], dim4[m], dim5[n]);
2413 return out;
2414}
2415
2419
2433Numeric interp(const ConstTensor6View& yi, const ConstTensor6View& iw,
2434 const Lagrange& dim0, const Lagrange& dim1, const Lagrange& dim2,
2435 const Lagrange& dim3, const Lagrange& dim4,
2436 const Lagrange& dim5);
2437
2451template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2452 std::size_t PolyOrder2, std::size_t PolyOrder3,
2453 std::size_t PolyOrder4, std::size_t PolyOrder5, class Tensor6Type>
2454constexpr Numeric interp(
2455 const Tensor6Type& yi,
2456 const FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2457 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>& iw,
2458 const FixedLagrange<PolyOrder0>& dim0,
2459 const FixedLagrange<PolyOrder1>& dim1,
2460 const FixedLagrange<PolyOrder2>& dim2,
2461 const FixedLagrange<PolyOrder3>& dim3,
2462 const FixedLagrange<PolyOrder4>& dim4,
2463 const FixedLagrange<PolyOrder5>& dim5) {
2464 Numeric out(0.0);
2465 const Index I = yi.nvitrines();
2466 const Index J = yi.nshelves();
2467 const Index K = yi.nbooks();
2468 const Index L = yi.npages();
2469 const Index M = yi.nrows();
2470 const Index N = yi.ncols();
2471 for (Index i = 0; i < dim0.size(); i++)
2472 for (Index j = 0; j < dim1.size(); j++)
2473 for (Index k = 0; k < dim2.size(); k++)
2474 for (Index l = 0; l < dim3.size(); l++)
2475 for (Index m = 0; m < dim4.size(); m++)
2476 for (Index n = 0; n < dim5.size(); n++)
2477 out += iw(i, j, k, l, m, n) *
2478 yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J),
2479 cycler(k + dim2.pos, K), cycler(l + dim3.pos, L),
2480 cycler(m + dim4.pos, M), cycler(n + dim5.pos, N));
2481 return out;
2482}
2483
2487
2502void reinterp(Tensor6View out, const ConstTensor6View& iy,
2503 const Grid<Tensor6, 6>& iw, const Array<Lagrange>& dim0,
2504 const Array<Lagrange>& dim1,
2505 const Array<Lagrange>& dim2,
2506 const Array<Lagrange>& dim3,
2507 const Array<Lagrange>& dim4,
2508 const Array<Lagrange>& dim5);
2509
2524 const Array<Lagrange>& dim0,
2525 const Array<Lagrange>& dim1,
2526 const Array<Lagrange>& dim2,
2527 const Array<Lagrange>& dim3,
2528 const Array<Lagrange>& dim4,
2529 const Array<Lagrange>& dim5);
2530
2544template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2545 std::size_t PolyOrder2, std::size_t PolyOrder3,
2546 std::size_t PolyOrder4, std::size_t PolyOrder5>
2548 const Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1,
2549 PolyOrder2 + 1, PolyOrder3 + 1,
2550 PolyOrder4 + 1, PolyOrder5 + 1>,
2551 6>& iw,
2552 const Array<FixedLagrange<PolyOrder0>>& dim0,
2553 const Array<FixedLagrange<PolyOrder1>>& dim1,
2554 const Array<FixedLagrange<PolyOrder2>>& dim2,
2555 const Array<FixedLagrange<PolyOrder3>>& dim3,
2556 const Array<FixedLagrange<PolyOrder4>>& dim4,
2557 const Array<FixedLagrange<PolyOrder5>>& dim5) {
2558 Tensor6 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
2559 dim5.size());
2560 for (std::size_t i = 0; i < dim0.size(); i++)
2561 for (std::size_t j = 0; j < dim1.size(); j++)
2562 for (std::size_t k = 0; k < dim2.size(); k++)
2563 for (std::size_t l = 0; l < dim3.size(); l++)
2564 for (std::size_t m = 0; m < dim4.size(); m++)
2565 for (std::size_t n = 0; n < dim5.size(); n++)
2566 out(i, j, k, l, m, n) =
2567 interp(iy, iw(i, j, k, l, m, n), dim0[i], dim1[j], dim2[k],
2568 dim3[l], dim4[m], dim5[n]);
2569 return out;
2570}
2571
2575
2579
2591Tensor7 interpweights(const Lagrange& dim0, const Lagrange& dim1,
2592 const Lagrange& dim2, const Lagrange& dim3,
2593 const Lagrange& dim4, const Lagrange& dim5,
2594 const Lagrange& dim6);
2595
2608 const Array<Lagrange>& dim1,
2609 const Array<Lagrange>& dim2,
2610 const Array<Lagrange>& dim3,
2611 const Array<Lagrange>& dim4,
2612 const Array<Lagrange>& dim5,
2613 const Array<Lagrange>& dim6);
2614
2626template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2627 std::size_t PolyOrder2, std::size_t PolyOrder3,
2628 std::size_t PolyOrder4, std::size_t PolyOrder5,
2629 std::size_t PolyOrder6>
2630constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2631 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1,
2632 PolyOrder6 + 1>
2634 const FixedLagrange<PolyOrder1>& dim1,
2635 const FixedLagrange<PolyOrder2>& dim2,
2636 const FixedLagrange<PolyOrder3>& dim3,
2637 const FixedLagrange<PolyOrder4>& dim4,
2638 const FixedLagrange<PolyOrder5>& dim5,
2639 const FixedLagrange<PolyOrder6>& dim6) {
2640 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2641 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1, PolyOrder6 + 1>
2642 out;
2643 for (Index i = 0; i < dim0.size(); i++)
2644 for (Index j = 0; j < dim1.size(); j++)
2645 for (Index k = 0; k < dim2.size(); k++)
2646 for (Index l = 0; l < dim3.size(); l++)
2647 for (Index m = 0; m < dim4.size(); m++)
2648 for (Index n = 0; n < dim5.size(); n++)
2649 for (Index o = 0; o < dim6.size(); o++)
2650 out(i, j, k, l, m, n, o) = dim0.lx[i] * dim1.lx[j] *
2651 dim2.lx[k] * dim3.lx[l] *
2652 dim4.lx[m] * dim5.lx[n] * dim6.lx[o];
2653 return out;
2654}
2655
2667template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2668 std::size_t PolyOrder2, std::size_t PolyOrder3,
2669 std::size_t PolyOrder4, std::size_t PolyOrder5,
2670 std::size_t PolyOrder6>
2671Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2672 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1, PolyOrder6 + 1>,
2673 7>
2675 const Array<FixedLagrange<PolyOrder1>>& dim1,
2676 const Array<FixedLagrange<PolyOrder2>>& dim2,
2677 const Array<FixedLagrange<PolyOrder3>>& dim3,
2678 const Array<FixedLagrange<PolyOrder4>>& dim4,
2679 const Array<FixedLagrange<PolyOrder5>>& dim5,
2680 const Array<FixedLagrange<PolyOrder6>>& dim6) {
2681 Grid<
2682 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2683 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1, PolyOrder6 + 1>,
2684 7>
2685 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
2686 dim5.size(), dim6.size());
2687 for (std::size_t i = 0; i < dim0.size(); i++)
2688 for (std::size_t j = 0; j < dim1.size(); j++)
2689 for (std::size_t k = 0; k < dim2.size(); k++)
2690 for (std::size_t l = 0; l < dim3.size(); l++)
2691 for (std::size_t m = 0; m < dim4.size(); m++)
2692 for (std::size_t n = 0; n < dim5.size(); n++)
2693 for (std::size_t o = 0; o < dim6.size(); o++)
2694 out(i, j, k, l, m, n, o) =
2695 interpweights(dim0[i], dim1[j], dim2[k], dim3[l], dim4[m],
2696 dim5[n], dim6[o]);
2697 return out;
2698}
2699
2703
2716Tensor7 dinterpweights(const Lagrange& dim0, const Lagrange& dim1,
2717 const Lagrange& dim2, const Lagrange& dim3,
2718 const Lagrange& dim4, const Lagrange& dim5,
2719 const Lagrange& dim6, Index dim);
2720
2734 const Array<Lagrange>& dim1,
2735 const Array<Lagrange>& dim2,
2736 const Array<Lagrange>& dim3,
2737 const Array<Lagrange>& dim4,
2738 const Array<Lagrange>& dim5,
2739 const Array<Lagrange>& dim6, Index dim);
2740
2753template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2754 std::size_t PolyOrder2, std::size_t PolyOrder3,
2755 std::size_t PolyOrder4, std::size_t PolyOrder5,
2756 std::size_t PolyOrder6, std::size_t DerivativeDim>
2757constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2758 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1,
2759 PolyOrder6 + 1>
2761 const FixedLagrange<PolyOrder1>& dim1,
2762 const FixedLagrange<PolyOrder2>& dim2,
2763 const FixedLagrange<PolyOrder3>& dim3,
2764 const FixedLagrange<PolyOrder4>& dim4,
2765 const FixedLagrange<PolyOrder5>& dim5,
2766 const FixedLagrange<PolyOrder6>& dim6) {
2767 static_assert(DerivativeDim < 7);
2768
2769 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2770 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1, PolyOrder6 + 1>
2771 out;
2772 for (Index i = 0; i < dim0.size(); i++)
2773 for (Index j = 0; j < dim1.size(); j++)
2774 for (Index k = 0; k < dim2.size(); k++)
2775 for (Index l = 0; l < dim3.size(); l++)
2776 for (Index m = 0; m < dim4.size(); m++)
2777 for (Index n = 0; n < dim5.size(); n++)
2778 for (Index o = 0; o < dim6.size(); o++)
2779 out(i, j, k, l, m, n, o) =
2780 (DerivativeDim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
2781 (DerivativeDim == 1 ? dim1.dlx[j] : dim1.lx[j]) *
2782 (DerivativeDim == 2 ? dim2.dlx[k] : dim2.lx[k]) *
2783 (DerivativeDim == 3 ? dim3.dlx[l] : dim3.lx[l]) *
2784 (DerivativeDim == 4 ? dim4.dlx[m] : dim4.lx[m]) *
2785 (DerivativeDim == 5 ? dim5.dlx[n] : dim5.lx[n]) *
2786 (DerivativeDim == 6 ? dim6.dlx[o] : dim6.lx[o]);
2787 return out;
2788}
2789
2802template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2803 std::size_t PolyOrder2, std::size_t PolyOrder3,
2804 std::size_t PolyOrder4, std::size_t PolyOrder5,
2805 std::size_t PolyOrder6, std::size_t DerivativeDim>
2806Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2807 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1, PolyOrder6 + 1>,
2808 7>
2810 const Array<FixedLagrange<PolyOrder1>>& dim1,
2811 const Array<FixedLagrange<PolyOrder2>>& dim2,
2812 const Array<FixedLagrange<PolyOrder3>>& dim3,
2813 const Array<FixedLagrange<PolyOrder4>>& dim4,
2814 const Array<FixedLagrange<PolyOrder5>>& dim5,
2815 const Array<FixedLagrange<PolyOrder6>>& dim6) {
2816 Grid<
2817 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2818 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1, PolyOrder6 + 1>,
2819 7>
2820 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
2821 dim5.size(), dim6.size());
2822 for (std::size_t i = 0; i < dim0.size(); i++)
2823 for (std::size_t j = 0; j < dim1.size(); j++)
2824 for (std::size_t k = 0; k < dim2.size(); k++)
2825 for (std::size_t l = 0; l < dim3.size(); l++)
2826 for (std::size_t m = 0; m < dim4.size(); m++)
2827 for (std::size_t n = 0; n < dim5.size(); n++)
2828 for (std::size_t o = 0; o < dim6.size(); o++)
2829 out(i, j, k, l, m, n, o) =
2830 dinterpweights<PolyOrder0, PolyOrder1, PolyOrder2,
2831 PolyOrder3, PolyOrder4, PolyOrder5,
2832 PolyOrder6, DerivativeDim>(
2833 dim0[i], dim1[j], dim2[k], dim3[l], dim4[m], dim5[n],
2834 dim6[o]);
2835 return out;
2836}
2837
2841
2856Numeric interp(const ConstTensor7View& yi, const ConstTensor7View& iw,
2857 const Lagrange& dim0, const Lagrange& dim1, const Lagrange& dim2,
2858 const Lagrange& dim3, const Lagrange& dim4, const Lagrange& dim5,
2859 const Lagrange& dim6);
2860
2875template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2876 std::size_t PolyOrder2, std::size_t PolyOrder3,
2877 std::size_t PolyOrder4, std::size_t PolyOrder5,
2878 std::size_t PolyOrder6, class Tensor7Type>
2879constexpr Numeric interp(
2880 const Tensor7Type& yi,
2881 const FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2882 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1,
2883 PolyOrder6 + 1>& iw,
2884 const FixedLagrange<PolyOrder0>& dim0,
2885 const FixedLagrange<PolyOrder1>& dim1,
2886 const FixedLagrange<PolyOrder2>& dim2,
2887 const FixedLagrange<PolyOrder3>& dim3,
2888 const FixedLagrange<PolyOrder4>& dim4,
2889 const FixedLagrange<PolyOrder5>& dim5,
2890 const FixedLagrange<PolyOrder6>& dim6) {
2891 Numeric out(0.0);
2892 const Index I = yi.nlibraries();
2893 const Index J = yi.nvitrines();
2894 const Index K = yi.nshelves();
2895 const Index L = yi.nbooks();
2896 const Index M = yi.npages();
2897 const Index N = yi.nrows();
2898 const Index O = yi.ncols();
2899 for (Index i = 0; i < dim0.size(); i++)
2900 for (Index j = 0; j < dim1.size(); j++)
2901 for (Index k = 0; k < dim2.size(); k++)
2902 for (Index l = 0; l < dim3.size(); l++)
2903 for (Index m = 0; m < dim4.size(); m++)
2904 for (Index n = 0; n < dim5.size(); n++)
2905 for (Index o = 0; o < dim6.size(); o++)
2906 out += iw(i, j, k, l, m, n, o) *
2907 yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J),
2908 cycler(k + dim2.pos, K), cycler(l + dim3.pos, L),
2909 cycler(m + dim4.pos, M), cycler(n + dim5.pos, N),
2910 cycler(o + dim6.pos, O));
2911 return out;
2912}
2913
2917
2933void reinterp(Tensor7View out, const ConstTensor7View& iy,
2934 const Grid<Tensor7, 7>& iw, const Array<Lagrange>& dim0,
2935 const Array<Lagrange>& dim1,
2936 const Array<Lagrange>& dim2,
2937 const Array<Lagrange>& dim3,
2938 const Array<Lagrange>& dim4,
2939 const Array<Lagrange>& dim5,
2940 const Array<Lagrange>& dim6);
2941
2957 const Array<Lagrange>& dim0,
2958 const Array<Lagrange>& dim1,
2959 const Array<Lagrange>& dim2,
2960 const Array<Lagrange>& dim3,
2961 const Array<Lagrange>& dim4,
2962 const Array<Lagrange>& dim5,
2963 const Array<Lagrange>& dim6);
2964
2979template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2980 std::size_t PolyOrder2, std::size_t PolyOrder3,
2981 std::size_t PolyOrder4, std::size_t PolyOrder5,
2982 std::size_t PolyOrder6>
2984 const ConstTensor7View& iy,
2985 const Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1,
2986 PolyOrder2 + 1, PolyOrder3 + 1, PolyOrder4 + 1,
2987 PolyOrder5 + 1, PolyOrder6 + 1>,
2988 7>& iw,
2989 const Array<FixedLagrange<PolyOrder0>>& dim0,
2990 const Array<FixedLagrange<PolyOrder1>>& dim1,
2991 const Array<FixedLagrange<PolyOrder2>>& dim2,
2992 const Array<FixedLagrange<PolyOrder3>>& dim3,
2993 const Array<FixedLagrange<PolyOrder4>>& dim4,
2994 const Array<FixedLagrange<PolyOrder5>>& dim5,
2995 const Array<FixedLagrange<PolyOrder6>>& dim6) {
2996 Tensor7 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
2997 dim5.size(), dim6.size());
2998 for (std::size_t i = 0; i < dim0.size(); i++)
2999 for (std::size_t j = 0; j < dim1.size(); j++)
3000 for (std::size_t k = 0; k < dim2.size(); k++)
3001 for (std::size_t l = 0; l < dim3.size(); l++)
3002 for (std::size_t m = 0; m < dim4.size(); m++)
3003 for (std::size_t n = 0; n < dim5.size(); n++)
3004 for (std::size_t o = 0; o < dim6.size(); o++)
3005 out(i, j, k, l, m, n, o) =
3006 interp(iy, iw(i, j, k, l, m, n, o), dim0[i], dim1[j],
3007 dim2[k], dim3[l], dim4[m], dim5[n], dim6[o]);
3008 return out;
3009}
3010} // namespace Interpolation
3011
3014
3017
3020
3023
3026
3028template <std::size_t N>
3030
3031#endif // interpolation_lagrange_h
A constant view of a Matrix.
Definition: matpackI.h:1050
A constant view of a Tensor3.
Definition: matpackIII.h:130
A constant view of a Tensor4.
Definition: matpackIV.h:131
A constant view of a Tensor5.
Definition: matpackV.h:141
A constant view of a Tensor6.
Definition: matpackVI.h:147
A constant view of a Tensor7.
Definition: matpackVII.h:145
A constant view of a Vector.
Definition: matpackI.h:517
Row-major fixed grid creation.
Definition: grids.h:132
Row-major grid creation.
Definition: grids.h:52
The MatrixView class.
Definition: matpackI.h:1169
The Matrix class.
Definition: matpackI.h:1270
The Tensor3View class.
Definition: matpackIII.h:244
The Tensor3 class.
Definition: matpackIII.h:344
The Tensor4View class.
Definition: matpackIV.h:290
The Tensor4 class.
Definition: matpackIV.h:427
The Tensor5View class.
Definition: matpackV.h:341
The Tensor5 class.
Definition: matpackV.h:514
The Tensor6View class.
Definition: matpackVI.h:630
The Tensor6 class.
Definition: matpackVI.h:1097
The Tensor7View class.
Definition: matpackVII.h:1301
The Tensor7 class.
Definition: matpackVII.h:2397
The VectorView class.
Definition: matpackI.h:658
The Vector class.
Definition: matpackI.h:908
Constants of physical expressions as constexpr.
Helper macros for debugging.
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
GridType
#define abs(x)
#define dx
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > MatrixType
Definition: matpackI.h:115
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
auto sind(T x) noexcept -> decltype(std::sin(deg2rad(x)))
Returns the sine of the deg2rad of the input.
Definition: constants.h:349
auto cosd(T x) noexcept -> decltype(std::cos(deg2rad(x)))
Returns the cosine of the deg2rad of the input.
Definition: constants.h:343
constexpr Numeric cyclic_clamp(const Numeric x, const std::pair< Numeric, Numeric > xlim) noexcept
constexpr Index cycler(const Index n, const Index N) noexcept
constexpr Index pos_finder(const Index pos0, const Numeric x, const SortedVectorType &xi, const Index polyorder, const bool cyclic, const bool ascending, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
constexpr double dl_dval(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const LagrangeVectorType &li, const Index j, const Index i, const std::pair< Numeric, Numeric > cycle) noexcept
Numeric interp(const ConstVectorView &yi, const ConstVectorView &iw, const Lagrange &dim0)
Array< FixedLagrange< PolyOrder > > FixedLagrangeVector(const UnsortedVectorType &xs, const SortedVectorType &xi, const bool do_derivs=false, const GridType type=GridType::Standard, const std::pair< Numeric, Numeric > cycle={-180, 180})
constexpr Index IMIN(const Index a, const Index b) noexcept
Return the minimum of two integer numbers.
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
constexpr bool is_ascending(const SortedVectorType &xi) noexcept
Checks whether the Sorted Vector is ascending or not by checking its first two elements.
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})
void reinterp(VectorView out, const ConstVectorView &iy, const Grid< Vector, 1 > &iw, const Array< Lagrange > &dim0)
constexpr Numeric dl(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const LagrangeVectorType &li, const Index j, const std::pair< Numeric, Numeric > cycle={-180, 180}) noexcept
Vector dinterpweights(const Lagrange &dim0)
constexpr bool full_cycle(const Numeric xo, const Numeric xn, const std::pair< Numeric, Numeric > xlim) noexcept
constexpr Numeric min_cyclic(const Numeric x, const std::pair< Numeric, Numeric > xlim) noexcept
ENUMCLASS(GridType, char, Standard, Cyclic, Log, Log10, Log2, SinDeg, SinRad, CosDeg, CosRad)
constexpr Index IMAX(const Index a, const Index b) noexcept
Return the maximum of two integer numbers.
constexpr Index start_pos_finder(const Numeric x, const SortedVectorType &xvec) noexcept
Vector interpweights(const Lagrange &dim0)
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: nonstd.h:5
constexpr T abs(T x) noexcept
Definition: nonstd.h:13
#define a
#define b
#define N
Definition: rng.cc:164
#define M
Definition: rng.cc:165
static constexpr Index size() noexcept
static constexpr std::array< Numeric, PolyOrder+1 > dlx_finder(const Index p0, const Numeric x, const SortedVectorType &xi, const std::array< Numeric, PolyOrder+1 > &li, const GridType type, const std::pair< Numeric, Numeric > cycle) noexcept
constexpr FixedLagrange & operator=(FixedLagrange &&l) noexcept
Enusre that the move operator exists.
friend std::ostream & operator<<(std::ostream &os, const FixedLagrange &l)
std::array< Numeric, PolyOrder+1 > dlx
static constexpr std::array< Numeric, PolyOrder+1 > lx_finder(const Index p0, const Numeric x, const SortedVectorType &xi, const GridType type, const std::pair< Numeric, Numeric > cycle) noexcept
constexpr FixedLagrange(FixedLagrange &&l) noexcept
Enusre that the move constructor exists.
std::array< Numeric, PolyOrder+1 > lx
constexpr FixedLagrange(const Index p0, const Numeric x, const SortedVectorType &xi, const bool do_derivs=false, const GridType type=GridType::Standard, const std::pair< Numeric, Numeric > cycle={-180, 180}) noexcept
A Lagrange interpolation computer.
static Array< Numeric > lx_finder(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const GridType type, const std::pair< Numeric, Numeric > cycle) noexcept
Lagrange() noexcept
Default constructor for zero-length elements.
Lagrange(const Index pos0, const Numeric x, const SortedVectorType &xi, const Index polyorder=1, const bool do_derivs=false, const GridType type=GridType::Standard, const std::pair< Numeric, Numeric > cycle={-180, 180}) noexcept
friend std::ostream & operator<<(std::ostream &os, const Lagrange &l)
Lagrange & operator=(const Lagrange &l)=default
Ensure that the copy operator exists.
Lagrange(const Lagrange &l)=default
Ensure that the copy constructor exists.
static Array< Numeric > dlx_finder(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Array< Numeric > &li, const GridType type, const std::pair< Numeric, Numeric > cycle) noexcept
Lagrange & operator=(Lagrange &&l) noexcept
Ensure that the move operator exists.
Lagrange(Lagrange &&l) noexcept
Ensure that the move constructor exists.
Index size() const noexcept