source: Sophya/trunk/Poubelle/DPC:FitsIOServer/NTools/integ.cc@ 1797

#include "integ.h"
#include "generalfit.h"

// A faire :
//   1 dans GLInteg, optimiser recalcul en utilisant des poids
//     calcules sur [0,1] et chgt de variable
//   2 dans GLInteg, distinction nStep = subdivisions de l'intervalle,
//     et on applique GL d'ordre "ordre" sur chaque subdivision, en
//     utilisant le (1)

//++
// Class        Integrator
// Lib          Outils++ 
// include      integ.h
//
//      Classe abstraite d'intégration numérique 1D.
//      On fournit une fonction double f(double) au constructeur, ou
//      une GeneralFunction avec des paramètres définis.
//      L'objet Integrator est convertible en valeur double qui est la valeur
//      de l'intégrale. Diverses méthodes permettent de choisir des options
//      de calcul, et ces méthodes retournent une référence sur l'objet, pour
//      permettre une notation chaînée.
//--
//++ 
// Links        Implementations
// TrpzInteg
// GLInteg
//--

//++
// Titre        Constructeurs
//--

//++
Integrator::Integrator()
//
//      Constructeur par défaut. L'objet n'est pas utilisable en l'état.
//--

: mFunc(NULL), mGFF(NULL), mGFFParm(NULL),
  mNStep(50), mDX(-1), mReqPrec(-1),
  mXMin(0), mXMax(1)
{}

//++
Integrator::Integrator(FUNC f, double xmin, double xmax)
//
//      Constructeur à partir de la fonction double->double, et des
//      bornes d'intégration.
//--
: mFunc(f), mGFF(NULL), mGFFParm(NULL),
  mNStep(50), mDX(-1), mReqPrec(-1),
  mXMin(xmin), mXMax(xmax)
{}

//++
Integrator::Integrator(FUNC f)
// 
//      Constructeur sans spécifier les bornes. Elles sont positionnées
//      à [0,1], et on pourra les modifier plus tard.
//--
: mFunc(f), mGFF(NULL), mGFFParm(NULL),
  mNStep(50), mDX(-1), mReqPrec(-1),
  mXMin(0), mXMax(1)
{}

//++
Integrator::Integrator(GeneralFunction* gff, double* par, double xmin, double xmax)
//
//      Constructeur à partir d'une GeneralFunction. La fonction doit être une
//      fonction de une variable, et on fournit les valeurs des paramètres ainsi
//      que les bornes d'intégration.
//--
: mFunc(NULL), mGFF(gff), mGFFParm(par),
  mNStep(50), mDX(-1), mReqPrec(-1),
  mXMin(xmin), mXMax(xmax)
{
  DBASSERT(gff->NVar() == 1);
}

//++
Integrator::Integrator(GeneralFunction* gff, double* par)
//
//      Constructeur à partir d'une GeneralFunction. La fonction doit être une
//      fonction de une variable, et on fournit les valeurs des paramètres. 
//      On ne spécifie pas les bornes. Elles sont positionnées
//      à [0,1], et on pourra les modifier plus tard.
//--
: mFunc(NULL), mGFF(gff), mGFFParm(par),
  mNStep(50), mDX(-1), mReqPrec(-1),
  mXMin(0), mXMax(1)
{
  DBASSERT(gff->NVar() == 1);
}

Integrator::~Integrator()
{}

//++
// Titre        Méthodes
//--

//++
Integrator& 
Integrator::NStep(int n)
//
//      Spécifie le nombre de pas pour l'intégration numérique.
//      La signification peut dépendre de la méthode d'intégration.
//--
{
  mNStep = n;
  mDX = mReqPrec = -1;
  StepsChanged();
  return *this;
}

//++
Integrator& 
Integrator::DX(double d)
//
//      Spécifie le pas d'intégration.
//      La signification peut dépendre de la méthode d'intégration.
//--
{
  mDX = d;
  mNStep = mReqPrec = -1;
  StepsChanged();
  return *this;
}

