integrate.hpp

Go to the documentation of this file.

#pragma once
/*  */

#include "gsl/gsl_integration.h"
#include <new>

namespace QQbar_threshold{


  template<class Integrand>
  inline
  double integrate(Integrand I, double x_min, double x_max);

  namespace detail{
    template<class Functor>
    inline
    double call_from_gsl(double t, void * p){
      Functor * f = static_cast<Functor *>(p);
      return (*f)(t);
    }
  }


  template<class Integrand>
  class integrator{
    static constexpr size_t workspace_size = 200;
    Integrand i;
    gsl_integration_cquad_workspace * space;
    double abserr;
    size_t nevals;
  public:

    integrator(Integrand in):
      i(in),
      space(gsl_integration_cquad_workspace_alloc(workspace_size)){
      if(space == nullptr) throw std::bad_alloc();
    }
    ~integrator(){
      gsl_integration_cquad_workspace_free(space);
    }

    double integrate(
    double a, double b,
    double epsabs = 10e-8, double epsrel = 10e-6
    ){
      double result;
      gsl_function F;
      F.function = &detail::call_from_gsl<Integrand>;
      F.params = &i;
      gsl_integration_cquad(
      &F, a, b, epsabs, epsrel, space, &result, &abserr, &nevals
      );
      return result;
    }
    double get_error() const {
      return abserr;
    }
    double get_num_evals() const {
      return nevals;
    }
  };


  template<class Integrand>
  inline
  integrator<Integrand> make_integrator(Integrand I){
    return integrator<Integrand>{I};
  }

  template<class Integrand>
  inline
  double integrate(Integrand I, double x_min, double x_max){
    return integrator<Integrand>{I}.integrate(x_min, x_max);
  }

}