check_input.cc

Go to the documentation of this file.

/* Copyright (C) 2002-2012
   Patrick Eriksson <patrick.eriksson@chalmers.se>
   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. */



/*===========================================================================
  === File description 
  ===========================================================================*/

/*===========================================================================
  === External declarations
  ===========================================================================*/

#include <cmath>
#include <stdexcept>
#include "auto_md.h"
#include "check_input.h"
#include "array.h"
#include "logic.h"
#include "gridded_fields.h"
#include <cfloat>

extern const Index GFIELD3_P_GRID;
extern const Index GFIELD3_LAT_GRID;
extern const Index GFIELD3_LON_GRID;





/*===========================================================================
  === Functions for Index
  ===========================================================================*/


void chk_if_bool( 
        const String&   x_name,
        const Index&    x )
{
  if ( !is_bool(x) )
    {
      ostringstream os;
      os << "The variable *" << x_name <<  "* must be a boolean (0 or 1).\n" 
         << "The present value of *"<< x_name <<  "* is " << x << ".";
      throw runtime_error( os.str() );
    }
}


void chk_if_in_range( 
        const String&   x_name,
        const Index&    x, 
        const Index&    x_low, 
        const Index&    x_high )
{
  if ( (x<x_low) || (x>x_high) )
    {
      ostringstream os;
      os << "The variable *" << x_name <<  "* must fulfill:\n"
         << "   " << x_low << " <= " << x_name << " <= " << x_high << "\n" 
         << "The present value of *"<< x_name <<  "* is " << x << ".";
      throw runtime_error( os.str() );
    }
}


void chk_if_increasing( 
        const String&       x_name,
        const ArrayOfIndex& x ) 
{
  if ( !is_increasing(x) )
    {
      ostringstream os;
      os << "The ArrayOfIndex *" << x_name <<  "* must have strictly\n"
         << "increasing values, but this is not the case.\n";
      os << "x = " << x << "\n";
      throw runtime_error( os.str() );
    }
}





/*===========================================================================
  === Functions for Numeric
  ===========================================================================*/


void chk_not_negative( 
        const String&    x_name,
        const Numeric&   x ) 
{
  if ( x < 0 )
    {
      ostringstream os;
      os << "The variable *" << x_name <<  "* must be >= 0.\n"
         << "The present value of *"<< x_name <<  "* is " << x << ".";
      throw runtime_error( os.str() );
    }
}




void chk_if_in_range( 
        const String&    x_name,
        const Numeric&   x, 
        const Numeric&   x_low, 
        const Numeric&   x_high )
{
  if ( (x<x_low) || (x>x_high) )
    {
      ostringstream os;
      os << "The variable *" << x_name <<  "* must fulfill:\n"
         << "   " << x_low << " <= " << x_name << " <= " << x_high << "\n" 
         << "The present value of *"<< x_name <<  "* is " << x << ".";
      throw runtime_error( os.str() );
    }
}





/*===========================================================================
  === Functions for Vector
  ===========================================================================*/


void chk_vector_length( 
        const String&      x_name,
        ConstVectorView    x,
        const Index&       l ) 
{
  if ( x.nelem() != l )
    {
      ostringstream os;
      os << "The vector *" << x_name <<  "* must have the length " << l 
         << ".\n" 
         << "The present length of *"<< x_name <<  "* is " << x.nelem() << ".";
      throw runtime_error( os.str() );
    }
}




void chk_vector_length( 
        const String&      x1_name,
        const String&      x2_name,
        ConstVectorView    x1, 
        ConstVectorView    x2 ) 
{
  if ( x1.nelem() != x2.nelem() )
    {
      ostringstream os;
      os << "The vectors *" << x1_name <<  "* and *" << x2_name 
         <<  "* must have the same length.\n"
         << "The length of *"<< x1_name <<  "* is " << x1.nelem() << ".\n"
         << "The length of *"<< x2_name <<  "* is " << x2.nelem() << ".";
      throw runtime_error( os.str() );
    }
}




void chk_if_increasing( 
        const String&      x_name,
        ConstVectorView    x ) 
{
  if ( !is_increasing(x) )
    {
      ostringstream os;
      os << "The vector *" << x_name <<  "* must have strictly\n"
         << "increasing values, but this is not the case.\n";
      os << "x = " << x << "\n";
      throw runtime_error( os.str() );
    }
}




void chk_if_decreasing( 
        const String&      x_name,
        ConstVectorView    x ) 
{
  if ( !is_decreasing(x) )
    {
      ostringstream os;
      os << "The vector *" << x_name <<  "* must have strictly\ndecreasing "
         << "values, but this is not the case.\n";
      throw runtime_error( os.str() );
    }
}




void chk_if_equal(
        const String&   x1_name,
        const String&   x2_name,
        ConstVectorView v1,
        ConstVectorView v2,
        Numeric         margin
        )
{
  chk_vector_length(x1_name, x2_name, v1, v2);

  for (Index i = 0; i<v1.nelem(); i++)
  {
    if (abs(v1[i] - v2[i]) > margin)
      {
        ostringstream os;
        os << "Vectors " << x1_name << " and " << x2_name 
           << " differ.\n"
           << x1_name << "[" << i << "]" << " = " << v1[i] << "\n"
           << x2_name << "[" << i << "]" << " = " << v2[i] << "\n"
           << "Difference should not exceed " << margin << "\n";
       throw runtime_error(os.str());
     }
  }
}





/*===========================================================================
  === Functions for Matrix
  ===========================================================================*/


void chk_matrix_ncols( 
        const String&      x_name,
        ConstMatrixView    x,
        const Index&       l ) 
{
  if ( x.ncols() != l )
    {
      ostringstream os;
      os << "The matrix *" << x_name <<  "* must have " << l << " columns,\n"
         << "but the number of columns is " << x.ncols() << ".";
      throw runtime_error( os.str() );
    }
}




