special_interp.cc

Go to the documentation of this file.

/* Copyright (C) 2002-2008 Stefan Buehler <sbuehler@ltu.se>
 
   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.
 
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */
 
#include <cmath>
#include <iostream>
#include <stdexcept>
#include "auto_md.h"
#include "check_input.h"
#include "math_funcs.h"
#include "messages.h"
#include "special_interp.h"
 
 
 
 
/*===========================================================================
  === Interpolation functions for atmospheric grids and fields
  ===========================================================================*/
 
 
void interp_atmfield_gp2itw( 
              Matrix&           itw, 
        const Index&            atmosphere_dim,
        const ArrayOfGridPos&   gp_p,
        const ArrayOfGridPos&   gp_lat,
        const ArrayOfGridPos&   gp_lon )
{
  const Index n = gp_p.nelem();
 
  if( atmosphere_dim == 1 )
    {
      itw.resize(n,2);
      interpweights( itw, gp_p );
    }
 
  else if( atmosphere_dim == 2 )
    {
      assert( gp_lat.nelem() == n );
      itw.resize(n,4);
      interpweights( itw, gp_p, gp_lat );
    }
 
  else if( atmosphere_dim == 3 )
    {
      assert( gp_lat.nelem() == n );
      assert( gp_lon.nelem() == n );
      itw.resize(n,8);
      interpweights( itw, gp_p, gp_lat, gp_lon );
    }
}
 
 
 
 
void interp_atmfield_by_itw( 
              VectorView        x, 
        const Index&            atmosphere_dim,
        ConstTensor3View        x_field,
        const ArrayOfGridPos&   gp_p,
        const ArrayOfGridPos&   gp_lat,
        const ArrayOfGridPos&   gp_lon,
        ConstMatrixView         itw )
{
  assert( x.nelem() == gp_p.nelem() );
 
  if( atmosphere_dim == 1 )
    { 
      assert( itw.ncols() == 2 );
      interp( x, itw, x_field(Range(joker),0,0), gp_p ); 
    }
 
  else if( atmosphere_dim == 2 )
    { 
      assert( itw.ncols() == 4 );
      interp( x, itw, x_field(Range(joker),Range(joker),0), gp_p, gp_lat ); 
    }
 
  else if( atmosphere_dim == 3 )
    {       
      assert( itw.ncols() == 8 );
      interp( x, itw, x_field, gp_p, gp_lat, gp_lon ); 
    }
}
 
 
 
 
void interp_atmfield_by_gp( 
              VectorView        x, 
        const Index&            atmosphere_dim,
        ConstTensor3View        x_field,
        const ArrayOfGridPos&   gp_p,
        const ArrayOfGridPos&   gp_lat,
        const ArrayOfGridPos&   gp_lon )
{
  Matrix itw;
 
  interp_atmfield_gp2itw( itw, atmosphere_dim, gp_p, gp_lat, gp_lon );
 
  interp_atmfield_by_itw( x, atmosphere_dim, x_field, gp_p, gp_lat, gp_lon, 
                                                                         itw );
}
 
 
 
 
Numeric interp_atmfield_by_gp( 
        const Index&            atmosphere_dim,
        ConstTensor3View        x_field,
        const GridPos&          gp_p,
        const GridPos&          gp_lat,
        const GridPos&          gp_lon )
{
  ArrayOfGridPos agp_p(1), agp_lat(0), agp_lon(0);
 
  gridpos_copy( agp_p[0], gp_p  ); 
 
  if( atmosphere_dim > 1 )
    {
      agp_lat.resize(1);
      gridpos_copy( agp_lat[0], gp_lat  );
    }
 
  if( atmosphere_dim > 2 )
    {
      agp_lon.resize(1);
      gridpos_copy( agp_lon[0], gp_lon  );
    }
 
  Vector x(1);
 
  interp_atmfield_by_gp( x, atmosphere_dim, x_field, agp_p, agp_lat, agp_lon );
 
  return x[0];
}
 
 
 
 
void interp_cloudfield_gp2itw( 
              Matrix&           itw, 
              ArrayOfGridPos&   gp_p_out,
              ArrayOfGridPos&   gp_lat_out,
              ArrayOfGridPos&   gp_lon_out,
        const ArrayOfGridPos&   gp_p_in,
        const ArrayOfGridPos&   gp_lat_in,
        const ArrayOfGridPos&   gp_lon_in,
        const Index&            atmosphere_dim,
        const ArrayOfIndex&     cloudbox_limits )
{
  const Index n = gp_p_in.nelem();
 
  // Shift grid positions to cloud box grids
  if( atmosphere_dim == 1 )
    {
      gp_p_out.resize(n);
      for (Index i = 0; i < n; i++ ) 
        { 
          gridpos_copy( gp_p_out[i], gp_p_in[i] ); 
          gp_p_out[i].idx -= cloudbox_limits[0]; 
        }
      const Index n1 = cloudbox_limits[1] - cloudbox_limits[0];
      gridpos_upperend_check( gp_p_out, n1 );
      //
      itw.resize(n,2);
      interpweights( itw, gp_p_out );
    }
  else if( atmosphere_dim == 2 )
    {
      gp_p_out.resize(n);
      gp_lat_out.resize(n);
      for (Index i = 0; i < n; i++ ) 
        { 
          gridpos_copy( gp_p_out[i], gp_p_in[i] ); 
          gridpos_copy( gp_lat_out[i], gp_lat_in[i] ); 
          gp_p_out[i].idx   -= cloudbox_limits[0]; 
          gp_lat_out[i].idx -= cloudbox_limits[2]; 
        }
      const Index n1 = cloudbox_limits[1] - cloudbox_limits[0];
      const Index n2 = cloudbox_limits[3] - cloudbox_limits[2];
      gridpos_upperend_check( gp_p_out,   n1 );
      gridpos_upperend_check( gp_lat_out, n2 );
      //
      itw.resize(n,4);
      interpweights( itw, gp_p_out, gp_lat_out );
    }
  else 
    {
      gp_p_out.resize(n);
      gp_lat_out.resize(n);
      gp_lon_out.resize(n);
      for (Index i = 0; i < n; i++ ) 
        { 
          gridpos_copy( gp_p_out[i], gp_p_in[i] ); 
          gridpos_copy( gp_lat_out[i], gp_lat_in[i] ); 
          gridpos_copy( gp_lon_out[i], gp_lon_in[i] ); 
          gp_p_out[i].idx   -= cloudbox_limits[0]; 
          gp_lat_out[i].idx -= cloudbox_limits[2]; 
          gp_lon_out[i].idx -= cloudbox_limits[4];
        }
      const Index n1 = cloudbox_limits[1] - cloudbox_limits[0];
      const Index n2 = cloudbox_limits[3] - cloudbox_limits[2];
      const Index n3 = cloudbox_limits[5] - cloudbox_limits[4];
      gridpos_upperend_check( gp_p_out,   n1 );
      gridpos_upperend_check( gp_lat_out, n2 );
      gridpos_upperend_check( gp_lon_out, n3 );
      //
      itw.resize(n,8);
      interpweights( itw, gp_p_out, gp_lat_out, gp_lon_out );
    }
}
 
 
 
 
void interp_atmsurface_gp2itw( 
              Matrix&           itw, 
        const Index&            atmosphere_dim,
        const ArrayOfGridPos&   gp_lat,
        const ArrayOfGridPos&   gp_lon )
{
  if( atmosphere_dim == 1 )
    {
      itw.resize(1,1);
      itw = 1;
    }
 
  else if( atmosphere_dim == 2 )
    {
      const Index n = gp_lat.nelem();
      itw.resize(n,2);
      interpweights( itw, gp_lat );
    }
 
  else if( atmosphere_dim == 3 )
    {
      const Index n = gp_lat.nelem();
      assert( n == gp_lon.nelem() );
      itw.resize(n,4);
      interpweights( itw, gp_lat, gp_lon );
    }
}
 
 
 
 
void interp_atmsurface_by_itw( 
              VectorView        x, 
        const Index&            atmosphere_dim,
        ConstMatrixView         x_surface,
        const ArrayOfGridPos&   gp_lat,
        const ArrayOfGridPos&   gp_lon,
        ConstMatrixView         itw )
{
  if( atmosphere_dim == 1 )
    { 
      assert( itw.ncols() == 1 );
      x = x_surface(0,0); 
    }
 
  else if( atmosphere_dim == 2 )
    { 
      assert( x.nelem() == gp_lat.nelem() );
      assert( itw.ncols() == 2 );
      interp( x, itw, x_surface(Range(joker),0), gp_lat ); 
    }
 
  else if( atmosphere_dim == 3 )
    { 
      assert( x.nelem() == gp_lat.nelem() );
      assert( itw.ncols() == 4 );
      interp( x, itw, x_surface, gp_lat, gp_lon ); 
    }
}
 
 
 
 
void interp_atmsurface_by_gp( 
              VectorView        x, 
        const Index&            atmosphere_dim,
        ConstMatrixView         x_surface,
        const ArrayOfGridPos&   gp_lat,
        const ArrayOfGridPos&   gp_lon )
{
  Matrix itw;
 
  interp_atmsurface_gp2itw( itw, atmosphere_dim, gp_lat, gp_lon );
 
  interp_atmsurface_by_itw( x, atmosphere_dim, x_surface, gp_lat, gp_lon, itw );
}
 
 
 
 
Numeric interp_atmsurface_by_gp( 
        const Index&       atmosphere_dim,
        ConstMatrixView    x_surface,
        const GridPos&     gp_lat,
        const GridPos&     gp_lon )
{
  ArrayOfGridPos agp_lat(0), agp_lon(0);
 
  if( atmosphere_dim > 1 )
    {
      agp_lat.resize(1);
      gridpos_copy( agp_lat[0], gp_lat  );
    }
 
  if( atmosphere_dim > 2 )
    {
      agp_lon.resize(1);
      gridpos_copy( agp_lon[0], gp_lon  );
    }
 
  Vector x(1);
 
  interp_atmsurface_by_gp( x, atmosphere_dim, x_surface, agp_lat, agp_lon );
 
  return x[0];
}
 
 
 
 
void interp_gfield3( 
                 Numeric&   value,
           const GriddedField3&   gfield3,
           const Index&     effective_dim,
           const Numeric&   x,
           const Numeric&   y,
           const Numeric&   z,
           const String&    dim0,
           const String&    dim1,
           const String&    dim2 )
{
  chk_if_in_range( "effective_dim", effective_dim, 1, 3 );
 
  ArrayOfGridPos gp0(1), gp1(1), gp2(1);
 
  // Check and grid position, dimension 0
  //
    {
      const String gridname = gfield3.get_grid_name(0);
      const ConstVectorView grid = gfield3.get_numeric_grid(0);
      //
      if( dim0 != gridname )
        {
          ostringstream os;
          os << "Wrong quantity found for grid of first dimension:\n"
             << "  expected quantity  : " << dim0 << "\n"
             << "  quantity in gfield : " << gridname << "\n";
          throw runtime_error( os.str() );
          }
      //
      chk_if_increasing( "first grid of gfield3", grid );
      //
      if( x < grid[0]  ||  x > last(grid) )
        {
          ostringstream os;
          os << "Interpolation outside covered range for first dimension:\n"
             << "  interpolation point : " << x << "\n"
             << "  gfield grid range   : [" << grid[0] << "," << last(grid) 
             << "]\n";
          throw runtime_error( os.str() );
        }       
      //
      gridpos( gp0, grid, Vector(1,x) );
    }
 
  // Check and grid position, dimension 2
  //
  if( effective_dim >= 2 )
    {
      const String gridname = gfield3.get_grid_name(1);
      const ConstVectorView grid = gfield3.get_numeric_grid(1);
      //
      if( dim1 != gridname )
        {
          ostringstream os;
          os << "Wrong quantity found for grid of second dimension:\n"
             << "  expected quantity  : " << dim1 << "\n"
             << "  quantity in gfield : " << gridname << "\n";
          throw runtime_error( os.str() );
          }
      //
      chk_if_increasing( "second grid of gfield3", grid );
      //
      if( y < grid[0]  ||  y > last(grid) )
        {
          ostringstream os;
          os << "Interpolation outside covered range for second dimension:\n"
             << "  interpolation point : " << y << "\n"
             << "  gfield grid range   : [" << grid[0] << "," << last(grid) 
             << "]\n";
          throw runtime_error( os.str() );
        }       
      //
      gridpos( gp1, grid, Vector(1,y) );
    }
 
  // Check and grid position, dimension 2
  //
  if( effective_dim >= 3 )
    {
      const String gridname = gfield3.get_grid_name(2);
      const ConstVectorView grid = gfield3.get_numeric_grid(2);
      //
      if( dim2 != gridname )
        {
          ostringstream os;
          os << "Wrong quantity found for grid of third dimension:\n"
             << "  expected quantity  : " << dim2 << "\n"
             << "  quantity in gfield : " << gridname << "\n";
          throw runtime_error( os.str() );
          }
      //
      chk_if_increasing( "third grid of gfield3", grid );
      //
      if( z < grid[0]  ||  z > last(grid) )
        {
          ostringstream os;
          os << "Interpolation outside covered range for second dimension:\n"
             << "  interpolation point : " << z << "\n"
             << "  gfield grid range   : [" << grid[0] << "," << last(grid) 
             << "]\n";
          throw runtime_error( os.str() );
        }       
      //
      gridpos( gp2, grid, Vector(1,z) );
    }
 
  // Perform interpolation
  //
  Vector result( 1 );
  //
  if( effective_dim == 1 )
    {
      if( gfield3.data.nrows() > 1  ||  gfield3.data.ncols() > 1 )
        {
          ostringstream os;
          os << "A 1D interpolation requested, but the provided gridded field "
             << "has an effective dimension of 2D or 3D.";
          throw runtime_error( os.str() );
        }       
        
      Matrix itw(1,2);
      interpweights( itw, gp0 );
      interp( result, itw, gfield3.data(joker,0,0), gp0 );
    }
  //
  else if( effective_dim == 2 )
    {
      if( gfield3.data.ncols() > 1 )
        {
          ostringstream os;
          os << "A 2D interpolation requested, but the provided gridded field "
             << "has an effective dimension of 3D.";
          throw runtime_error( os.str() );
        }       
 
      Matrix itw(1,4);
      interpweights( itw, gp0, gp1 );
      interp( result, itw, gfield3.data(joker,joker,0), gp0, gp1 );
    }
  //
  else if( effective_dim == 3 )
    {
      Matrix itw(1,8);
      interpweights( itw, gp0, gp1, gp2 );
      interp( result, itw, gfield3.data, gp0, gp1, gp2 );
    }
 
  value = result[0];
}
 
 
 
 
 
 
/*===========================================================================
  === Conversion altitudes / pressure
  ===========================================================================*/
 
 
