// Copyright (c) 1997-1999 Matra Datavision
// Copyright (c) 1999-2012 OPEN CASCADE SAS
//
// The content of this file is subject to the Open CASCADE Technology Public
// License Version 6.5 (the "License"). You may not use the content of this file
// except in compliance with the License. Please obtain a copy of the License
// at http://www.opencascade.org and read it completely before using this file.
//
// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
//
// The Original Code and all software distributed under the License is
// distributed on an "AS IS" basis, without warranty of any kind, and the
// Initial Developer hereby disclaims all such warranties, including without
// limitation, any warranties of merchantability, fitness for a particular
// purpose or non-infringement. Please see the License for the specific terms
// and conditions governing the rights and limitations under the License.


/*
Par Gauss le calcul d'une integrale simple se transforme en sommation des 
valeurs de la fonction donnee aux <Order> points de Gauss affectee des poids
de Gauss.

Les points et poids de Gauss sont stockes dans GaussPoints.cxx.
Les points sont compris entre les valeurs -1 et +1, ce qui necessite un 
changement de variable pour les faire varier dans l'intervalle [Lower, Upper].


On veut calculer Integrale( f(u)* du) entre a et b.


Etapes du calcul:

1- calcul de la fonction au ieme point de Gauss (apres changement de variable).

2- multiplication de cette valeur par le ieme poids de Gauss.

3- sommation de toutes ces valeurs.

4- retour a l'intervalle [Lower, Upper] de notre integrale.

*/

//#ifndef DEB
#define No_Standard_RangeError
#define No_Standard_OutOfRange
#define No_Standard_DimensionError
//#endif

#include <math_GaussSingleIntegration.ixx>

#include <math.hxx>
#include <math_Vector.hxx>
#include <math_Function.hxx>

math_GaussSingleIntegration::math_GaussSingleIntegration() : Done(Standard_False)  
{
}

math_GaussSingleIntegration::
                 math_GaussSingleIntegration(math_Function& F,
					     const Standard_Real Lower,
					     const Standard_Real Upper,
					     const Standard_Integer Order)
{
  Standard_Integer theOrder = Min(math::GaussPointsMax(), Order);
  Perform(F, Lower, Upper, theOrder);
}

math_GaussSingleIntegration::
                 math_GaussSingleIntegration(math_Function& F,
					     const Standard_Real Lower,
					     const Standard_Real Upper,
					     const Standard_Integer Order,
					     const Standard_Real Tol)
{
  Standard_Integer theOrder = Min(math::GaussPointsMax(), Order);

  const Standard_Integer IterMax = 13;   // Max number of iteration
  Standard_Integer NIter = 1;            // current number of iteration
  Standard_Integer NbInterval = 1;       // current number of subintervals
  Standard_Real    dU,OldLen,Len;

  Perform(F, Lower, Upper, theOrder);
  Len = Val;
  do {
    OldLen = Len;
    Len = 0.;
    NbInterval *= 2;
    dU = (Upper-Lower)/NbInterval;    
    for (Standard_Integer i=1; i<=NbInterval; i++) {
      Perform(F, Lower+(i-1)*dU, Lower+i*dU, theOrder);
      if (!Done) return;
      Len += Val;
    }
    NIter++;
  }    
  while (fabs(OldLen-Len) > Tol && NIter <= IterMax);

  Val = Len;
}

void math_GaussSingleIntegration::Perform(math_Function& F,
					  const Standard_Real Lower,
					  const Standard_Real Upper,
					  const Standard_Integer Order)
{
  Standard_Real xr, xm, dx;
  Standard_Integer j;
  Standard_Real F1, F2;
  Standard_Boolean Ok1;
  math_Vector GaussP(1, Order);
  math_Vector GaussW(1, Order);
  Done = Standard_False;

//Recuperation des points de Gauss dans le fichier GaussPoints.
  math::GaussPoints(Order,GaussP);
  math::GaussWeights(Order,GaussW);

// Calcul de l'integrale aux points de Gauss :

// Changement de variable pour la mise a l'echelle [Lower, Upper] :
  xm = 0.5*(Upper + Lower);
  xr = 0.5*(Upper - Lower);
  Val = 0.;

  Standard_Integer ind = Order/2, ind1 = (Order+1)/2;
  if(ind1 > ind) { // odder case
    Ok1 = F.Value(xm, Val);
    if (!Ok1) return;
    Val *= GaussW(ind1);
  }
// Sommation sur tous les points de Gauss: avec utilisation de la symetrie.
  for (j = 1; j <= ind; j++) {
    dx = xr*GaussP(j);
    Ok1 = F.Value(xm-dx, F1);
    if(!Ok1) return;
    Ok1 = F.Value(xm+dx, F2);
    if(!Ok1) return;
    // Multiplication par les poids de Gauss.
    Standard_Real FT = F1+F2;
    Val += GaussW(j)*FT;  
  }
  // Mise a l'echelle de l'intervalle [Lower, Upper]
  Val *= xr;
  Done = Standard_True;
}

void math_GaussSingleIntegration::Dump(Standard_OStream& o) const {

  o <<"math_GaussSingleIntegration ";
   if (Done) {
     o << " Status = Done \n";
     o << "Integration Value = " << Val<<"\n";
   }
   else {
     o << "Status = not Done \n";
   }
}