math_funcs.cc

Go to the documentation of this file.

/* Copyright (C) 2000, 2001 Patrick Eriksson <patrick@rss.chalmers.se>
                            Stefan Buehler   <sbuehler@uni-bremen.de>

   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 <time.h>
#include <math.h>
#include <stdexcept>
#include "arts.h"
#include "math_funcs.h"
#include "make_vector.h"
#include "array.h"



//
Numeric first( ConstVectorView x )
{
  return x[0]; 
}

Numeric last( ConstVectorView x )
{
  return x[x.nelem()-1]; 
}



//   Logical functions

//
bool any( const ArrayOfIndex& x ) 
{
  for ( Index i=0; i<x.nelem(); i++ ) {
    if ( x[i] )
      return true;
  }
  return false;
}




//   Functions to generate vectors

//
void linspace(                      
              Vector&     x,           
              const Numeric     start,    
              const Numeric     stop,        
              const Numeric     step )
{
  Index n = (Index) floor( (stop-start)/step ) + 1;
  if ( n<1 )
    n=1;
  x.resize(n);
  for ( Index i=0; i<n; i++ )
    x[i] = start + i*step;
}

void nlinspace(         
               Vector&     x, 
               const Numeric     start,     
               const Numeric     stop,        
               const Index       n )
{
  assert( 1<n );                // Number of points must be greatere 1.
  x.resize(n);
  Numeric step = (stop-start)/(n-1) ;
  for ( Index i=0; i<n; i++ )
    x[i] = start + i*step;
}

//
void nlogspace(         
               Vector&     x, 
               const Numeric     start,     
               const Numeric     stop,        
               const Index         n )
{
  // Number of points must be greatere 1:
  assert( 1<n );        
  // Only positive numbers are allowed for start and stop:
  assert( 0<start );
  assert( 0<stop );

  x.resize(n);
  Numeric a = log(start);
  Numeric step = (log(stop)-a)/(n-1);
  x[0] = start;
  for ( Index i=1; i<n-1; i++ )
    x[i] = exp(a + i*step);
  x[n-1] = stop;
}

Vector nlogspace(  
                 const Numeric start, 
                 const Numeric stop,  
                 const Index     n )
{
  Vector x;
  nlogspace( x, start, stop, n );
  return x; 
}                     




//   Interpolation routines

Index interp_check( ConstVectorView  x,        
                    ConstVectorView  xi,
                    const Index   n_y )
{
  const Index  n  = x.nelem();
  const Index  ni = xi.nelem();
  Index  order=1;          // flag for pos. or neg. order, assume pos.

  // Determine the order, -1=decreasing and 1=increasing
  if ( x[0] > x[n-1] )
    order = -1;

  if ( n < 2 )
    throw runtime_error("Vector length for interpolation must be >= 2.");

  if ( n != n_y ) 
    throw runtime_error("Sizes of input data to interpolation do not match.");

  // Make lower and upper bounds tolerate the new point to be half a
  // grid distance outside the original grid.
  const Numeric lower_bound = x[0]   - 0.5*(x[1]-x[0]);
  const Numeric upper_bound = x[n-1] + 0.5*(x[n-1]-x[n-2]);

  if ( (order*xi[0]<order*lower_bound) || (order*xi[ni-1]>order*upper_bound) ) 
    {
      ostringstream os;
      os << "Interpolation points must be not more than\n"
         << "half a grid spacing outside the original range.\n"
         << "Int.:  xi[0] = " << xi[0] << ", xi[ni-1] = " << xi[ni-1] << '\n'
         << "Orig.: x[0]  = " << x[0]  << ", x[n-1]   = " << x[n-1];
      throw runtime_error(os.str());
    }

  for ( Index i=0; i<n-1; i++ )
  {
    if ( order*x[i+1] < order*x[i] ) 
      throw runtime_error("Original interpolation grid must be ordered");
  }

  for ( Index i=0; i<ni-1; i++ )
  {
    if ( order*xi[i+1] < order*xi[i] ) 
      throw runtime_error("Interpolation points must be ordered");
  }

  return order;
}

void interp_lin_vector( VectorView       yi,
                        ConstVectorView  x, 
                        ConstVectorView  y, 
                        ConstVectorView  xi )
{
  // Check grids and get order of grids
  Index order = interp_check( x, xi, y.nelem() ); 

  Index        i, j=0;
  const Index  n=xi.nelem();
  const Index  nx=x.nelem();
  Numeric      w;

  // Check that output vector has the right size:
  assert( n==yi.nelem() ); 

  for ( i=0; i<n; i++ )
  {
    // Find right place:
    while ( j<nx-2 && order*x[j+1]<order*xi[i] )
      {
        ++j;
      }
    // j should now point to the place in the old grid below the
    // interpolation point. If the interpolation point is below the
    // original grid, j=0. If the interpolation point is above the
    // original grid, j=nx-2.
    
    assert( x[j+1]!=x[j] );
    w = (xi[i]-x[j]) / (x[j+1]-x[j]);
    // This expression should also be correct for points just outside
    // the original grid. In that case w can be negative or larger
    // than 1.
    yi[i] = y[j] + w * (y[j+1]-y[j]); 
  }
}      

void interp_lin_matrix(    
                       MatrixView        Yi,
                       ConstVectorView   x, 
                       ConstMatrixView   Y, 
                       ConstVectorView   xi )
{
  // Check grids and get order of grids
  Index order = interp_check( x, xi, Y.ncols() ); 

  Index j=0;
  const Index n=xi.nelem(), nrow=Y.nrows(), nx=x.nelem();
  Numeric w;

  assert( nrow == Yi.nrows() );
  assert( n    == Yi.ncols() );

  for (Index i=0; i<n; i++ )
  {
    // Find right place:
    while ( j<nx-2 && order*x[j+1]<order*xi[i] )
      {
        ++j;
      }
    // j should now point to the place in the old grid below the
    // interpolation point. If the interpolation point is below the
    // original grid, j=0. If the interpolation point is above the
    // original grid, j=nx-2.

    assert( x[j+1]!=x[j] );
    w = (xi[i]-x[j]) / (x[j+1]-x[j]);
    // This expression should also be correct for points just outside
    // the original grid. In that case w can be negative or larger
    // than 1.
    for( Index k=0; k<nrow; k++ )
    {
      Yi(k,i) = Y(k,j) + w * (Y(k,j+1)-Y(k,j));
    }
  }
}        

Numeric interp_lin(
                   ConstVectorView  x, 
                   ConstVectorView  y, 
                   const Numeric  xi )
{
  Vector Yi(1);
  MakeVector Xi(xi);   // MakeVector is a special kind of vector that
                       // can be initialized explicitly with one or
                       // more arguments of type Numeric. 

  interp_lin_vector( Yi, x, y, Xi );
  return Yi[0];
}        




//   Check of function input

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



//
void check_if_in_range( 
   const Numeric& x_low, 
   const Numeric& x_high, 
   const Numeric& x, 
   const String&  x_name ) 
{
  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 check_lengths( const Vector& x1, const String& x1_name,
                    const Vector& x2, const String& x2_name ) 
{
  if ( x1.nelem() != x2.nelem() )
  {
    ostringstream os;
    os << "The vectors *" << x1_name <<  "* and *" << x2_name << "*\n"
       << "must have the same lengths. \nThe lengths are: \n"
       << x1_name << ": " << x1.nelem() << "\n"
       << x2_name << ": " << x2.nelem();
    throw runtime_error( os.str() );
  }
}



//
void check_length_nrow( const Vector& x, const String& x_name,
                        const Matrix& A, const String& A_name ) 
{
  if ( x.nelem() != A.nrows() )
  {
    ostringstream os;
    os << "The length of vector *" << x_name <<  "* must be the\n"
       << "same as the number of rows of *" << A_name << "*.\n"
       << "The length of *" << x_name <<  "* is " << x.nelem() << ".\n"
       << "The number of rows of *" << A_name <<  "* is " << A.nrows() 
       << ".";
    throw runtime_error( os.str() );
  }
}



//
void check_length_ncol( const Vector& x, const String& x_name,
                        const Matrix& A, const String& A_name ) 
{
  if ( x.nelem() != A.ncols() )
  {
    ostringstream os;
    os << "The length of vector *" << x_name <<  "* must be the\n"
       << "same as the number of columns of *" << A_name << "*.\n"
       << "The length of *" << x_name <<  "* is " << x.nelem() << ".\n"
       << "The number of columns of *" << A_name <<  "* is " << A.ncols()
       << ".";
    throw runtime_error( os.str() );
  }
}

//
void check_ncol_nrow( const Matrix& A, const String& A_name,
                      const Matrix& B, const String& B_name ) 
{
  if ( A.ncols() != B.nrows() )
  {
    ostringstream os;
    os << "The number of columns of *" << A_name << "* must be the\n"
       << "same as the number of rows of *" << B_name << "*."
       << "The number of columns of *" << A_name <<  "* is " << A.ncols()
       << ".\n"
       << "The number of rows of *" << B_name <<  "* is " << B.nrows()
       << ".";
    throw runtime_error( os.str() );
  }
}

---