//++
Integrator& 
Integrator::ReqPrec(double p)
//
//      Spécifie la précision souhaitée. 
//      La signification peut dépendre de la méthode d'intégration.
//      Non disponible dans toutes les méthodes d'intégration.
//--
{
  DBASSERT( !"Pas encore implemente !");
  mReqPrec = p;
  mDX = mNStep = -1;
  StepsChanged();
  return *this;
}

//++
Integrator& 
Integrator::Limits(double xmin, double xmax)
//
//      Spécifie les bornes de l'intégrale.
//--
{
  mXMin = xmin;
  mXMax = xmax;
  LimitsChanged();
  return *this;
}

//++
Integrator& 
Integrator::Func(FUNC f)
//
//      Spécifie la fonction à intégrer, sous forme double f(double).
//--
{
  mFunc = f;
  mGFFParm = NULL;
  mFunc = NULL;
  FuncChanged();
  return *this;
}

//++
Integrator& 
Integrator::Func(GeneralFunction* gff, double* par)
//
//      Spécifie la fonction à intégrer, sous forme de GeneralFunction
//      à une variable, et les paramètres sont fournis.
//--
{
  mGFF = gff;
  mGFFParm = par;
  mFunc = NULL;
  DBASSERT(gff->NVar() == 1);
  FuncChanged();
  return *this;
}

double
Integrator::FVal(double x) const
{
  DBASSERT( mFunc || (mGFF && mGFFParm) );
  return mFunc ? (*mFunc)(x) : mGFF->Value(&x, mGFFParm); 
}

double
Integrator::ValueBetween(double xmin, double xmax)
{
  Limits(xmin,xmax);
  return Value();
}

//++
// Class        TrpzInteg
// Lib          Outils++ 
// include      integ.h
//
//      Classe d'intégration par la méthode des trapèzes.
//      Voir Integrator pour les méthodes. Le nombre de pas
//      est le nombre de trapèze, le pas d'intégration est
//      la largeur des trapèzez. Impossible de demander une
//      précision.
//
//--
//++
// Links        Parents
// Integrator
//--

//++
// Titre        Constructeurs
//      Voir Integrator pour les détails.
//--

//++
TrpzInteg::TrpzInteg()
//
//--
{}

//++
TrpzInteg::TrpzInteg(FUNC f, double xmin, double xmax)
//
//--
: Integrator(f, xmin, xmax)
{}

//++
TrpzInteg::TrpzInteg(FUNC f)
//
//--
: Integrator(f)
{}

//++
TrpzInteg::TrpzInteg(GeneralFunction* gff, double* par, double xmin, double xmax)
//
//--
: Integrator(gff, par, xmin, xmax)
{}

//++
TrpzInteg::TrpzInteg(GeneralFunction* gff, double* par)
//
//--
: Integrator(gff, par)
{}

TrpzInteg::~TrpzInteg()
{}

double
TrpzInteg::Value()
{
  double dx = mDX;
  double nstep = mNStep;

  if (dx > 0)
    nstep = (mXMax - mXMin)/dx;

  if (nstep <= 0)
    nstep = 10;

  dx = (mXMax - mXMin) / nstep;

  double s = (FVal(mXMin) + FVal(mXMax))/2;
  double x = dx + mXMin;
  for (int i=1; i<nstep; i++, x += dx)
    s += FVal(x);
    
  return s * dx;
}



//++
// Class        GLInteg
// Lib          Outils++ 
// include      integ.h
//
//      Classe d'intégration par la méthode de Gauss-Legendre.
//      Voir Integrator pour les méthodes. 
//      Pour le moment, nstep est l'ordre de la méthode.
//      Il est prévu un jour de spécifier l'ordre, et que NStep
//      découpe en intervalles sur chacun desquels on applique GL.
//      Le principe de la méthode est de calculer les valeurs de la 
//      fonction aux zéros des polynomes de Legendre. Avec les poids
//      qui vont bien, GL d'ordre n est exacte pour des polynomes de 
//      degré <= 2n+1.
//      Impossible de demander une précision donnée.
//
//--
//++
// Links        Parents
// Integrator
//--