void chk_matrix_nrows( 
        const String&      x_name,
        ConstMatrixView    x,
        const Index&       l ) 
{
  if ( x.nrows() != l )
    {
      ostringstream os;
      os << "The matrix *" << x_name <<  "* must have " << l << " rows,\n"
         << "but the number of rows is " << x.nrows() << ".";
      throw runtime_error( os.str() );
    }
}





/*===========================================================================
  === Functions for Tensors
  ===========================================================================*/


void chk_size( const String&    x_name,
               ConstVectorView  x,
               const Index&     c ) 
{
  if ( !is_size(x,c) )
    {
      ostringstream os;
      os << "The object *" << x_name
         << "* does not have the right size.\n"
         << "Dimension should be:"
         << " " << c 
         << ",\nbut it is:          "
         << " " << x.nelem()      
         << ".";
      throw runtime_error( os.str() );
    }
}


void chk_size( const String&    x_name,
               ConstMatrixView  x,
               const Index&     r,
               const Index&     c ) 
{
  if ( !is_size(x,r,c) )
    {
      ostringstream os;
      os << "The object *" << x_name
         << "* does not have the right size.\n"
         << "Dimensions should be:"
         << " " << r 
         << " " << c 
         << ",\nbut they are:         "
         << " " << x.nrows()      
         << " " << x.ncols()      
         << ".";
      throw runtime_error( os.str() );
    }
}


void chk_size( const String&    x_name,
               ConstTensor3View x,
               const Index&     p,
               const Index&     r,
               const Index&     c ) 
{
  if ( !is_size(x,p,r,c) )
    {
      ostringstream os;
      os << "The object *" << x_name
         << "* does not have the right size.\n"
         << "Dimensions should be:"
         << " " << p 
         << " " << r 
         << " " << c 
         << ",\nbut they are:         "
         << " " << x.npages()     
         << " " << x.nrows()      
         << " " << x.ncols()      
         << ".";
      throw runtime_error( os.str() );
    }
}


void chk_size( const String&    x_name,
               ConstTensor4View x,
               const Index&     b,
               const Index&     p,
               const Index&     r,
               const Index&     c ) 
{
  if ( !is_size(x,b,p,r,c) )
    {
      ostringstream os;
      os << "The object *" << x_name
         << "* does not have the right size.\n"
         << "Dimensions should be:"
         << " " << b 
         << " " << p 
         << " " << r 
         << " " << c 
         << ",\nbut they are:         "
         << " " << x.nbooks()     
         << " " << x.npages()     
         << " " << x.nrows()      
         << " " << x.ncols()      
         << ".";
      throw runtime_error( os.str() );
    }
}


void chk_size( const String&    x_name,
               ConstTensor5View x,
               const Index&     s,
               const Index&     b,
               const Index&     p,
               const Index&     r,
               const Index&     c ) 
{
  if ( !is_size(x,s,b,p,r,c) )
    {
      ostringstream os;
      os << "The object *" << x_name
         << "* does not have the right size.\n"
         << "Dimensions should be:"
         << " " << s 
         << " " << b 
         << " " << p 
         << " " << r 
         << " " << c 
         << ",\nbut they are:         "
         << " " << x.nshelves()   
         << " " << x.nbooks()     
         << " " << x.npages()     
         << " " << x.nrows()      
         << " " << x.ncols()      
         << ".";
      throw runtime_error( os.str() );
    }
}


void chk_size( const String&    x_name,
               ConstTensor6View x,
               const Index&     v,
               const Index&     s,
               const Index&     b,
               const Index&     p,
               const Index&     r,
               const Index&     c ) 
{
  if ( !is_size(x,v,s,b,p,r,c) )
    {
      ostringstream os;
      os << "The object *" << x_name
         << "* does not have the right size.\n"
         << "Dimensions should be:"
         << " " << v 
         << " " << s 
         << " " << b 
         << " " << p 
         << " " << r 
         << " " << c 
         << ",\nbut they are:         "
         << " " << x.nvitrines()  
         << " " << x.nshelves()   
         << " " << x.nbooks()     
         << " " << x.npages()     
         << " " << x.nrows()      
         << " " << x.ncols()      
         << ".";
      throw runtime_error( os.str() );
    }
}


void chk_size( const String&    x_name,
               ConstTensor7View x,
               const Index&     l,
               const Index&     v,
               const Index&     s,
               const Index&     b,
               const Index&     p,
               const Index&     r,
               const Index&     c ) 
{
  if ( !is_size(x,l,v,s,b,p,r,c) )
    {
      ostringstream os;
      os << "The object *" << x_name
         << "* does not have the right size.\n"
         << "Dimensions should be:"
         << " " << l 
         << " " << v 
         << " " << s 
         << " " << b 
         << " " << p 
         << " " << r 
         << " " << c 
         << ",\nbut they are:         "
         << " " << x.nlibraries() 
         << " " << x.nvitrines()  
         << " " << x.nshelves()   
         << " " << x.nbooks()     
         << " " << x.npages()     
         << " " << x.nrows()      
         << " " << x.ncols()      
         << ".";
      throw runtime_error( os.str() );
    }
}






/*===========================================================================
  === Functions for Agendas
  ===========================================================================*/


void chk_not_empty( 
        const String&      x_name,
        const Agenda&      x ) 
{
  if( x.nelem() == 0 )
    {
      ostringstream os;
      os << "The agenda *" << x_name << "* is empty.\nIt is not allowed \n"
         << "that an empty agenda that is actually used.\n"
         << "Empty agendas are only created of methods setting dummy values \n"
         << "to variables.";
      throw runtime_error( os.str() );
    }
}







/*===========================================================================
  === Functions for interpolation grids
  ===========================================================================*/


