interpolation.cc

Go to the documentation of this file.

/* Copyright (C) 2002-2012 Stefan Buehler <sbuehler@ltu.se>
 
   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.
 
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */
 
#include "interpolation.h"
#include <cmath>
#include <iostream>
#include "array.h"
#include "check_input.h"
#include "logic.h"
 
// File-global constants:
 
 
DEBUG_ONLY(const Numeric sum_check_epsilon = 1e-6;)
 
 
const Numeric FD_TOL = 1.5e-3;
 
 
#define LOOPIT(x) for (const Numeric* x = &t##x.fd[1]; x >= &t##x.fd[0]; --x)
 
 
ostream& operator<<(ostream& os, const GridPos& gp) {
  os << gp.idx << " " << gp.fd[0] << " " << gp.fd[1] << "\n";
  return os;
}
 
 
void gridpos(ArrayOfGridPos& gp,
             ConstVectorView old_grid,
             ConstVectorView new_grid,
             const Numeric& extpolfac) {
  const Index n_old = old_grid.nelem();
  const Index n_new = new_grid.nelem();
 
  // Assert that gp has the right size:
  assert(is_size(gp, n_new));
 
  // Assert, that the old grid has more than one element
  assert(1 < n_old);
 
  // This function hast two parts, depending on whether old_grid is
  // sorted in ascending or descending order. Maybe that's not too
  // elegant, but it's the most efficient, because in this way there
  // is no additional runtime overhead to handle both cases.
 
  // We use only the first two elements to decide how the grid is
  // sorted. (The rest of the grid is checked later by an assert.)
  // If both values are the same, we still assume the grid is
  // ascending. However, this will lead to an assertion fail later on,
  // because the grid has to be strictly sorted.
  bool ascending = (old_grid[0] <= old_grid[1]);
 
  if (ascending) {
    // So, old_grid should always be sorted in strictly ascending order.
    // (No duplicate values.)
    // Assert that this is so. This may depend on user input, however,
    // inside this elementary function is not the place to check for
    // that. There must be runtime checks on higher levels to insure
    // that all grids are sorted. The assertion here is just as a last
    // safety check.
    assert(is_increasing(old_grid));
 
    // Limits of extrapolation.
    const Numeric og_min =
        old_grid[0] - extpolfac * (old_grid[1] - old_grid[0]);
    const Numeric og_max =
        old_grid[n_old - 1] +
        extpolfac * (old_grid[n_old - 1] - old_grid[n_old - 2]);
 
    // We will make no firm assumptions about the new grid. But the case
    // that we have in mind is that it is either also sorted, or at
    // least partially sorted, for example like this:
    // 5 4 3 2 2.5 3 4
    // This kind of sequence should be typical if we interpolate
    // atmospheric fields along a limb propagation path.
 
    // Let's get some idea where the first point in the new grid is,
    // relative to the old grid. We use linear interpolation between the
    // maximum and the minimum of the old grid for this. (Taking
    // into account the small allowed extrapolation.)
    Numeric frac = (new_grid[0] - og_min) / (og_max - og_min);
 
    // We are not checking if the value of frac is reasonable,
    // because there is another assertion below to catch the
    // consequences.
 
    // Initialize current_position:
    Index current_position = (Index)rint(frac * (Numeric)(n_old - 2));
 
    // The above statement should satisfy
    // 0 <= current_position <= n_old-2
    // Assert that this is indeed the case.
    //       cout << "frac             = " << frac             << "\n";
    //       cout << "current_position = " << current_position << "\n";
    assert(0 <= current_position);
    assert(current_position <= n_old - 2);
 
    // The variables lower and upper are used to remember the value of
    // the old grid at current_position and one above current_position:
    Numeric lower = old_grid[current_position];
    Numeric upper = old_grid[current_position + 1];
 
    // Loop over all points in the new grid:
    for (Index i_new = 0; i_new < n_new; ++i_new) {
      // Get a reference to the current element of gp:
      GridPos& tgp = gp[i_new];
      // And on the current value of the new grid:
      const Numeric tng = new_grid[i_new];
 
      // Verify that the new grid is within the limits of
      // extrapolation that we have defined by og_min and og_max:
      assert(og_min <= tng);
      assert(tng <= og_max);
 
      //           cout << "lower / tng / upper = "
      //                << lower << " / "
      //                << tng   << " / "
      //                << upper << "\n";
 
      // Is current_position too high?
      // (The current_position>0 condition is there so that the position
      // stays 0 for extrapolation.)
      if (tng < lower && current_position > 0) {
        do {
          --current_position;
          lower = old_grid[current_position];
        } while (tng < lower && current_position > 0);
 
        upper = old_grid[current_position + 1];
 
        tgp.idx = current_position;
        tgp.fd[0] = (tng - lower) / (upper - lower);
        tgp.fd[1] = 1.0 - tgp.fd[0];
      } else {
        // Is it too low?
        // (The current_position<n_old condition is there so
        // that uppers stays n_old-1 for extrapolation.)
        if (tng >= upper && current_position < n_old - 2) {
          do {
            ++current_position;
            upper = old_grid[current_position + 1];
          } while (tng >= upper && current_position < n_old - 2);
 
          lower = old_grid[current_position];
 
          tgp.idx = current_position;
          tgp.fd[0] = (tng - lower) / (upper - lower);
          tgp.fd[1] = 1.0 - tgp.fd[0];
        } else {
          // None of the other two conditions were true. That means:
          // lower <= tng < upper. The current_position is still
          // valid.
 
          // SAB 2010-04-28: Note that if a new grid point is exactly on
          // top of an old grid point, then you are now guaranteed to get
          // fd[0] = 0 and fd[1] = 1. (Except at the upper grid end.)
 
          tgp.idx = current_position;
          tgp.fd[0] = (tng - lower) / (upper - lower);
          tgp.fd[1] = 1.0 - tgp.fd[0];
        }
      }
 
      //          cout << tgp.idx << " " << tgp.fd[0] << " " << tgp.fd[1] << endl;
 
      // Safety check to ensure the above promise:
      assert(old_grid[tgp.idx] <= tng || tgp.idx == 0);
    }
  } else  //   if (ascending)
  {
    // Now we are in the "descending old grid" part. We do exactly
    // the same as in the other part, just accounting for the
    // different order of things. Comments here refer only to
    // interesting differences from the ascending case. See that
    // case for more general comments.
 
    // This time ensure strictly descending order:
    assert(is_decreasing(old_grid));
 
    // The max is now the first point, the min the last point!
    // I think the sign is right here...
    const Numeric og_max =
        old_grid[0] - extpolfac * (old_grid[1] - old_grid[0]);
    const Numeric og_min =
        old_grid[n_old - 1] +
        extpolfac * (old_grid[n_old - 1] - old_grid[n_old - 2]);
 
    // We have to take 1- here, because we are starting from the
    // high end.
    Numeric frac = 1 - (new_grid[0] - og_min) / (og_max - og_min);
 
    // We are not checking if the value of frac is reasonable,
    // because there is another assertion below to catch the
    // consequences.
 
    Index current_position = (Index)rint(frac * (Numeric)(n_old - 2));
 
    // The above statement should satisfy
    // 0 <= current_position <= n_old-2
    // Assert that this is indeed the case.
    //       cout << "frac             = " << frac             << "\n";
    //       cout << "current_position = " << current_position << "\n";
    assert(0 <= current_position);
    assert(current_position <= n_old - 2);
 
    // Note, that old_grid[lower] has a higher numerical value than
    // old_grid[upper]!
    Numeric lower = old_grid[current_position];
    Numeric upper = old_grid[current_position + 1];
 
    for (Index i_new = 0; i_new < n_new; ++i_new) {
      GridPos& tgp = gp[i_new];
      const Numeric tng = new_grid[i_new];
 
      // Verify that the new grid is within the limits of
      // extrapolation that we have defined by og_min and og_max:
      assert(og_min <= tng);
      assert(tng <= og_max);
 
      //           cout << "lower / tng / upper = "
      //                << lower << " / "
      //                << tng   << " / "
      //                << upper << "\n";
 
      // Is current_position too high? (Sign of comparison changed
      // compared to ascending case!)
      if (tng > lower && current_position > 0) {
        do {
          --current_position;
          lower = old_grid[current_position];
        } while (tng > lower && current_position > 0);
 
        upper = old_grid[current_position + 1];
 
        tgp.idx = current_position;
        tgp.fd[0] = (tng - lower) / (upper - lower);
        tgp.fd[1] = 1.0 - tgp.fd[0];
      } else {
        // Is it too low? (Sign of comparison changed
        // compared to ascending case!)
        if (tng <= upper && current_position < n_old - 2) {
          do {
            ++current_position;
            upper = old_grid[current_position + 1];
          } while (tng <= upper && current_position < n_old - 2);
 
          lower = old_grid[current_position];
 
          tgp.idx = current_position;
          tgp.fd[0] = (tng - lower) / (upper - lower);
          tgp.fd[1] = 1.0 - tgp.fd[0];
        } else {
          // None of the other two conditions were true. That means:
          // lower >= tng > upper. The current_position is still
          // valid. (Note that upper and lower have switched
          // place compared to the ascending case.)
 
          // SAB 2010-04-28: Note that if a new grid point is exactly on
          // top of an old grid point, then you are now guaranteed to get
          // fd[0] = 0 and fd[1] = 1. (Except at the upper grid end.)
 
          tgp.idx = current_position;
          tgp.fd[0] = (tng - lower) / (upper - lower);
          tgp.fd[1] = 1.0 - tgp.fd[0];
        }
      }
 
      // Safety check to ensure the above promise:
      assert(old_grid[tgp.idx] >= tng || tgp.idx == 0);
    }
  }
}
 
 
void gridpos(GridPos& gp,
             ConstVectorView old_grid,
             const Numeric& new_grid,
             const Numeric& extpolfac) {
  ArrayOfGridPos agp(1);
  Vector v(1, new_grid);
  gridpos(agp, old_grid, v, extpolfac);
  gridpos_copy(gp, agp[0]);
}
 
 
void gridpos_1to1(ArrayOfGridPos& gp, ConstVectorView grid) {
  const Index n = grid.nelem();
  gp.resize(n);
 
  for (Index i = 0; i < n - 1; i++) {
    gp[i].idx = i;
    gp[i].fd[0] = 0;
    gp[i].fd[1] = 1;
  }
  gp[n - 1].idx = n - 2;
  gp[n - 1].fd[0] = 1;
  gp[n - 1].fd[1] = 0;
}
 
 
void gridpos_copy(GridPos& gp_new, const GridPos& gp_old) {
  gp_new.idx = gp_old.idx;
  gp_new.fd[0] = gp_old.fd[0];
  gp_new.fd[1] = gp_old.fd[1];
}
 
 
Numeric fractional_gp(const GridPos& gp) {
  return (Numeric(gp.idx) + gp.fd[0]);
}
 
 
void gridpos_check_fd(GridPos& gp) {
  // Catch values that "must" be wrong
  assert(gp.fd[0] > -FD_TOL);
  assert(gp.fd[0] < 1.0 + FD_TOL);
  assert(gp.fd[1] > -FD_TOL);
  assert(gp.fd[1] < 1.0 + FD_TOL);
 
  if (gp.fd[0] < 0.0) {
    gp.fd[0] = 0.0;
    gp.fd[1] = 1.0;
  } else if (gp.fd[0] > 1.0) {
    gp.fd[0] = 1.0;
    gp.fd[1] = 0.0;
  }
 
  if (gp.fd[1] < 0.0) {
    gp.fd[1] = 0.0;
    gp.fd[0] = 1.0;
  } else if (gp.fd[1] > 1.0) {
    gp.fd[1] = 1.0;
    gp.fd[0] = 0.0;
  }
}
 
 
void gridpos_force_end_fd(GridPos& gp, const Index& n) {
  assert(gp.idx >= 0);
 
  // If fd=1, shift to grid index above
  if (gp.fd[0] > 0.5) {
    gp.idx += 1;
  }
  gp.fd[0] = 0;
  gp.fd[1] = 1;
 
  assert(gp.idx < n);
 
  // End of complete grid range must be handled separately
  if (gp.idx == n - 1) {
    gp.idx -= 1;
    gp.fd[0] = 1;
    gp.fd[1] = 0;
  }
}
 
 
void gridpos_upperend_check(GridPos& gp, const Index& ie) {
  if (gp.idx == ie) {
    assert(gp.fd[0] < 0.005);  // To capture obviously bad cases
    gp.idx -= 1;
    gp.fd[0] = 1.0;
    gp.fd[1] = 0.0;
  }
}
 
 
void gridpos_upperend_check(ArrayOfGridPos& gp, const Index& ie) {
  for (Index i = 0; i < gp.nelem(); i++) {
    if (gp[i].idx == ie) {
      assert(gp[i].fd[0] < 0.005);  // To capture obviously bad cases
      gp[i].idx -= 1;
      gp[i].fd[0] = 1.0;
      gp[i].fd[1] = 0.0;
    }
  }
}
 
 
void gp4length1grid(ArrayOfGridPos& gp) {
  for (Index i = 0; i < gp.nelem(); i++) {
    gp[i].idx = 0;
    gp[i].fd[0] = 0;
    gp[i].fd[1] = 1;
  }
}
 
 
bool is_gridpos_at_index_i(const GridPos& gp,
                           const Index& i,
                           const bool& strict) {
  if (strict) {
    // Assume that gridpos_force_end_fd has been used. The expression 0==0
    // should be safer than 1==1.
    if (gp.idx == i && gp.fd[0] == 0) {
      return true;
    } else if (gp.idx == i - 1 && gp.fd[1] == 0) {
      return true;
    }
  } else {
    if (gp.idx == i && gp.fd[0] < FD_TOL) {
      return true;
    } else if (gp.idx == i - 1 && gp.fd[1] < FD_TOL) {
      return true;
    }
  }
  return false;
}
 
 
Index gridpos2gridrange(const GridPos& gp, const bool& upwards) {
  assert(gp.fd[0] >= 0);
  assert(gp.fd[1] >= 0);
 
  // Not at a grid point
  if (gp.fd[0] > 0 && gp.fd[1] > 0) {
    return gp.idx;
  }
 
  // Fractional distance 0
  else if (gp.fd[0] == 0) {
    if (upwards)
      return gp.idx;
    else
      return gp.idx - 1;
  }
 
  // Fractional distance 1
  else {
    if (upwards)
      return gp.idx + 1;
    else
      return gp.idx;
  }
}
 
