Integrator.h

#ifndef QNANO_NEW_INTEGRATOR_DEFINED_H
#define QNANO_NEW_INTEGRATOR_DEFINED_H

#include "tools/Integrable_Function.h"
#include "tools/List_Class.h"

class Integrator_Discretization_Parameter{
public: 
  double lim_a, lim_b;
  size_t N;
};
class Integrator_Discretization: public List_Class<Integrator_Discretization_Parameter>{
public: 

  Integrator_Discretization(size_t D, double lim_a, double lim_b, size_t N){
    list.resize(D);
    for(size_t i=0; i<D; i++){
      list[i].lim_a=lim_a;
      list[i].lim_b=lim_b;
      list[i].N=N;
    }
  }
};

class Integrator{
public:
  virtual double integrate(const Integrable_Function& ifct, const Integrator_Discretization& discr)const=0;
};

template<int D> class Integrator_Riemann: public Integrator{
public:
  virtual double integrate(const Integrable_Function& ifct, const Integrator_Discretization& discr, std::vector<double> &x)const{

    if(discr.size() < D){
      std::cerr<<"Integrator_Riemann: discr.size()="<<discr.size()<<" < D="<<D<<std::endl;
      exit(1);
    }
    if(x.size() <D){
      std::cerr<<"Integrator_Riemann: x.size()="<<x.size()<<" < D="<<D<<std::endl;
      exit(1);
    }

    
    double h=(discr[D-1].lim_b-discr[D-1].lim_a)/discr[D-1].N;
    double res=0.;
    for(size_t i=0; i<discr[D-1].N; i++){
      x[D-1] = discr[D-1].lim_a + i * h;
      res += h * ( Integrator_Riemann<D-1>().integrate(ifct, discr, x) );
    }
    return res;
  }
  virtual double integrate(const Integrable_Function& ifct, const Integrator_Discretization& discr)const{
    std::vector<double> x(D,0);
    return integrate(ifct, discr, x);
  }
};

template<> class Integrator_Riemann<0>{
public:
  double integrate(const Integrable_Function& ifct, const Integrator_Discretization& discr, std::vector<double> &x)const{
    return ifct.evaluate(x);
  }
};

template<int D> class Integrator_Simpson: public Integrator{
public:
  virtual double integrate(const Integrable_Function& ifct, const Integrator_Discretization& discr, std::vector<double> &x)const{

    if(discr.size() < D){
      std::cerr<<"Integrator_Riemann: discr.size()="<<discr.size()<<" < D="<<D<<std::endl;
      exit(1);
    }
    if(x.size() <D){
      std::cerr<<"Integrator_Riemann: x.size()="<<x.size()<<" < D="<<D<<std::endl;
      exit(1);
    }

    double h=(discr[D-1].lim_b-discr[D-1].lim_a)/discr[D-1].N;
    double res=0.;

    x[D-1]=discr[D-1].lim_a;
    res+=h/3*Integrator_Simpson<D-1>().integrate(ifct, discr, x);

    x[D-1]=discr[D-1].lim_b;
    res+=h/3*Integrator_Simpson<D-1>().integrate(ifct, discr, x);

    double pref=4./3.*h;
    for(size_t i=1; i<discr[D-1].N; i+=2){
      x[D-1]=discr[D-1].lim_a+h*i;
      res+=pref*Integrator_Simpson<D-1>().integrate(ifct, discr, x);
    }
    pref=2./3.*h;
    for(size_t i=2; i<discr[D-1].N-1; i+=2){
      x[D-1]=discr[D-1].lim_a+h*i;
      res+=pref*Integrator_Simpson<D-1>().integrate(ifct, discr, x);
    }
    return res;
  }

  virtual double integrate(const Integrable_Function& ifct, const Integrator_Discretization& discr)const{
    std::vector<double> x(D,0);
    return integrate(ifct, discr, x);
  }
};
template <> class Integrator_Simpson<0>{
public:
  double integrate(const Integrable_Function& ifct, const Integrator_Discretization& discr, std::vector<double> &x) const{
    return ifct.evaluate(x);
  }
};





#endif