void chk_interpolation_grids_loose(Index&          ing_min,
                                   Index&          ing_max,
                                   const String&   which_interpolation,
                                   ConstVectorView old_grid,
                                   ConstVectorView new_grid,
                                   ConstVectorView data,
                                   const Index     order)
{
    chk_interpolation_grids_loose_no_data_check(ing_min, ing_max, which_interpolation,
                                                old_grid, new_grid, order);

    chk_interpolation_grids_loose_check_data(ing_min, ing_max, which_interpolation,
                                             old_grid, new_grid, data);
}



void chk_interpolation_grids_loose_no_data_check(Index&          ing_min,
                                                 Index&          ing_max,
                                                 const String&   which_interpolation,
                                                 ConstVectorView old_grid,
                                                 ConstVectorView new_grid,
                                                 const Index     order)
{
  const Index n_old = old_grid.nelem();
  
  if (!new_grid.nelem()) throw runtime_error("The new grid is not allowed to be empty.");

  ostringstream os;
  os << "There is a problem with the grids for the\n"
  << "following interpolation: " << which_interpolation << ".\n";
  
  // Old grid must have at least order+1 elements:
  if (n_old < order+1)
  {
    os << "The original grid must have at least " << order+1 << " elements.";
    throw runtime_error( os.str() );
  }
  
  // Decide whether we have an ascending or descending grid:
  const bool ascending = ( old_grid[0] <= old_grid[1] );
  
  // Minimum and maximum allowed value from old grid. (Will include
  // extrapolation tolerance.)
  Numeric og_min, og_max;
  
  ing_min = 0;
  ing_max = new_grid.nelem()-1;
  if (ascending)  
  {
    // Old grid must be strictly increasing (no duplicate values.)
    if ( !is_increasing(old_grid) )
    {
      os << "The original grid must be strictly sorted\n"
      << "(no duplicate values). Yours is:\n"
      << old_grid << ".";
      throw runtime_error( os.str() );
    }
    
    // Limits of extrapolation. 
    og_min = old_grid[0];
    og_max = old_grid[n_old-1];
  }
  else
  {
    // Old grid must be strictly decreasing (no duplicate values.)
    if ( !is_decreasing(old_grid) )
    {
      os << "The original grid must be strictly sorted\n"
      << "(no duplicate values). Yours is:\n"
      << old_grid << ".";
      throw runtime_error( os.str() );
    }
    
    // The max is now the first point, the min the last point!
    og_max = old_grid[0];
    og_min = old_grid[n_old-1];
  }
  
  // Min and max of new grid:
  const Numeric ng_min = min(new_grid);
  const Numeric ng_max = max(new_grid);
  
  // If new grid is not inside old grid, determine the indexes of the range
  // that is.
  
  const Index iog_min = 0;
  const Index iog_max = old_grid.nelem()-1;

  ing_min = 0;
  ing_max = new_grid.nelem()-1;

  if (ascending)
  {
      if (ng_max > og_max)
      {
          while (ing_max > 0 && new_grid[ing_max] > old_grid[iog_max])
              ing_max--;
      }

      if (ng_min < og_min)
      {
          while (ing_min < new_grid.nelem()-1 && new_grid[ing_min] < old_grid[iog_min])
              ing_min++;
      }
  }
  else
  {
      if (ng_min < og_min)
      {
          while (ing_max > 0 && new_grid[ing_max] < old_grid[iog_max])
              ing_max--;
      }

      if (ng_max > og_max)
      {
          while (ing_min < new_grid.nelem()-1 && new_grid[ing_min] > old_grid[iog_min])
              ing_min++;
      }
  }
}



void chk_interpolation_pgrids_loose_no_data_check(Index&          ing_min,
                                                  Index&          ing_max,
                                                  const String&   which_interpolation,
                                                  ConstVectorView old_pgrid,
                                                  ConstVectorView new_pgrid,
                                                  const Index     order)
{
    // 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 );

    chk_interpolation_grids_loose_no_data_check(ing_min, ing_max,
                                                which_interpolation,
                                                logold, lognew,
                                                order);
}



void chk_interpolation_grids_loose_check_data(Index&          ing_min,
                                              Index&          ing_max,
                                              const String&   which_interpolation,
                                              ConstVectorView old_grid,
                                              ConstVectorView new_grid,
                                              ConstVectorView data)
{
  if (!new_grid.nelem()) throw runtime_error("The new grid is not allowed to be empty.");

  ostringstream os;
  os << "There is a problem with the grids for the\n"
     << "following interpolation: " << which_interpolation << ".\n";

  // Decide whether we have an ascending or descending grid:
  const bool ascending = ( old_grid[0] <= old_grid[1] );

  // If new grid is not inside old grid, determine the indexes of the range
  // that is.

  const Index iog_min = ascending?old_grid.nelem()-1:0;
  const Index iog_max = ascending?0:old_grid.nelem()-1;

  if (ing_min > 0 && data[iog_min] != 0)
  {
    os << "\nThe new grid is not fully inside the original grid.\n"
    << "This is allowed if the corresponding boundary value of raw data is 0.\n"
    << "New grid point: " << new_grid[ing_min] << "\n"
    << "Old grid point: " << old_grid[iog_min] << "\n"
    << "Boundary value: " << data[iog_min];
    throw runtime_error(os.str());
  }
  
  if (ing_max < new_grid.nelem()-1 && data[iog_max] != 0)
  {
    os << "\nThe the new grid is not fully inside the original grid.\n"
    << "This is allowed if the corresponding boundary value of raw data is 0.\n"
    << "New grid point: " << new_grid[ing_max] << "\n"
    << "Old grid point: " << old_grid[iog_max] << "\n"
    << "Boundary value: " << data[iog_max];
    throw runtime_error(os.str());
  }
}