//                      Red Interpolation
 
 
void interpweights(VectorView itw, const GridPos& tc) {
  assert(is_size(itw, 2));  // We must store 2 interpolation
                            // weights.
 
  // Interpolation weights are stored in this order (l=lower
  // u=upper, c=column):
  // 1. l-c
  // 2. u-c
 
  Index iti = 0;
 
  // We could use a straight-forward for loop here:
  //
  //       for ( Index c=1; c>=0; --c )
  //    {
  //      ti[iti] = tc.fd[c];
  //      ++iti;
  //    }
  //
  // However, it is slightly faster to use pointers (I tried it,
  // the speed gain is about 20%). So we should write something
  // like:
  //
  //       for ( const Numeric* c=&tc.fd[1]; c>=&tc.fd[0]; --c )
  //    {
  //      ti[iti] = *c;
  //      ++iti;
  //    }
  //
  // For higher dimensions we have to nest these loops. To avoid
  // typos and safe typing, I use the LOOPIT macro, which expands
  // to the for loop above. Note: NO SEMICOLON AFTER THE LOOPIT
  // COMMAND!
 
  LOOPIT(c) {
    itw.get(iti) = *c;
    ++iti;
  }
}
 
 
void interpweights(VectorView itw, const GridPos& tr, const GridPos& tc) {
  assert(is_size(itw, 4));  // We must store 4 interpolation
                            // weights.
  Index iti = 0;
 
  LOOPIT(r)
  LOOPIT(c) {
    itw.get(iti) = (*r) * (*c);
    ++iti;
  }
}
 
 
void interpweights(VectorView itw,
                   const GridPos& tp,
                   const GridPos& tr,
                   const GridPos& tc) {
  assert(is_size(itw, 8));  // We must store 8 interpolation
                            // weights.
  Index iti = 0;
 
  LOOPIT(p)
  LOOPIT(r)
  LOOPIT(c) {
    itw.get(iti) = (*p) * (*r) * (*c);
    ++iti;
  }
}
 
 
void interpweights(VectorView itw,
                   const GridPos& tb,
                   const GridPos& tp,
                   const GridPos& tr,
                   const GridPos& tc) {
  assert(is_size(itw, 16));  // We must store 16 interpolation
                             // weights.
  Index iti = 0;
 
  LOOPIT(b)
  LOOPIT(p)
  LOOPIT(r)
  LOOPIT(c) {
    itw.get(iti) = (*b) * (*p) * (*r) * (*c);
    ++iti;
  }
}
 
 
void interpweights(VectorView itw,
                   const GridPos& ts,
                   const GridPos& tb,
                   const GridPos& tp,
                   const GridPos& tr,
                   const GridPos& tc) {
  assert(is_size(itw, 32));  // We must store 32 interpolation
                             // weights.
  Index iti = 0;
 
  LOOPIT(s)
  LOOPIT(b)
  LOOPIT(p)
  LOOPIT(r)
  LOOPIT(c) {
    itw.get(iti) = (*s) * (*b) * (*p) * (*r) * (*c);
    ++iti;
  }
}
 
 
void interpweights(VectorView itw,
                   const GridPos& tv,
                   const GridPos& ts,
                   const GridPos& tb,
                   const GridPos& tp,
                   const GridPos& tr,
                   const GridPos& tc) {
  assert(is_size(itw, 64));  // We must store 64 interpolation
                             // weights.
  Index iti = 0;
 
  LOOPIT(v)
  LOOPIT(s)
  LOOPIT(b)
  LOOPIT(p)
  LOOPIT(r)
  LOOPIT(c) {
    itw.get(iti) = (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
    ++iti;
  }
}
 
 
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos& tc) {
  assert(is_size(itw, 2));  // We need 2 interpolation
                            // weights.
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one.
  assert(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
 
  // To store interpolated value:
  Numeric tia = 0;
 
  Index iti = 0;
  for (Index c = 0; c < 2; ++c) {
    tia += a.get(tc.idx + c) * itw.get(iti);
    ++iti;
  }
 
  return tia;
}
 
 
Numeric interp(ConstVectorView itw,
               ConstMatrixView a,
               const GridPos& tr,
               const GridPos& tc) {
  assert(is_size(itw, 4));  // We need 4 interpolation
                            // weights.
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one.
  assert(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
 
  // To store interpolated value:
  Numeric tia = 0;
 
  Index iti = 0;
  for (Index r = 0; r < 2; ++r)
    for (Index c = 0; c < 2; ++c) {
      tia += a.get(tr.idx + r, tc.idx + c) * itw.get(iti);
      ++iti;
    }
 
  return tia;
}
 
 
Numeric interp(ConstVectorView itw,
               ConstTensor3View a,
               const GridPos& tp,
               const GridPos& tr,
               const GridPos& tc) {
  assert(is_size(itw, 8));  // We need 8 interpolation
                            // weights.
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one.
  assert(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
 
  // To store interpolated value:
  Numeric tia = 0;
 
  Index iti = 0;
  for (Index p = 0; p < 2; ++p)
    for (Index r = 0; r < 2; ++r)
      for (Index c = 0; c < 2; ++c) {
        tia += a.get(tp.idx + p, tr.idx + r, tc.idx + c) * itw.get(iti);
        ++iti;
      }
 
  return tia;
}
 
 
Numeric interp(ConstVectorView itw,
               ConstTensor4View a,
               const GridPos& tb,
               const GridPos& tp,
               const GridPos& tr,
               const GridPos& tc) {
  assert(is_size(itw, 16));  // We need 16 interpolation
                             // weights.
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one.
  assert(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
 
  // To store interpolated value:
  Numeric tia = 0;
 
  Index iti = 0;
  for (Index b = 0; b < 2; ++b)
    for (Index p = 0; p < 2; ++p)
      for (Index r = 0; r < 2; ++r)
        for (Index c = 0; c < 2; ++c) {
          tia += a.get(tb.idx + b, tp.idx + p, tr.idx + r, tc.idx + c) *
                 itw.get(iti);
          ++iti;
        }
 
  return tia;
}
 
 
Numeric interp(ConstVectorView itw,
               ConstTensor5View a,
               const GridPos& ts,
               const GridPos& tb,
               const GridPos& tp,
               const GridPos& tr,
               const GridPos& tc) {
  assert(is_size(itw, 32));  // We need 32 interpolation
                             // weights.
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one.
  assert(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
 
  // To store interpolated value:
  Numeric tia = 0;
 
  Index iti = 0;
  for (Index s = 0; s < 2; ++s)
    for (Index b = 0; b < 2; ++b)
      for (Index p = 0; p < 2; ++p)
        for (Index r = 0; r < 2; ++r)
          for (Index c = 0; c < 2; ++c) {
            tia += a.get(ts.idx + s,
                         tb.idx + b,
                         tp.idx + p,
                         tr.idx + r,
                         tc.idx + c) *
                   itw.get(iti);
            ++iti;
          }
 
  return tia;
}
 
 
Numeric interp(ConstVectorView itw,
               ConstTensor6View a,
               const GridPos& tv,
               const GridPos& ts,
               const GridPos& tb,
               const GridPos& tp,
               const GridPos& tr,
               const GridPos& tc) {
  assert(is_size(itw, 64));  // We need 64 interpolation
                             // weights.
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one.
  assert(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
 
  // To store interpolated value:
  Numeric tia = 0;
 
  Index iti = 0;
  for (Index v = 0; v < 2; ++v)
    for (Index s = 0; s < 2; ++s)
      for (Index b = 0; b < 2; ++b)
        for (Index p = 0; p < 2; ++p)
          for (Index r = 0; r < 2; ++r)
            for (Index c = 0; c < 2; ++c) {
              tia += a.get(tv.idx + v,
                           ts.idx + s,
                           tb.idx + b,
                           tp.idx + p,
                           tr.idx + r,
                           tc.idx + c) *
                     itw.get(iti);
              ++iti;
            }
 
  return tia;
}
 
//                      Blue interpolation
 
 
void interpweights(MatrixView itw, const ArrayOfGridPos& cgp) {
  Index n = cgp.nelem();
  assert(is_size(itw, n, 2));  // We must store 2 interpolation
                               // weights for each position.
 
  // We have to loop all the points in the sequence:
  for (Index i = 0; i < n; ++i) {
    // Current grid positions:
    const GridPos& tc = cgp[i];
 
    // Interpolation weights are stored in this order (l=lower
    // u=upper, c=column):
    // 1. l-c
    // 2. u-c
 
    Index iti = 0;
 
    // We could use a straight-forward for loop here:
    //
    //       for ( Index c=1; c>=0; --c )
    //        {
    //          ti[iti] = tc.fd[c];
    //          ++iti;
    //        }
    //
    // However, it is slightly faster to use pointers (I tried it,
    // the speed gain is about 20%). So we should write something
    // like:
    //
    //       for ( const Numeric* c=&tc.fd[1]; c>=&tc.fd[0]; --c )
    //        {
    //          ti[iti] = *c;
    //          ++iti;
    //        }
    //
    // For higher dimensions we have to nest these loops. To avoid
    // typos and safe typing, I use the LOOPIT macro, which expands
    // to the for loop above. Note: NO SEMICOLON AFTER THE LOOPIT
    // COMMAND!
 
    LOOPIT(c) {
      itw.get(i, iti) = *c;
      ++iti;
    }
  }
}
 
 
void interpweights(MatrixView itw,
                   const ArrayOfGridPos& rgp,
                   const ArrayOfGridPos& cgp) {
  Index n = cgp.nelem();
  assert(is_size(rgp, n));     // rgp must have same size as cgp.
  assert(is_size(itw, n, 4));  // We must store 4 interpolation
                               // weights for each position.
 
  // We have to loop all the points in the sequence:
  for (Index i = 0; i < n; ++i) {
    // Current grid positions:
    const GridPos& tr = rgp[i];
    const GridPos& tc = cgp[i];
 
    // Interpolation weights are stored in this order (l=lower
    // u=upper, r=row, c=column):
    // 1. l-r l-c
    // 2. l-r u-c
    // 3. u-r l-c
    // 4. u-r u-c
 
    Index iti = 0;
 
    LOOPIT(r)
    LOOPIT(c) {
      itw.get(i, iti) = (*r) * (*c);
      ++iti;
    }
  }
}
 
 
void interpweights(MatrixView itw,
                   const ArrayOfGridPos& pgp,
                   const ArrayOfGridPos& rgp,
                   const ArrayOfGridPos& cgp) {
  Index n = cgp.nelem();
  assert(is_size(pgp, n));     // pgp must have same size as cgp.
  assert(is_size(rgp, n));     // rgp must have same size as cgp.
  assert(is_size(itw, n, 8));  // We must store 8 interpolation
                               // weights for each position.
 
  // We have to loop all the points in the sequence:
  for (Index i = 0; i < n; ++i) {
    // Current grid positions:
    const GridPos& tp = pgp[i];
    const GridPos& tr = rgp[i];
    const GridPos& tc = cgp[i];
 
    Index iti = 0;
    LOOPIT(p)
    LOOPIT(r)
    LOOPIT(c) {
      itw.get(i, iti) = (*p) * (*r) * (*c);
      ++iti;
    }
  }
}
 
 
void interpweights(MatrixView itw,
                   const ArrayOfGridPos& bgp,
                   const ArrayOfGridPos& pgp,
                   const ArrayOfGridPos& rgp,
                   const ArrayOfGridPos& cgp) {
  Index n = cgp.nelem();
  assert(is_size(bgp, n));      // bgp must have same size as cgp.
  assert(is_size(pgp, n));      // pgp must have same size as cgp.
  assert(is_size(rgp, n));      // rgp must have same size as cgp.
  assert(is_size(itw, n, 16));  // We must store 16 interpolation
                                // weights for each position.
 
  // We have to loop all the points in the sequence:
  for (Index i = 0; i < n; ++i) {
    // Current grid positions:
    const GridPos& tb = bgp[i];
    const GridPos& tp = pgp[i];
    const GridPos& tr = rgp[i];
    const GridPos& tc = cgp[i];
 
    Index iti = 0;
    LOOPIT(b)
    LOOPIT(p)
    LOOPIT(r)
    LOOPIT(c) {
      itw.get(i, iti) = (*b) * (*p) * (*r) * (*c);
      ++iti;
    }
  }
}
 
 
void interpweights(MatrixView itw,
                   const ArrayOfGridPos& sgp,
                   const ArrayOfGridPos& bgp,
                   const ArrayOfGridPos& pgp,
                   const ArrayOfGridPos& rgp,
                   const ArrayOfGridPos& cgp) {
  Index n = cgp.nelem();
  assert(is_size(sgp, n));      // sgp must have same size as cgp.
  assert(is_size(bgp, n));      // bgp must have same size as cgp.
  assert(is_size(pgp, n));      // pgp must have same size as cgp.
  assert(is_size(rgp, n));      // rgp must have same size as cgp.
  assert(is_size(itw, n, 32));  // We must store 32 interpolation
                                // weights for each position.
 
  // We have to loop all the points in the sequence:
  for (Index i = 0; i < n; ++i) {
    // Current grid positions:
    const GridPos& ts = sgp[i];
    const GridPos& tb = bgp[i];
    const GridPos& tp = pgp[i];
    const GridPos& tr = rgp[i];
    const GridPos& tc = cgp[i];
 
    Index iti = 0;
    LOOPIT(s)
    LOOPIT(b)
    LOOPIT(p)
    LOOPIT(r)
    LOOPIT(c) {
      itw.get(i, iti) = (*s) * (*b) * (*p) * (*r) * (*c);
      ++iti;
    }
  }
}
 
 
void interpweights(MatrixView itw,
                   const ArrayOfGridPos& vgp,
                   const ArrayOfGridPos& sgp,
                   const ArrayOfGridPos& bgp,
                   const ArrayOfGridPos& pgp,
                   const ArrayOfGridPos& rgp,
                   const ArrayOfGridPos& cgp) {
  Index n = cgp.nelem();
  assert(is_size(vgp, n));      // vgp must have same size as cgp.
  assert(is_size(sgp, n));      // sgp must have same size as cgp.
  assert(is_size(bgp, n));      // bgp must have same size as cgp.
  assert(is_size(pgp, n));      // pgp must have same size as cgp.
  assert(is_size(rgp, n));      // rgp must have same size as cgp.
  assert(is_size(itw, n, 64));  // We must store 64 interpolation
                                // weights for each position.
 
  // We have to loop all the points in the sequence:
  for (Index i = 0; i < n; ++i) {
    // Current grid positions:
    const GridPos& tv = vgp[i];
    const GridPos& ts = sgp[i];
    const GridPos& tb = bgp[i];
    const GridPos& tp = pgp[i];
    const GridPos& tr = rgp[i];
    const GridPos& tc = cgp[i];
 
    Index iti = 0;
    LOOPIT(v)
    LOOPIT(s)
    LOOPIT(b)
    LOOPIT(p)
    LOOPIT(r)
    LOOPIT(c) {
      itw.get(i, iti) = (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
      ++iti;
    }
  }
}
 
 
void interp(VectorView ia,
            ConstMatrixView itw,
            ConstVectorView a,
            const ArrayOfGridPos& cgp) {
  Index n = cgp.nelem();
  assert(is_size(ia, n));      //  ia must have same size as cgp.
  assert(is_size(itw, n, 2));  // We need 2 interpolation
                               // weights for each position.
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one. We
  // only check the first element.
  assert(
      is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
 
  // We have to loop all the points in the sequence:
  for (Index i = 0; i < n; ++i) {
    // Current grid positions:
    const GridPos& tc = cgp[i];
 
    // Get handle to current element of output vector and initialize
    // it to zero:
    Numeric& tia = ia[i];
    tia = 0;
 
    Index iti = 0;
    for (Index c = 0; c < 2; ++c) {
      assert(tc.idx + c < a.nelem());  // Temporary !?
      tia += a.get(tc.idx + c) * itw.get(i, iti);
      ++iti;
    }
  }
}
 
 
void interp(VectorView ia,
            ConstMatrixView itw,
            ConstMatrixView a,
            const ArrayOfGridPos& rgp,
            const ArrayOfGridPos& cgp) {
  Index n = cgp.nelem();
  assert(is_size(ia, n));      //  ia must have same size as cgp.
  assert(is_size(rgp, n));     // rgp must have same size as cgp.
  assert(is_size(itw, n, 4));  // We need 4 interpolation
                               // weights for each position.
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one. We
  // only check the first element.
  assert(
      is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
 
  // We have to loop all the points in the sequence:
  for (Index i = 0; i < n; ++i) {
    // Current grid positions:
    const GridPos& tr = rgp[i];
    const GridPos& tc = cgp[i];
 
    // Get handle to current element of output vector and initialize
    // it to zero:
    Numeric& tia = ia[i];
    tia = 0;
 
    Index iti = 0;
    for (Index r = 0; r < 2; ++r)
      for (Index c = 0; c < 2; ++c) {
        tia += a.get(tr.idx + r, tc.idx + c) * itw.get(i, iti);
        ++iti;
      }
  }
}
 
 
void interp(VectorView ia,
            ConstMatrixView itw,
            ConstTensor3View a,
            const ArrayOfGridPos& pgp,
            const ArrayOfGridPos& rgp,
            const ArrayOfGridPos& cgp) {
  Index n = cgp.nelem();
  assert(is_size(ia, n));      //  ia must have same size as cgp.
  assert(is_size(pgp, n));     // pgp must have same size as cgp.
  assert(is_size(rgp, n));     // rgp must have same size as cgp.
  assert(is_size(itw, n, 8));  // We need 8 interpolation
                               // weights for each position.
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one. We
  // only check the first element.
  assert(
      is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
 
  // We have to loop all the points in the sequence:
  for (Index i = 0; i < n; ++i) {
    // Current grid positions:
    const GridPos& tp = pgp[i];
    const GridPos& tr = rgp[i];
    const GridPos& tc = cgp[i];
 
    // Get handle to current element of output vector and initialize
    // it to zero:
    Numeric& tia = ia[i];
    tia = 0;
 
    Index iti = 0;
    for (Index p = 0; p < 2; ++p)
      for (Index r = 0; r < 2; ++r)
        for (Index c = 0; c < 2; ++c) {
          tia += a.get(tp.idx + p, tr.idx + r, tc.idx + c) * itw.get(i, iti);
          ++iti;
        }
  }
}
 
 
void interp(VectorView ia,
            ConstMatrixView itw,
            ConstTensor4View a,
            const ArrayOfGridPos& bgp,
            const ArrayOfGridPos& pgp,
            const ArrayOfGridPos& rgp,
            const ArrayOfGridPos& cgp) {
  Index n = cgp.nelem();
  assert(is_size(ia, n));       //  ia must have same size as cgp.
  assert(is_size(bgp, n));      // bgp must have same size as cgp.
  assert(is_size(pgp, n));      // pgp must have same size as cgp.
  assert(is_size(rgp, n));      // rgp must have same size as cgp.
  assert(is_size(itw, n, 16));  // We need 16 interpolation
                                // weights for each position.
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one. We
  // only check the first element.
  assert(
      is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
 
  // We have to loop all the points in the sequence:
  for (Index i = 0; i < n; ++i) {
    // Current grid positions:
    const GridPos& tb = bgp[i];
    const GridPos& tp = pgp[i];
    const GridPos& tr = rgp[i];
    const GridPos& tc = cgp[i];
 
    // Get handle to current element of output vector and initialize
    // it to zero:
    Numeric& tia = ia[i];
    tia = 0;
 
    Index iti = 0;
    for (Index b = 0; b < 2; ++b)
      for (Index p = 0; p < 2; ++p)
        for (Index r = 0; r < 2; ++r)
          for (Index c = 0; c < 2; ++c) {
            tia += a.get(tb.idx + b, tp.idx + p, tr.idx + r, tc.idx + c) *
                   itw.get(i, iti);
            ++iti;
          }
  }
}
 
 
void interp(VectorView ia,
            ConstMatrixView itw,
            ConstTensor5View a,
            const ArrayOfGridPos& sgp,
            const ArrayOfGridPos& bgp,
            const ArrayOfGridPos& pgp,
            const ArrayOfGridPos& rgp,
            const ArrayOfGridPos& cgp) {
  Index n = cgp.nelem();
  assert(is_size(ia, n));       //  ia must have same size as cgp.
  assert(is_size(sgp, n));      // sgp must have same size as cgp.
  assert(is_size(bgp, n));      // bgp must have same size as cgp.
  assert(is_size(pgp, n));      // pgp must have same size as cgp.
  assert(is_size(rgp, n));      // rgp must have same size as cgp.
  assert(is_size(itw, n, 32));  // We need 32 interpolation
                                // weights for each position.
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one. We
  // only check the first element.
  assert(
      is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
 
  // We have to loop all the points in the sequence:
  for (Index i = 0; i < n; ++i) {
    // Current grid positions:
    const GridPos& ts = sgp[i];
    const GridPos& tb = bgp[i];
    const GridPos& tp = pgp[i];
    const GridPos& tr = rgp[i];
    const GridPos& tc = cgp[i];
 
    // Get handle to current element of output vector and initialize
    // it to zero:
    Numeric& tia = ia[i];
    tia = 0;
 
    Index iti = 0;
    for (Index s = 0; s < 2; ++s)
      for (Index b = 0; b < 2; ++b)
        for (Index p = 0; p < 2; ++p)
          for (Index r = 0; r < 2; ++r)
            for (Index c = 0; c < 2; ++c) {
              tia += a.get(ts.idx + s,
                           tb.idx + b,
                           tp.idx + p,
                           tr.idx + r,
                           tc.idx + c) *
                     itw.get(i, iti);
              ++iti;
            }
  }
}
 
 
void interp(VectorView ia,
            ConstMatrixView itw,
            ConstTensor6View a,
            const ArrayOfGridPos& vgp,
            const ArrayOfGridPos& sgp,
            const ArrayOfGridPos& bgp,
            const ArrayOfGridPos& pgp,
            const ArrayOfGridPos& rgp,
            const ArrayOfGridPos& cgp) {
  Index n = cgp.nelem();
  assert(is_size(ia, n));       //  ia must have same size as cgp.
  assert(is_size(vgp, n));      // vgp must have same size as cgp.
  assert(is_size(sgp, n));      // sgp must have same size as cgp.
  assert(is_size(bgp, n));      // bgp must have same size as cgp.
  assert(is_size(pgp, n));      // pgp must have same size as cgp.
  assert(is_size(rgp, n));      // rgp must have same size as cgp.
  assert(is_size(itw, n, 64));  // We need 64 interpolation
                                // weights for each position.
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one. We
  // only check the first element.
  assert(
      is_same_within_epsilon(itw(0, Range(joker)).sum(), 1, sum_check_epsilon));
 
  // We have to loop all the points in the sequence:
  for (Index i = 0; i < n; ++i) {
    // Current grid positions:
    const GridPos& tv = vgp[i];
    const GridPos& ts = sgp[i];
    const GridPos& tb = bgp[i];
    const GridPos& tp = pgp[i];
    const GridPos& tr = rgp[i];
    const GridPos& tc = cgp[i];
 
    // Get handle to current element of output vector and initialize
    // it to zero:
    Numeric& tia = ia[i];
    tia = 0;
 
    Index iti = 0;
    for (Index v = 0; v < 2; ++v)
      for (Index s = 0; s < 2; ++s)
        for (Index b = 0; b < 2; ++b)
          for (Index p = 0; p < 2; ++p)
            for (Index r = 0; r < 2; ++r)
              for (Index c = 0; c < 2; ++c) {
                tia += a.get(tv.idx + v,
                             ts.idx + s,
                             tb.idx + b,
                             tp.idx + p,
                             tr.idx + r,
                             tc.idx + c) *
                       itw.get(i, iti);
                ++iti;
              }
  }
}
 
//                      Green interpolation
 
 
void interpweights(Tensor3View itw,
                   const ArrayOfGridPos& rgp,
                   const ArrayOfGridPos& cgp) {
  Index nr = rgp.nelem();
  Index nc = cgp.nelem();
  assert(is_size(itw, nr, nc, 4));  // We must store 4 interpolation
                                    // weights for each position.
 
  // We have to loop all the points in the new grid:
  for (Index ir = 0; ir < nr; ++ir) {
    // Current grid position:
    const GridPos& tr = rgp[ir];
 
    for (Index ic = 0; ic < nc; ++ic) {
      // Current grid position:
      const GridPos& tc = cgp[ic];
 
      // Interpolation weights are stored in this order (l=lower
      // u=upper, r=row, c=column):
      // 1. l-r l-c
      // 2. l-r u-c
      // 3. u-r l-c
      // 4. u-r u-c
 
      Index iti = 0;
 
      LOOPIT(r)
      LOOPIT(c) {
        itw.get(ir, ic, iti) = (*r) * (*c);
        ++iti;
      }
    }
  }
}
 
 
void interpweights(Tensor4View itw,
                   const ArrayOfGridPos& pgp,
                   const ArrayOfGridPos& rgp,
                   const ArrayOfGridPos& cgp) {
  Index np = pgp.nelem();
  Index nr = rgp.nelem();
  Index nc = cgp.nelem();
  // We must store 8 interpolation weights for each position:
  assert(is_size(itw, np, nr, nc, 8));
 
  // We have to loop all the points in the new grid:
  for (Index ip = 0; ip < np; ++ip) {
    const GridPos& tp = pgp[ip];
    for (Index ir = 0; ir < nr; ++ir) {
      const GridPos& tr = rgp[ir];
      for (Index ic = 0; ic < nc; ++ic) {
        const GridPos& tc = cgp[ic];
 
        Index iti = 0;
 
        LOOPIT(p)
        LOOPIT(r)
        LOOPIT(c) {
          itw.get(ip, ir, ic, iti) = (*p) * (*r) * (*c);
          ++iti;
        }
      }
    }
  }
}
 
 
void interpweights(Tensor5View itw,
                   const ArrayOfGridPos& bgp,
                   const ArrayOfGridPos& pgp,
                   const ArrayOfGridPos& rgp,
                   const ArrayOfGridPos& cgp) {
  Index nb = bgp.nelem();
  Index np = pgp.nelem();
  Index nr = rgp.nelem();
  Index nc = cgp.nelem();
  // We must store 16 interpolation weights for each position:
  assert(is_size(itw, nb, np, nr, nc, 16));
 
  // We have to loop all the points in the new grid:
  for (Index ib = 0; ib < nb; ++ib) {
    const GridPos& tb = bgp[ib];
    for (Index ip = 0; ip < np; ++ip) {
      const GridPos& tp = pgp[ip];
      for (Index ir = 0; ir < nr; ++ir) {
        const GridPos& tr = rgp[ir];
        for (Index ic = 0; ic < nc; ++ic) {
          const GridPos& tc = cgp[ic];
 
          Index iti = 0;
 
          LOOPIT(b)
          LOOPIT(p)
          LOOPIT(r)
          LOOPIT(c) {
            itw.get(ib, ip, ir, ic, iti) = (*b) * (*p) * (*r) * (*c);
            ++iti;
          }
        }
      }
    }
  }
}
 
 
void interpweights(Tensor6View itw,
                   const ArrayOfGridPos& sgp,
                   const ArrayOfGridPos& bgp,
                   const ArrayOfGridPos& pgp,
                   const ArrayOfGridPos& rgp,
                   const ArrayOfGridPos& cgp) {
  Index ns = sgp.nelem();
  Index nb = bgp.nelem();
  Index np = pgp.nelem();
  Index nr = rgp.nelem();
  Index nc = cgp.nelem();
  // We must store 32 interpolation weights for each position:
  assert(is_size(itw, ns, nb, np, nr, nc, 32));
 
  // We have to loop all the points in the new grid:
  for (Index is = 0; is < ns; ++is) {
    const GridPos& ts = sgp[is];
    for (Index ib = 0; ib < nb; ++ib) {
      const GridPos& tb = bgp[ib];
      for (Index ip = 0; ip < np; ++ip) {
        const GridPos& tp = pgp[ip];
        for (Index ir = 0; ir < nr; ++ir) {
          const GridPos& tr = rgp[ir];
          for (Index ic = 0; ic < nc; ++ic) {
            const GridPos& tc = cgp[ic];
 
            Index iti = 0;
 
            LOOPIT(s)
            LOOPIT(b)
            LOOPIT(p)
            LOOPIT(r)
            LOOPIT(c) {
              itw.get(is, ib, ip, ir, ic, iti) =
                  (*s) * (*b) * (*p) * (*r) * (*c);
              ++iti;
            }
          }
        }
      }
    }
  }
}
 
 
void interpweights(Tensor7View itw,
                   const ArrayOfGridPos& vgp,
                   const ArrayOfGridPos& sgp,
                   const ArrayOfGridPos& bgp,
                   const ArrayOfGridPos& pgp,
                   const ArrayOfGridPos& rgp,
                   const ArrayOfGridPos& cgp) {
  Index nv = vgp.nelem();
  Index ns = sgp.nelem();
  Index nb = bgp.nelem();
  Index np = pgp.nelem();
  Index nr = rgp.nelem();
  Index nc = cgp.nelem();
  // We must store 64 interpolation weights for each position:
  assert(is_size(itw, nv, ns, nb, np, nr, nc, 64));
 
  // We have to loop all the points in the new grid:
  for (Index iv = 0; iv < nv; ++iv) {
    const GridPos& tv = vgp[iv];
    for (Index is = 0; is < ns; ++is) {
      const GridPos& ts = sgp[is];
      for (Index ib = 0; ib < nb; ++ib) {
        const GridPos& tb = bgp[ib];
        for (Index ip = 0; ip < np; ++ip) {
          const GridPos& tp = pgp[ip];
          for (Index ir = 0; ir < nr; ++ir) {
            const GridPos& tr = rgp[ir];
            for (Index ic = 0; ic < nc; ++ic) {
              const GridPos& tc = cgp[ic];
 
              Index iti = 0;
 
              LOOPIT(v)
              LOOPIT(s)
              LOOPIT(b)
              LOOPIT(p)
              LOOPIT(r)
              LOOPIT(c) {
                itw.get(iv, is, ib, ip, ir, ic, iti) =
                    (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
                ++iti;
              }
            }
          }
        }
      }
    }
  }
}
 
 
void interp(MatrixView ia,
            ConstTensor3View itw,
            ConstMatrixView a,
            const ArrayOfGridPos& rgp,
            const ArrayOfGridPos& cgp) {
  Index nr = rgp.nelem();
  Index nc = cgp.nelem();
  assert(is_size(ia, nr, nc));
  assert(is_size(itw, nr, nc, 4));  // We need 4 interpolation
                                    // weights for each position.
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one. We
  // only check the first element.
  assert(is_same_within_epsilon(
      itw(0, 0, Range(joker)).sum(), 1, sum_check_epsilon));
 
  // We have to loop all the points in the new grid:
  for (Index ir = 0; ir < nr; ++ir) {
    // Current grid position:
    const GridPos& tr = rgp[ir];
 
    for (Index ic = 0; ic < nc; ++ic) {
      // Current grid position:
      const GridPos& tc = cgp[ic];
 
      // Get handle to current element of output tensor and initialize
      // it to zero:
      Numeric& tia = ia(ir, ic);
      tia = 0;
 
      Index iti = 0;
      for (Index r = 0; r < 2; ++r)
        for (Index c = 0; c < 2; ++c) {
          assert(tr.idx + r < a.nrows());  // Temporary !?
          assert(tc.idx + c < a.ncols());  // Temporary !?
          tia += a.get(tr.idx + r, tc.idx + c) * itw.get(ir, ic, iti);
          ++iti;
        }
    }
  }
}
 
 
void interp(Tensor3View ia,
            ConstTensor4View itw,
            ConstTensor3View a,
            const ArrayOfGridPos& pgp,
            const ArrayOfGridPos& rgp,
            const ArrayOfGridPos& cgp) {
  Index np = pgp.nelem();
  Index nr = rgp.nelem();
  Index nc = cgp.nelem();
  assert(is_size(ia, np, nr, nc));
  assert(is_size(itw, np, nr, nc, 8));
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one. We
  // only check the first element.
  assert(is_same_within_epsilon(
      itw(0, 0, 0, Range(joker)).sum(), 1, sum_check_epsilon));
 
  // We have to loop all the points in the new grid:
  for (Index ip = 0; ip < np; ++ip) {
    const GridPos& tp = pgp[ip];
    for (Index ir = 0; ir < nr; ++ir) {
      const GridPos& tr = rgp[ir];
      for (Index ic = 0; ic < nc; ++ic) {
        // Current grid position:
        const GridPos& tc = cgp[ic];
 
        // Get handle to current element of output tensor and
        // initialize it to zero:
        Numeric& tia = ia(ip, ir, ic);
        tia = 0;
 
        Index iti = 0;
        for (Index p = 0; p < 2; ++p)
          for (Index r = 0; r < 2; ++r)
            for (Index c = 0; c < 2; ++c) {
              assert(tp.idx + p < a.npages());  // Temporary !?
              assert(tr.idx + r < a.nrows());   // Temporary !?
              assert(tc.idx + c < a.ncols());   // Temporary !?
              tia += a.get(tp.idx + p, tr.idx + r, tc.idx + c) *
                     itw.get(ip, ir, ic, iti);
              ++iti;
            }
      }
    }
  }
}
 
 
void interp(Tensor4View ia,
            ConstTensor5View itw,
            ConstTensor4View a,
            const ArrayOfGridPos& bgp,
            const ArrayOfGridPos& pgp,
            const ArrayOfGridPos& rgp,
            const ArrayOfGridPos& cgp) {
  Index nb = bgp.nelem();
  Index np = pgp.nelem();
  Index nr = rgp.nelem();
  Index nc = cgp.nelem();
  assert(is_size(ia, nb, np, nr, nc));
  assert(is_size(itw, nb, np, nr, nc, 16));
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one. We
  // only check the first element.
  assert(is_same_within_epsilon(
      itw(0, 0, 0, 0, Range(joker)).sum(), 1, sum_check_epsilon));
 
  // We have to loop all the points in the new grid:
  for (Index ib = 0; ib < nb; ++ib) {
    const GridPos& tb = bgp[ib];
    for (Index ip = 0; ip < np; ++ip) {
      const GridPos& tp = pgp[ip];
      for (Index ir = 0; ir < nr; ++ir) {
        const GridPos& tr = rgp[ir];
        for (Index ic = 0; ic < nc; ++ic) {
          // Current grid position:
          const GridPos& tc = cgp[ic];
 
          // Get handle to current element of output tensor and
          // initialize it to zero:
          Numeric& tia = ia(ib, ip, ir, ic);
          tia = 0;
 
          Index iti = 0;
          for (Index b = 0; b < 2; ++b)
            for (Index p = 0; p < 2; ++p)
              for (Index r = 0; r < 2; ++r)
                for (Index c = 0; c < 2; ++c) {
                  tia += a.get(tb.idx + b, tp.idx + p, tr.idx + r, tc.idx + c) *
                         itw.get(ib, ip, ir, ic, iti);
                  ++iti;
                }
        }
      }
    }
  }
}
 
 
void interp(Tensor5View ia,
            ConstTensor6View itw,
            ConstTensor5View a,
            const ArrayOfGridPos& sgp,
            const ArrayOfGridPos& bgp,
            const ArrayOfGridPos& pgp,
            const ArrayOfGridPos& rgp,
            const ArrayOfGridPos& cgp) {
  Index ns = sgp.nelem();
  Index nb = bgp.nelem();
  Index np = pgp.nelem();
  Index nr = rgp.nelem();
  Index nc = cgp.nelem();
  assert(is_size(ia, ns, nb, np, nr, nc));
  assert(is_size(itw, ns, nb, np, nr, nc, 32));
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one. We
  // only check the first element.
  assert(is_same_within_epsilon(
      itw(0, 0, 0, 0, 0, Range(joker)).sum(), 1, sum_check_epsilon));
 
  // We have to loop all the points in the new grid:
  for (Index is = 0; is < ns; ++is) {
    const GridPos& ts = sgp[is];
    for (Index ib = 0; ib < nb; ++ib) {
      const GridPos& tb = bgp[ib];
      for (Index ip = 0; ip < np; ++ip) {
        const GridPos& tp = pgp[ip];
        for (Index ir = 0; ir < nr; ++ir) {
          const GridPos& tr = rgp[ir];
          for (Index ic = 0; ic < nc; ++ic) {
            // Current grid position:
            const GridPos& tc = cgp[ic];
 
            // Get handle to current element of output tensor and
            // initialize it to zero:
            Numeric& tia = ia(is, ib, ip, ir, ic);
            tia = 0;
 
            Index iti = 0;
            for (Index s = 0; s < 2; ++s)
              for (Index b = 0; b < 2; ++b)
                for (Index p = 0; p < 2; ++p)
                  for (Index r = 0; r < 2; ++r)
                    for (Index c = 0; c < 2; ++c) {
                      tia += a.get(ts.idx + s,
                                   tb.idx + b,
                                   tp.idx + p,
                                   tr.idx + r,
                                   tc.idx + c) *
                             itw.get(is, ib, ip, ir, ic, iti);
                      ++iti;
                    }
          }
        }
      }
    }
  }
}
 
 
void interp(Tensor6View ia,
            ConstTensor7View itw,
            ConstTensor6View a,
            const ArrayOfGridPos& vgp,
            const ArrayOfGridPos& sgp,
            const ArrayOfGridPos& bgp,
            const ArrayOfGridPos& pgp,
            const ArrayOfGridPos& rgp,
            const ArrayOfGridPos& cgp) {
  Index nv = vgp.nelem();
  Index ns = sgp.nelem();
  Index nb = bgp.nelem();
  Index np = pgp.nelem();
  Index nr = rgp.nelem();
  Index nc = cgp.nelem();
  assert(is_size(ia, nv, ns, nb, np, nr, nc));
  assert(is_size(itw, nv, ns, nb, np, nr, nc, 64));
 
  // Check that interpolation weights are valid. The sum of all
  // weights (last dimension) must always be approximately one. We
  // only check the first element.
  assert(is_same_within_epsilon(
      itw(0, 0, 0, 0, 0, 0, Range(joker)).sum(), 1, sum_check_epsilon));
 
  // We have to loop all the points in the new grid:
  for (Index iv = 0; iv < nv; ++iv) {
    const GridPos& tv = vgp[iv];
    for (Index is = 0; is < ns; ++is) {
      const GridPos& ts = sgp[is];
      for (Index ib = 0; ib < nb; ++ib) {
        const GridPos& tb = bgp[ib];
        for (Index ip = 0; ip < np; ++ip) {
          const GridPos& tp = pgp[ip];
          for (Index ir = 0; ir < nr; ++ir) {
            const GridPos& tr = rgp[ir];
            for (Index ic = 0; ic < nc; ++ic) {
              // Current grid position:
              const GridPos& tc = cgp[ic];
 
              // Get handle to current element of output tensor and
              // initialize it to zero:
              Numeric& tia = ia(iv, is, ib, ip, ir, ic);
              tia = 0;
 
              Index iti = 0;
              for (Index v = 0; v < 2; ++v)
                for (Index s = 0; s < 2; ++s)
                  for (Index b = 0; b < 2; ++b)
                    for (Index p = 0; p < 2; ++p)
                      for (Index r = 0; r < 2; ++r)
                        for (Index c = 0; c < 2; ++c) {
                          tia += a.get(tv.idx + v,
                                       ts.idx + s,
                                       tb.idx + b,
                                       tp.idx + p,
                                       tr.idx + r,
                                       tc.idx + c) *
                                 itw.get(iv, is, ib, ip, ir, ic, iti);
                          ++iti;
                        }
            }
          }
        }
      }
    }
  }
}
 
 
Numeric interp_poly(ConstVectorView x,
                    ConstVectorView y,
                    const Numeric& x_i,
                    const GridPos& gp) {
  Index N_x = x.nelem();
 
  assert(N_x == y.nelem());
  assert(N_x > 2);
 
  Vector xa(4), ya(4);
  Numeric y_int;
  Numeric dy_int;
  y_int = 0.;
 
  // 1 - polynomial interpolation (3 points) with grid position search
  // 2 - polynomial interpolation (3 points) without grid position search
  // 3 - polynomial interpolation (4 points)
 
  Index interp_method = 1;
 
  if (interp_method == 1) {
    // Pick out three points for interpolation
    if ((gp.fd[0] <= 0.5 && gp.idx > 0) || gp.idx == N_x - 2) {
      xa[0] = x[gp.idx - 1];
      xa[1] = x[gp.idx];
      xa[2] = x[gp.idx + 1];
 
      ya[0] = y[gp.idx - 1];
      ya[1] = y[gp.idx];
      ya[2] = y[gp.idx + 1];
    }
 
    else if ((gp.fd[0] > 0.5 && gp.idx < N_x - 2) || gp.idx == 0) {
      xa[0] = x[gp.idx];
      xa[1] = x[gp.idx + 1];
      xa[2] = x[gp.idx + 2];
 
      ya[0] = y[gp.idx];
      ya[1] = y[gp.idx + 1];
      ya[2] = y[gp.idx + 2];
    }
 
    else if (gp.idx == N_x - 1) {
      xa[0] = x[N_x - 2];
      xa[1] = x[N_x - 1];
      xa[2] = x[N_x];
 
      ya[0] = y[N_x - 2];
      ya[1] = y[N_x - 1];
      ya[2] = y[N_x];
    } else {
      assert(false);
      arts_exit();
    }
 
    polint(y_int, dy_int, xa, ya, 3, x_i);
 
  }
 
  else if (interp_method == 2) {
    if (gp.idx == 0) {
      xa[0] = x[gp.idx];
      xa[1] = x[gp.idx + 1];
      xa[2] = x[gp.idx + 2];
 
      ya[0] = y[gp.idx];
      ya[1] = y[gp.idx + 1];
      ya[2] = y[gp.idx + 2];
    } else if (gp.idx == N_x - 1) {
      xa[0] = x[gp.idx - 2];
      xa[1] = x[gp.idx - 1];
      xa[2] = x[gp.idx];
 
      ya[0] = y[gp.idx - 2];
      ya[1] = y[gp.idx - 1];
      ya[2] = y[gp.idx];
    } else {
      xa[0] = x[gp.idx - 1];
      xa[1] = x[gp.idx];
      xa[2] = x[gp.idx + 1];
 
      ya[0] = y[gp.idx - 1];
      ya[1] = y[gp.idx];
      ya[2] = y[gp.idx + 1];
    }
 
    // Polynominal interpolation, n = 3
    polint(y_int, dy_int, xa, ya, 3, x_i);
  }
 
  else if (interp_method == 3) {
    // Take 4 points
    if (gp.idx == 0) {
      xa[0] = -x[gp.idx + 1];
      xa[1] = x[gp.idx + 0];
      xa[2] = x[gp.idx + 1];
      xa[3] = x[gp.idx + 2];
 
      ya[0] = y[gp.idx + 1];
      ya[1] = y[gp.idx + 0];
      ya[2] = y[gp.idx + 1];
      ya[3] = y[gp.idx + 2];
    } else if (gp.idx == N_x - 1) {
      xa[0] = x[gp.idx - 1];
      xa[1] = x[gp.idx - 0];
      xa[2] = 2 * x[gp.idx] - x[gp.idx - 1];
      xa[3] = 2 * x[gp.idx] - x[gp.idx - 2];
 
      ya[0] = y[gp.idx - 1];
      ya[1] = y[gp.idx - 0];
      ya[2] = y[gp.idx - 1];
      ya[3] = y[gp.idx - 2];
    } else if (gp.idx == N_x - 2) {
      xa[0] = x[gp.idx - 2];
      xa[1] = x[gp.idx - 1];
      xa[2] = x[gp.idx];
      xa[3] = x[gp.idx + 1];
 
      ya[0] = y[gp.idx - 2];
      ya[1] = y[gp.idx - 1];
      ya[2] = y[gp.idx];
      ya[3] = y[gp.idx + 1];
    } else {
      xa[0] = x[gp.idx - 1];
      xa[1] = x[gp.idx];
      xa[2] = x[gp.idx + 1];
      xa[3] = x[gp.idx + 2];
 
      ya[0] = y[gp.idx - 1];
      ya[1] = y[gp.idx];
      ya[2] = y[gp.idx + 1];
      ya[3] = y[gp.idx + 2];
    }
    // Polinominal interpolation, n = 4
    polint(y_int, dy_int, xa, ya, 4, x_i);
  }
 
  return y_int;
}
 
 
void polint(Numeric& y_int,
            Numeric& dy_int,
            ConstVectorView xa,
            ConstVectorView ya,
            const Index& n,
            const Numeric& x) {
  Index ns = 1;
  Numeric den, dif, dift, ho, hp, w;
 
  dif = abs(x - xa[0]);
 
  Vector c(n);
  Vector d(n);
 
  // Find the index of the closest table entry
  for (Index i = 0; i < n; i++) {
    if ((dift = abs(x - xa[i])) < dif) {
      ns = i;
      dif = dift;
    }
    // Initialize c and d
    c[i] = ya[i];
    d[i] = ya[i];
  }
  // Initial approximation to y
  y_int = ya[ns--];
 
  for (Index m = 1; m < n; m++) {
    for (Index i = 0; i < n - m; i++) {
      ho = xa[i] - x;
      hp = xa[i + m] - x;
      w = c[i + 1] - d[i];
      den = ho - hp;
      // This error occurs when two input xa's are identical.
      assert(den != 0.);
      den = w / den;
      d[i] = hp * den;
      c[i] = ho * den;
    }
    y_int += (dy_int = (2 * (ns + 1) < (n - m) ? c[ns + 1] : d[ns--]));
  }
}