source: Sophya/trunk/SophyaLib/NTools/integ.h@ 1516

// This may look like C code, but it is really -*- C++ -*-
// Simple integrator 
//        E. Aubourg        1996-...
// DAPNIA/SPP (Saclay) / CEA    LAL - IN2P3/CNRS  (Orsay)

#ifndef INTEG_H_SEEN
#define INTEG_H_SEEN

#include "machdefs.h"
#include "pexceptions.h"
#include <set>

namespace SOPHYA {

class GeneralFunction;

class Integrator {
public:
  typedef double(*FUNC)(double);

  Integrator();
  Integrator(FUNC, double xmin, double xmax);
  Integrator(FUNC);
  Integrator(GeneralFunction*, double* par, double xmin, double xmax);
  Integrator(GeneralFunction*, double* par);
  virtual ~Integrator();

  virtual double Value() = 0;
  virtual double ValueBetween(double xmin, double xmax);
  virtual operator double() {return Value();}

  virtual Integrator& NStep(int);
  virtual Integrator& DX(double);
  virtual Integrator& ReqPrec(double);   // NOT YET !
  virtual Integrator& Limits(double xmin, double xmax);
  virtual Integrator& Func(FUNC);
  virtual Integrator& Func(GeneralFunction*, double* par);

protected:
  FUNC             mFunc;
  GeneralFunction* mGFF;
  double*          mGFFParm;

  int              mNStep;   // -1 si pas fixe
  double           mDX;
  double           mReqPrec;
  double           mXMin, mXMax;
  
  double           FVal(double x) const;
  virtual void     LimitsChanged() {}
  virtual void     StepsChanged() {}
  virtual void     FuncChanged() {}
};



class TrpzInteg : public Integrator {
public:
  TrpzInteg();
  TrpzInteg(FUNC, double xmin, double xmax);
  TrpzInteg(FUNC);
  TrpzInteg(GeneralFunction*, double* par, double xmin, double xmax);
  TrpzInteg(GeneralFunction*, double* par);

  ~TrpzInteg();

  virtual double Value();
};

class GLInteg : public Integrator {
public:
  GLInteg();
  GLInteg(FUNC, double xmin, double xmax);
  GLInteg(FUNC);
  GLInteg(GeneralFunction*, double* par, double xmin, double xmax);
  GLInteg(GeneralFunction*, double* par);

  ~GLInteg();

  virtual double Value();
  virtual GLInteg& AddBound(double x);
  virtual GLInteg& SetOrder(int order);
protected:

  double* mXPos;     // integration entre 0 et 1
  double* mWeights;  // integration entre 0 et 1
  int     mOrder;    // ordre GL
  set<double> mBounds; // limites des intervalles d'integration
  
  void ComputeWeights();
  void LimitsChanged();
  void StepsChanged();
  void InvalWeights();
  void ComputeBounds();
};

} // Fin du namespace

#endif