void chk_interpolation_grids(const String&   which_interpolation,
                             ConstVectorView old_grid,
                             ConstVectorView new_grid,
                             const Index     order,
                             const Numeric&  extpolfac,
                             const bool      islog)
{
  const Index n_old = old_grid.nelem();

  if (!new_grid.nelem()) throw runtime_error("The new grid is not allowed to be empty.");


  // Old grid must have at least order+1 elements:
  if (n_old < order+1)
    {
      ostringstream os;
      os << "There is a problem with the grids for the\n"
         << "following interpolation: " << which_interpolation << ".\n"
         << "The original grid must have at least " << order+1 << " elements.";
      throw runtime_error( os.str() );
    }
  
  // Decide whether we have an ascending or descending grid:
  const bool ascending = ( old_grid[0] <= old_grid[1] );

  // Minimum and maximum allowed value from old grid. (Will include
  // extrapolation tolerance.)
  Numeric og_min, og_max;

  if (ascending)  
    {
      // Old grid must be strictly increasing (no duplicate values.)
      if ( !is_increasing(old_grid) )
        {
          ostringstream os;
          os << "There is a problem with the grids for the\n"
             << "following interpolation: " << which_interpolation << ".\n"
             << "The original grid must be strictly sorted\n"
             << "(no duplicate values). Yours is:\n"
             << old_grid << ".";
          throw runtime_error( os.str() );
        }

      // Limits of extrapolation. 
      og_min = old_grid[0] - 
               extpolfac * ( old_grid[1] - old_grid[0] );
      og_max = old_grid[n_old-1] + 
               extpolfac * ( old_grid[n_old-1] - old_grid[n_old-2] );
    }
  else
    {
      // Old grid must be strictly decreasing (no duplicate values.)
      if ( !is_decreasing(old_grid) )
        {
          ostringstream os;
          os << "There is a problem with the grids for the\n"
             << "following interpolation: " << which_interpolation << ".\n"
             << "The original grid must be strictly sorted\n"
             << "(no duplicate values). Yours is:\n"
             << old_grid << ".";
          throw runtime_error( os.str() );
        }

      // The max is now the first point, the min the last point!
      // I think the sign is right here...
      og_max = old_grid[0] - 
               extpolfac * ( old_grid[1] - old_grid[0] );
      og_min = old_grid[n_old-1] + 
               extpolfac * ( old_grid[n_old-1] - old_grid[n_old-2] );
    }
  
  // Min and max of new grid:
  const Numeric ng_min = min(new_grid);
  const Numeric ng_max = max(new_grid);

  // New grid must be inside old grid (plus extpolfac).
  // (Values right on the edge (ng_min==og_min) are still allowed.)

  if (ng_min < og_min)
    {
      ostringstream os;
      os << "There is a problem with the grids for the\n"
         << "following interpolation: " << which_interpolation << ".\n"
         << "The minimum of the new grid must be inside\n"
         << "the original grid. (We allow a bit of extrapolation,\n"
         << "but not so much).\n"
         << "Minimum of original grid:           " << min(old_grid);
      if (islog) os << " (" << exp(min(old_grid)) << ")";
      os << "\nMinimum allowed value for new grid: " << og_min;
      if (islog) os << " (" << exp(og_min) << ")";
      os << "\nActual minimum of new grid:         " << ng_min;
      if (islog) os << " (" << exp(ng_min) << ")";
      throw runtime_error( os.str() );
    }

  if (ng_max > og_max)
    {
      ostringstream os;
      os << "There is a problem with the grids for the\n"
         << "following interpolation: " << which_interpolation << ".\n"
         << "The maximum of the new grid must be inside\n"
         << "the original grid. (We allow a bit of extrapolation,\n"
         << "but not so much).\n"
         << "Maximum of original grid:           " << max(old_grid);
      if (islog) os << " (" << exp(max(old_grid)) << ")";
      os << "\nMaximum allowed value for new grid: " << og_max;
      if (islog) os << " (" << exp(og_max) << ")";
      os << "\nActual maximum of new grid:         " << ng_max;
      if (islog) os << " (" << exp(ng_max) << ")";
      throw runtime_error( os.str() );
    }

  // If we get here, than everything should be fine.
}



void chk_interpolation_grids(const String&   which_interpolation,
                             ConstVectorView old_grid,
                             const Numeric&  new_grid,
                             const Index     order,
                             const Numeric&  extpolfac )
{
  const Vector v(1, new_grid);
  chk_interpolation_grids(which_interpolation,
                          old_grid,
                          v,
                          order,
                          extpolfac );
}



void chk_interpolation_pgrids(const String&   which_interpolation,
                              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 );

    chk_interpolation_grids(which_interpolation, logold, lognew, order, extpolfac, true);
}





/*===========================================================================
  === Functions related to atmospheric and surface grids and fields.
  ===========================================================================*/


void chk_atm_grids( 
        const Index&      dim,
        ConstVectorView   p_grid,
        ConstVectorView   lat_grid,
        ConstVectorView   lon_grid )
{
  // p_grid
  if( p_grid.nelem() < 2 )
    throw runtime_error( "The length of *p_grid* must be >= 2." );
  chk_if_decreasing( "p_grid", p_grid );

  // lat_grid
  if( dim == 1 )
    {
      if( lat_grid.nelem() > 0 )
        throw runtime_error("For dim=1, the length of *lat_grid* must be 0.");
    }
  else
    {
      if( lat_grid.nelem() < 2 )
        throw runtime_error(
                         "For dim>1, the length of *lat_grid* must be >= 2.");
      chk_if_increasing( "lat_grid", lat_grid );
    }

  // lon_grid
  if( dim < 3 )
    { 
      if( lon_grid.nelem() > 0 )
        throw runtime_error("For dim<3, the length of *lon_grid* must be 0.");
    }
  else
    {
      if( lon_grid.nelem() < 2 )
        throw runtime_error(
                         "For dim=3, the length of *lon_grid* must be >= 2.");
      chk_if_increasing( "lon_grid", lon_grid );
    }

  // Check that latitude and longitude grids are inside OK ranges for 3D
  if( dim == 3 )
    {
      if( lat_grid[0] < -90 )
        throw runtime_error( 
                  "The latitude grid cannot extend below -90 degrees for 3D" );
      if( lat_grid[lat_grid.nelem() - 1] > 90 )
        throw runtime_error( 
                  "The latitude grid cannot extend above 90 degrees for 3D" );
      if( lon_grid[0] < -360 )
        throw runtime_error( 
                "No longitude (in lon_grid) can be below -360 degrees." );
      if( lon_grid[lon_grid.nelem() - 1] > 360 )
        throw runtime_error( 
                "No longitude (in lon_grid) can be above 360 degrees." );
      if( lon_grid[lon_grid.nelem() - 1]-lon_grid[0] > 360 )
        throw runtime_error( 
         "The longitude grid is not allowed to cover more than 360 degrees." );
    }
}