void itw2p(
              VectorView       p_values,
        ConstVectorView        p_grid,
        const ArrayOfGridPos&  gp,
        ConstMatrixView        itw )
{
  assert( itw.ncols() == 2 );
  assert( p_values.nelem() == itw.nrows() );
 
  // Local variable to store log of the pressure grid:
  Vector logpgrid( p_grid.nelem() );
 
  transform( logpgrid, log, p_grid );
  
  interp( p_values, itw, logpgrid, gp ); 
 
  transform( p_values, exp, p_values );
}
 
 
 
 
void p2gridpos( ArrayOfGridPos& gp,
                ConstVectorView old_pgrid,
                ConstVectorView new_pgrid,
                const Numeric&  extpolfac )
{
  // Local variable to store log of the pressure grids
  Vector logold( old_pgrid.nelem() );
  Vector lognew( new_pgrid.nelem() );
 
  transform( logold, log, old_pgrid );
  transform( lognew, log, new_pgrid );
 
  gridpos( gp, logold, lognew, extpolfac );
}
 
 
void p2gridpos_poly( ArrayOfGridPosPoly& gp,
                     ConstVectorView old_pgrid,
                     ConstVectorView new_pgrid,
                     const Index     order,
                     const Numeric&  extpolfac )
{
  // Local variable to store log of the pressure grids
  Vector logold( old_pgrid.nelem() );
  Vector lognew( new_pgrid.nelem() );
  
  transform( logold, log, old_pgrid );
  transform( lognew, log, new_pgrid );
  
  gridpos_poly( gp, logold, lognew, order, extpolfac );
}
 
 
/*===========================================================================
  === Various functions connected to the atmospheric fields
  ===========================================================================*/
 
 
