85#define LOOPIT(x) for (const Numeric* x = &t##x.fd[1]; x >= &t##x.fd[0]; --x)
97 os << gp.
idx <<
" " << gp.
fd[0] <<
" " << gp.
fd[1] <<
"\n";
179 bool ascending = (old_grid[0] <= old_grid[1]);
193 old_grid[0] - extpolfac * (old_grid[1] - old_grid[0]);
195 old_grid[n_old - 1] +
196 extpolfac * (old_grid[n_old - 1] - old_grid[n_old - 2]);
209 Numeric frac = (new_grid[0] - og_min) / (og_max - og_min);
228 Numeric lower = old_grid[current_position];
229 Numeric upper = old_grid[current_position + 1];
232 for (
Index i_new = 0; i_new < n_new; ++i_new) {
236 const Numeric tng = new_grid[i_new];
251 if (tng < lower && current_position > 0) {
254 lower = old_grid[current_position];
255 }
while (tng < lower && current_position > 0);
257 upper = old_grid[current_position + 1];
259 tgp.
idx = current_position;
260 tgp.
fd[0] = (tng - lower) / (upper - lower);
261 tgp.
fd[1] = 1.0 - tgp.
fd[0];
266 if (tng >= upper && current_position < n_old - 2) {
269 upper = old_grid[current_position + 1];
270 }
while (tng >= upper && current_position < n_old - 2);
272 lower = old_grid[current_position];
274 tgp.
idx = current_position;
275 tgp.
fd[0] = (tng - lower) / (upper - lower);
276 tgp.
fd[1] = 1.0 - tgp.
fd[0];
286 tgp.
idx = current_position;
287 tgp.
fd[0] = (tng - lower) / (upper - lower);
288 tgp.
fd[1] = 1.0 - tgp.
fd[0];
311 old_grid[0] - extpolfac * (old_grid[1] - old_grid[0]);
313 old_grid[n_old - 1] +
314 extpolfac * (old_grid[n_old - 1] - old_grid[n_old - 2]);
318 Numeric frac = 1 - (new_grid[0] - og_min) / (og_max - og_min);
336 Numeric lower = old_grid[current_position];
337 Numeric upper = old_grid[current_position + 1];
339 for (
Index i_new = 0; i_new < n_new; ++i_new) {
341 const Numeric tng = new_grid[i_new];
355 if (tng > lower && current_position > 0) {
358 lower = old_grid[current_position];
359 }
while (tng > lower && current_position > 0);
361 upper = old_grid[current_position + 1];
363 tgp.
idx = current_position;
364 tgp.
fd[0] = (tng - lower) / (upper - lower);
365 tgp.
fd[1] = 1.0 - tgp.
fd[0];
369 if (tng <= upper && current_position < n_old - 2) {
372 upper = old_grid[current_position + 1];
373 }
while (tng <= upper && current_position < n_old - 2);
375 lower = old_grid[current_position];
377 tgp.
idx = current_position;
378 tgp.
fd[0] = (tng - lower) / (upper - lower);
379 tgp.
fd[1] = 1.0 - tgp.
fd[0];
390 tgp.
idx = current_position;
391 tgp.
fd[0] = (tng - lower) / (upper - lower);
392 tgp.
fd[1] = 1.0 - tgp.
fd[0];
430 gridpos(agp, old_grid,
v, extpolfac);
453 for (
Index i = 0; i < n - 1; i++) {
458 gp[n - 1].idx = n - 2;
475 gp_new.
fd[0] = gp_old.
fd[0];
476 gp_new.
fd[1] = gp_old.
fd[1];
516 if (gp.
fd[0] < 0.0) {
519 }
else if (gp.
fd[0] > 1.0) {
524 if (gp.
fd[1] < 0.0) {
527 }
else if (gp.
fd[1] > 1.0) {
559 if (gp.
fd[0] > 0.5) {
568 if (gp.
idx == n - 1) {
614 if (gp[i].idx == ie) {
659 const bool& strict) {
663 if (gp.
idx == i && gp.
fd[0] == 0) {
665 }
else if (gp.
idx == i - 1 && gp.
fd[1] == 0) {
702 if (gp.
fd[0] > 0 && gp.
fd[1] > 0) {
707 else if (gp.
fd[0] == 0) {
803 itw.
get(iti) = (*r) * (*c);
835 itw.
get(iti) = (*p) * (*r) * (*c);
870 itw.
get(iti) = (*b) * (*p) * (*r) * (*c);
908 itw.
get(iti) = (*s) * (*b) * (*p) * (*r) * (*c);
949 itw.
get(iti) = (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
983 tia +=
a.get(tc.
idx +
c) * itw.
get(iti);
1022 for (
Index r = 0; r < 2; ++r)
1024 tia +=
a.get(tr.
idx + r, tc.
idx +
c) * itw.
get(iti);
1065 for (
Index p = 0; p < 2; ++p)
1066 for (
Index r = 0; r < 2; ++r)
1112 for (
Index p = 0; p < 2; ++p)
1113 for (
Index r = 0; r < 2; ++r)
1161 for (
Index s = 0; s < 2; ++s)
1163 for (
Index p = 0; p < 2; ++p)
1164 for (
Index r = 0; r < 2; ++r)
1166 tia +=
a.get(ts.
idx + s,
1219 for (
Index s = 0; s < 2; ++s)
1221 for (
Index p = 0; p < 2; ++p)
1222 for (
Index r = 0; r < 2; ++r)
1224 tia +=
a.get(tv.
idx +
v,
1263 for (
Index i = 0; i < n; ++i) {
1298 itw.
get(i, iti) = *
c;
1335 for (
Index i = 0; i < n; ++i) {
1351 itw.
get(i, iti) = (*r) * (*c);
1391 for (
Index i = 0; i < n; ++i) {
1401 itw.
get(i, iti) = (*p) * (*r) * (*c);
1444 for (
Index i = 0; i < n; ++i) {
1456 itw.
get(i, iti) = (*b) * (*p) * (*r) * (*c);
1502 for (
Index i = 0; i < n; ++i) {
1516 itw.
get(i, iti) = (*s) * (*b) * (*p) * (*r) * (*c);
1565 for (
Index i = 0; i < n; ++i) {
1581 itw.
get(i, iti) = (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
1620 for (
Index i = 0; i < n; ++i) {
1632 tia +=
a.get(tc.
idx +
c) * itw.
get(i, iti);
1678 for (
Index i = 0; i < n; ++i) {
1689 for (
Index r = 0; r < 2; ++r)
1691 tia +=
a.get(tr.
idx + r, tc.
idx +
c) * itw.
get(i, iti);
1740 for (
Index i = 0; i < n; ++i) {
1752 for (
Index p = 0; p < 2; ++p)
1753 for (
Index r = 0; r < 2; ++r)
1755 tia +=
a.get(tp.
idx + p, tr.
idx + r, tc.
idx +
c) * itw.
get(i, iti);
1807 for (
Index i = 0; i < n; ++i) {
1821 for (
Index p = 0; p < 2; ++p)
1822 for (
Index r = 0; r < 2; ++r)
1880 for (
Index i = 0; i < n; ++i) {
1894 for (
Index s = 0; s < 2; ++s)
1896 for (
Index p = 0; p < 2; ++p)
1897 for (
Index r = 0; r < 2; ++r)
1899 tia +=
a.get(ts.
idx + s,
1962 for (
Index i = 0; i < n; ++i) {
1978 for (
Index s = 0; s < 2; ++s)
1980 for (
Index p = 0; p < 2; ++p)
1981 for (
Index r = 0; r < 2; ++r)
1983 tia +=
a.get(tv.
idx +
v,
2030 for (
Index ir = 0; ir < nr; ++ir) {
2034 for (
Index ic = 0; ic < nc; ++ic) {
2049 itw.
get(ir, ic, iti) = (*r) * (*c);
2090 for (
Index ip = 0; ip < np; ++ip) {
2092 for (
Index ir = 0; ir < nr; ++ir) {
2094 for (
Index ic = 0; ic < nc; ++ic) {
2102 itw.
get(ip, ir, ic, iti) = (*p) * (*r) * (*c);
2147 for (
Index ib = 0; ib < nb; ++ib) {
2149 for (
Index ip = 0; ip < np; ++ip) {
2151 for (
Index ir = 0; ir < nr; ++ir) {
2153 for (
Index ic = 0; ic < nc; ++ic) {
2162 itw.
get(ib, ip, ir, ic, iti) = (*b) * (*p) * (*r) * (*c);
2211 for (
Index is = 0; is < ns; ++is) {
2213 for (
Index ib = 0; ib < nb; ++ib) {
2215 for (
Index ip = 0; ip < np; ++ip) {
2217 for (
Index ir = 0; ir < nr; ++ir) {
2219 for (
Index ic = 0; ic < nc; ++ic) {
2229 itw.
get(is, ib, ip, ir, ic, iti) =
2230 (*s) * (*b) * (*p) * (*r) * (*c);
2283 for (
Index iv = 0; iv < nv; ++iv) {
2285 for (
Index is = 0; is < ns; ++is) {
2287 for (
Index ib = 0; ib < nb; ++ib) {
2289 for (
Index ip = 0; ip < np; ++ip) {
2291 for (
Index ir = 0; ir < nr; ++ir) {
2293 for (
Index ic = 0; ic < nc; ++ic) {
2304 itw.
get(iv, is, ib, ip, ir, ic, iti) =
2305 (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
2353 itw(0, 0,
Range(
joker)).sum(), 1, sum_check_epsilon));
2356 for (
Index ir = 0; ir < nr; ++ir) {
2360 for (
Index ic = 0; ic < nc; ++ic) {
2370 for (
Index r = 0; r < 2; ++r)
2374 tia +=
a.get(tr.
idx + r, tc.
idx +
c) * itw.
get(ir, ic, iti);
2420 itw(0, 0, 0,
Range(
joker)).sum(), 1, sum_check_epsilon));
2423 for (
Index ip = 0; ip < np; ++ip) {
2425 for (
Index ir = 0; ir < nr; ++ir) {
2427 for (
Index ic = 0; ic < nc; ++ic) {
2433 Numeric& tia = ia(ip, ir, ic);
2437 for (
Index p = 0; p < 2; ++p)
2438 for (
Index r = 0; r < 2; ++r)
2443 tia +=
a.get(tp.
idx + p, tr.
idx + r, tc.
idx +
c) *
2444 itw.
get(ip, ir, ic, iti);
2494 itw(0, 0, 0, 0,
Range(
joker)).sum(), 1, sum_check_epsilon));
2497 for (
Index ib = 0; ib < nb; ++ib) {
2499 for (
Index ip = 0; ip < np; ++ip) {
2501 for (
Index ir = 0; ir < nr; ++ir) {
2503 for (
Index ic = 0; ic < nc; ++ic) {
2509 Numeric& tia = ia(ib, ip, ir, ic);
2514 for (
Index p = 0; p < 2; ++p)
2515 for (
Index r = 0; r < 2; ++r)
2518 itw.
get(ib, ip, ir, ic, iti);
2572 itw(0, 0, 0, 0, 0,
Range(
joker)).sum(), 1, sum_check_epsilon));
2575 for (
Index is = 0; is < ns; ++is) {
2577 for (
Index ib = 0; ib < nb; ++ib) {
2579 for (
Index ip = 0; ip < np; ++ip) {
2581 for (
Index ir = 0; ir < nr; ++ir) {
2583 for (
Index ic = 0; ic < nc; ++ic) {
2589 Numeric& tia = ia(is, ib, ip, ir, ic);
2593 for (
Index s = 0; s < 2; ++s)
2595 for (
Index p = 0; p < 2; ++p)
2596 for (
Index r = 0; r < 2; ++r)
2598 tia +=
a.get(ts.
idx + s,
2603 itw.
get(is, ib, ip, ir, ic, iti);
2661 itw(0, 0, 0, 0, 0, 0,
Range(
joker)).sum(), 1, sum_check_epsilon));
2664 for (
Index iv = 0; iv < nv; ++iv) {
2666 for (
Index is = 0; is < ns; ++is) {
2668 for (
Index ib = 0; ib < nb; ++ib) {
2670 for (
Index ip = 0; ip < np; ++ip) {
2672 for (
Index ir = 0; ir < nr; ++ir) {
2674 for (
Index ic = 0; ic < nc; ++ic) {
2680 Numeric& tia = ia(iv, is, ib, ip, ir, ic);
2685 for (
Index s = 0; s < 2; ++s)
2687 for (
Index p = 0; p < 2; ++p)
2688 for (
Index r = 0; r < 2; ++r)
2690 tia +=
a.get(tv.
idx +
v,
2696 itw.
get(iv, is, ib, ip, ir, ic, iti);
2737 dif =
abs(x - xa[0]);
2743 for (
Index i = 0; i < n; i++) {
2744 if ((dift =
abs(x - xa[i])) < dif) {
2755 for (
Index m = 1; m < n; m++) {
2756 for (
Index i = 0; i < n - m; i++) {
2759 w =
c[i + 1] -
d[i];
2767 y_int += (dy_int = (2 * (ns + 1) < (n - m) ?
c[ns + 1] :
d[ns--]));
2805 Index interp_method = 1;
2807 if (interp_method == 1) {
2809 if ((gp.
fd[0] <= 0.5 && gp.
idx > 0) || gp.
idx == N_x - 2) {
2810 xa[0] = x[gp.
idx - 1];
2812 xa[2] = x[gp.
idx + 1];
2814 ya[0] = y[gp.
idx - 1];
2816 ya[2] = y[gp.
idx + 1];
2819 else if ((gp.
fd[0] > 0.5 && gp.
idx < N_x - 2) || gp.
idx == 0) {
2821 xa[1] = x[gp.
idx + 1];
2822 xa[2] = x[gp.
idx + 2];
2825 ya[1] = y[gp.
idx + 1];
2826 ya[2] = y[gp.
idx + 2];
2829 else if (gp.
idx == N_x - 1) {
2842 polint(y_int, dy_int, xa, ya, 3, x_i);
2846 else if (interp_method == 2) {
2849 xa[1] = x[gp.
idx + 1];
2850 xa[2] = x[gp.
idx + 2];
2853 ya[1] = y[gp.
idx + 1];
2854 ya[2] = y[gp.
idx + 2];
2855 }
else if (gp.
idx == N_x - 1) {
2856 xa[0] = x[gp.
idx - 2];
2857 xa[1] = x[gp.
idx - 1];
2860 ya[0] = y[gp.
idx - 2];
2861 ya[1] = y[gp.
idx - 1];
2864 xa[0] = x[gp.
idx - 1];
2866 xa[2] = x[gp.
idx + 1];
2868 ya[0] = y[gp.
idx - 1];
2870 ya[2] = y[gp.
idx + 1];
2874 polint(y_int, dy_int, xa, ya, 3, x_i);
2877 else if (interp_method == 3) {
2880 xa[0] = -x[gp.
idx + 1];
2881 xa[1] = x[gp.
idx + 0];
2882 xa[2] = x[gp.
idx + 1];
2883 xa[3] = x[gp.
idx + 2];
2885 ya[0] = y[gp.
idx + 1];
2886 ya[1] = y[gp.
idx + 0];
2887 ya[2] = y[gp.
idx + 1];
2888 ya[3] = y[gp.
idx + 2];
2889 }
else if (gp.
idx == N_x - 1) {
2890 xa[0] = x[gp.
idx - 1];
2891 xa[1] = x[gp.
idx - 0];
2892 xa[2] = 2 * x[gp.
idx] - x[gp.
idx - 1];
2893 xa[3] = 2 * x[gp.
idx] - x[gp.
idx - 2];
2895 ya[0] = y[gp.
idx - 1];
2896 ya[1] = y[gp.
idx - 0];
2897 ya[2] = y[gp.
idx - 1];
2898 ya[3] = y[gp.
idx - 2];
2899 }
else if (gp.
idx == N_x - 2) {
2900 xa[0] = x[gp.
idx - 2];
2901 xa[1] = x[gp.
idx - 1];
2903 xa[3] = x[gp.
idx + 1];
2905 ya[0] = y[gp.
idx - 2];
2906 ya[1] = y[gp.
idx - 1];
2908 ya[3] = y[gp.
idx + 1];
2910 xa[0] = x[gp.
idx - 1];
2912 xa[2] = x[gp.
idx + 1];
2913 xa[3] = x[gp.
idx + 2];
2915 ya[0] = y[gp.
idx - 1];
2917 ya[2] = y[gp.
idx + 1];
2918 ya[3] = y[gp.
idx + 2];
2921 polint(y_int, dy_int, xa, ya, 4, x_i);
This file contains the definition of Array.
void arts_exit(int status)
This is the exit function of ARTS.
Index nelem() const ARTS_NOEXCEPT
A constant view of a Matrix.
Numeric get(Index r, Index c) const ARTS_NOEXCEPT
Get element implementation without assertions.
A constant view of a Tensor3.
Numeric get(Index p, Index r, Index c) const
Get element implementation without assertions.
A constant view of a Tensor4.
Numeric get(Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
A constant view of a Tensor5.
Numeric get(Index s, Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
A constant view of a Tensor6.
Numeric get(Index v, Index s, Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
A constant view of a Tensor7.
Numeric get(Index l, Index v, Index s, Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
A constant view of a Vector.
Numeric get(Index n) const ARTS_NOEXCEPT
Get element implementation without assertions.
Numeric sum() const ARTS_NOEXCEPT
The sum of all elements of a Vector.
Index nelem() const noexcept
Returns the number of elements.
Numeric & get(Index r, Index c) ARTS_NOEXCEPT
Get element implementation without assertions.
Numeric & get(Index p, Index r, Index c)
Get element implementation without assertions.
Numeric & get(Index b, Index p, Index r, Index c)
Get element implementation without assertions.
Numeric & get(Index s, Index b, Index p, Index r, Index c)
Get element implementation without assertions.
Numeric & get(Index v, Index s, Index b, Index p, Index r, Index c)
Get element implementation without assertions.
Numeric & get(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Get element implementation without assertions.
Numeric & get(Index n) ARTS_NOEXCEPT
Get element implementation without assertions.
#define ARTS_ASSERT(condition,...)
Index gridpos2gridrange(const GridPos &gp, const bool &upwards)
gridpos2gridrange
void gridpos_force_end_fd(GridPos &gp, const Index &n)
gridpos_force_end_fd
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void gridpos_check_fd(GridPos &gp)
gridpos_check_fd
void gridpos_upperend_check(GridPos &gp, const Index &ie)
gridpos_upperend_check
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
#define LOOPIT(x)
Macro for interpolation weight loops.
Numeric interp_poly(ConstVectorView x, ConstVectorView y, const Numeric &x_i, const GridPos &gp)
Polynomial interpolation.
const Numeric FD_TOL
The maximum difference from 1 that we allow for a sum check.
ostream & operator<<(ostream &os, const GridPos &gp)
Output operator for GridPos.
void polint(Numeric &y_int, Numeric &dy_int, ConstVectorView xa, ConstVectorView ya, const Index &n, const Numeric &x)
Polynomial interpolation.
bool is_gridpos_at_index_i(const GridPos &gp, const Index &i, const bool &strict)
is_gridpos_at_index_i
void gp4length1grid(ArrayOfGridPos &gp)
Grid position matching a grid of length 1.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
void gridpos_copy(GridPos &gp_new, const GridPos &gp_old)
gridpos_copy
void gridpos_1to1(ArrayOfGridPos &gp, ConstVectorView grid)
gridpos_1to1
Numeric fractional_gp(const GridPos &gp)
fractional_gp
Header file for interpolation.cc.
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
bool is_decreasing(ConstVectorView x)
Checks if a vector is sorted in reversed order and is strictly decreasing.
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
bool is_same_within_epsilon(const Numeric &a, const Numeric &b, const Numeric &epsilon)
Check, if two numbers agree within a given epsilon.
Header file for logic.cc.
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
NUMERIC Numeric
The type to use for all floating point numbers.
INDEX Index
The type to use for all integer numbers and indices.
Structure to store a grid position.
std::array< Numeric, 2 > fd