void chk_atm_field( 
        const String&     x_name,
        ConstTensor3View  x, 
        const Index&      dim,
        ConstVectorView   p_grid,
        ConstVectorView   lat_grid,
        ConstVectorView   lon_grid,
        const bool&       chk_lat90)
{
  // It is assumed that grids OK-ed through chk_atm_grids
  Index npages=p_grid.nelem(), nrows=1, ncols=1;
  if( dim > 1 )
    nrows = lat_grid.nelem();
  if( dim > 2 )
    ncols = lon_grid.nelem();
  if( x.ncols()!=ncols || x.nrows()!=nrows || x.npages()!=npages ) 
    {
      ostringstream os;
      os << "The atmospheric field *" << x_name <<  "* has wrong size.\n"
         << "Expected size is " << npages << " x " << nrows << " x " 
         << ncols << ", while actual size is " << x.npages() << " x " 
         << x.nrows() << " x " << x.ncols() << ".";
      throw runtime_error( os.str() );
    }

  // NaNs are not allowed
  for( Index ip=0; ip<npages; ip++ )
    { for( Index ir=0; ir<nrows; ir++ )
        { for( Index ic=0; ic<ncols; ic++ )
            {
              if( isnan( x(ip,ir,ic) ) )
                {
                  ostringstream os;
                  os << "The variable *" << x_name <<  "* contains one or "
                     << "several NaNs. This is not allowed!";
                  throw runtime_error( os.str() );
                }
            }
        }
    }

  // Special 3D checks:
  if( dim == 3  )
    {
      // If all lons are covered, check if cyclic
      if( is_lon_cyclic(lon_grid) )
        {
          const Index ic = ncols-1;
          for( Index ip=0; ip<npages; ip++ )
            {
              for( Index ir=0; ir<nrows; ir++ )
                {
                  if( !is_same_within_epsilon(x(ip,ir,ic),x(ip,ir,0),4*DBL_EPSILON) )
                    {
                      ostringstream os;
                      os << "The variable *" << x_name <<  "* covers 360 "
                         << "degrees in the longitude direction, but the field "
                         << "seems to deviate between first and last longitude "
                         << "point. The field must be \"cyclic\".\n"
                         << "Difference: " << setprecision(16) << x(ip,ir,ic)- x(ip,ir,0) << "\n"
                         << "Epsilon   : " << 4*DBL_EPSILON * max(x(ip,ir,ic),x(ip,ir,0));
                      throw runtime_error( os.str() );
                    }
                }
            }
        }

      chk_if_bool("chk_lat90", chk_lat90);
      if( chk_lat90 )
        {
          // No variation at the South pole!
          if( lat_grid[0] == -90 )
            {
              for( Index ip=0; ip<npages; ip++ )
                {
                  for( Index ic=1; ic<ncols; ic++ )
                    {
                      if( !is_same_within_epsilon(x(ip,0,ic),x(ip,0,ic-1),2*DBL_EPSILON) )
                        {
                          ostringstream os;
                          os << "The variable *" << x_name <<  "* covers the South\n"
                             << "pole. The data corresponding to the pole can not\n"
                             << "vary with longitude, but this appears to be the\n"
                             << "case.";
/*                             << "case: at " << ip << ".th  pressure it has val\n"
                             << x(ip,0,ic-1) << ", but " << x(ip,0,ic)
                             << " (i.e., a diff of " << fabs(x(ip,0,ic)-x(ip,0,ic-1))
                             << ") at " << ic-1 << "th and " << ic << "th longitudes!\n";
*/
                          throw runtime_error( os.str() );
                        }
                    }
                }
            }
          // No variation at the North pole!
          if( lat_grid[nrows-1] == 90 )
            {
              const Index ir = nrows-1;
              for( Index ip=0; ip<npages; ip++ )
                {
                  for( Index ic=1; ic<ncols; ic++ )
                    {
                      if( !is_same_within_epsilon(x(ip,ir,ic),x(ip,ir,ic-1),2*DBL_EPSILON) )
                        {
                          ostringstream os;
                          os << "The variable *" << x_name <<  "* covers the North\n"
                             << "pole. The data corresponding to the pole can not\n"
                             << "vary with longitude, but this appears to be the "
                             << "case.";
/*                             << "case: at " << ip << ".th  pressure it has val\n"
                             << x(ip,ir,ic-1) << ", but " << x(ip,ir,ic)
                             << " (i.e., a diff of " << fabs(x(ip,ir,ic)-x(ip,ir,ic-1))
                             << ") at " << ic-1 << "th and " << ic << "th longitudes!\n";
*/
                          throw runtime_error( os.str() );
                        }
                    }
                }
            }
        }
    }
}




