special_interp.cc

Go to the documentation of this file.

/* Copyright (C) 2002-2012 Stefan Buehler <sbuehler@ltu.se>

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */

#include <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( 
              VectorView      itw, 
              GridPos&        gp_p_out,
              GridPos&        gp_lat_out,
              GridPos&        gp_lon_out,
        const GridPos&        gp_p_in,
        const GridPos&        gp_lat_in,
        const GridPos&        gp_lon_in,
        const Index&          atmosphere_dim,
        const ArrayOfIndex&   cloudbox_limits )
{
  // Shift grid positions to cloud box grids
  if( atmosphere_dim == 1 )
    {
      gridpos_copy( gp_p_out, gp_p_in ); 
      gp_p_out.idx -= cloudbox_limits[0]; 
      gridpos_upperend_check( gp_p_out, cloudbox_limits[1]-cloudbox_limits[0] );
      assert(itw.nelem() == 2 );
      interpweights( itw, gp_p_out );
    }
  else if( atmosphere_dim == 2 )
    {
      gridpos_copy( gp_p_out, gp_p_in ); 
      gridpos_copy( gp_lat_out, gp_lat_in ); 
      gp_p_out.idx   -= cloudbox_limits[0]; 
      gp_lat_out.idx -= cloudbox_limits[2]; 
      gridpos_upperend_check( gp_p_out,  cloudbox_limits[1]-cloudbox_limits[0]);
      gridpos_upperend_check( gp_lat_out,cloudbox_limits[3]-cloudbox_limits[2]);
      assert(itw.nelem() == 4 );
      interpweights( itw, gp_p_out, gp_lat_out );
    }
  else 
    {
      gridpos_copy( gp_p_out, gp_p_in ); 
      gridpos_copy( gp_lat_out, gp_lat_in ); 
      gridpos_copy( gp_lon_out, gp_lon_in ); 
      gp_p_out.idx   -= cloudbox_limits[0]; 
      gp_lat_out.idx -= cloudbox_limits[2]; 
      gp_lon_out.idx -= cloudbox_limits[4];
      gridpos_upperend_check( gp_p_out,  cloudbox_limits[1]-cloudbox_limits[0]);
      gridpos_upperend_check( gp_lat_out,cloudbox_limits[3]-cloudbox_limits[2]);
      gridpos_upperend_check( gp_lon_out,cloudbox_limits[5]-cloudbox_limits[4]);
      assert(itw.nelem() == 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];
}





/*===========================================================================
  === 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 );
}




void rte_pos2gridpos(
         GridPos&     gp_p,
         GridPos&     gp_lat,
         GridPos&     gp_lon,
   const Index&       atmosphere_dim,
   ConstVectorView    p_grid,
   ConstVectorView    lat_grid,
   ConstVectorView    lon_grid,
   ConstTensor3View   z_field,
   ConstVectorView    rte_pos )
{
  chk_rte_pos( atmosphere_dim, rte_pos );

  if( atmosphere_dim == 1 )
    { 
      chk_interpolation_grids( "Altitude interpolation", z_field(joker,0,0), 
                                                                  rte_pos[0] );
      gridpos( gp_p, z_field(joker,0,0), rte_pos[0] ); 
    }
  else
    {
      // Determine z at lat/lon (z_grid) by blue interpolation
      const Index np = p_grid.nelem();
      Vector z_grid( np );
      ArrayOfGridPos agp_z, agp_lat(np);
      //
      gridpos_1to1( agp_z, p_grid );
      //
      chk_interpolation_grids( "Latitude interpolation", lat_grid, rte_pos[1] );
      gridpos( gp_lat, lat_grid, rte_pos[1] );

      if( atmosphere_dim == 2 )
        {
          for( Index i=0; i<np; i++ )
            { agp_lat[i] = gp_lat; }
          Matrix itw( np, 4 );
          interpweights( itw, agp_z, agp_lat );
          interp( z_grid, itw, z_field(joker,joker,0), agp_z, agp_lat );
        }
      else
        {
          chk_interpolation_grids( "Longitude interpolation", lon_grid, 
                                                                  rte_pos[2] );
          gridpos( gp_lon, lon_grid, rte_pos[2] );
          ArrayOfGridPos agp_lon(np);
          for( Index i=0; i<np; i++ )
            { 
              agp_lat[i] = gp_lat;  
              agp_lon[i] = gp_lon; 
            }
          Matrix itw( np, 8 );
          interpweights( itw, agp_z, agp_lat, agp_lon );
          interp( z_grid, itw, z_field, agp_z, agp_lat, agp_lon );

        }

      // And use z_grid to get gp_p (gp_al and gp_lon determined above)
      chk_interpolation_grids( "Altitude interpolation", z_grid, rte_pos[0] );
      gridpos( gp_p, z_grid, rte_pos[0] );
    }
}




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);
}



/*===========================================================================
  === Interpolation of GriddedFields
  ===========================================================================*/


void complex_n_interp(
         MatrixView       n_real,
         MatrixView       n_imag,
   const GriddedField3&   complex_n,
   const String&          varname,
   ConstVectorView        f_grid,
   ConstVectorView        t_grid )
{
  // Set expected order of grids
  Index gfield_fID    = 0;
  Index gfield_tID    = 1;
  Index gfield_compID = 2;

  // Check of complex_n
  //
  complex_n.checksize_strict();
  //
  chk_griddedfield_gridname( complex_n, gfield_fID,    "Frequency"   );
  chk_griddedfield_gridname( complex_n, gfield_tID,    "Temperature" );
  chk_griddedfield_gridname( complex_n, gfield_compID, "Complex"     );
  //
  if( complex_n.data.ncols() != 2 )
    {
      ostringstream os;
      os << "The data in *" << varname << "* must have exactly two pages. One page "
         << "each\nfor the real and imaginary part of the complex refractive index.";
    } 

  // Frequency and temperature grid sizes
  const Index nf_in = complex_n.data.npages();
  const Index nt_in = complex_n.data.nrows(); 
  const Index nf_out = f_grid.nelem();
  const Index nt_out = t_grid.nelem();

  //Assert size of output variables
  assert( n_real.nrows() == nf_out  &&  n_real.ncols() == nt_out );
  assert( n_imag.nrows() == nf_out  &&  n_imag.ncols() == nt_out );

  const Vector& f_grid_in = complex_n.get_numeric_grid(gfield_fID);
  const Vector& t_grid_in = complex_n.get_numeric_grid(gfield_tID);

  // Expand/interpolate in frequency dimension
  //
  Matrix nrf(nf_out,nt_in), nif(nf_out,nt_in);
  //
  if( nf_in == 1 )
    {
      for( Index i=0; i<nf_out; i++ )
        { 
          nrf(i,joker) = complex_n.data(0,joker,0); 
          nif(i,joker) = complex_n.data(0,joker,1); 
        }
    }
  else
    {
      chk_interpolation_grids( "Frequency interpolation", f_grid_in, f_grid );
      //
      ArrayOfGridPos gp( nf_out );
      Matrix         itw( nf_out, 2 );
      gridpos( gp, f_grid_in, f_grid );
      interpweights( itw, gp );
      for( Index i=0; i<nt_in; i++ )
        { 
          interp( nrf(joker,i), itw, complex_n.data(joker,i,0), gp );
          interp( nif(joker,i), itw, complex_n.data(joker,i,1), gp );
        }
    }

  // Expand/interpolate in temperature dimension
  //
  if( nt_in == 1 )
    {
      for( Index i=0; i<nt_out; i++ )
        { 
          n_real(joker,i) = nrf(joker,0); 
          n_imag(joker,i) = nif(joker,0); 
        }
    }
  else
    {
      chk_interpolation_grids( "Temperature interpolation", t_grid_in, t_grid );
      //
      ArrayOfGridPos gp( nt_out );
      Matrix         itw( nt_out, 2 );
      gridpos( gp, t_grid_in, t_grid );
      interpweights( itw, gp );
      for( Index i=0; i<nf_out; i++ )
        { 
          interp( n_real(i,joker), itw, nrf(i,joker), gp );
          interp( n_imag(i,joker), itw, nif(i,joker), gp );
        }
    }
}