/*
    C++ Mesh Generation Library
    Copyright (c) 2000echet <eric at bechet dot ca>

    This file is part of the C++ Mesh Generation Library.

    See the NOTICE & LICENSE files for conditions.
*/
//---------------------------------------------------------------------------
#ifndef INTEGRATORS_H
#define INTEGRATORS_H
//---------------------------------------------------------------------------


#include "mesh_const.h"

///integrateur stupide a 1 point d'integration : F(b)-F(a)=f((a+b)/2)*((b-a))
template <class F> class Super_Simple_Integrator
{

///
  F* Func;

public:

///
  Super_Simple_Integrator ( F& function )
  {
    Func=&function;
  }

///
  double operator () ( double a, double b )
  {
    return ( ( *Func ) ( ( a+b ) /2 ) * ( b-a ) );
  }
};


///integrateur d'ordre 1 � 2 points d'integration : F(b)-F(a)=(f(a)+f(b))*((b-a)/2)
template <class F> class Simple_Integrator
{

///
  F* Func;

public:

///
  Simple_Integrator ( F& function )
  {
    Func=&function;
  }

///
  double operator () ( double a, double b )
  {
    return ( ( ( *Func ) ( a ) + ( *Func ) ( b ) ) * ( b-a ) /2. );
  }
};


/// integrateur de simpson (2nd ordre a n points d'integration)
template<class F> class Simpson_Integrator
{

///
  F* Func;

public:

///
  Simpson_Integrator ( F& function )
  {
    Func=&function;
  }

///
  double operator () ( double a, double b, int n )
  {

    n*=2;
    double h = ( b-a ) /n;
    double sum1 = 0.0;
    double sum2 = 0.0;
    double x;
    int i;

    x=a+2.0*h;

    for ( i=1; i<=n/2-1; ++i )
    {
      sum1+= ( *Func ) ( x );
      x+=2.0*h;
    }

    x=a+h;

    for ( i=1; i<=n/2; ++i )
    {
      sum2+= ( *Func ) ( x );
      x+=2.0*h;
    }

    return h/3.0* ( ( *Func ) ( a ) +2.0*sum1+4.0*sum2+ ( *Func ) ( b ) );
  }
};




#endif // INTEGRATORS_H