void chk_atm_field( 
        const String&   x_name,
        ConstTensor4View  x, 
        const Index&    dim,
        const Index&    nspecies,
        ConstVectorView p_grid,
        ConstVectorView lat_grid,
        ConstVectorView lon_grid )
{
  const Index nbooks=nspecies;
  // 
  if( nbooks == 0 )
    {
      if( x.nbooks() )
        {
          ostringstream os;
          os << "The atmospheric field *" << x_name <<  "* should be empty.\n";
          throw runtime_error( os.str() );
        }
      else
        { return; }
    }

  Index npages=p_grid.nelem(), nrows=1, ncols=1;
  if( dim > 1 )
    nrows = lat_grid.nelem();
  if( dim > 2 )
    ncols = lon_grid.nelem();

  if( x.ncols()!=ncols || x.nrows()!=nrows || x.npages()!=npages ||
      x.nbooks()!=nbooks ) 
    {
      ostringstream os;
      os << "The atmospheric field *" << x_name <<  "* has wrong size.\n"
         << "Expected size is "
         << nbooks << " x " << npages << " x "
         << nrows << " x " << ncols << ",\n"
         << "while actual size is "
         << x.nbooks() << " x " << x.npages() << " x "
         << x.nrows() << " x " << x.ncols() << ".";
      throw runtime_error( os.str() );
    }

  // NaNs are not allowed
  for( Index ib=0; ib<nbooks; ib++ )
    { for( Index ip=0; ip<npages; ip++ )
      { for( Index ir=0; ir<nrows; ir++ )
          { for( Index ic=0; ic<ncols; ic++ )
              {
                if( isnan( x(ib,ip,ir,ic) ) )
                  {
                    ostringstream os;
                    os << "The variable *" << x_name <<  "* contains one or "
                       << "several NaNs. This is not allowed!";
                    throw runtime_error( os.str() );
                  }
              }
          }
      }
    }

  // Special 3D checks:
  if( dim == 3  )
    {
      // If all lons are covered, check if cyclic
      if( (lon_grid[ncols-1]-lon_grid[0]) == 360 )
        {
          const Index ic = ncols-1;
          for( Index is=0; is<nspecies; is++ )
            {
          for( Index ip=0; ip<npages; ip++ )
            {
              for( Index ir=0; ir<nrows; ir++ )
                {
                  if( !is_same_within_epsilon(x(is,ip,ir,ic),x(is,ip,ir,0),2*DBL_EPSILON) )
                    {
                      ostringstream os;
                      os << "The variable *" << x_name <<  "* covers 360 "
                         << "degrees in the longitude direction, but at least "
                         << "one field seems to deviate between first and last "
                         << "longitude point. The field must be \"cyclic\". "
                         << "This was found for field with index " 
                         << is <<" (0-based).";
                      throw runtime_error( os.str() );
                    }
                }
            }
            }
        }
      // No variation at the South pole!
      if( lat_grid[0] == -90 )
        {
          for( Index is=0; is<nspecies; is++ )
            {
          for( Index ip=0; ip<npages; ip++ )
            {
              for( Index ic=1; ic<ncols; ic++ )
                {
                  if( !is_same_within_epsilon(x(is,ip,0,ic),x(is,ip,0,ic-1),2*DBL_EPSILON) )
                    {
                      ostringstream os;
                      os << "The variable *" << x_name <<  "* covers the South "
                         << "pole. The data corresponding to the pole can not "
                         << "vary with longitude, but this appears to be the "
                         << "case. This was found for field with index " 
                         << is <<" (0-based).";
                      throw runtime_error( os.str() );
                    }
                }
            }
            }
        }
      // No variation at the North pole!
      if( lat_grid[nrows-1] == 90 )
        {
          const Index ir = nrows-1;
          for( Index is=0; is<nspecies; is++ )
            {
          for( Index ip=0; ip<npages; ip++ )
            {
              for( Index ic=1; ic<ncols; ic++ )
                {
                  if( !is_same_within_epsilon(x(is,ip,ir,ic),x(is,ip,ir,ic-1),2*DBL_EPSILON) )
                    {
                      ostringstream os;
                      os << "The variable *" << x_name <<  "* covers the North "
                         << "pole. The data corresponding to the pole can not "
                         << "vary with longitude, but this appears to be the "
                         << "case. This was found for field with index " 
                         << is <<" (0-based).";
                      throw runtime_error( os.str() );
                    }
                }
            }
            }
        }
    }
}