//++
// Titre        Constructeurs
//      Voir Integrator pour les détails.
//--




//++
GLInteg::GLInteg()
//
//--
: mXPos(NULL), mWeights(NULL)
{}


//++
GLInteg::GLInteg(FUNC f, double xmin, double xmax)
//
//--
: Integrator(f, xmin, xmax), mXPos(NULL), mWeights(NULL), mOrder(8)
{
  NStep(1);
}

//++
GLInteg::GLInteg(FUNC f)
//
//--
: Integrator(f), mXPos(NULL), mWeights(NULL), mOrder(8)
{
  NStep(1);
}

//++
GLInteg::GLInteg(GeneralFunction* gff, double* par, double xmin, double xmax)
//
//--
: Integrator(gff, par, xmin, xmax), mXPos(NULL), mWeights(NULL), mOrder(8)
{
  NStep(1);
}

//++
GLInteg::GLInteg(GeneralFunction* gff, double* par)
//
//--
: Integrator(gff, par), mXPos(NULL), mWeights(NULL), mOrder(8)
{
  NStep(1);
}

GLInteg::~GLInteg()
{
  delete[] mXPos;
  delete[] mWeights;
  mXPos = mWeights = NULL;
}

void
GLInteg::InvalWeights()
{
  delete[] mXPos;
  delete[] mWeights;
  mXPos = mWeights = NULL;
}


void
GLInteg::LimitsChanged()
{
  ComputeBounds();
}

void
GLInteg::StepsChanged()
{
  ComputeBounds();
}

GLInteg&
GLInteg::AddBound(double x)
{
  if (x<=mXMin || x>=mXMax) THROW(rangeCheckErr); 
  // On introduira les classes d'exections apres reflexion et de maniere systematique (Rz+cmv)
  // if (x<=mXMin || x>=mXMax) throw range_error("GLInteg::AddBound  bound outside interval");
  mBounds.insert(x);
  return *this;
}


GLInteg&
GLInteg::SetOrder(int order)
{
  mOrder = order;
  InvalWeights();
  return *this;
}



double
GLInteg::Value()
{
  if (!mXPos) ComputeWeights();
  double s = 0;
  set<double>::iterator i=mBounds.begin();
  set<double>::iterator j=i; j++;
  while (j != mBounds.end()) {
    double s1 = 0;
    double x1 = *i;
    double x2 = *j;
    for (int k=0; k<mOrder; k++)
       s1 += mWeights[k] * FVal(x1 + (x2-x1)*mXPos[k]);
    s += s1*(x2-x1);
    i++; j++;
  }
  return s;
}

void
GLInteg::ComputeBounds()
{
  mBounds.erase(mBounds.begin(), mBounds.end());
  for (int i=0; i<=mNStep; i++) 
    mBounds.insert(mXMin + (mXMax-mXMin)*i/mNStep);
}

void
GLInteg::ComputeWeights()
{
   const double EPS_gauleg = 3.0e-11;

   int m=(mOrder+1)/2;
   const double xxMin = 0;
   const double xxMax = 1;
   
   mXPos    = new double[mOrder];
   mWeights = new double[mOrder];
   
   double xm=0.5*(xxMax+xxMin);
   double xl=0.5*(xxMax-xxMin);
   for (int i=1;i<=m;i++)  {
      double z=cos(3.141592654*(i-0.25)/(mOrder+0.5));
      double z1, pp;
      do {
         double p1=1.0;
         double p2=0.0;
         for (int j=1;j<=mOrder;j++) {
            double p3=p2;
            p2=p1;
            p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
         }
         pp=mOrder*(z*p1-p2)/(z*z-1.0);
         z1=z;
         z=z1-p1/pp;
      } while (fabs(z-z1) > EPS_gauleg);
      mXPos[i-1]         = xm-xl*z;
      mXPos[mOrder-i]    = xm+xl*z;
      mWeights[i-1]      = 2.0*xl/((1.0-z*z)*pp*pp);
      mWeights[mOrder-i] = mWeights[i-1];
   }

}