interpolation_poly.cc

Go to the documentation of this file.

/* Copyright (C) 2012 Stefan Buehler <sbuehler(at)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 <iostream>
#include <cmath>
#include "interpolation_poly.h"
#include "interpolation.h"
#include "logic.h"
#include "check_input.h"
#include "arts_omp.h"



Index IMAX(Index a, Index b)
{
  return a>b ? a:b;
}


Index IMIN(Index a, Index b)
{
  return a<b ? a:b;
}

// File-global constants:


const Numeric sum_check_epsilon = 1e-6;


void gridpos_poly(ArrayOfGridPosPoly& gp,
                  ConstVectorView old_grid,
                  ConstVectorView new_grid,
                  const Index     order,
                  const Numeric&  extpolfac)
{
  // Number of points used in the interpolation (order + 1):
  Index m=order+1;

  const Index n_old = old_grid.nelem();
  const Index n_new = new_grid.nelem();

  // Since we need m interpolation points, the old grid must have a
  // least m elements. 
  assert(n_old>=m);

  // Consistently with gridpos, the array size of gp has to be set
  // outside. Here, we only assert that it is correct:
  assert( is_size(gp,n_new) );
  
  // First call the traditional gridpos to find the grid positions:
  ArrayOfGridPos gp_trad(n_new);
  if (n_old>1)
    {
      gridpos( gp_trad, old_grid, new_grid, extpolfac );
    }
  else if (n_old==1)
    {
      // This case is not handled by traditional gridpos, but is ok for zero
      // order interpolation, so we handle it explicitly here. Since there is
      // only 1 point in the old grid, it is nearest neighbour to all
      // interpolation points.
      for (Index i=0; i<n_new; ++i) {
          gp_trad[i].idx = 0;
          gp_trad[i].fd[0] = 0;
          gp_trad[i].fd[1] = 1;
      }
    }
  else
    {
      // We should never be here.
      assert(false);
    }
    
  for (Index s=0; s<n_new; ++s)
    {
                   
      // Here we calculate the index of the first of the range of
      // points used for interpolation. For linear interpolation this
      // is identical to j. The idea for this expression is from
      // Numerical Receipes (Chapter 3, section "after the hunt"), but
      // there it is for 1-based arrays.
      Index k;
      if (m!=1)
        {
          k = IMIN(IMAX(gp_trad[s].idx-(m-1)/2, 0),
                     n_old-m);
        }
      else
        {
          // The above formula for k is not valid for m==1
          // (nearest neighbour interpolation).
          if (gp_trad[s].fd[0]<=0.5)
              k=gp_trad[s].idx;
          else
              k=gp_trad[s].idx+1;

          // It is a matter of definition what we do with the exact fd==0.5 case.
          // Here I arbitrarily decided to stick with the "left" point (smaller
          // index). I believe this is consistent with the behaviour for m=3,
          // where 2 points on the left and 1 point on the right is used. (So,
          // we always prefer the left side.)
        }

      //      cout << "m: "<< m << ", k: " << k << endl;


      // Make gp[s].idx and gp[s].w the right size:
      gp[s].idx.resize(m);
      gp[s].w.resize(m);
      
      // Calculate w for each interpolation point. In the linear case
      // these are just the fractional distances to each interpolation
      // point. The w here correspond exactly to the terms in front of
      // the yi in Numerical Recipes, 2nd edition, section 3.1,
      // eq. 3.1.1.
      for (Index i=0; i<m; ++i)
        {
          gp[s].idx[i] = k+i;

          //  Numerical Recipes, 2nd edition, section 3.1, eq. 3.1.1.

          // Numerator:
          Numeric num = 1;
          for (Index j=0; j<m; ++j)
            if (j!=i)
              num *= new_grid[s] - old_grid[k+j];
      
          // Denominator:
          Numeric denom = 1;
          for (Index j=0; j<m; ++j)
            if (j!=i)
              denom *= old_grid[k+i] - old_grid[k+j];

          gp[s].w[i] = num / denom;
        }

      // Debugging: Test if sum of all w is 1, as it should be:
//       Numeric testsum = 0;
//       for (Index i=0; i<m; ++i) testsum += gp[s].w[i];
//       cout << "Testsum = " << testsum << endl;        
      
    }
}



void gridpos_poly( GridPosPoly&    gp,
                   ConstVectorView old_grid,
                   const Numeric&  new_grid,
                   const Index     order,
                   const Numeric&  extpolfac )
{
  ArrayOfGridPosPoly  agp(1);
  gridpos_poly( agp, old_grid, new_grid, order, extpolfac );
  gp = agp[0];  
}



void gridpos_poly_longitudinal(const String&   error_msg,
                               ArrayOfGridPosPoly& gp,
                               ConstVectorView old_grid,
                               ConstVectorView new_grid,
                               const Index     order,
                               const Numeric&  extpolfac)
{
    bool global_lons = is_lon_cyclic(old_grid);

    if (global_lons)
        gridpos_poly_cyclic_longitudinal(gp, old_grid, new_grid, order, extpolfac);
    else
    {
        if (old_grid[0] >= new_grid[new_grid.nelem()-1])
        {
            Vector shifted_old_grid = old_grid;
            shifted_old_grid -= 360;
            chk_interpolation_grids(error_msg, shifted_old_grid, new_grid, order, extpolfac);

            gridpos_poly(gp, shifted_old_grid, new_grid, order, extpolfac);
        }
        else if (old_grid[old_grid.nelem()-1] <= new_grid[0])
        {
            Vector shifted_old_grid = old_grid;
            shifted_old_grid += 360;
            chk_interpolation_grids(error_msg, shifted_old_grid, new_grid, order, extpolfac);
            gridpos_poly(gp, shifted_old_grid, new_grid, order, extpolfac);
        }
        else
        {
            chk_interpolation_grids(error_msg, old_grid, new_grid, order, extpolfac);
            gridpos_poly(gp, old_grid, new_grid, order, extpolfac);
        }
    }
}



void gridpos_poly_cyclic_longitudinal(ArrayOfGridPosPoly& gp,
                                      ConstVectorView old_grid,
                                      ConstVectorView new_grid,
                                      const Index     order,
                                      const Numeric&  extpolfac)
{
    assert(is_lon_cyclic(old_grid));

    Index new_n = old_grid.nelem() - 1;
    ConstVectorView old_grid_n1 = old_grid[Range(0, new_n)];
    Range r_left = Range(0, new_n);
    Range r_right = Range(new_n*2, new_n);

    Vector large_grid(new_n*3);

    large_grid[r_left] = old_grid_n1;
    large_grid[r_left] -= 360.;

    large_grid[Range(new_n, new_n)] = old_grid_n1;

    large_grid[r_right] = old_grid_n1;
    large_grid[r_right] += 360.;

    gridpos_poly(gp, large_grid, new_grid, order, extpolfac);

    for (ArrayOfGridPosPoly::iterator itgp = gp.begin(); itgp != gp.end(); itgp++)
        for (ArrayOfIndex::iterator iti = itgp->idx.begin(); iti != itgp->idx.end(); iti++)
            *iti %= new_n;
}




#define LOOPW(x) for ( ConstIterator1D x = t##x##begin; x!=t##x##end; ++x)

#define CACHEW(x) const ConstIterator1D t##x##begin = t##x.w.begin(); \
                  const ConstIterator1D t##x##end = t##x.w.end();


#define LOOPIDX(x) for (ArrayOfIndex::const_iterator x=t##x##begin; x!=t##x##end; ++x)

#define CACHEIDX(x) const ArrayOfIndex::const_iterator t##x##begin = t##x.idx.begin(); \
                    const ArrayOfIndex::const_iterator t##x##end = t##x.idx.end();



ostream& operator<<(ostream& os, const GridPosPoly& gp)
{
  os << "idx: " << gp.idx << "\n";
  os << "w:   " << gp.w   << "\n";

//   cout << "Test iterator:\n";
//   for (ArrayOfIndex::const_iterator x=gp.idx.begin(); x!=gp.idx.end(); ++x)
//     cout << *x << ":";
//   cout << "\n";

  return os;
}






//                      Red Interpolation


void interpweights( VectorView itw,
                    const GridPosPoly& tc )
{
  assert(is_size(itw,tc.w.nelem()));       

  // Interpolation weights are stored in this order (l=lower
  // u=upper, c=column):
  // 1. l-c
  // 2. u-c

  Index iti = 0;

  CACHEW(c)
  LOOPW(c)
    {
      itw.get(iti) = *c;
      ++iti;
    }
}


void interpweights( VectorView itw,
                    const GridPosPoly& tr,
                    const GridPosPoly& tc )
{
  assert(is_size(itw,
                 tr.w.nelem()*
                 tc.w.nelem())); 
  Index iti = 0;

  CACHEW(r)
  CACHEW(c)
  LOOPW(r)
  LOOPW(c)
    {
      itw.get(iti) = (*r) * (*c);
      ++iti;
    }
}


void interpweights( VectorView itw,
                    const GridPosPoly& tp,
                    const GridPosPoly& tr,
                    const GridPosPoly& tc )
{
  assert(is_size(itw,
                 tp.w.nelem()*
                 tr.w.nelem()*
                 tc.w.nelem()));       
                               
  Index iti = 0;

  CACHEW(p)
  CACHEW(r)
  CACHEW(c)
  LOOPW(p)
  LOOPW(r)
  LOOPW(c)
    {
      itw.get(iti) = (*p) * (*r) * (*c);
      ++iti;
    }
}


void interpweights( VectorView itw,
                    const GridPosPoly& tb,
                    const GridPosPoly& tp,
                    const GridPosPoly& tr,
                    const GridPosPoly& tc )
{
  assert(is_size(itw,
                 tb.w.nelem()*
                 tp.w.nelem()*
                 tr.w.nelem()*
                 tc.w.nelem()));      
                                
  Index iti = 0;

  CACHEW(b)
  CACHEW(p)
  CACHEW(r)
  CACHEW(c)
  LOOPW(b)
  LOOPW(p)
  LOOPW(r)
  LOOPW(c)
    {
      itw.get(iti) = (*b) * (*p) * (*r) * (*c);
      ++iti;
    }
}


void interpweights( VectorView itw,
                    const GridPosPoly& ts,
                    const GridPosPoly& tb,
                    const GridPosPoly& tp,
                    const GridPosPoly& tr,
                    const GridPosPoly& tc )
{
  assert(is_size(itw,
                 ts.w.nelem()*
                 tb.w.nelem()*
                 tp.w.nelem()*
                 tr.w.nelem()*
                 tc.w.nelem()));      
                                
  Index iti = 0;

  CACHEW(s)
  CACHEW(b)
  CACHEW(p)
  CACHEW(r)
  CACHEW(c)
  LOOPW(s)
  LOOPW(b)
  LOOPW(p)
  LOOPW(r)
  LOOPW(c)
    {
      itw.get(iti) = (*s) * (*b) * (*p) * (*r) * (*c);
      ++iti;
    }
}


void interpweights( VectorView itw,
                    const GridPosPoly& tv,
                    const GridPosPoly& ts,
                    const GridPosPoly& tb,
                    const GridPosPoly& tp,
                    const GridPosPoly& tr,
                    const GridPosPoly& tc )
{
  assert(is_size(itw,
                 tv.w.nelem()*
                 ts.w.nelem()*
                 tb.w.nelem()*
                 tp.w.nelem()*
                 tr.w.nelem()*
                 tc.w.nelem()));      
                                
  Index iti = 0;

  CACHEW(v)
  CACHEW(s)
  CACHEW(b)
  CACHEW(p)
  CACHEW(r)
  CACHEW(c)
  LOOPW(v)
  LOOPW(s)
  LOOPW(b)
  LOOPW(p)
  LOOPW(r)
  LOOPW(c)
    {
      itw.get(iti) = (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
      ++iti;
    }
}


Numeric interp( ConstVectorView itw,
                ConstVectorView a,    
                const GridPosPoly&  tc )
{
  assert(is_size(itw,tc.w.nelem()));       

  // 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;
  CACHEIDX(c)
  LOOPIDX(c)
    {
      tia += a.get(*c) * itw.get(iti);
      ++iti;
    }

  return tia;
}


Numeric interp( ConstVectorView  itw,
                ConstMatrixView  a,    
                const GridPosPoly&   tr,
                const GridPosPoly&   tc )
{
  assert(is_size(itw,
                 tr.w.nelem()*
                 tc.w.nelem()));       

  // 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;
  CACHEIDX(r)
  CACHEIDX(c)
  LOOPIDX(r)
    LOOPIDX(c)
      {
        tia += a.get(*r,
                     *c) * itw.get(iti);
        ++iti;
      }

  return tia;
}


Numeric interp( ConstVectorView  itw,
                ConstTensor3View a,    
                const GridPosPoly&   tp,
                const GridPosPoly&   tr,
                const GridPosPoly&   tc )
{
  assert(is_size(itw,
                 tp.w.nelem()*
                 tr.w.nelem()*
                 tc.w.nelem()));      

  // 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;
  CACHEIDX(p)
  CACHEIDX(r)
  CACHEIDX(c)
  LOOPIDX(p)
  LOOPIDX(r)
  LOOPIDX(c)
        {
          tia += a.get(*p,
                       *r,
                       *c) * itw.get(iti);
          ++iti;
        }

  return tia;
}


Numeric interp( ConstVectorView  itw,
                ConstTensor4View a,    
                const GridPosPoly&   tb,
                const GridPosPoly&   tp,
                const GridPosPoly&   tr,
                const GridPosPoly&   tc )
{
  assert(is_size(itw,
                 tb.w.nelem()*
                 tp.w.nelem()*
                 tr.w.nelem()*
                 tc.w.nelem()));      

  // 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;
  CACHEIDX(b)
  CACHEIDX(p)
  CACHEIDX(r)
  CACHEIDX(c)
  LOOPIDX(b)
  LOOPIDX(p)
  LOOPIDX(r)
  LOOPIDX(c)
          {
            tia += a.get(*b,
                         *p,
                         *r,
                         *c) * itw.get(iti);
            ++iti;
          }

  return tia;
}


Numeric interp( ConstVectorView  itw,
                ConstTensor5View a,    
                const GridPosPoly&   ts,
                const GridPosPoly&   tb,
                const GridPosPoly&   tp,
                const GridPosPoly&   tr,
                const GridPosPoly&   tc )
{
  assert(is_size(itw,
                 ts.w.nelem()*
                 tb.w.nelem()*
                 tp.w.nelem()*
                 tr.w.nelem()*
                 tc.w.nelem()));      

  // 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;
  CACHEIDX(s)
  CACHEIDX(b)
  CACHEIDX(p)
  CACHEIDX(r)
  CACHEIDX(c)
  LOOPIDX(s)
  LOOPIDX(b)
  LOOPIDX(p)
  LOOPIDX(r)
  LOOPIDX(c)
            {
              tia += a.get(*s,
                           *b,
                           *p,
                           *r,
                           *c) * itw.get(iti);
              ++iti;
            }

  return tia;
}


Numeric interp( ConstVectorView  itw,
                ConstTensor6View a,    
                const GridPosPoly&   tv,
                const GridPosPoly&   ts,
                const GridPosPoly&   tb,
                const GridPosPoly&   tp,
                const GridPosPoly&   tr,
                const GridPosPoly&   tc )
{
  assert(is_size(itw,
                 tv.w.nelem()*
                 ts.w.nelem()*
                 tb.w.nelem()*
                 tp.w.nelem()*
                 tr.w.nelem()*
                 tc.w.nelem()));      

  // 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;
  CACHEIDX(v)
  CACHEIDX(s)
  CACHEIDX(b)
  CACHEIDX(p)
  CACHEIDX(r)
  CACHEIDX(c)
  LOOPIDX(v)
  LOOPIDX(s)
  LOOPIDX(b)
  LOOPIDX(p)
  LOOPIDX(r)
  LOOPIDX(c)
              {
                tia += a.get(*v,
                             *s,
                             *b,
                             *p,
                             *r,
                             *c) * itw.get(iti);
                ++iti;
              }

  return tia;
}



//                      Blue interpolation


void interpweights( MatrixView itw,
                    const ArrayOfGridPosPoly& cgp )
{
  Index n = cgp.nelem();
  assert(is_size(itw,n,
                 cgp[0].w.nelem()));     
                                

  // We have to loop all the points in the sequence:
  for ( Index i=0; i<n; ++i )
    {
      // Current grid positions:
      const GridPosPoly& 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 LOOPW macro, which expands
      // to the for loop above. Note: NO SEMICOLON AFTER THE LOOPW
      // COMMAND! 

      CACHEW(c)
      LOOPW(c)
        {
          itw.get(i,iti) = *c;
          ++iti;
        }
    }
}


void interpweights( MatrixView itw,
                    const ArrayOfGridPosPoly& rgp,
                    const ArrayOfGridPosPoly& cgp )
{
  Index n = cgp.nelem();
  assert(is_size(rgp,n));       // rgp must have same size as cgp.
  assert(is_size(itw,n,
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));     
                                

  // We have to loop all the points in the sequence:
  for ( Index i=0; i<n; ++i )
    {
      // Current grid positions:
      const GridPosPoly& tr = rgp[i];
      const GridPosPoly& 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;

      CACHEW(r)
      CACHEW(c)
      LOOPW(r)
      LOOPW(c)
          {
            itw.get(i,iti) = (*r) * (*c);
            ++iti;
          }
    }
}


void interpweights( MatrixView itw,
                    const ArrayOfGridPosPoly& pgp,
                    const ArrayOfGridPosPoly& rgp,
                    const ArrayOfGridPosPoly& 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,
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));     
                                

  // We have to loop all the points in the sequence:
  for ( Index i=0; i<n; ++i )
    {
      // Current grid positions:
      const GridPosPoly& tp = pgp[i];
      const GridPosPoly& tr = rgp[i];
      const GridPosPoly& tc = cgp[i];

      Index iti = 0;
      CACHEW(p)
      CACHEW(r)
      CACHEW(c)
      LOOPW(p)
      LOOPW(r)
      LOOPW(c)
        {
          itw.get(i,iti) = (*p) * (*r) * (*c);
          ++iti;
        }
    }
}


void interpweights( MatrixView itw,
                    const ArrayOfGridPosPoly& bgp,
                    const ArrayOfGridPosPoly& pgp,
                    const ArrayOfGridPosPoly& rgp,
                    const ArrayOfGridPosPoly& 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,
                 bgp[0].w.nelem()*
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));    
                                

  // We have to loop all the points in the sequence:
  for ( Index i=0; i<n; ++i )
    {
      // Current grid positions:
      const GridPosPoly& tb = bgp[i];
      const GridPosPoly& tp = pgp[i];
      const GridPosPoly& tr = rgp[i];
      const GridPosPoly& tc = cgp[i];

      Index iti = 0;
      CACHEW(b)
      CACHEW(p)
      CACHEW(r)
      CACHEW(c)
      LOOPW(b)
      LOOPW(p)
      LOOPW(r)
      LOOPW(c)
        {
          itw.get(i,iti) = (*b) * (*p) * (*r) * (*c);
          ++iti;
        }
    }
}


void interpweights( MatrixView itw,
                    const ArrayOfGridPosPoly& sgp,
                    const ArrayOfGridPosPoly& bgp,
                    const ArrayOfGridPosPoly& pgp,
                    const ArrayOfGridPosPoly& rgp,
                    const ArrayOfGridPosPoly& 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,
                 sgp[0].w.nelem()*
                 bgp[0].w.nelem()*
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));    
                                

  // We have to loop all the points in the sequence:
  for ( Index i=0; i<n; ++i )
    {
      // Current grid positions:
      const GridPosPoly& ts = sgp[i];
      const GridPosPoly& tb = bgp[i];
      const GridPosPoly& tp = pgp[i];
      const GridPosPoly& tr = rgp[i];
      const GridPosPoly& tc = cgp[i];

      Index iti = 0;
      CACHEW(s)
      CACHEW(b)
      CACHEW(p)
      CACHEW(r)
      CACHEW(c)
      LOOPW(s)
      LOOPW(b)
      LOOPW(p)
      LOOPW(r)
      LOOPW(c)
        {
          itw.get(i,iti) = (*s) * (*b) * (*p) * (*r) * (*c);
          ++iti;
        }
    }
}


void interpweights( MatrixView itw,
                    const ArrayOfGridPosPoly& vgp,
                    const ArrayOfGridPosPoly& sgp,
                    const ArrayOfGridPosPoly& bgp,
                    const ArrayOfGridPosPoly& pgp,
                    const ArrayOfGridPosPoly& rgp,
                    const ArrayOfGridPosPoly& 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,
                 vgp[0].w.nelem()*
                 sgp[0].w.nelem()*
                 bgp[0].w.nelem()*
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));    
                                

  // We have to loop all the points in the sequence:
  for ( Index i=0; i<n; ++i )
    {
      // Current grid positions:
      const GridPosPoly& tv = vgp[i];
      const GridPosPoly& ts = sgp[i];
      const GridPosPoly& tb = bgp[i];
      const GridPosPoly& tp = pgp[i];
      const GridPosPoly& tr = rgp[i];
      const GridPosPoly& tc = cgp[i];

      Index iti = 0;
      CACHEW(v)
      CACHEW(s)
      CACHEW(b)
      CACHEW(p)
      CACHEW(r)
      CACHEW(c)
      LOOPW(v)
      LOOPW(s)
      LOOPW(b)
      LOOPW(p)
      LOOPW(r)
      LOOPW(c)
        {
          itw.get(i,iti) = (*v) * (*s) * (*b) * (*p) * (*r) * (*c);
          ++iti;
        }
    }
}


void interp( VectorView            ia,
             ConstMatrixView       itw,
             ConstVectorView       a,    
             const ArrayOfGridPosPoly& cgp)
{
  Index n = cgp.nelem();
  assert(is_size(ia,n));        //  ia must have same size as cgp.
  assert(is_size(itw,n,
                 cgp[0].w.nelem()));     

  // 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 GridPosPoly& 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;
      CACHEIDX(c)
      LOOPIDX(c)
        {
          tia += a.get(*c) * itw.get(i,iti);
          ++iti;
        }
    }
}


void interp( VectorView            ia,
             ConstMatrixView       itw,
             ConstMatrixView       a,    
             const ArrayOfGridPosPoly& rgp,
             const ArrayOfGridPosPoly& 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,
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));     
  
  // 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 GridPosPoly& tr = rgp[i];
      const GridPosPoly& 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;
      CACHEIDX(r)
      CACHEIDX(c)
      LOOPIDX(r)
        LOOPIDX(c)
          {
            tia += a.get(*r,
                         *c) * itw.get(i,iti);
            ++iti;
          }
    }
}


void interp( VectorView            ia,
             ConstMatrixView       itw,
             ConstTensor3View      a,    
             const ArrayOfGridPosPoly& pgp,
             const ArrayOfGridPosPoly& rgp,
             const ArrayOfGridPosPoly& 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,
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));     
  
  // 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 GridPosPoly& tp = pgp[i];
      const GridPosPoly& tr = rgp[i];
      const GridPosPoly& 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;
      CACHEIDX(p)
      CACHEIDX(r)
      CACHEIDX(c)
      LOOPIDX(p)
      LOOPIDX(r)
      LOOPIDX(c)
            {
              tia += a.get(*p,
                           *r,
                           *c) * itw.get(i,iti);
              ++iti;
            }
    }
}


void interp( VectorView            ia,
             ConstMatrixView       itw,
             ConstTensor4View      a,    
             const ArrayOfGridPosPoly& bgp,
             const ArrayOfGridPosPoly& pgp,
             const ArrayOfGridPosPoly& rgp,
             const ArrayOfGridPosPoly& 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,
                 bgp[0].w.nelem()*
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));    
  
  // 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 GridPosPoly& tb = bgp[i];
      const GridPosPoly& tp = pgp[i];
      const GridPosPoly& tr = rgp[i];
      const GridPosPoly& 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;
      CACHEIDX(b)
      CACHEIDX(p)
      CACHEIDX(r)
      CACHEIDX(c)
      LOOPIDX(b)
      LOOPIDX(p)
      LOOPIDX(r)
      LOOPIDX(c)
              {
                tia += a.get(*b,
                             *p,
                             *r,
                             *c) * itw.get(i,iti);
                ++iti;
              }
    }
}


void interp( VectorView            ia,
             ConstMatrixView       itw,
             ConstTensor5View      a,    
             const ArrayOfGridPosPoly& sgp,
             const ArrayOfGridPosPoly& bgp,
             const ArrayOfGridPosPoly& pgp,
             const ArrayOfGridPosPoly& rgp,
             const ArrayOfGridPosPoly& 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,
                 sgp[0].w.nelem()*
                 bgp[0].w.nelem()*
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));    
  
  // 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 GridPosPoly& ts = sgp[i];
      const GridPosPoly& tb = bgp[i];
      const GridPosPoly& tp = pgp[i];
      const GridPosPoly& tr = rgp[i];
      const GridPosPoly& 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;
      CACHEIDX(s)
      CACHEIDX(b)
      CACHEIDX(p)
      CACHEIDX(r)
      CACHEIDX(c)
      LOOPIDX(s)
      LOOPIDX(b)
      LOOPIDX(p)
      LOOPIDX(r)
      LOOPIDX(c)
                {
                  tia += a.get(*s,
                               *b,
                               *p,
                               *r,
                               *c) * itw.get(i,iti);
                  ++iti;
                }
    }
}


void interp( VectorView            ia,
             ConstMatrixView       itw,
             ConstTensor6View      a,    
             const ArrayOfGridPosPoly& vgp,
             const ArrayOfGridPosPoly& sgp,
             const ArrayOfGridPosPoly& bgp,
             const ArrayOfGridPosPoly& pgp,
             const ArrayOfGridPosPoly& rgp,
             const ArrayOfGridPosPoly& 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,
                 vgp[0].w.nelem()*
                 sgp[0].w.nelem()*
                 bgp[0].w.nelem()*
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));    
  
  // 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 GridPosPoly& tv = vgp[i];
      const GridPosPoly& ts = sgp[i];
      const GridPosPoly& tb = bgp[i];
      const GridPosPoly& tp = pgp[i];
      const GridPosPoly& tr = rgp[i];
      const GridPosPoly& 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;
      CACHEIDX(v)
      CACHEIDX(s)
      CACHEIDX(b)
      CACHEIDX(p)
      CACHEIDX(r)
      CACHEIDX(c)
      LOOPIDX(v)
      LOOPIDX(s)
      LOOPIDX(b)
      LOOPIDX(p)
      LOOPIDX(r)
      LOOPIDX(c)
                  {
                    tia += a.get(*v,
                                 *s,
                                 *b,
                                 *p,
                                 *r,
                                 *c) * itw.get(i,iti);
                    ++iti;
                  }
    }
}

//                      Green interpolation


void interpweights( Tensor3View itw,
                    const ArrayOfGridPosPoly& rgp,
                    const ArrayOfGridPosPoly& cgp )
{
  Index nr = rgp.nelem();
  Index nc = cgp.nelem();
  assert(is_size(itw,nr,nc,
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem())); 
                                

  // We have to loop all the points in the new grid:
  for ( Index ir=0; ir<nr; ++ir )
    {
      // Current grid position:
      const GridPosPoly& tr = rgp[ir];
      CACHEW(r)

      for ( Index ic=0; ic<nc; ++ic )
        {
          // Current grid position:
          const GridPosPoly& tc = cgp[ic];
          CACHEW(c)

          // 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;

          LOOPW(r)
            LOOPW(c)
            {
              itw.get(ir,ic,iti) = (*r) * (*c);
              ++iti;
            }
        }
    }
}


void interpweights( Tensor4View itw,
                    const ArrayOfGridPosPoly& pgp,
                    const ArrayOfGridPosPoly& rgp,
                    const ArrayOfGridPosPoly& cgp )
{
  Index np = pgp.nelem();
  Index nr = rgp.nelem();
  Index nc = cgp.nelem();

  assert(is_size(itw,np,nr,nc,
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));      

  // We have to loop all the points in the new grid:
  for ( Index ip=0; ip<np; ++ip )
    {
      const GridPosPoly& tp = pgp[ip];
      CACHEW(p)
      for ( Index ir=0; ir<nr; ++ir )
        {
          const GridPosPoly& tr = rgp[ir];
          CACHEW(r)
          for ( Index ic=0; ic<nc; ++ic )
            {
              const GridPosPoly& tc = cgp[ic];
              CACHEW(c)

              Index iti = 0;

              LOOPW(p)
                LOOPW(r)
                LOOPW(c)
                {
                  itw.get(ip,ir,ic,iti) =
                    (*p) * (*r) * (*c);
                  ++iti;
                }
            }
        }
    }
}


void interpweights( Tensor5View itw,
                    const ArrayOfGridPosPoly& bgp,
                    const ArrayOfGridPosPoly& pgp,
                    const ArrayOfGridPosPoly& rgp,
                    const ArrayOfGridPosPoly& cgp )
{
  Index nb = bgp.nelem();
  Index np = pgp.nelem();
  Index nr = rgp.nelem();
  Index nc = cgp.nelem();

  assert(is_size(itw,nb,np,nr,nc,
                 bgp[0].w.nelem()*
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));  

  // We have to loop all the points in the new grid:
  for ( Index ib=0; ib<nb; ++ib )
    {
      const GridPosPoly& tb = bgp[ib];
      CACHEW(b)
      for ( Index ip=0; ip<np; ++ip )
        {
          const GridPosPoly& tp = pgp[ip];
          CACHEW(p)
          for ( Index ir=0; ir<nr; ++ir )
            {
              const GridPosPoly& tr = rgp[ir];
              CACHEW(r)
              for ( Index ic=0; ic<nc; ++ic )
                {
                  const GridPosPoly& tc = cgp[ic];
                  CACHEW(c)

                  Index iti = 0;

                  LOOPW(b)
                    LOOPW(p)
                    LOOPW(r)
                    LOOPW(c)
                    {
                      itw.get(ib,ip,ir,ic,iti) =
                        (*b) * (*p) * (*r) * (*c);
                      ++iti;
                    }
                }
            }
        }
    }
}


void interpweights( Tensor6View itw,
                    const ArrayOfGridPosPoly& sgp,
                    const ArrayOfGridPosPoly& bgp,
                    const ArrayOfGridPosPoly& pgp,
                    const ArrayOfGridPosPoly& rgp,
                    const ArrayOfGridPosPoly& cgp )
{
  Index ns = sgp.nelem();
  Index nb = bgp.nelem();
  Index np = pgp.nelem();
  Index nr = rgp.nelem();
  Index nc = cgp.nelem();

  assert(is_size(itw,ns,nb,np,nr,nc,
                 sgp[0].w.nelem()*
                 bgp[0].w.nelem()*
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));       

  // We have to loop all the points in the new grid:
  for ( Index is=0; is<ns; ++is )
    {
      const GridPosPoly& ts = sgp[is];
      CACHEW(s)
      for ( Index ib=0; ib<nb; ++ib )
        {
          const GridPosPoly& tb = bgp[ib];
          CACHEW(b)
          for ( Index ip=0; ip<np; ++ip )
            {
              const GridPosPoly& tp = pgp[ip];
              CACHEW(p)
              for ( Index ir=0; ir<nr; ++ir )
                {
                  const GridPosPoly& tr = rgp[ir];
                  CACHEW(r)
                  for ( Index ic=0; ic<nc; ++ic )
                    {
                      const GridPosPoly& tc = cgp[ic];
                      CACHEW(c)

                      Index iti = 0;

                      LOOPW(s)
                        LOOPW(b)
                        LOOPW(p)
                        LOOPW(r)
                        LOOPW(c)
                        {
                          itw.get(is,ib,ip,ir,ic,iti) =
                            (*s) * (*b) * (*p) * (*r) * (*c);
                          ++iti;
                        }
                    }
                }
            }
        }
    }
}


void interpweights( Tensor7View itw,
                    const ArrayOfGridPosPoly& vgp,
                    const ArrayOfGridPosPoly& sgp,
                    const ArrayOfGridPosPoly& bgp,
                    const ArrayOfGridPosPoly& pgp,
                    const ArrayOfGridPosPoly& rgp,
                    const ArrayOfGridPosPoly& 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(itw,nv,ns,nb,np,nr,nc,
                 vgp[0].w.nelem()*
                 sgp[0].w.nelem()*
                 bgp[0].w.nelem()*
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));    

  // We have to loop all the points in the new grid:
  for ( Index iv=0; iv<nv; ++iv )
    {
      const GridPosPoly& tv = vgp[iv];
      CACHEW(v)
      for ( Index is=0; is<ns; ++is )
        {
          const GridPosPoly& ts = sgp[is];
          CACHEW(s)
          for ( Index ib=0; ib<nb; ++ib )
            {
              const GridPosPoly& tb = bgp[ib];
              CACHEW(b)
              for ( Index ip=0; ip<np; ++ip )
                {
                  const GridPosPoly& tp = pgp[ip];
                  CACHEW(p)
                  for ( Index ir=0; ir<nr; ++ir )
                    {
                      const GridPosPoly& tr = rgp[ir];
                      CACHEW(r)
                      for ( Index ic=0; ic<nc; ++ic )
                        {
                          const GridPosPoly& tc = cgp[ic];
                          CACHEW(c)

                          Index iti = 0;

                          LOOPW(v)
                            LOOPW(s)
                            LOOPW(b)
                            LOOPW(p)
                            LOOPW(r)
                            LOOPW(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 ArrayOfGridPosPoly& rgp,
             const ArrayOfGridPosPoly& cgp)
{
  Index nr = rgp.nelem();
  Index nc = cgp.nelem();
  assert(is_size(ia,nr,nc));    
  assert(is_size(itw,nr,nc,
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem())); 

  // 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 GridPosPoly& tr = rgp[ir];
      CACHEIDX(r)

      for ( Index ic=0; ic<nc; ++ic )
        {
          // Current grid position:
          const GridPosPoly& tc = cgp[ic];
          CACHEIDX(c)

          // Get handle to current element of output tensor and initialize
          // it to zero:
          Numeric& tia = ia(ir,ic);
          tia = 0;

          Index iti = 0;
          LOOPIDX(r)
          LOOPIDX(c)
            {
              tia += a.get(*r,
                           *c) * itw.get(ir,ic,iti);
              ++iti;
            }
        }
    }
}


void interp( Tensor3View           ia,
             ConstTensor4View      itw,
             ConstTensor3View      a,   
             const ArrayOfGridPosPoly& pgp,
             const ArrayOfGridPosPoly& rgp,
             const ArrayOfGridPosPoly& 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,
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));

  // 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 GridPosPoly& tp = pgp[ip];
      CACHEIDX(p);
      for ( Index ir=0; ir<nr; ++ir )
        {
          const GridPosPoly& tr = rgp[ir];
          CACHEIDX(r);
          for ( Index ic=0; ic<nc; ++ic )
            {
              // Current grid position:
              const GridPosPoly& tc = cgp[ic];
              CACHEIDX(c);

              // Get handle to current element of output tensor and
              // initialize it to zero:
              Numeric& tia = ia(ip,ir,ic);
              tia = 0;

              Index iti = 0;
              LOOPIDX(p)
              LOOPIDX(r)
              LOOPIDX(c)
                    {
                      tia += a.get(*p,
                                   *r,
                                   *c) * itw.get(ip,ir,ic,
                                                 iti);
                      ++iti;
                    }
            }
        }
    }
}


void interp( Tensor4View           ia,
             ConstTensor5View      itw,
             ConstTensor4View      a,   
             const ArrayOfGridPosPoly& bgp,
             const ArrayOfGridPosPoly& pgp,
             const ArrayOfGridPosPoly& rgp,
             const ArrayOfGridPosPoly& 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,
                 bgp[0].w.nelem()*
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));

  // 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 GridPosPoly& tb = bgp[ib];
      CACHEIDX(b);
      for ( Index ip=0; ip<np; ++ip )
        {
          const GridPosPoly& tp = pgp[ip];
          CACHEIDX(p);
          for ( Index ir=0; ir<nr; ++ir )
            {
              const GridPosPoly& tr = rgp[ir];
              CACHEIDX(r);
              for ( Index ic=0; ic<nc; ++ic )
                {
                  // Current grid position:
                  const GridPosPoly& tc = cgp[ic];
                  CACHEIDX(c);

                  // 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;
                  LOOPIDX(b)
                  LOOPIDX(p)
                  LOOPIDX(r)
                  LOOPIDX(c)
                          {
                            tia += a.get(*b,
                                         *p,
                                         *r,
                                         *c) * itw.get(ib,ip,ir,ic,
                                                       iti);
                            ++iti;
                          }
                }
            }
        }
    }
}


void interp( Tensor5View           ia,
             ConstTensor6View      itw,
             ConstTensor5View      a,   
             const ArrayOfGridPosPoly& sgp,
             const ArrayOfGridPosPoly& bgp,
             const ArrayOfGridPosPoly& pgp,
             const ArrayOfGridPosPoly& rgp,
             const ArrayOfGridPosPoly& 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,
                 sgp[0].w.nelem()*
                 bgp[0].w.nelem()*
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));

  // 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 GridPosPoly& ts = sgp[is];
      CACHEIDX(s);
      for ( Index ib=0; ib<nb; ++ib )
        {
          const GridPosPoly& tb = bgp[ib];
          CACHEIDX(b);
          for ( Index ip=0; ip<np; ++ip )
            {
              const GridPosPoly& tp = pgp[ip];
              CACHEIDX(p);
              for ( Index ir=0; ir<nr; ++ir )
                {
                  const GridPosPoly& tr = rgp[ir];
                  CACHEIDX(r);
                  for ( Index ic=0; ic<nc; ++ic )
                    {
                      // Current grid position:
                      const GridPosPoly& tc = cgp[ic];
                      CACHEIDX(c);

                      // 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;
                      LOOPIDX(s)
                      LOOPIDX(b)
                      LOOPIDX(p)
                      LOOPIDX(r)
                      LOOPIDX(c)
                                {
                                  tia += a.get(*s,
                                               *b,
                                               *p,
                                               *r,
                                               *c) * itw.get(is,ib,ip,ir,ic,
                                                             iti);
                                  ++iti;
                                }
                    }
                }
            }
        }
    }
}


void interp( Tensor6View           ia,
             ConstTensor7View      itw,
             ConstTensor6View      a,   
             const ArrayOfGridPosPoly& vgp,
             const ArrayOfGridPosPoly& sgp,
             const ArrayOfGridPosPoly& bgp,
             const ArrayOfGridPosPoly& pgp,
             const ArrayOfGridPosPoly& rgp,
             const ArrayOfGridPosPoly& 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,
                 vgp[0].w.nelem()*
                 sgp[0].w.nelem()*
                 bgp[0].w.nelem()*
                 pgp[0].w.nelem()*
                 rgp[0].w.nelem()*
                 cgp[0].w.nelem()));

  // 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 GridPosPoly& tv = vgp[iv];
      CACHEIDX(v);
      for ( Index is=0; is<ns; ++is )
        {
          const GridPosPoly& ts = sgp[is];
          CACHEIDX(s);
          for ( Index ib=0; ib<nb; ++ib )
            {
              const GridPosPoly& tb = bgp[ib];
              CACHEIDX(b);
              for ( Index ip=0; ip<np; ++ip )
                {
                  const GridPosPoly& tp = pgp[ip];
                  CACHEIDX(p);
                  for ( Index ir=0; ir<nr; ++ir )
                    {
                      const GridPosPoly& tr = rgp[ir];
                      CACHEIDX(r);
                      for ( Index ic=0; ic<nc; ++ic )
                        {
                          // Current grid position:
                          const GridPosPoly& tc = cgp[ic];
                          CACHEIDX(c);

                          // 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;
                          LOOPIDX(v)
                          LOOPIDX(s)
                          LOOPIDX(b)
                          LOOPIDX(p)
                          LOOPIDX(r)
                          LOOPIDX(c)
                                      {
                                        tia += a.get(*v,
                                                     *s,
                                                     *b,
                                                     *p,
                                                     *r,
                                                     *c) * itw.get(iv,is,ib,ip,ir,ic,
                                                                   iti);
                                        ++iti;
                                      }
                        }
                    }
                }
            }
        }
    }
}