fp_gsl_integrate.cc

Go to the documentation of this file.

#include "fp_gsl_integrate.h"
#include "fp_exception.h"
#include <boost/foreach.hpp>
using namespace FullPhysics;

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

GslIntegrate::GslIntegrate(int Max_intervals)
  : max_ninterval(1000)
{
  w = gsl_integration_workspace_alloc(Max_intervals);
}

GslIntegrate::~GslIntegrate()
{
  gsl_integration_workspace_free(w);
}


//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

double GslIntegrate::integrate
(const boost::function<double (double)>& F, double xmin, double xmax, 
 double eps_abs, double eps_rel, int key) const
{
  double res, error_est;
  integrate_err_est(F, xmin, xmax,res, error_est, eps_abs, eps_rel, key);
  return res;
}

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

double GslIntegrate::integrate
(const boost::function<double (double)>& F, double xmin, double xmax, 
 const std::vector<double>& breakpoints,
 double eps_abs, double eps_rel, int key) const
{
  double res, error_est;
  integrate_err_est(F, xmin, xmax, breakpoints, res, error_est, 
                    eps_abs, eps_rel, key);
  return res;
}

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

double gsl_integrate_integrand_wrapper(double x, void* params)
{
  const boost::function<double (double)>* f = 
    (const boost::function<double (double)>*) params;
  return (*f)(x);
}

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

void GslIntegrate::integrate_err_est
(const boost::function<double (double)>& F,
 double xmin, double xmax, double &Res, double& Error_est,
 double eps_abs, double eps_rel,
 int key) const
{
  gsl_function gf;
  gf.function = &gsl_integrate_integrand_wrapper;
  gf.params = (void *) &F;
  int status = gsl_integration_qags(&gf, xmin, xmax, eps_abs, eps_rel,
                                    max_ninterval, w, &Res, &Error_est);
  gsl_check(status);
}

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

void GslIntegrate::integrate_err_est
(const boost::function<double (double)>& F,
 double xmin, double xmax, const std::vector<double>& breakpoints,
 double &Res, double& Error_est,
 double eps_abs, double eps_rel,
 int key) const
{
  gsl_function gf;
  gf.function = &gsl_integrate_integrand_wrapper;
  gf.params = (void *) &F;
  std::vector<double> bp;
  bp.push_back(xmin);
  BOOST_FOREACH(double x, breakpoints)
    if(x > xmin && x < xmax)
      bp.push_back(x);
  bp.push_back(xmax);
  std::sort(bp.begin(), bp.end());
  int status = gsl_integration_qagp(&gf, &bp[0], bp.size(),
                                    eps_abs, eps_rel,
                                    max_ninterval, w, &Res, &Error_est);
  gsl_check(status);
}