b2integration1d.H

//------------------------------------------------------------------------
// b2integration1d.H --
//
//
// written by Mathias Doreille
//
// Copyright (c) 2007-2012 SMR Engineering & Development SA
//               2502 Bienne, Switzerland
//
// All Rights Reserved.  Proprietary source code.  The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR.
//------------------------------------------------------------------------
 
#ifndef B2INTEGRATION1D_H_
#define B2INTEGRATION1D_H_
 
#include "b2exception.H"
 
typedef double(quadpack_function_t)(const double*);
 
extern "C" void FC_GLOBAL(dqagse, DQAGSE)(
      quadpack_function_t f, double* a, double* b, double* epsabs, double* epsrel, int* limit,
      double* result, double* abserr, int* neval, int* ier, double* alist, double* blist,
      double* rlist, double* elist, int* iord, int* last);
 
namespace b2000 { namespace quadpack {
 
template <typename F>
double simpson_integration(const F& f, double a, double b, int n) {
    if (n % 2 != 0) { n += 1; }
    double r = (f(a) + f(b)) / 2;
    const double h = (b - a) / n;
    a += h;
    for (int i = 1; i != n; ++i, a += h) { r += (1 + i % 2) * f(a); }
    return 2 * r * h / 3;
}
 
template <typename F>
class FOToF {
public:
    static void set(const F* f) { f_static = f; }
    static double f(const double* x) { return (*f_static)(*x); }
 
private:
    // Put it in thread local variable for a thread safe implementation.
    static const F* f_static;
};
 
template <typename F>
const F* FOToF<F>::f_static;
 
class Exception : public b2000::Exception {
public:
    Exception(int ier) {
        switch (ier) {
            case 1:
                set_msg("Maximum number of subdivisions allowed has been achieved");
                return;
            case 2:
                set_msg(
                      "The occurrence of roundoff error is detected, "
                      "which prevents the requested tolerance from being achieved.");
                return;
            case 3:
                set_msg(
                      "Extremely bad integrand behaviour occurs at some "
                      "points of the integration interval.");
                return;
            case 4:
                set_msg("The algorithm does not converge.");
                return;
            case 5:
                set_msg("The integral is probably divergent, or slowly convergent.");
                return;
            case 6:
                set_msg("The input is invali.");
                return;
            default:
                set_msg("The integral is probably divergent, or slowly convergent.");
                return;
        }
    }
};
 
template <typename F>
double qags(
      const F& f, double a, double b, double epsabs, double epsrel, int subintevals_limit,
      double& abserr, int& nb_subinterval) {
    double res;
    int neval;
    int ier;
    double* work = new double[subintevals_limit * 4];
    int* iwork = new int[subintevals_limit];
    FOToF<F>::set(&f);
    FC_GLOBAL(dqagse, DQAGSE)
    (FOToF<F>::f, &a, &b, &epsabs, &epsrel, &subintevals_limit, &res, &abserr, &neval, &ier, work,
     work + subintevals_limit, work + 2 * subintevals_limit, work + 3 * subintevals_limit, iwork,
     &nb_subinterval);
    nb_subinterval = neval;
    delete[] work;
    delete[] iwork;
    if (ier > 1) { Exception(ier) << THROW; }
    return res;
}
 
}}  // namespace b2000::quadpack
 
#endif /*B2INTEGRATION1D_H_*/