interpolation_poly.cc

Go to the documentation of this file.

/* Copyright (C) 2008 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 "interpolation_poly.h"
#include "interpolation.h"
#include "logic.h"

// These two macros here are from Numerical Receipes. They give the
// maxium and minimum of two Index arguments.
static Index imaxarg1, imaxarg2;
#define IMAX(a,b) (imaxarg1=(a), imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\
                   (imaxarg1) : (imaxarg2))

static Index iminarg1, iminarg2;
#define IMIN(a,b) (iminarg1=(a), iminarg2=(b),(iminarg1) < (iminarg2) ?\
                   (iminarg1) : (iminarg2))

// 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);
  gridpos( gp_trad, old_grid, new_grid, extpolfac );


  // This loop is parallelized. I thought this could help a bit, since
  // we spend some time waiting for new memory being allocated by
  // resize. Dunno if the benefit is very large, though.

#pragma omp parallel 
#pragma omp for 
  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 is is for 1-based arrays.
      Index k = IMIN(IMAX(gp_trad[s].idx-(m-1)/2, 0),
                     n_old-m);

      //      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 fron 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);
  Vector          v( 1, new_grid );
  gridpos_poly( agp, old_grid, v, order, extpolfac );
  gp = agp[0];  
}






// FIXME: Below here is copied code from interpolation.cc that must be adapted.


//#define LOOPW(x) for ( const Numeric* x=&t##x.fd[1]; x>=&t##x.fd[0]; --x )
//#define LOOPW(x) for ( Index x=0; x<t##x.w.nelem(); ++x )
#define LOOPW(x) for ( ConstIterator1D x=t##x.w.begin(); x!=t##x.w.end(); ++x )


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



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;

  LOOPW(c)
    {
      itw[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;

  LOOPW(r)
  LOOPW(c)
    {
      itw[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;

  LOOPW(p)
  LOOPW(r)
  LOOPW(c)
    {
      itw[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;

  LOOPW(b)
  LOOPW(p)
  LOOPW(r)
  LOOPW(c)
    {
      itw[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;

  LOOPW(s)
  LOOPW(b)
  LOOPW(p)
  LOOPW(r)
  LOOPW(c)
    {
      itw[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;

  LOOPW(v)
  LOOPW(s)
  LOOPW(b)
  LOOPW(p)
  LOOPW(r)
  LOOPW(c)
    {
      itw[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;
  LOOPIDX(c)
    {
      tia += a[*c] * itw[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;
  LOOPIDX(r)
    LOOPIDX(c)
      {
        tia += a(*r,
                 *c) * itw[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;
  LOOPIDX(p)
  LOOPIDX(r)
  LOOPIDX(c)
        {
          tia += a(*p,
                   *r,
                   *c) * itw[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;
  LOOPIDX(b)
  LOOPIDX(p)
  LOOPIDX(r)
  LOOPIDX(c)
          {
            tia += a(*b,
                     *p,
                     *r,
                     *c) * itw[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;
  LOOPIDX(s)
  LOOPIDX(b)
  LOOPIDX(p)
  LOOPIDX(r)
  LOOPIDX(c)
            {
              tia += a(*s,
                       *b,
                       *p,
                       *r,
                       *c) * itw[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;
  LOOPIDX(v)
  LOOPIDX(s)
  LOOPIDX(b)
  LOOPIDX(p)
  LOOPIDX(r)
  LOOPIDX(c)
              {
                tia += a(*v,
                         *s,
                         *b,
                         *p,
                         *r,
                         *c) * itw[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! 

      LOOPW(c)
        {
          itw(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;

      LOOPW(r)
      LOOPW(c)
          {
            itw(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;
      LOOPW(p)
      LOOPW(r)
      LOOPW(c)
        {
          itw(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;
      LOOPW(b)
      LOOPW(p)
      LOOPW(r)
      LOOPW(c)
        {
          itw(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;
      LOOPW(s)
      LOOPW(b)
      LOOPW(p)
      LOOPW(r)
      LOOPW(c)
        {
          itw(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;
      LOOPW(v)
      LOOPW(s)
      LOOPW(b)
      LOOPW(p)
      LOOPW(r)
      LOOPW(c)
        {
          itw(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;
      LOOPIDX(c)
        {
          tia += a[*c] * itw(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;
      LOOPIDX(r)
        LOOPIDX(c)
          {
            tia += a(*r,
                     *c) * itw(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;
      LOOPIDX(p)
      LOOPIDX(r)
      LOOPIDX(c)
            {
              tia += a(*p,
                       *r,
                       *c) * itw(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;
      LOOPIDX(b)
      LOOPIDX(p)
      LOOPIDX(r)
      LOOPIDX(c)
              {
                tia += a(*b,
                         *p,
                         *r,
                         *c) * itw(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;
      LOOPIDX(s)
      LOOPIDX(b)
      LOOPIDX(p)
      LOOPIDX(r)
      LOOPIDX(c)
                {
                  tia += a(*s,
                           *b,
                           *p,
                           *r,
                           *c) * itw(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;
      LOOPIDX(v)
      LOOPIDX(s)
      LOOPIDX(b)
      LOOPIDX(p)
      LOOPIDX(r)
      LOOPIDX(c)
                  {
                    tia += a(*v,
                             *s,
                             *b,
                             *p,
                             *r,
                             *c) * itw(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];

      for ( Index ic=0; ic<nc; ++ic )
        {
          // Current grid position:
          const GridPosPoly& 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;

          LOOPW(r)
            LOOPW(c)
            {
              itw(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];
      for ( Index ir=0; ir<nr; ++ir )
        {
          const GridPosPoly& tr = rgp[ir];
          for ( Index ic=0; ic<nc; ++ic )
            {
              const GridPosPoly& tc = cgp[ic];

              Index iti = 0;

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

                  Index iti = 0;

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

                      Index iti = 0;

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

                          Index iti = 0;

                          LOOPW(v)
                            LOOPW(s)
                            LOOPW(b)
                            LOOPW(p)
                            LOOPW(r)
                            LOOPW(c)
                            {
                              itw(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];

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

---