void z_at_lat_2d(
             VectorView   z,
        ConstVectorView   p_grid,
// FIXME only used in assertion
#ifndef NDEBUG
        ConstVectorView   lat_grid,
#else
        ConstVectorView,
#endif
        ConstMatrixView   z_field,
        const GridPos&    gp_lat )
{
  const Index   np = p_grid.nelem();
 
  assert( z.nelem() == np );
  assert( z_field.nrows() == np );
  assert( z_field.ncols() == lat_grid.nelem() );
 
  Matrix           z_matrix(np,1);
  ArrayOfGridPos   gp_z(np), agp_lat(1);
  Tensor3          itw(np,1,4);
 
  gridpos_copy( agp_lat[0], gp_lat );
  gridpos( gp_z, p_grid, p_grid );
  interpweights( itw, gp_z, agp_lat );
  interp( z_matrix, itw, z_field, gp_z, agp_lat );
 
  z = z_matrix(Range(joker),0);
}
 
 
 
 
void z_at_latlon(
             VectorView    z,
        ConstVectorView    p_grid,
//FIXME only used in assertion
#ifndef NDEBUG
        ConstVectorView    lat_grid,
        ConstVectorView    lon_grid,
#else
        ConstVectorView,
        ConstVectorView,
#endif
        ConstTensor3View   z_field,
        const GridPos&     gp_lat,
        const GridPos&     gp_lon )
{
  const Index   np = p_grid.nelem();
 
  assert( z.nelem() == np );
  assert( z_field.npages() == np );
  assert( z_field.nrows() == lat_grid.nelem() );
  assert( z_field.ncols() == lon_grid.nelem() );
 
  Tensor3          z_tensor(np,1,1);
  ArrayOfGridPos   agp_z(np), agp_lat(1), agp_lon(1);
  Tensor4          itw(np,1,1,8);
 
  gridpos_copy( agp_lat[0], gp_lat );
  gridpos_copy( agp_lon[0], gp_lon );
  gridpos( agp_z, p_grid, p_grid );
  interpweights( itw, agp_z, agp_lat, agp_lon );
 
  interp( z_tensor, itw, z_field, agp_z, agp_lat, agp_lon );
 
  z = z_tensor(Range(joker),0,0);
}