ARTS  2.3.1285(git:92a29ea9-dirty)
interpolation.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2012 Stefan Buehler <sbuehler@ltu.se>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  later version.
7 
8  This program is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  GNU General Public License for more details.
12 
13  You should have received a copy of the GNU General Public License
14  along with this program; if not, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA. */
17 
40 #include "interpolation.h"
41 #include <cmath>
42 #include <iostream>
43 #include "array.h"
44 #include "check_input.h"
45 #include "logic.h"
46 
47 // File-global constants:
48 
50 
61 DEBUG_ONLY(const Numeric sum_check_epsilon = 1e-6;)
62 
64 
73 const Numeric FD_TOL = 1.5e-3;
74 
76 
85 #define LOOPIT(x) for (const Numeric* x = &t##x.fd[1]; x >= &t##x.fd[0]; --x)
86 
88 
96 ostream& operator<<(ostream& os, const GridPos& gp) {
97  os << gp.idx << " " << gp.fd[0] << " " << gp.fd[1] << "\n";
98  return os;
99 }
100 
102 
157  ConstVectorView old_grid,
158  ConstVectorView new_grid,
159  const Numeric& extpolfac) {
160  const Index n_old = old_grid.nelem();
161  const Index n_new = new_grid.nelem();
162 
163  // Assert that gp has the right size:
164  assert(is_size(gp, n_new));
165 
166  // Assert, that the old grid has more than one element
167  assert(1 < n_old);
168 
169  // This function hast two parts, depending on whether old_grid is
170  // sorted in ascending or descending order. Maybe that's not too
171  // elegant, but it's the most efficient, because in this way there
172  // is no additional runtime overhead to handle both cases.
173 
174  // We use only the first two elements to decide how the grid is
175  // sorted. (The rest of the grid is checked later by an assert.)
176  // If both values are the same, we still assume the grid is
177  // ascending. However, this will lead to an assertion fail later on,
178  // because the grid has to be strictly sorted.
179  bool ascending = (old_grid[0] <= old_grid[1]);
180 
181  if (ascending) {
182  // So, old_grid should always be sorted in strictly ascending order.
183  // (No duplicate values.)
184  // Assert that this is so. This may depend on user input, however,
185  // inside this elementary function is not the place to check for
186  // that. There must be runtime checks on higher levels to insure
187  // that all grids are sorted. The assertion here is just as a last
188  // safety check.
189  assert(is_increasing(old_grid));
190 
191  // Limits of extrapolation.
192  const Numeric og_min =
193  old_grid[0] - extpolfac * (old_grid[1] - old_grid[0]);
194  const Numeric og_max =
195  old_grid[n_old - 1] +
196  extpolfac * (old_grid[n_old - 1] - old_grid[n_old - 2]);
197 
198  // We will make no firm assumptions about the new grid. But the case
199  // that we have in mind is that it is either also sorted, or at
200  // least partially sorted, for example like this:
201  // 5 4 3 2 2.5 3 4
202  // This kind of sequence should be typical if we interpolate
203  // atmospheric fields along a limb propagation path.
204 
205  // Let's get some idea where the first point in the new grid is,
206  // relative to the old grid. We use linear interpolation between the
207  // maximum and the minimum of the old grid for this. (Taking
208  // into account the small allowed extrapolation.)
209  Numeric frac = (new_grid[0] - og_min) / (og_max - og_min);
210 
211  // We are not checking if the value of frac is reasonable,
212  // because there is another assertion below to catch the
213  // consequences.
214 
215  // Initialize current_position:
216  Index current_position = (Index)rint(frac * (Numeric)(n_old - 2));
217 
218  // The above statement should satisfy
219  // 0 <= current_position <= n_old-2
220  // Assert that this is indeed the case.
221  // cout << "frac = " << frac << "\n";
222  // cout << "current_position = " << current_position << "\n";
223  assert(0 <= current_position);
224  assert(current_position <= n_old - 2);
225 
226  // The variables lower and upper are used to remember the value of
227  // the old grid at current_position and one above current_position:
228  Numeric lower = old_grid[current_position];
229  Numeric upper = old_grid[current_position + 1];
230 
231  // Loop over all points in the new grid:
232  for (Index i_new = 0; i_new < n_new; ++i_new) {
233  // Get a reference to the current element of gp:
234  GridPos& tgp = gp[i_new];
235  // And on the current value of the new grid:
236  const Numeric tng = new_grid[i_new];
237 
238  // Verify that the new grid is within the limits of
239  // extrapolation that we have defined by og_min and og_max:
240  assert(og_min <= tng);
241  assert(tng <= og_max);
242 
243  // cout << "lower / tng / upper = "
244  // << lower << " / "
245  // << tng << " / "
246  // << upper << "\n";
247 
248  // Is current_position too high?
249  // (The current_position>0 condition is there so that the position
250  // stays 0 for extrapolation.)
251  if (tng < lower && current_position > 0) {
252  do {
253  --current_position;
254  lower = old_grid[current_position];
255  } while (tng < lower && current_position > 0);
256 
257  upper = old_grid[current_position + 1];
258 
259  tgp.idx = current_position;
260  tgp.fd[0] = (tng - lower) / (upper - lower);
261  tgp.fd[1] = 1.0 - tgp.fd[0];
262  } else {
263  // Is it too low?
264  // (The current_position<n_old condition is there so
265  // that uppers stays n_old-1 for extrapolation.)
266  if (tng >= upper && current_position < n_old - 2) {
267  do {
268  ++current_position;
269  upper = old_grid[current_position + 1];
270  } while (tng >= upper && current_position < n_old - 2);
271 
272  lower = old_grid[current_position];
273 
274  tgp.idx = current_position;
275  tgp.fd[0] = (tng - lower) / (upper - lower);
276  tgp.fd[1] = 1.0 - tgp.fd[0];
277  } else {
278  // None of the other two conditions were true. That means:
279  // lower <= tng < upper. The current_position is still
280  // valid.
281 
282  // SAB 2010-04-28: Note that if a new grid point is exactly on
283  // top of an old grid point, then you are now guaranteed to get
284  // fd[0] = 0 and fd[1] = 1. (Except at the upper grid end.)
285 
286  tgp.idx = current_position;
287  tgp.fd[0] = (tng - lower) / (upper - lower);
288  tgp.fd[1] = 1.0 - tgp.fd[0];
289  }
290  }
291 
292  // cout << tgp.idx << " " << tgp.fd[0] << " " << tgp.fd[1] << endl;
293 
294  // Safety check to ensure the above promise:
295  assert(old_grid[tgp.idx] <= tng || tgp.idx == 0);
296  }
297  } else // if (ascending)
298  {
299  // Now we are in the "descending old grid" part. We do exactly
300  // the same as in the other part, just accounting for the
301  // different order of things. Comments here refer only to
302  // interesting differences from the ascending case. See that
303  // case for more general comments.
304 
305  // This time ensure strictly descending order:
306  assert(is_decreasing(old_grid));
307 
308  // The max is now the first point, the min the last point!
309  // I think the sign is right here...
310  const Numeric og_max =
311  old_grid[0] - extpolfac * (old_grid[1] - old_grid[0]);
312  const Numeric og_min =
313  old_grid[n_old - 1] +
314  extpolfac * (old_grid[n_old - 1] - old_grid[n_old - 2]);
315 
316  // We have to take 1- here, because we are starting from the
317  // high end.
318  Numeric frac = 1 - (new_grid[0] - og_min) / (og_max - og_min);
319 
320  // We are not checking if the value of frac is reasonable,
321  // because there is another assertion below to catch the
322  // consequences.
323 
324  Index current_position = (Index)rint(frac * (Numeric)(n_old - 2));
325 
326  // The above statement should satisfy
327  // 0 <= current_position <= n_old-2
328  // Assert that this is indeed the case.
329  // cout << "frac = " << frac << "\n";
330  // cout << "current_position = " << current_position << "\n";
331  assert(0 <= current_position);
332  assert(current_position <= n_old - 2);
333 
334  // Note, that old_grid[lower] has a higher numerical value than
335  // old_grid[upper]!
336  Numeric lower = old_grid[current_position];
337  Numeric upper = old_grid[current_position + 1];
338 
339  for (Index i_new = 0; i_new < n_new; ++i_new) {
340  GridPos& tgp = gp[i_new];
341  const Numeric tng = new_grid[i_new];
342 
343  // Verify that the new grid is within the limits of
344  // extrapolation that we have defined by og_min and og_max:
345  assert(og_min <= tng);
346  assert(tng <= og_max);
347 
348  // cout << "lower / tng / upper = "
349  // << lower << " / "
350  // << tng << " / "
351  // << upper << "\n";
352 
353  // Is current_position too high? (Sign of comparison changed
354  // compared to ascending case!)
355  if (tng > lower && current_position > 0) {
356  do {
357  --current_position;
358  lower = old_grid[current_position];
359  } while (tng > lower && current_position > 0);
360 
361  upper = old_grid[current_position + 1];
362 
363  tgp.idx = current_position;
364  tgp.fd[0] = (tng - lower) / (upper - lower);
365  tgp.fd[1] = 1.0 - tgp.fd[0];
366  } else {
367  // Is it too low? (Sign of comparison changed
368  // compared to ascending case!)
369  if (tng <= upper && current_position < n_old - 2) {
370  do {
371  ++current_position;
372  upper = old_grid[current_position + 1];
373  } while (tng <= upper && current_position < n_old - 2);
374 
375  lower = old_grid[current_position];
376 
377  tgp.idx = current_position;
378  tgp.fd[0] = (tng - lower) / (upper - lower);
379  tgp.fd[1] = 1.0 - tgp.fd[0];
380  } else {
381  // None of the other two conditions were true. That means:
382  // lower >= tng > upper. The current_position is still
383  // valid. (Note that upper and lower have switched
384  // place compared to the ascending case.)
385 
386  // SAB 2010-04-28: Note that if a new grid point is exactly on
387  // top of an old grid point, then you are now guaranteed to get
388  // fd[0] = 0 and fd[1] = 1. (Except at the upper grid end.)
389 
390  tgp.idx = current_position;
391  tgp.fd[0] = (tng - lower) / (upper - lower);
392  tgp.fd[1] = 1.0 - tgp.fd[0];
393  }
394  }
395 
396  // Safety check to ensure the above promise:
397  assert(old_grid[tgp.idx] >= tng || tgp.idx == 0);
398  }
399  }
400 }
401 
403 
424 void gridpos(GridPos& gp,
425  ConstVectorView old_grid,
426  const Numeric& new_grid,
427  const Numeric& extpolfac) {
428  ArrayOfGridPos agp(1);
429  Vector v(1, new_grid);
430  gridpos(agp, old_grid, v, extpolfac);
431  gridpos_copy(gp, agp[0]);
432 }
433 
435 
450  const Index n = grid.nelem();
451  gp.resize(n);
452 
453  for (Index i = 0; i < n - 1; i++) {
454  gp[i].idx = i;
455  gp[i].fd[0] = 0;
456  gp[i].fd[1] = 1;
457  }
458  gp[n - 1].idx = n - 2;
459  gp[n - 1].fd[0] = 1;
460  gp[n - 1].fd[1] = 0;
461 }
462 
464 
473 void gridpos_copy(GridPos& gp_new, const GridPos& gp_old) {
474  gp_new.idx = gp_old.idx;
475  gp_new.fd[0] = gp_old.fd[0];
476  gp_new.fd[1] = gp_old.fd[1];
477 }
478 
480 
493  return (Numeric(gp.idx) + gp.fd[0]);
494 }
495 
497 
510  // Catch values that "must" be wrong
511  assert(gp.fd[0] > -FD_TOL);
512  assert(gp.fd[0] < 1.0 + FD_TOL);
513  assert(gp.fd[1] > -FD_TOL);
514  assert(gp.fd[1] < 1.0 + FD_TOL);
515 
516  if (gp.fd[0] < 0.0) {
517  gp.fd[0] = 0.0;
518  gp.fd[1] = 1.0;
519  } else if (gp.fd[0] > 1.0) {
520  gp.fd[0] = 1.0;
521  gp.fd[1] = 0.0;
522  }
523 
524  if (gp.fd[1] < 0.0) {
525  gp.fd[1] = 0.0;
526  gp.fd[0] = 1.0;
527  } else if (gp.fd[1] > 1.0) {
528  gp.fd[1] = 1.0;
529  gp.fd[0] = 0.0;
530  }
531 }
532 
534 
555 void gridpos_force_end_fd(GridPos& gp, const Index& n) {
556  assert(gp.idx >= 0);
557 
558  // If fd=1, shift to grid index above
559  if (gp.fd[0] > 0.5) {
560  gp.idx += 1;
561  }
562  gp.fd[0] = 0;
563  gp.fd[1] = 1;
564 
565  assert(gp.idx < n);
566 
567  // End of complete grid range must be handled separately
568  if (gp.idx == n - 1) {
569  gp.idx -= 1;
570  gp.fd[0] = 1;
571  gp.fd[1] = 0;
572  }
573 }
574 
576 
589 void gridpos_upperend_check(GridPos& gp, const Index& ie) {
590  if (gp.idx == ie) {
591  assert(gp.fd[0] < 0.005); // To capture obviously bad cases
592  gp.idx -= 1;
593  gp.fd[0] = 1.0;
594  gp.fd[1] = 0.0;
595  }
596 }
597 
599 
613  for (Index i = 0; i < gp.nelem(); i++) {
614  if (gp[i].idx == ie) {
615  assert(gp[i].fd[0] < 0.005); // To capture obviously bad cases
616  gp[i].idx -= 1;
617  gp[i].fd[0] = 1.0;
618  gp[i].fd[1] = 0.0;
619  }
620  }
621 }
622 
624 
637  for (Index i = 0; i < gp.nelem(); i++) {
638  gp[i].idx = 0;
639  gp[i].fd[0] = 0;
640  gp[i].fd[1] = 1;
641  }
642 }
643 
645 
658  const Index& i,
659  const bool& strict) {
660  if (strict) {
661  // Assume that gridpos_force_end_fd has been used. The expression 0==0
662  // should be safer than 1==1.
663  if (gp.idx == i && gp.fd[0] == 0) {
664  return true;
665  } else if (gp.idx == i - 1 && gp.fd[1] == 0) {
666  return true;
667  }
668  } else {
669  if (gp.idx == i && gp.fd[0] < FD_TOL) {
670  return true;
671  } else if (gp.idx == i - 1 && gp.fd[1] < FD_TOL) {
672  return true;
673  }
674  }
675  return false;
676 }
677 
679 
697 Index gridpos2gridrange(const GridPos& gp, const bool& upwards) {
698  assert(gp.fd[0] >= 0);
699  assert(gp.fd[1] >= 0);
700 
701  // Not at a grid point
702  if (gp.fd[0] > 0 && gp.fd[1] > 0) {
703  return gp.idx;
704  }
705 
706  // Fractional distance 0
707  else if (gp.fd[0] == 0) {
708  if (upwards)
709  return gp.idx;
710  else
711  return gp.idx - 1;
712  }
713 
714  // Fractional distance 1
715  else {
716  if (upwards)
717  return gp.idx + 1;
718  else
719  return gp.idx;
720  }
721 }
722 
724 // Red Interpolation
726 
728 
741 void interpweights(VectorView itw, const GridPos& tc) {
742  assert(is_size(itw, 2)); // We must store 2 interpolation
743  // weights.
744 
745  // Interpolation weights are stored in this order (l=lower
746  // u=upper, c=column):
747  // 1. l-c
748  // 2. u-c
749 
750  Index iti = 0;
751 
752  // We could use a straight-forward for loop here:
753  //
754  // for ( Index c=1; c>=0; --c )
755  // {
756  // ti[iti] = tc.fd[c];
757  // ++iti;
758  // }
759  //
760  // However, it is slightly faster to use pointers (I tried it,
761  // the speed gain is about 20%). So we should write something
762  // like:
763  //
764  // for ( const Numeric* c=&tc.fd[1]; c>=&tc.fd[0]; --c )
765  // {
766  // ti[iti] = *c;
767  // ++iti;
768  // }
769  //
770  // For higher dimensions we have to nest these loops. To avoid
771  // typos and safe typing, I use the LOOPIT macro, which expands
772  // to the for loop above. Note: NO SEMICOLON AFTER THE LOOPIT
773  // COMMAND!
774 
775  LOOPIT(c) {
776  itw.get(iti) = *c;
777  ++iti;
778  }
779 }
780 
782 
796 void interpweights(VectorView itw, const GridPos& tr, const GridPos& tc) {
797  assert(is_size(itw, 4)); // We must store 4 interpolation
798  // weights.
799  Index iti = 0;
800 
801  LOOPIT(r)
802  LOOPIT(c) {
803  itw.get(iti) = (*r) * (*c);
804  ++iti;
805  }
806 }
807 
809 
825  const GridPos& tp,
826  const GridPos& tr,
827  const GridPos& tc) {
828  assert(is_size(itw, 8)); // We must store 8 interpolation
829  // weights.
830  Index iti = 0;
831 
832  LOOPIT(p)
833  LOOPIT(r)
834  LOOPIT(c) {
835  itw.get(iti) = (*p) * (*r) * (*c);
836  ++iti;
837  }
838 }
839 
841 
858  const GridPos& tb,
859  const GridPos& tp,
860  const GridPos& tr,
861  const GridPos& tc) {
862  assert(is_size(itw, 16)); // We must store 16 interpolation
863  // weights.
864  Index iti = 0;
865 
866  LOOPIT(b)
867  LOOPIT(p)
868  LOOPIT(r)
869  LOOPIT(c) {
870  itw.get(iti) = (*b) * (*p) * (*r) * (*c);
871  ++iti;
872  }
873 }
874 
876 
894  const GridPos& ts,
895  const GridPos& tb,
896  const GridPos& tp,
897  const GridPos& tr,
898  const GridPos& tc) {
899  assert(is_size(itw, 32)); // We must store 32 interpolation
900  // weights.
901  Index iti = 0;
902 
903  LOOPIT(s)
904  LOOPIT(b)
905  LOOPIT(p)
906  LOOPIT(r)
907  LOOPIT(c) {
908  itw.get(iti) = (*s) * (*b) * (*p) * (*r) * (*c);
909  ++iti;
910  }
911 }
912 
914 
933  const GridPos& tv,
934  const GridPos& ts,
935  const GridPos& tb,
936  const GridPos& tp,
937  const GridPos& tr,
938  const GridPos& tc) {
939  assert(is_size(itw, 64)); // We must store 64 interpolation
940  // weights.
941  Index iti = 0;
942 
943  LOOPIT(v)
944  LOOPIT(s)
945  LOOPIT(b)
946  LOOPIT(p)
947  LOOPIT(r)
948  LOOPIT(c) {
949  itw.get(iti) = (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
950  ++iti;
951  }
952 }
953 
955 
971  assert(is_size(itw, 2)); // We need 2 interpolation
972  // weights.
973 
974  // Check that interpolation weights are valid. The sum of all
975  // weights (last dimension) must always be approximately one.
976  assert(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
977 
978  // To store interpolated value:
979  Numeric tia = 0;
980 
981  Index iti = 0;
982  for (Index c = 0; c < 2; ++c) {
983  tia += a.get(tc.idx + c) * itw.get(iti);
984  ++iti;
985  }
986 
987  return tia;
988 }
989 
991 
1008  ConstMatrixView a,
1009  const GridPos& tr,
1010  const GridPos& tc) {
1011  assert(is_size(itw, 4)); // We need 4 interpolation
1012  // weights.
1013 
1014  // Check that interpolation weights are valid. The sum of all
1015  // weights (last dimension) must always be approximately one.
1016  assert(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
1017 
1018  // To store interpolated value:
1019  Numeric tia = 0;
1020 
1021  Index iti = 0;
1022  for (Index r = 0; r < 2; ++r)
1023  for (Index c = 0; c < 2; ++c) {
1024  tia += a.get(tr.idx + r, tc.idx + c) * itw.get(iti);
1025  ++iti;
1026  }
1027 
1028  return tia;
1029 }
1030 
1032 
1050  ConstTensor3View a,
1051  const GridPos& tp,
1052  const GridPos& tr,
1053  const GridPos& tc) {
1054  assert(is_size(itw, 8)); // We need 8 interpolation
1055  // weights.
1056 
1057  // Check that interpolation weights are valid. The sum of all
1058  // weights (last dimension) must always be approximately one.
1059  assert(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
1060 
1061  // To store interpolated value:
1062  Numeric tia = 0;
1063 
1064  Index iti = 0;
1065  for (Index p = 0; p < 2; ++p)
1066  for (Index r = 0; r < 2; ++r)
1067  for (Index c = 0; c < 2; ++c) {
1068  tia += a.get(tp.idx + p, tr.idx + r, tc.idx + c) * itw.get(iti);
1069  ++iti;
1070  }
1071 
1072  return tia;
1073 }
1074 
1076 
1095  ConstTensor4View a,
1096  const GridPos& tb,
1097  const GridPos& tp,
1098  const GridPos& tr,
1099  const GridPos& tc) {
1100  assert(is_size(itw, 16)); // We need 16 interpolation
1101  // weights.
1102 
1103  // Check that interpolation weights are valid. The sum of all
1104  // weights (last dimension) must always be approximately one.
1105  assert(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
1106 
1107  // To store interpolated value:
1108  Numeric tia = 0;
1109 
1110  Index iti = 0;
1111  for (Index b = 0; b < 2; ++b)
1112  for (Index p = 0; p < 2; ++p)
1113  for (Index r = 0; r < 2; ++r)
1114  for (Index c = 0; c < 2; ++c) {
1115  tia += a.get(tb.idx + b, tp.idx + p, tr.idx + r, tc.idx + c) *
1116  itw.get(iti);
1117  ++iti;
1118  }
1119 
1120  return tia;
1121 }
1122 
1124 
1144  ConstTensor5View a,
1145  const GridPos& ts,
1146  const GridPos& tb,
1147  const GridPos& tp,
1148  const GridPos& tr,
1149  const GridPos& tc) {
1150  assert(is_size(itw, 32)); // We need 32 interpolation
1151  // weights.
1152 
1153  // Check that interpolation weights are valid. The sum of all
1154  // weights (last dimension) must always be approximately one.
1155  assert(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
1156 
1157  // To store interpolated value:
1158  Numeric tia = 0;
1159 
1160  Index iti = 0;
1161  for (Index s = 0; s < 2; ++s)
1162  for (Index b = 0; b < 2; ++b)
1163  for (Index p = 0; p < 2; ++p)
1164  for (Index r = 0; r < 2; ++r)
1165  for (Index c = 0; c < 2; ++c) {
1166  tia += a.get(ts.idx + s,
1167  tb.idx + b,
1168  tp.idx + p,
1169  tr.idx + r,
1170  tc.idx + c) *
1171  itw.get(iti);
1172  ++iti;
1173  }
1174 
1175  return tia;
1176 }
1177 
1179 
1200  ConstTensor6View a,
1201  const GridPos& tv,
1202  const GridPos& ts,
1203  const GridPos& tb,
1204  const GridPos& tp,
1205  const GridPos& tr,
1206  const GridPos& tc) {
1207  assert(is_size(itw, 64)); // We need 64 interpolation
1208  // weights.
1209 
1210  // Check that interpolation weights are valid. The sum of all
1211  // weights (last dimension) must always be approximately one.
1212  assert(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
1213 
1214  // To store interpolated value:
1215  Numeric tia = 0;
1216 
1217  Index iti = 0;
1218  for (Index v = 0; v < 2; ++v)
1219  for (Index s = 0; s < 2; ++s)
1220  for (Index b = 0; b < 2; ++b)
1221  for (Index p = 0; p < 2; ++p)
1222  for (Index r = 0; r < 2; ++r)
1223  for (Index c = 0; c < 2; ++c) {
1224  tia += a.get(tv.idx + v,
1225  ts.idx + s,
1226  tb.idx + b,
1227  tp.idx + p,
1228  tr.idx + r,
1229  tc.idx + c) *
1230  itw.get(iti);
1231  ++iti;
1232  }
1233 
1234  return tia;
1235 }
1236 
1238 // Blue interpolation
1240 
1242 
1257 void interpweights(MatrixView itw, const ArrayOfGridPos& cgp) {
1258  Index n = cgp.nelem();
1259  assert(is_size(itw, n, 2)); // We must store 2 interpolation
1260  // weights for each position.
1261 
1262  // We have to loop all the points in the sequence:
1263  for (Index i = 0; i < n; ++i) {
1264  // Current grid positions:
1265  const GridPos& tc = cgp[i];
1266 
1267  // Interpolation weights are stored in this order (l=lower
1268  // u=upper, c=column):
1269  // 1. l-c
1270  // 2. u-c
1271 
1272  Index iti = 0;
1273 
1274  // We could use a straight-forward for loop here:
1275  //
1276  // for ( Index c=1; c>=0; --c )
1277  // {
1278  // ti[iti] = tc.fd[c];
1279  // ++iti;
1280  // }
1281  //
1282  // However, it is slightly faster to use pointers (I tried it,
1283  // the speed gain is about 20%). So we should write something
1284  // like:
1285  //
1286  // for ( const Numeric* c=&tc.fd[1]; c>=&tc.fd[0]; --c )
1287  // {
1288  // ti[iti] = *c;
1289  // ++iti;
1290  // }
1291  //
1292  // For higher dimensions we have to nest these loops. To avoid
1293  // typos and safe typing, I use the LOOPIT macro, which expands
1294  // to the for loop above. Note: NO SEMICOLON AFTER THE LOOPIT
1295  // COMMAND!
1296 
1297  LOOPIT(c) {
1298  itw.get(i, iti) = *c;
1299  ++iti;
1300  }
1301  }
1302 }
1303 
1305 
1327  const ArrayOfGridPos& rgp,
1328  const ArrayOfGridPos& cgp) {
1329  Index n = cgp.nelem();
1330  assert(is_size(rgp, n)); // rgp must have same size as cgp.
1331  assert(is_size(itw, n, 4)); // We must store 4 interpolation
1332  // weights for each position.
1333 
1334  // We have to loop all the points in the sequence:
1335  for (Index i = 0; i < n; ++i) {
1336  // Current grid positions:
1337  const GridPos& tr = rgp[i];
1338  const GridPos& tc = cgp[i];
1339 
1340  // Interpolation weights are stored in this order (l=lower
1341  // u=upper, r=row, c=column):
1342  // 1. l-r l-c
1343  // 2. l-r u-c
1344  // 3. u-r l-c
1345  // 4. u-r u-c
1346 
1347  Index iti = 0;
1348 
1349  LOOPIT(r)
1350  LOOPIT(c) {
1351  itw.get(i, iti) = (*r) * (*c);
1352  ++iti;
1353  }
1354  }
1355 }
1356 
1358 
1381  const ArrayOfGridPos& pgp,
1382  const ArrayOfGridPos& rgp,
1383  const ArrayOfGridPos& cgp) {
1384  Index n = cgp.nelem();
1385  assert(is_size(pgp, n)); // pgp must have same size as cgp.
1386  assert(is_size(rgp, n)); // rgp must have same size as cgp.
1387  assert(is_size(itw, n, 8)); // We must store 8 interpolation
1388  // weights for each position.
1389 
1390  // We have to loop all the points in the sequence:
1391  for (Index i = 0; i < n; ++i) {
1392  // Current grid positions:
1393  const GridPos& tp = pgp[i];
1394  const GridPos& tr = rgp[i];
1395  const GridPos& tc = cgp[i];
1396 
1397  Index iti = 0;
1398  LOOPIT(p)
1399  LOOPIT(r)
1400  LOOPIT(c) {
1401  itw.get(i, iti) = (*p) * (*r) * (*c);
1402  ++iti;
1403  }
1404  }
1405 }
1406 
1408 
1432  const ArrayOfGridPos& bgp,
1433  const ArrayOfGridPos& pgp,
1434  const ArrayOfGridPos& rgp,
1435  const ArrayOfGridPos& cgp) {
1436  Index n = cgp.nelem();
1437  assert(is_size(bgp, n)); // bgp must have same size as cgp.
1438  assert(is_size(pgp, n)); // pgp must have same size as cgp.
1439  assert(is_size(rgp, n)); // rgp must have same size as cgp.
1440  assert(is_size(itw, n, 16)); // We must store 16 interpolation
1441  // weights for each position.
1442 
1443  // We have to loop all the points in the sequence:
1444  for (Index i = 0; i < n; ++i) {
1445  // Current grid positions:
1446  const GridPos& tb = bgp[i];
1447  const GridPos& tp = pgp[i];
1448  const GridPos& tr = rgp[i];
1449  const GridPos& tc = cgp[i];
1450 
1451  Index iti = 0;
1452  LOOPIT(b)
1453  LOOPIT(p)
1454  LOOPIT(r)
1455  LOOPIT(c) {
1456  itw.get(i, iti) = (*b) * (*p) * (*r) * (*c);
1457  ++iti;
1458  }
1459  }
1460 }
1461 
1463 
1488  const ArrayOfGridPos& sgp,
1489  const ArrayOfGridPos& bgp,
1490  const ArrayOfGridPos& pgp,
1491  const ArrayOfGridPos& rgp,
1492  const ArrayOfGridPos& cgp) {
1493  Index n = cgp.nelem();
1494  assert(is_size(sgp, n)); // sgp must have same size as cgp.
1495  assert(is_size(bgp, n)); // bgp must have same size as cgp.
1496  assert(is_size(pgp, n)); // pgp must have same size as cgp.
1497  assert(is_size(rgp, n)); // rgp must have same size as cgp.
1498  assert(is_size(itw, n, 32)); // We must store 32 interpolation
1499  // weights for each position.
1500 
1501  // We have to loop all the points in the sequence:
1502  for (Index i = 0; i < n; ++i) {
1503  // Current grid positions:
1504  const GridPos& ts = sgp[i];
1505  const GridPos& tb = bgp[i];
1506  const GridPos& tp = pgp[i];
1507  const GridPos& tr = rgp[i];
1508  const GridPos& tc = cgp[i];
1509 
1510  Index iti = 0;
1511  LOOPIT(s)
1512  LOOPIT(b)
1513  LOOPIT(p)
1514  LOOPIT(r)
1515  LOOPIT(c) {
1516  itw.get(i, iti) = (*s) * (*b) * (*p) * (*r) * (*c);
1517  ++iti;
1518  }
1519  }
1520 }
1521 
1523 
1549  const ArrayOfGridPos& vgp,
1550  const ArrayOfGridPos& sgp,
1551  const ArrayOfGridPos& bgp,
1552  const ArrayOfGridPos& pgp,
1553  const ArrayOfGridPos& rgp,
1554  const ArrayOfGridPos& cgp) {
1555  Index n = cgp.nelem();
1556  assert(is_size(vgp, n)); // vgp must have same size as cgp.
1557  assert(is_size(sgp, n)); // sgp must have same size as cgp.
1558  assert(is_size(bgp, n)); // bgp must have same size as cgp.
1559  assert(is_size(pgp, n)); // pgp must have same size as cgp.
1560  assert(is_size(rgp, n)); // rgp must have same size as cgp.
1561  assert(is_size(itw, n, 64)); // We must store 64 interpolation
1562  // weights for each position.
1563 
1564  // We have to loop all the points in the sequence:
1565  for (Index i = 0; i < n; ++i) {
1566  // Current grid positions:
1567  const GridPos& tv = vgp[i];
1568  const GridPos& ts = sgp[i];
1569  const GridPos& tb = bgp[i];
1570  const GridPos& tp = pgp[i];
1571  const GridPos& tr = rgp[i];
1572  const GridPos& tc = cgp[i];
1573 
1574  Index iti = 0;
1575  LOOPIT(v)
1576  LOOPIT(s)
1577  LOOPIT(b)
1578  LOOPIT(p)
1579  LOOPIT(r)
1580  LOOPIT(c) {
1581  itw.get(i, iti) = (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
1582  ++iti;
1583  }
1584  }
1585 }
1586 
1588 
1605  ConstMatrixView itw,
1606  ConstVectorView a,
1607  const ArrayOfGridPos& cgp) {
1608  Index n = cgp.nelem();
1609  assert(is_size(ia, n)); // ia must have same size as cgp.
1610  assert(is_size(itw, n, 2)); // We need 2 interpolation
1611  // weights for each position.
1612 
1613  // Check that interpolation weights are valid. The sum of all
1614  // weights (last dimension) must always be approximately one. We
1615  // only check the first element.
1616  assert(
1617  is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
1618 
1619  // We have to loop all the points in the sequence:
1620  for (Index i = 0; i < n; ++i) {
1621  // Current grid positions:
1622  const GridPos& tc = cgp[i];
1623 
1624  // Get handle to current element of output vector and initialize
1625  // it to zero:
1626  Numeric& tia = ia[i];
1627  tia = 0;
1628 
1629  Index iti = 0;
1630  for (Index c = 0; c < 2; ++c) {
1631  assert(tc.idx + c < a.nelem()); // Temporary !?
1632  tia += a.get(tc.idx + c) * itw.get(i, iti);
1633  ++iti;
1634  }
1635  }
1636 }
1637 
1639 
1661  ConstMatrixView itw,
1662  ConstMatrixView a,
1663  const ArrayOfGridPos& rgp,
1664  const ArrayOfGridPos& cgp) {
1665  Index n = cgp.nelem();
1666  assert(is_size(ia, n)); // ia must have same size as cgp.
1667  assert(is_size(rgp, n)); // rgp must have same size as cgp.
1668  assert(is_size(itw, n, 4)); // We need 4 interpolation
1669  // weights for each position.
1670 
1671  // Check that interpolation weights are valid. The sum of all
1672  // weights (last dimension) must always be approximately one. We
1673  // only check the first element.
1674  assert(
1675  is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
1676 
1677  // We have to loop all the points in the sequence:
1678  for (Index i = 0; i < n; ++i) {
1679  // Current grid positions:
1680  const GridPos& tr = rgp[i];
1681  const GridPos& tc = cgp[i];
1682 
1683  // Get handle to current element of output vector and initialize
1684  // it to zero:
1685  Numeric& tia = ia[i];
1686  tia = 0;
1687 
1688  Index iti = 0;
1689  for (Index r = 0; r < 2; ++r)
1690  for (Index c = 0; c < 2; ++c) {
1691  tia += a.get(tr.idx + r, tc.idx + c) * itw.get(i, iti);
1692  ++iti;
1693  }
1694  }
1695 }
1696 
1698 
1721  ConstMatrixView itw,
1722  ConstTensor3View a,
1723  const ArrayOfGridPos& pgp,
1724  const ArrayOfGridPos& rgp,
1725  const ArrayOfGridPos& cgp) {
1726  Index n = cgp.nelem();
1727  assert(is_size(ia, n)); // ia must have same size as cgp.
1728  assert(is_size(pgp, n)); // pgp must have same size as cgp.
1729  assert(is_size(rgp, n)); // rgp must have same size as cgp.
1730  assert(is_size(itw, n, 8)); // We need 8 interpolation
1731  // weights for each position.
1732 
1733  // Check that interpolation weights are valid. The sum of all
1734  // weights (last dimension) must always be approximately one. We
1735  // only check the first element.
1736  assert(
1737  is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
1738 
1739  // We have to loop all the points in the sequence:
1740  for (Index i = 0; i < n; ++i) {
1741  // Current grid positions:
1742  const GridPos& tp = pgp[i];
1743  const GridPos& tr = rgp[i];
1744  const GridPos& tc = cgp[i];
1745 
1746  // Get handle to current element of output vector and initialize
1747  // it to zero:
1748  Numeric& tia = ia[i];
1749  tia = 0;
1750 
1751  Index iti = 0;
1752  for (Index p = 0; p < 2; ++p)
1753  for (Index r = 0; r < 2; ++r)
1754  for (Index c = 0; c < 2; ++c) {
1755  tia += a.get(tp.idx + p, tr.idx + r, tc.idx + c) * itw.get(i, iti);
1756  ++iti;
1757  }
1758  }
1759 }
1760 
1762 
1786  ConstMatrixView itw,
1787  ConstTensor4View a,
1788  const ArrayOfGridPos& bgp,
1789  const ArrayOfGridPos& pgp,
1790  const ArrayOfGridPos& rgp,
1791  const ArrayOfGridPos& cgp) {
1792  Index n = cgp.nelem();
1793  assert(is_size(ia, n)); // ia must have same size as cgp.
1794  assert(is_size(bgp, n)); // bgp must have same size as cgp.
1795  assert(is_size(pgp, n)); // pgp must have same size as cgp.
1796  assert(is_size(rgp, n)); // rgp must have same size as cgp.
1797  assert(is_size(itw, n, 16)); // We need 16 interpolation
1798  // weights for each position.
1799 
1800  // Check that interpolation weights are valid. The sum of all
1801  // weights (last dimension) must always be approximately one. We
1802  // only check the first element.
1803  assert(
1804  is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
1805 
1806  // We have to loop all the points in the sequence:
1807  for (Index i = 0; i < n; ++i) {
1808  // Current grid positions:
1809  const GridPos& tb = bgp[i];
1810  const GridPos& tp = pgp[i];
1811  const GridPos& tr = rgp[i];
1812  const GridPos& tc = cgp[i];
1813 
1814  // Get handle to current element of output vector and initialize
1815  // it to zero:
1816  Numeric& tia = ia[i];
1817  tia = 0;
1818 
1819  Index iti = 0;
1820  for (Index b = 0; b < 2; ++b)
1821  for (Index p = 0; p < 2; ++p)
1822  for (Index r = 0; r < 2; ++r)
1823  for (Index c = 0; c < 2; ++c) {
1824  tia += a.get(tb.idx + b, tp.idx + p, tr.idx + r, tc.idx + c) *
1825  itw.get(i, iti);
1826  ++iti;
1827  }
1828  }
1829 }
1830 
1832 
1857  ConstMatrixView itw,
1858  ConstTensor5View a,
1859  const ArrayOfGridPos& sgp,
1860  const ArrayOfGridPos& bgp,
1861  const ArrayOfGridPos& pgp,
1862  const ArrayOfGridPos& rgp,
1863  const ArrayOfGridPos& cgp) {
1864  Index n = cgp.nelem();
1865  assert(is_size(ia, n)); // ia must have same size as cgp.
1866  assert(is_size(sgp, n)); // sgp must have same size as cgp.
1867  assert(is_size(bgp, n)); // bgp must have same size as cgp.
1868  assert(is_size(pgp, n)); // pgp must have same size as cgp.
1869  assert(is_size(rgp, n)); // rgp must have same size as cgp.
1870  assert(is_size(itw, n, 32)); // We need 32 interpolation
1871  // weights for each position.
1872 
1873  // Check that interpolation weights are valid. The sum of all
1874  // weights (last dimension) must always be approximately one. We
1875  // only check the first element.
1876  assert(
1877  is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
1878 
1879  // We have to loop all the points in the sequence:
1880  for (Index i = 0; i < n; ++i) {
1881  // Current grid positions:
1882  const GridPos& ts = sgp[i];
1883  const GridPos& tb = bgp[i];
1884  const GridPos& tp = pgp[i];
1885  const GridPos& tr = rgp[i];
1886  const GridPos& tc = cgp[i];
1887 
1888  // Get handle to current element of output vector and initialize
1889  // it to zero:
1890  Numeric& tia = ia[i];
1891  tia = 0;
1892 
1893  Index iti = 0;
1894  for (Index s = 0; s < 2; ++s)
1895  for (Index b = 0; b < 2; ++b)
1896  for (Index p = 0; p < 2; ++p)
1897  for (Index r = 0; r < 2; ++r)
1898  for (Index c = 0; c < 2; ++c) {
1899  tia += a.get(ts.idx + s,
1900  tb.idx + b,
1901  tp.idx + p,
1902  tr.idx + r,
1903  tc.idx + c) *
1904  itw.get(i, iti);
1905  ++iti;
1906  }
1907  }
1908 }
1909 
1911 
1937  ConstMatrixView itw,
1938  ConstTensor6View a,
1939  const ArrayOfGridPos& vgp,
1940  const ArrayOfGridPos& sgp,
1941  const ArrayOfGridPos& bgp,
1942  const ArrayOfGridPos& pgp,
1943  const ArrayOfGridPos& rgp,
1944  const ArrayOfGridPos& cgp) {
1945  Index n = cgp.nelem();
1946  assert(is_size(ia, n)); // ia must have same size as cgp.
1947  assert(is_size(vgp, n)); // vgp must have same size as cgp.
1948  assert(is_size(sgp, n)); // sgp must have same size as cgp.
1949  assert(is_size(bgp, n)); // bgp must have same size as cgp.
1950  assert(is_size(pgp, n)); // pgp must have same size as cgp.
1951  assert(is_size(rgp, n)); // rgp must have same size as cgp.
1952  assert(is_size(itw, n, 64)); // We need 64 interpolation
1953  // weights for each position.
1954 
1955  // Check that interpolation weights are valid. The sum of all
1956  // weights (last dimension) must always be approximately one. We
1957  // only check the first element.
1958  assert(
1959  is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
1960 
1961  // We have to loop all the points in the sequence:
1962  for (Index i = 0; i < n; ++i) {
1963  // Current grid positions:
1964  const GridPos& tv = vgp[i];
1965  const GridPos& ts = sgp[i];
1966  const GridPos& tb = bgp[i];
1967  const GridPos& tp = pgp[i];
1968  const GridPos& tr = rgp[i];
1969  const GridPos& tc = cgp[i];
1970 
1971  // Get handle to current element of output vector and initialize
1972  // it to zero:
1973  Numeric& tia = ia[i];
1974  tia = 0;
1975 
1976  Index iti = 0;
1977  for (Index v = 0; v < 2; ++v)
1978  for (Index s = 0; s < 2; ++s)
1979  for (Index b = 0; b < 2; ++b)
1980  for (Index p = 0; p < 2; ++p)
1981  for (Index r = 0; r < 2; ++r)
1982  for (Index c = 0; c < 2; ++c) {
1983  tia += a.get(tv.idx + v,
1984  ts.idx + s,
1985  tb.idx + b,
1986  tp.idx + p,
1987  tr.idx + r,
1988  tc.idx + c) *
1989  itw.get(i, iti);
1990  ++iti;
1991  }
1992  }
1993 }
1994 
1996 // Green interpolation
1998 
2000 
2022  const ArrayOfGridPos& rgp,
2023  const ArrayOfGridPos& cgp) {
2024  Index nr = rgp.nelem();
2025  Index nc = cgp.nelem();
2026  assert(is_size(itw, nr, nc, 4)); // We must store 4 interpolation
2027  // weights for each position.
2028 
2029  // We have to loop all the points in the new grid:
2030  for (Index ir = 0; ir < nr; ++ir) {
2031  // Current grid position:
2032  const GridPos& tr = rgp[ir];
2033 
2034  for (Index ic = 0; ic < nc; ++ic) {
2035  // Current grid position:
2036  const GridPos& tc = cgp[ic];
2037 
2038  // Interpolation weights are stored in this order (l=lower
2039  // u=upper, r=row, c=column):
2040  // 1. l-r l-c
2041  // 2. l-r u-c
2042  // 3. u-r l-c
2043  // 4. u-r u-c
2044 
2045  Index iti = 0;
2046 
2047  LOOPIT(r)
2048  LOOPIT(c) {
2049  itw.get(ir, ic, iti) = (*r) * (*c);
2050  ++iti;
2051  }
2052  }
2053  }
2054 }
2055 
2057 
2080  const ArrayOfGridPos& pgp,
2081  const ArrayOfGridPos& rgp,
2082  const ArrayOfGridPos& cgp) {
2083  Index np = pgp.nelem();
2084  Index nr = rgp.nelem();
2085  Index nc = cgp.nelem();
2086  // We must store 8 interpolation weights for each position:
2087  assert(is_size(itw, np, nr, nc, 8));
2088 
2089  // We have to loop all the points in the new grid:
2090  for (Index ip = 0; ip < np; ++ip) {
2091  const GridPos& tp = pgp[ip];
2092  for (Index ir = 0; ir < nr; ++ir) {
2093  const GridPos& tr = rgp[ir];
2094  for (Index ic = 0; ic < nc; ++ic) {
2095  const GridPos& tc = cgp[ic];
2096 
2097  Index iti = 0;
2098 
2099  LOOPIT(p)
2100  LOOPIT(r)
2101  LOOPIT(c) {
2102  itw.get(ip, ir, ic, iti) = (*p) * (*r) * (*c);
2103  ++iti;
2104  }
2105  }
2106  }
2107  }
2108 }
2109 
2111 
2135  const ArrayOfGridPos& bgp,
2136  const ArrayOfGridPos& pgp,
2137  const ArrayOfGridPos& rgp,
2138  const ArrayOfGridPos& cgp) {
2139  Index nb = bgp.nelem();
2140  Index np = pgp.nelem();
2141  Index nr = rgp.nelem();
2142  Index nc = cgp.nelem();
2143  // We must store 16 interpolation weights for each position:
2144  assert(is_size(itw, nb, np, nr, nc, 16));
2145 
2146  // We have to loop all the points in the new grid:
2147  for (Index ib = 0; ib < nb; ++ib) {
2148  const GridPos& tb = bgp[ib];
2149  for (Index ip = 0; ip < np; ++ip) {
2150  const GridPos& tp = pgp[ip];
2151  for (Index ir = 0; ir < nr; ++ir) {
2152  const GridPos& tr = rgp[ir];
2153  for (Index ic = 0; ic < nc; ++ic) {
2154  const GridPos& tc = cgp[ic];
2155 
2156  Index iti = 0;
2157 
2158  LOOPIT(b)
2159  LOOPIT(p)
2160  LOOPIT(r)
2161  LOOPIT(c) {
2162  itw.get(ib, ip, ir, ic, iti) = (*b) * (*p) * (*r) * (*c);
2163  ++iti;
2164  }
2165  }
2166  }
2167  }
2168  }
2169 }
2170 
2172 
2197  const ArrayOfGridPos& sgp,
2198  const ArrayOfGridPos& bgp,
2199  const ArrayOfGridPos& pgp,
2200  const ArrayOfGridPos& rgp,
2201  const ArrayOfGridPos& cgp) {
2202  Index ns = sgp.nelem();
2203  Index nb = bgp.nelem();
2204  Index np = pgp.nelem();
2205  Index nr = rgp.nelem();
2206  Index nc = cgp.nelem();
2207  // We must store 32 interpolation weights for each position:
2208  assert(is_size(itw, ns, nb, np, nr, nc, 32));
2209 
2210  // We have to loop all the points in the new grid:
2211  for (Index is = 0; is < ns; ++is) {
2212  const GridPos& ts = sgp[is];
2213  for (Index ib = 0; ib < nb; ++ib) {
2214  const GridPos& tb = bgp[ib];
2215  for (Index ip = 0; ip < np; ++ip) {
2216  const GridPos& tp = pgp[ip];
2217  for (Index ir = 0; ir < nr; ++ir) {
2218  const GridPos& tr = rgp[ir];
2219  for (Index ic = 0; ic < nc; ++ic) {
2220  const GridPos& tc = cgp[ic];
2221 
2222  Index iti = 0;
2223 
2224  LOOPIT(s)
2225  LOOPIT(b)
2226  LOOPIT(p)
2227  LOOPIT(r)
2228  LOOPIT(c) {
2229  itw.get(is, ib, ip, ir, ic, iti) =
2230  (*s) * (*b) * (*p) * (*r) * (*c);
2231  ++iti;
2232  }
2233  }
2234  }
2235  }
2236  }
2237  }
2238 }
2239 
2241 
2267  const ArrayOfGridPos& vgp,
2268  const ArrayOfGridPos& sgp,
2269  const ArrayOfGridPos& bgp,
2270  const ArrayOfGridPos& pgp,
2271  const ArrayOfGridPos& rgp,
2272  const ArrayOfGridPos& cgp) {
2273  Index nv = vgp.nelem();
2274  Index ns = sgp.nelem();
2275  Index nb = bgp.nelem();
2276  Index np = pgp.nelem();
2277  Index nr = rgp.nelem();
2278  Index nc = cgp.nelem();
2279  // We must store 64 interpolation weights for each position:
2280  assert(is_size(itw, nv, ns, nb, np, nr, nc, 64));
2281 
2282  // We have to loop all the points in the new grid:
2283  for (Index iv = 0; iv < nv; ++iv) {
2284  const GridPos& tv = vgp[iv];
2285  for (Index is = 0; is < ns; ++is) {
2286  const GridPos& ts = sgp[is];
2287  for (Index ib = 0; ib < nb; ++ib) {
2288  const GridPos& tb = bgp[ib];
2289  for (Index ip = 0; ip < np; ++ip) {
2290  const GridPos& tp = pgp[ip];
2291  for (Index ir = 0; ir < nr; ++ir) {
2292  const GridPos& tr = rgp[ir];
2293  for (Index ic = 0; ic < nc; ++ic) {
2294  const GridPos& tc = cgp[ic];
2295 
2296  Index iti = 0;
2297 
2298  LOOPIT(v)
2299  LOOPIT(s)
2300  LOOPIT(b)
2301  LOOPIT(p)
2302  LOOPIT(r)
2303  LOOPIT(c) {
2304  itw.get(iv, is, ib, ip, ir, ic, iti) =
2305  (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
2306  ++iti;
2307  }
2308  }
2309  }
2310  }
2311  }
2312  }
2313  }
2314 }
2315 
2317 
2339  ConstTensor3View itw,
2340  ConstMatrixView a,
2341  const ArrayOfGridPos& rgp,
2342  const ArrayOfGridPos& cgp) {
2343  Index nr = rgp.nelem();
2344  Index nc = cgp.nelem();
2345  assert(is_size(ia, nr, nc));
2346  assert(is_size(itw, nr, nc, 4)); // We need 4 interpolation
2347  // weights for each position.
2348 
2349  // Check that interpolation weights are valid. The sum of all
2350  // weights (last dimension) must always be approximately one. We
2351  // only check the first element.
2352  assert(is_same_within_epsilon(
2353  itw(0, 0, Range(joker)).sum(), 1, sum_check_epsilon));
2354 
2355  // We have to loop all the points in the new grid:
2356  for (Index ir = 0; ir < nr; ++ir) {
2357  // Current grid position:
2358  const GridPos& tr = rgp[ir];
2359 
2360  for (Index ic = 0; ic < nc; ++ic) {
2361  // Current grid position:
2362  const GridPos& tc = cgp[ic];
2363 
2364  // Get handle to current element of output tensor and initialize
2365  // it to zero:
2366  Numeric& tia = ia(ir, ic);
2367  tia = 0;
2368 
2369  Index iti = 0;
2370  for (Index r = 0; r < 2; ++r)
2371  for (Index c = 0; c < 2; ++c) {
2372  assert(tr.idx + r < a.nrows()); // Temporary !?
2373  assert(tc.idx + c < a.ncols()); // Temporary !?
2374  tia += a.get(tr.idx + r, tc.idx + c) * itw.get(ir, ic, iti);
2375  ++iti;
2376  }
2377  }
2378  }
2379 }
2380 
2382 
2405  ConstTensor4View itw,
2406  ConstTensor3View a,
2407  const ArrayOfGridPos& pgp,
2408  const ArrayOfGridPos& rgp,
2409  const ArrayOfGridPos& cgp) {
2410  Index np = pgp.nelem();
2411  Index nr = rgp.nelem();
2412  Index nc = cgp.nelem();
2413  assert(is_size(ia, np, nr, nc));
2414  assert(is_size(itw, np, nr, nc, 8));
2415 
2416  // Check that interpolation weights are valid. The sum of all
2417  // weights (last dimension) must always be approximately one. We
2418  // only check the first element.
2419  assert(is_same_within_epsilon(
2420  itw(0, 0, 0, Range(joker)).sum(), 1, sum_check_epsilon));
2421 
2422  // We have to loop all the points in the new grid:
2423  for (Index ip = 0; ip < np; ++ip) {
2424  const GridPos& tp = pgp[ip];
2425  for (Index ir = 0; ir < nr; ++ir) {
2426  const GridPos& tr = rgp[ir];
2427  for (Index ic = 0; ic < nc; ++ic) {
2428  // Current grid position:
2429  const GridPos& tc = cgp[ic];
2430 
2431  // Get handle to current element of output tensor and
2432  // initialize it to zero:
2433  Numeric& tia = ia(ip, ir, ic);
2434  tia = 0;
2435 
2436  Index iti = 0;
2437  for (Index p = 0; p < 2; ++p)
2438  for (Index r = 0; r < 2; ++r)
2439  for (Index c = 0; c < 2; ++c) {
2440  assert(tp.idx + p < a.npages()); // Temporary !?
2441  assert(tr.idx + r < a.nrows()); // Temporary !?
2442  assert(tc.idx + c < a.ncols()); // Temporary !?
2443  tia += a.get(tp.idx + p, tr.idx + r, tc.idx + c) *
2444  itw.get(ip, ir, ic, iti);
2445  ++iti;
2446  }
2447  }
2448  }
2449  }
2450 }
2451 
2453 
2477  ConstTensor5View itw,
2478  ConstTensor4View a,
2479  const ArrayOfGridPos& bgp,
2480  const ArrayOfGridPos& pgp,
2481  const ArrayOfGridPos& rgp,
2482  const ArrayOfGridPos& cgp) {
2483  Index nb = bgp.nelem();
2484  Index np = pgp.nelem();
2485  Index nr = rgp.nelem();
2486  Index nc = cgp.nelem();
2487  assert(is_size(ia, nb, np, nr, nc));
2488  assert(is_size(itw, nb, np, nr, nc, 16));
2489 
2490  // Check that interpolation weights are valid. The sum of all
2491  // weights (last dimension) must always be approximately one. We
2492  // only check the first element.
2493  assert(is_same_within_epsilon(
2494  itw(0, 0, 0, 0, Range(joker)).sum(), 1, sum_check_epsilon));
2495 
2496  // We have to loop all the points in the new grid:
2497  for (Index ib = 0; ib < nb; ++ib) {
2498  const GridPos& tb = bgp[ib];
2499  for (Index ip = 0; ip < np; ++ip) {
2500  const GridPos& tp = pgp[ip];
2501  for (Index ir = 0; ir < nr; ++ir) {
2502  const GridPos& tr = rgp[ir];
2503  for (Index ic = 0; ic < nc; ++ic) {
2504  // Current grid position:
2505  const GridPos& tc = cgp[ic];
2506 
2507  // Get handle to current element of output tensor and
2508  // initialize it to zero:
2509  Numeric& tia = ia(ib, ip, ir, ic);
2510  tia = 0;
2511 
2512  Index iti = 0;
2513  for (Index b = 0; b < 2; ++b)
2514  for (Index p = 0; p < 2; ++p)
2515  for (Index r = 0; r < 2; ++r)
2516  for (Index c = 0; c < 2; ++c) {
2517  tia += a.get(tb.idx + b, tp.idx + p, tr.idx + r, tc.idx + c) *
2518  itw.get(ib, ip, ir, ic, iti);
2519  ++iti;
2520  }
2521  }
2522  }
2523  }
2524  }
2525 }
2526 
2528 
2553  ConstTensor6View itw,
2554  ConstTensor5View a,
2555  const ArrayOfGridPos& sgp,
2556  const ArrayOfGridPos& bgp,
2557  const ArrayOfGridPos& pgp,
2558  const ArrayOfGridPos& rgp,
2559  const ArrayOfGridPos& cgp) {
2560  Index ns = sgp.nelem();
2561  Index nb = bgp.nelem();
2562  Index np = pgp.nelem();
2563  Index nr = rgp.nelem();
2564  Index nc = cgp.nelem();
2565  assert(is_size(ia, ns, nb, np, nr, nc));
2566  assert(is_size(itw, ns, nb, np, nr, nc, 32));
2567 
2568  // Check that interpolation weights are valid. The sum of all
2569  // weights (last dimension) must always be approximately one. We
2570  // only check the first element.
2571  assert(is_same_within_epsilon(
2572  itw(0, 0, 0, 0, 0, Range(joker)).sum(), 1, sum_check_epsilon));
2573 
2574  // We have to loop all the points in the new grid:
2575  for (Index is = 0; is < ns; ++is) {
2576  const GridPos& ts = sgp[is];
2577  for (Index ib = 0; ib < nb; ++ib) {
2578  const GridPos& tb = bgp[ib];
2579  for (Index ip = 0; ip < np; ++ip) {
2580  const GridPos& tp = pgp[ip];
2581  for (Index ir = 0; ir < nr; ++ir) {
2582  const GridPos& tr = rgp[ir];
2583  for (Index ic = 0; ic < nc; ++ic) {
2584  // Current grid position:
2585  const GridPos& tc = cgp[ic];
2586 
2587  // Get handle to current element of output tensor and
2588  // initialize it to zero:
2589  Numeric& tia = ia(is, ib, ip, ir, ic);
2590  tia = 0;
2591 
2592  Index iti = 0;
2593  for (Index s = 0; s < 2; ++s)
2594  for (Index b = 0; b < 2; ++b)
2595  for (Index p = 0; p < 2; ++p)
2596  for (Index r = 0; r < 2; ++r)
2597  for (Index c = 0; c < 2; ++c) {
2598  tia += a.get(ts.idx + s,
2599  tb.idx + b,
2600  tp.idx + p,
2601  tr.idx + r,
2602  tc.idx + c) *
2603  itw.get(is, ib, ip, ir, ic, iti);
2604  ++iti;
2605  }
2606  }
2607  }
2608  }
2609  }
2610  }
2611 }
2612 
2614 
2640  ConstTensor7View itw,
2641  ConstTensor6View a,
2642  const ArrayOfGridPos& vgp,
2643  const ArrayOfGridPos& sgp,
2644  const ArrayOfGridPos& bgp,
2645  const ArrayOfGridPos& pgp,
2646  const ArrayOfGridPos& rgp,
2647  const ArrayOfGridPos& cgp) {
2648  Index nv = vgp.nelem();
2649  Index ns = sgp.nelem();
2650  Index nb = bgp.nelem();
2651  Index np = pgp.nelem();
2652  Index nr = rgp.nelem();
2653  Index nc = cgp.nelem();
2654  assert(is_size(ia, nv, ns, nb, np, nr, nc));
2655  assert(is_size(itw, nv, ns, nb, np, nr, nc, 64));
2656 
2657  // Check that interpolation weights are valid. The sum of all
2658  // weights (last dimension) must always be approximately one. We
2659  // only check the first element.
2660  assert(is_same_within_epsilon(
2661  itw(0, 0, 0, 0, 0, 0, Range(joker)).sum(), 1, sum_check_epsilon));
2662 
2663  // We have to loop all the points in the new grid:
2664  for (Index iv = 0; iv < nv; ++iv) {
2665  const GridPos& tv = vgp[iv];
2666  for (Index is = 0; is < ns; ++is) {
2667  const GridPos& ts = sgp[is];
2668  for (Index ib = 0; ib < nb; ++ib) {
2669  const GridPos& tb = bgp[ib];
2670  for (Index ip = 0; ip < np; ++ip) {
2671  const GridPos& tp = pgp[ip];
2672  for (Index ir = 0; ir < nr; ++ir) {
2673  const GridPos& tr = rgp[ir];
2674  for (Index ic = 0; ic < nc; ++ic) {
2675  // Current grid position:
2676  const GridPos& tc = cgp[ic];
2677 
2678  // Get handle to current element of output tensor and
2679  // initialize it to zero:
2680  Numeric& tia = ia(iv, is, ib, ip, ir, ic);
2681  tia = 0;
2682 
2683  Index iti = 0;
2684  for (Index v = 0; v < 2; ++v)
2685  for (Index s = 0; s < 2; ++s)
2686  for (Index b = 0; b < 2; ++b)
2687  for (Index p = 0; p < 2; ++p)
2688  for (Index r = 0; r < 2; ++r)
2689  for (Index c = 0; c < 2; ++c) {
2690  tia += a.get(tv.idx + v,
2691  ts.idx + s,
2692  tb.idx + b,
2693  tp.idx + p,
2694  tr.idx + r,
2695  tc.idx + c) *
2696  itw.get(iv, is, ib, ip, ir, ic, iti);
2697  ++iti;
2698  }
2699  }
2700  }
2701  }
2702  }
2703  }
2704  }
2705 }
2706 
2708 
2724  ConstVectorView y,
2725  const Numeric& x_i,
2726  const GridPos& gp) {
2727  Index N_x = x.nelem();
2728 
2729  assert(N_x == y.nelem());
2730  assert(N_x > 2);
2731 
2732  Vector xa(4), ya(4);
2733  Numeric y_int;
2734  Numeric dy_int;
2735  y_int = 0.;
2736 
2737  // 1 - polynomial interpolation (3 points) with grid position search
2738  // 2 - polynomial interpolation (3 points) without grid position search
2739  // 3 - polynomial interpolation (4 points)
2740 
2741  Index interp_method = 1;
2742 
2743  if (interp_method == 1) {
2744  // Pick out three points for interpolation
2745  if ((gp.fd[0] <= 0.5 && gp.idx > 0) || gp.idx == N_x - 2) {
2746  xa[0] = x[gp.idx - 1];
2747  xa[1] = x[gp.idx];
2748  xa[2] = x[gp.idx + 1];
2749 
2750  ya[0] = y[gp.idx - 1];
2751  ya[1] = y[gp.idx];
2752  ya[2] = y[gp.idx + 1];
2753  }
2754 
2755  else if ((gp.fd[0] > 0.5 && gp.idx < N_x - 2) || gp.idx == 0) {
2756  xa[0] = x[gp.idx];
2757  xa[1] = x[gp.idx + 1];
2758  xa[2] = x[gp.idx + 2];
2759 
2760  ya[0] = y[gp.idx];
2761  ya[1] = y[gp.idx + 1];
2762  ya[2] = y[gp.idx + 2];
2763  }
2764 
2765  else if (gp.idx == N_x - 1) {
2766  xa[0] = x[N_x - 2];
2767  xa[1] = x[N_x - 1];
2768  xa[2] = x[N_x];
2769 
2770  ya[0] = y[N_x - 2];
2771  ya[1] = y[N_x - 1];
2772  ya[2] = y[N_x];
2773  } else {
2774  assert(false);
2775  arts_exit();
2776  }
2777 
2778  polint(y_int, dy_int, xa, ya, 3, x_i);
2779 
2780  }
2781 
2782  else if (interp_method == 2) {
2783  if (gp.idx == 0) {
2784  xa[0] = x[gp.idx];
2785  xa[1] = x[gp.idx + 1];
2786  xa[2] = x[gp.idx + 2];
2787 
2788  ya[0] = y[gp.idx];
2789  ya[1] = y[gp.idx + 1];
2790  ya[2] = y[gp.idx + 2];
2791  } else if (gp.idx == N_x - 1) {
2792  xa[0] = x[gp.idx - 2];
2793  xa[1] = x[gp.idx - 1];
2794  xa[2] = x[gp.idx];
2795 
2796  ya[0] = y[gp.idx - 2];
2797  ya[1] = y[gp.idx - 1];
2798  ya[2] = y[gp.idx];
2799  } else {
2800  xa[0] = x[gp.idx - 1];
2801  xa[1] = x[gp.idx];
2802  xa[2] = x[gp.idx + 1];
2803 
2804  ya[0] = y[gp.idx - 1];
2805  ya[1] = y[gp.idx];
2806  ya[2] = y[gp.idx + 1];
2807  }
2808 
2809  // Polynominal interpolation, n = 3
2810  polint(y_int, dy_int, xa, ya, 3, x_i);
2811  }
2812 
2813  else if (interp_method == 3) {
2814  // Take 4 points
2815  if (gp.idx == 0) {
2816  xa[0] = -x[gp.idx + 1];
2817  xa[1] = x[gp.idx + 0];
2818  xa[2] = x[gp.idx + 1];
2819  xa[3] = x[gp.idx + 2];
2820 
2821  ya[0] = y[gp.idx + 1];
2822  ya[1] = y[gp.idx + 0];
2823  ya[2] = y[gp.idx + 1];
2824  ya[3] = y[gp.idx + 2];
2825  } else if (gp.idx == N_x - 1) {
2826  xa[0] = x[gp.idx - 1];
2827  xa[1] = x[gp.idx - 0];
2828  xa[2] = 2 * x[gp.idx] - x[gp.idx - 1];
2829  xa[3] = 2 * x[gp.idx] - x[gp.idx - 2];
2830 
2831  ya[0] = y[gp.idx - 1];
2832  ya[1] = y[gp.idx - 0];
2833  ya[2] = y[gp.idx - 1];
2834  ya[3] = y[gp.idx - 2];
2835  } else if (gp.idx == N_x - 2) {
2836  xa[0] = x[gp.idx - 2];
2837  xa[1] = x[gp.idx - 1];
2838  xa[2] = x[gp.idx];
2839  xa[3] = x[gp.idx + 1];
2840 
2841  ya[0] = y[gp.idx - 2];
2842  ya[1] = y[gp.idx - 1];
2843  ya[2] = y[gp.idx];
2844  ya[3] = y[gp.idx + 1];
2845  } else {
2846  xa[0] = x[gp.idx - 1];
2847  xa[1] = x[gp.idx];
2848  xa[2] = x[gp.idx + 1];
2849  xa[3] = x[gp.idx + 2];
2850 
2851  ya[0] = y[gp.idx - 1];
2852  ya[1] = y[gp.idx];
2853  ya[2] = y[gp.idx + 1];
2854  ya[3] = y[gp.idx + 2];
2855  }
2856  // Polinominal interpolation, n = 4
2857  polint(y_int, dy_int, xa, ya, 4, x_i);
2858  }
2859 
2860  return y_int;
2861 }
2862 
2864 
2883 void polint(Numeric& y_int,
2884  Numeric& dy_int,
2885  ConstVectorView xa,
2886  ConstVectorView ya,
2887  const Index& n,
2888  const Numeric& x) {
2889  Index ns = 1;
2890  Numeric den, dif, dift, ho, hp, w;
2891 
2892  dif = abs(x - xa[0]);
2893 
2894  Vector c(n);
2895  Vector d(n);
2896 
2897  // Find the index of the closest table entry
2898  for (Index i = 0; i < n; i++) {
2899  if ((dift = abs(x - xa[i])) < dif) {
2900  ns = i;
2901  dif = dift;
2902  }
2903  // Initialize c and d
2904  c[i] = ya[i];
2905  d[i] = ya[i];
2906  }
2907  // Initial approximation to y
2908  y_int = ya[ns--];
2909 
2910  for (Index m = 1; m < n; m++) {
2911  for (Index i = 0; i < n - m; i++) {
2912  ho = xa[i] - x;
2913  hp = xa[i + m] - x;
2914  w = c[i + 1] - d[i];
2915  den = ho - hp;
2916  // This error occurs when two input xa's are identical.
2917  assert(den != 0.);
2918  den = w / den;
2919  d[i] = hp * den;
2920  c[i] = ho * den;
2921  }
2922  y_int += (dy_int = (2 * (ns + 1) < (n - m) ? c[ns + 1] : d[ns--]));
2923  }
2924 }
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
Numeric interp_poly(ConstVectorView x, ConstVectorView y, const Numeric &x_i, const GridPos &gp)
Polynomial interpolation.
The VectorView class.
Definition: matpackI.h:610
#define ns
The Tensor4View class.
Definition: matpackIV.h:284
Index nelem() const
Number of elements.
Definition: array.h:195
void arts_exit(int status)
This is the exit function of ARTS.
Definition: arts.cc:42
A constant view of a Tensor7.
Definition: matpackVII.h:147
The Tensor7View class.
Definition: matpackVII.h:1286
void gridpos_upperend_check(GridPos &gp, const Index &ie)
gridpos_upperend_check
Numeric & get(Index n)
Get element implementation without assertions.
Definition: matpackI.h:638
Numeric get(Index v, Index s, Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackVI.h:550
The Vector class.
Definition: matpackI.h:860
#define abs(x)
The MatrixView class.
Definition: matpackI.h:1093
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
A constant view of a Tensor6.
Definition: matpackVI.h:149
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:215
The range class.
Definition: matpackI.h:160
Numeric get(Index r, Index c) const
Get element implementation without assertions.
Definition: matpackI.h:1009
Numeric fractional_gp(const GridPos &gp)
fractional_gp
Numeric & get(Index v, Index s, Index b, Index p, Index r, Index c)
Get element implementation without assertions.
Definition: matpackVI.h:1015
cmplx FADDEEVA() w(cmplx z, double relerr)
Definition: Faddeeva.cc:680
The Tensor6View class.
Definition: matpackVI.h:621
Header file for interpolation.cc.
Index gridpos2gridrange(const GridPos &gp, const bool &upwards)
gridpos2gridrange
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:147
bool is_gridpos_at_index_i(const GridPos &gp, const Index &i, const bool &strict)
is_gridpos_at_index_i
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
Structure to store a grid position.
Definition: interpolation.h:73
A constant view of a Tensor4.
Definition: matpackIV.h:133
void polint(Numeric &y_int, Numeric &dy_int, ConstVectorView xa, ConstVectorView ya, const Index &n, const Numeric &x)
Polynomial interpolation.
This file contains the definition of Array.
void gridpos_check_fd(GridPos &gp)
gridpos_check_fd
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
#define DEBUG_ONLY(...)
Definition: debug.h:36
ostream & operator<<(ostream &os, const GridPos &gp)
Output operator for GridPos.
bool is_same_within_epsilon(const Numeric &a, const Numeric &b, const Numeric &epsilon)
Check, if two numbers agree within a given epsilon.
Definition: logic.cc:351
Numeric fd[2]
Definition: interpolation.h:75
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
Numeric sum() const
The sum of all elements of a Vector.
Definition: matpackI.cc:53
The Tensor3View class.
Definition: matpackIII.h:239
Numeric get(Index n) const
Get element implementation without assertions.
Definition: matpackI.h:514
Numeric & get(Index r, Index c)
Get element implementation without assertions.
Definition: matpackI.h:1118
#define LOOPIT(x)
Macro for interpolation weight loops.
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
A constant view of a Tensor5.
Definition: matpackV.h:143
const Joker joker
Index idx
Definition: interpolation.h:74
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Numeric & get(Index p, Index r, Index c)
Get element implementation without assertions.
Definition: matpackIII.h:275
void gridpos_copy(GridPos &gp_new, const GridPos &gp_old)
gridpos_copy
The Tensor5View class.
Definition: matpackV.h:333
Header file for logic.cc.
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
This can be used to make arrays out of anything.
Definition: array.h:40
Numeric & get(Index b, Index p, Index r, Index c)
Get element implementation without assertions.
Definition: matpackIV.h:348
Numeric get(Index l, Index v, Index s, Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackVII.h:1211
A constant view of a Tensor3.
Definition: matpackIII.h:132
A constant view of a Vector.
Definition: matpackI.h:476
Numeric & get(Index s, Index b, Index p, Index r, Index c)
Get element implementation without assertions.
Definition: matpackV.h:431
A constant view of a Matrix.
Definition: matpackI.h:982
Numeric & get(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Get element implementation without assertions.
Definition: matpackVII.h:2308
Numeric get(Index s, Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackV.h:265
bool is_decreasing(ConstVectorView x)
Checks if a vector is sorted in reversed order and is strictly decreasing.
Definition: logic.cc:252
void gridpos_1to1(ArrayOfGridPos &gp, ConstVectorView grid)
gridpos_1to1
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:79
const Numeric FD_TOL
The maximum difference from 1 that we allow for a sum check.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
Numeric get(Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackIII.h:180
void gridpos_force_end_fd(GridPos &gp, const Index &n)
gridpos_force_end_fd
Numeric get(Index b, Index p, Index r, Index c) const
Get element implementation without assertions.
Definition: matpackIV.h:220
void gp4length1grid(ArrayOfGridPos &gp)
Grid position matching a grid of length 1.