void chk_atm_vecfield_lat90( 
        const String&     x1_name,
        ConstTensor3View  x1, 
        const String&     x2_name,
        ConstTensor3View  x2, 
        const Index&      dim,
        ConstVectorView   lat_grid,
        const Numeric&    threshold)
{
  // It is assumed that grids OK-ed through chk_atm_grids and fields
  // individually OK-ed.

  // We only need to check 3D cases. Else there is no variation in longitude
  // anyways.
  if( dim == 3  )
    {
      Index npages= x1.npages();
      Index nrows = x1.nrows();
      Index ncols = x1.ncols();

      // For safety check that both fields have identical dimensions 
      if( x2.ncols()!=ncols || x2.nrows()!=nrows || x2.npages()!=npages ) 
        {
          ostringstream os;
          os << "The atmospheric fields *" << x1_name <<  "* and *"
             << x2_name <<  "* do not match in size.\n"
             << "*" << x1_name <<  "*'s size is " << npages << " x " << nrows
             << " x " << ncols << ", while *" << x1_name <<  "*'s size is "
             << x2.npages() << " x " << x2.nrows() << " x " << x2.ncols() << ".";
          throw runtime_error( os.str() );
        }

      // redefine ratio deviation threshold of vector lengths to ratio of
      // squared vector lengths, cause don't want to calc squareroot everytime.
      // (val1**2/val2**2 - 1) / (val1/val2 - 1) = val1/val2 + 1
      // and with val1~val2                      = 2
      // i.e., (val1**2/val2**2 - 1) ~ 2 * (val1/val2 - 1)
      //
      // with val1/1 = sqrt(vec1/2), hence val1/2**2 = vec1/2
      //       (vec1/vec2 - 1) ~ 2 * (sqrt(vec1)/sqrt(vec2) - 1)
      //       (sqrt(vec1)/sqrt(vec2) - 1) ~ (vec1/vec2 - 1) / 2
      // 
      // we want to check: sqrt(vec1)/sqrt(vec2) - 1. < threshold
      // i.e., with the above,
      //       (vec1/vec2 - 1) / 2 < threshold
      //       (vec1/vec2 - 1) < threshold*2
      Numeric th = threshold * 2.;

      // No variation at the South pole!
      if( lat_grid[0] == -90 )
        {
          Numeric vec1, vec2;
          for( Index ip=0; ip<npages; ip++ )
            {
              for( Index ic=1; ic<ncols; ic++ )
                {
                  vec1 = x1(ip,0,ic)*x1(ip,0,ic) + x2(ip,0,ic)*x2(ip,0,ic);
                  vec2 = x1(ip,0,ic-1)*x1(ip,0,ic-1) + x2(ip,0,ic-1)*x2(ip,0,ic-1);
                  if( fabs( vec1/vec2-1. ) > th )
                    {
                      ostringstream os;
                      os << "The variables *" << x1_name <<  "* and *" << x2_name
                         << "* are assumed\n"
                         << "to be two horizontal components of a vector field.\n"
                         << "At the pole, the data (here: the total length of\n"
                         << "the horizontal vector) can NOT vary with longitude,\n"
                         << "but this appears to be the case on the South pole.\n"
                         << "The threshold is " << threshold << ", but the actual\n"
                         << "deviation at pressure level " << ip << " and longitude\n"
                         << "points " << ic-1 << " and " << ic << " is "
                         << sqrt(vec1)/sqrt(vec2)-1.;
                      throw runtime_error( os.str() );
                    }
                }
            }
        }
      // No variation at the North pole!
      if( lat_grid[nrows-1] == 90 )
        {
          Numeric vec1, vec2;
          const Index ir = nrows-1;
          for( Index ip=0; ip<npages; ip++ )
            {
              for( Index ic=1; ic<ncols; ic++ )
                {
                  vec1 = x1(ip,ir,ic)*x1(ip,ir,ic) + x2(ip,ir,ic)*x2(ip,ir,ic);
                  vec2 = x1(ip,ir,ic-1)*x1(ip,ir,ic-1) + x2(ip,ir,ic-1)*x2(ip,ir,ic-1);
                  if( fabs( vec1/vec2-1. ) > th )
                    {
                      ostringstream os;
                      os << "The variables *" << x1_name <<  "* and *" << x2_name
                         << "* are assumed\n"
                         << "to be two horizontal components of a vector field.\n"
                         << "At the pole, the data (here: the total length of\n"
                         << "the horizontal vector) can NOT vary with longitude,\n"
                         << "but this appears to be the case on the North pole.\n"
                         << "The threshold is " << threshold << ", but the actual\n"
                         << "deviation at pressure level " << ip << " and longitude\n"
                         << "points " << ic-1 << " and " << ic << " is "
                         << sqrt(vec1)/sqrt(vec2)-1.;
                      throw runtime_error( os.str() );
                    }
                }
            }
        }
    }
}




void chk_latlon_true(
   const Index&      atmosphere_dim,
   ConstVectorView   lat_grid,
   ConstVectorView   lat_true,
   ConstVectorView   lon_true )
{
  if( atmosphere_dim == 1 )
    {
      if( lat_true.nelem()!=1  ||  lon_true.nelem()!=1 )
        { 
          throw runtime_error( "For 1D, the method requires that *lat_true* "
                                              "and *lon_true* have length 1." );
        }
    }
  //
  else if( atmosphere_dim == 2 ) 
    {
      if( lat_true.nelem() != lat_grid.nelem()   ||  
          lon_true.nelem() != lat_grid.nelem() )
        { 
          throw runtime_error( "For 2D, the method requires that *lat_true* "
                         "and *lon_true* have the same length as *lat_grid*." );
        }
    }
}




void chk_atm_surface( 
        const String&     x_name,
        const Matrix&     x, 
        const Index&      dim,
        ConstVectorView   lat_grid,
        ConstVectorView   lon_grid )
{
  Index ncols=1, nrows=1;
  if( dim > 1 )
    nrows = lat_grid.nelem();
  if( dim > 2 )
    ncols = lon_grid.nelem();
  if( x.ncols()!=ncols || x.nrows()!=nrows ) 
    {
      ostringstream os;
      os << "The surface variable *" << x_name <<  "* has wrong size.\n"
         << "Expected size is " << nrows << " x " << ncols << ","
         << " while actual size is " << x.nrows() << " x " << x.ncols() << ".";
      throw runtime_error( os.str() );
    }

  // Special 3D checks:
  if( dim == 3  )
    {
      // If all lons are covered, check if cyclic
      if( (lon_grid[ncols-1]-lon_grid[0]) == 360 )
        {
          const Index ic = ncols-1;
          for( Index ir=0; ir<nrows; ir++ )
            {
              if( !is_same_within_epsilon(x(ir,ic),x(ir,0),2*DBL_EPSILON) )
                {
                  ostringstream os;
                  os << "The variable *" << x_name <<  "* covers 360 "
                     << "degrees in the longitude direction, but the field "
                     << "seems to deviate between first and last longitude "
                     << "point. The field must be \"cyclic\".";
                  throw runtime_error( os.str() );
                }
            }
        }

      // No variation at the South pole!
      if( lat_grid[0] == -90 )
        {
          for( Index ic=1; ic<ncols; ic++ )
            {
              if( !is_same_within_epsilon(x(0,ic),x(0,ic-1),2*DBL_EPSILON) )
                {
                  ostringstream os;
                  os << "The variable *" << x_name <<  "* covers the South "
                     << "pole. The data corresponding to the pole can not "
                     << "vary with longitude, but this appears to be the "
                     << "case.";
                  throw runtime_error( os.str() );
                }
            }
        }
      // No variation at the North pole!
      if( lat_grid[nrows-1] == 90 )
        {
          const Index ir = nrows-1;
          for( Index ic=1; ic<ncols; ic++ )
            {
              if( !is_same_within_epsilon(x(ir,ic),x(ir,ic-1),2*DBL_EPSILON) )
                {
                  ostringstream os;
                  os << "The variable *" << x_name <<  "* covers the North "
                     << "pole. The data corresponding to the pole can not "
                     << "vary with longitude, but this appears to be the "
                     << "case.";
                  throw runtime_error( os.str() );
                }
            }
        }
    }
}





