math_funcs.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 <iostream>
#include <cmath>
#include <stdexcept>
#include "array.h"
#include "math_funcs.h"
#include "logic.h"
#include "mystring.h"
extern const Numeric DEG2RAD;
extern const Numeric PI;



/*****************************************************************************
 *** The functions (in alphabetical order)
 *****************************************************************************/


Numeric fac(const Index n)
{
  Numeric sum;

  if (n == 0) return (1.0);

  sum = 1.0;
  for (Index i = 1; i <= n; i++)
    sum *= Numeric(i);

  return(sum);
}



Index integer_div( const Index& x, const Index& y )
{
  assert( is_multiple( x, y ) );
  return x/y;
}




Numeric LagrangeInterpol4( ConstVectorView x,
                           ConstVectorView y,
                           const Numeric a)
{
  // lowermost grid spacing on x-axis
  const Numeric Dlimit = 1.00000e-15;

  // Check that dimensions of x and y vector agree
  const Index n_x = x.nelem();
  const Index n_y = y.nelem();
  if ( (n_x != 4) || (n_y != 4) )
    {
      ostringstream os;
      os << "The vectors x and y must all have the same length of 4 elements!\n"
        << "Actual lengths:\n"
        << "x:" << n_x << ", " << "y:" << n_y << ".";
      throw runtime_error(os.str());
    }

  // assure that x1 =< a < x2 holds
  if ( (a < x[1]) || (a > x[2]) )
    {
      ostringstream os;
      os << "LagrangeInterpol4: the relation x[1] =< a < x[2] is not satisfied. " 
         << "No interpolation can be calculated.\n";
      throw runtime_error(os.str());
    };

  // calculate the Lagrange polynomial coefficients for a polynomial of the order of 3
  Numeric b[4];
  for (Index i=0 ; i < 4 ; ++i)
    {
      b[i] = 1.000e0;
      for (Index k=0 ; k < 4 ; ++k)
        {
          if ( (k != i) && (fabs(x[i]-x[k]) > Dlimit) )  
            b[i] = b[i] * ( (a-x[k]) / (x[i]-x[k]) );
        };
    };

  Numeric ya = 0.000e0;
  for (Index i=0 ; i < n_x ; ++i) ya = ya + b[i]*y[i];

  return ya;
}





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




Index last( const ArrayOfIndex& x )
{
  assert( x.nelem() > 0 );
  return x[x.nelem()-1]; 
}




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 + (double)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)/((double)n-1) ;
  for ( Index i=0; i<n-1; i++ )
    x[i] = start + (double)i*step;
  x[n-1] = stop;
}




void nlogspace(         
               Vector&     x, 
               const Numeric     start,     
               const Numeric     stop,        
               const Index         n )
{
  // Number of points must be greater than 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)/((double)n-1);
  x[0] = start;
  for ( Index i=1; i<n-1; i++ )
    x[i] = exp(a + (double)i*step);
  x[n-1] = stop;
}



Numeric AngIntegrate_trapezoid(ConstMatrixView Integrand,
                               ConstVectorView za_grid,
                               ConstVectorView aa_grid)
{

  Index n = za_grid.nelem();
  Index m = aa_grid.nelem();
  Vector res1(n);
  assert (is_size(Integrand, n, m));
  
  for (Index i = 0; i < n ; ++i)
    {
      res1[i] = 0.0;
      
      for (Index j = 0; j < m - 1; ++j)
        {
          res1[i] +=  0.5 * DEG2RAD * (Integrand(i, j) + Integrand(i, j + 1)) *
            (aa_grid[j + 1] - aa_grid[j]) * sin(za_grid[i] * DEG2RAD);
        }
    }
  Numeric res = 0.0;
  for (Index i = 0; i < n - 1; ++i)
    {
      res += 0.5 * DEG2RAD * (res1[i] + res1[i + 1]) * 
        (za_grid[i + 1] - za_grid[i]);
    }
  
  return res;
}



Numeric AngIntegrate_trapezoid_opti(ConstMatrixView Integrand,
                                    ConstVectorView za_grid,
                                    ConstVectorView aa_grid,
                                    ConstVectorView grid_stepsize)
{
  Numeric res = 0;
  if ((grid_stepsize[0] > 0) && (grid_stepsize[1] > 0))
    {
      Index n = za_grid.nelem();
      Index m = aa_grid.nelem();
      Numeric stepsize_za = grid_stepsize[0];
      Numeric stepsize_aa = grid_stepsize[1];
      Vector res1(n);
      assert (is_size(Integrand, n, m));

      Numeric temp = 0.0;
      
      for (Index i = 0; i < n ; ++i)
        {
          temp = Integrand(i, 0);
          for (Index j = 1; j < m - 1; j++)
            {
              temp += Integrand(i, j) * 2;
            }
          temp += Integrand(i, m-1);
          temp *= 0.5 * DEG2RAD * stepsize_aa * sin(za_grid[i] * DEG2RAD);
          res1[i] = temp;
        }

      res = res1[0];
      for (Index i = 1; i < n - 1; i++)
        {
          res += res1[i] * 2;
        }
      res += res1[n-1];
      res *= 0.5 * DEG2RAD * stepsize_za;
    }
  else
    {
      res = AngIntegrate_trapezoid(Integrand, za_grid, aa_grid);
    }

  return res;
}



Numeric AngIntegrate_trapezoid(ConstVectorView Integrand,
                               ConstVectorView za_grid)
{

  Index n = za_grid.nelem();
  assert (is_size(Integrand, n));
  
  Numeric res = 0.0;
  for (Index i = 0; i < n - 1; ++i)
    {
      // in this place 0.5 * 2 * PI is calculated:
      res += PI * DEG2RAD * (Integrand[i]* sin(za_grid[i] * DEG2RAD) 
                             + Integrand[i + 1] * sin(za_grid[i + 1] * DEG2RAD))
        * (za_grid[i + 1] - za_grid[i]);
    }
  
  return res;
}





Numeric sign( const Numeric& x )
{
  if( x < 0 )
    return -1.0;
  else if( x == 0 )
    return 0.0;
  else
    return 1.0;
}



Numeric gamma_func(Numeric xx)
{
  //double lgamma(double xx);
  
  Numeric gam;
  Index i;
  
  
  if (xx > 0.0) {
    if (xx == (int)xx) {
      gam = 1.0;               // use factorial
      for (i=2;i<xx;i++) {
        gam *= (Numeric)i;
      }
    }    
    else {       
            return exp(lgamma_func(xx));
    }
  } else {
    ostringstream os;
    os << "Argument is zero or negative."
    << "Gamma function can not be calculated.\n";
    throw runtime_error(os.str());
  }
  
  
  return gam;
}



Numeric lgamma_func(Numeric xx)
{
  
  Numeric x,y,tmp,ser;
  static const Numeric cof[6] = {
    76.18009172947146, -86.50532032941677, 24.01409824083091,
    -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5
  };
  
  if (xx > 0.0)
  {
    y=x=xx;
    tmp = x+5.5;
    tmp -= (x+0.5)*log(tmp);
    ser = 1.000000000190015;
    for (Index j=0;j<=5;j++) ser += cof[j]/++y;
    return -tmp+log(2.5066282746310005*ser/x);
  }
  else
  {
    ostringstream os;
    os << "Argument is zero or negative.\n"
    << "log Gamma function can not be calculated.\n";
    throw runtime_error(os.str());
  }
}




void unitl( Vector& x )
{
  assert( x.nelem() > 0 );
 
  const Numeric l = sqrt(x*x);
  for(Index i=0; i<x.nelem(); i++ )
    x[i] /= l;
}