/*===========================================================================
  === Functions related to sensor variables.
  ===========================================================================*/


void chk_rte_pos( 
        const Index&      atmosphere_dim,
        ConstVectorView   rte_pos,
        const bool&       is_rte_pos2 )

{  
  String vname = "*rte_pos*";
  if( is_rte_pos2 )
    { vname = "*rte_pos2*"; }

  if( atmosphere_dim == 1 )
    {
      if( !is_rte_pos2 )
        {
          if( rte_pos.nelem() != 1 )
            {
              ostringstream os;
              os << "For 1D, " << vname << " must have length 1.";
              throw runtime_error(os.str());
            }
        }
      else
        {
          if( rte_pos.nelem() != 2 )
            {
              ostringstream os;
              os << "For 1D, " << vname << " must have length 2.";
              throw runtime_error(os.str());
            }
          if( rte_pos[1] < -180  ||  rte_pos[1] > 180 )
            {
              ostringstream os;
              os << "For 1D, the latitude in " << vname << " must be in the "
                 << "range [-180,180].";
              throw runtime_error(os.str());
            }
        }
    }
  else if( atmosphere_dim == 2 )
    { 
      if( rte_pos.nelem() != 2 )
        {
          ostringstream os;
          os << "For 2D, " << vname << " must have length 2.";
          throw runtime_error(os.str());
        }
    }
  else
    {
      if( rte_pos.nelem() != 3 )
        {
          ostringstream os;
          os << "For 3D, " << vname << " must have length 3.";
          throw runtime_error(os.str());
        }
      if( rte_pos[1] < -90  ||  rte_pos[1] > 90 )
        {
          ostringstream os;
          os << "The (3D) latitude in " << vname << " must be in the "
             << "range [-90,90].";
          throw runtime_error(os.str());
        }
      if( rte_pos[2] < -360  ||  rte_pos[2] > 360 )
        {
          ostringstream os;
          os << "The longitude in " << vname << " must be in the "
             << "range [-360,360].";
          throw runtime_error(os.str());
        } 
    }
}




void chk_rte_los( 
        const Index&      atmosphere_dim,
        ConstVectorView   rte_los )

{  
  if( atmosphere_dim == 1 )
    {
      if( rte_los.nelem() != 1 )
        { throw runtime_error( "For 1D, *rte_los* must have length 1." ); }
      if( rte_los[0] < 0  ||  rte_los[0] > 180 )
        { throw runtime_error( "For 1D, the zenith angle of *rte_los* must "
                               "be in the range [0,180]." ); }
    }
  else if( atmosphere_dim == 2 )
    { 
      if( rte_los.nelem() != 1 )
        { throw runtime_error( "For 2D, *rte_los* must have length 1." ); }
      if( rte_los[0] < -180  ||  rte_los[0] > 180 )
        { throw runtime_error( "For 2D, the zenith angle of *rte_los* must "
                               "be in the range [-180,180]." ); }
    }
  else
    {
      if( rte_los.nelem() != 2 )
        { throw runtime_error( "For 3D, *rte_los* must have length 2." ); }
      if( rte_los[0] < 0  ||  rte_los[0] > 180 )
        { throw runtime_error( "For 3D, the zenith angle of *rte_los* must "
                               "be in the range [0,180]." ); }
      if( rte_los[1] < -180  ||  rte_los[1] > 180 )
        { throw runtime_error( "For 3D, the azimuth angle of *rte_los* must "
                               "be in the range [-180,180]." ); }
    }
}


/*===========================================================================
 === Functions related to GriddedFields.
 ===========================================================================*/


void chk_griddedfield_gridname(const GriddedField& gf,
                               const Index gridindex,
                               const String& gridname)
{
  if (gf.get_dim()-1 < gridindex)
  {
    ostringstream os;
    os << "Grid index " << gridindex << " exceeds dimension of GriddedField";
    if (gf.get_name().nelem()) os << " \"" << gf.get_name() << "\"";
    throw runtime_error(os.str());
  }

  String gfgridnameupper = gf.get_grid_name(gridindex);
  gfgridnameupper.toupper();

  String gridnameupper = gridname;
  gridnameupper.toupper();

  if (gfgridnameupper != gridnameupper)
  {
    ostringstream os;
    os << "Name of grid " << gridindex << " in GriddedField";
    if (gf.get_name().nelem()) os << " \"" << gf.get_name() << "\"";
    os << " is \"" << gf.get_grid_name(gridindex) << "\".\n"
    << "The expected name is \"" << gridname << "\".";
    throw runtime_error(os.str());
  }
}





/*===========================================================================
 === Functions checking agendas
 ===========================================================================*/


bool chk_if_std_blackbody_agenda( 
         Workspace&   ws,
   const Agenda&      blackbody_radiation_agenda )
{
  Vector btest;
  blackbody_radiation_agendaExecute( ws, btest, 300, Vector(1,100e9),
                                     blackbody_radiation_agenda );
  if( abs( btest[0] - 9.1435e-16) > 1e-19 )
    { return false; }
  else
    { return true; }
}