source: Sophya/trunk/Cosmo/SimLSS/integfunc.cc@ 3160

#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include "pexceptions.h"

#include "integfunc.h"

static unsigned short IntegrateFunc_GlOrder = 0;
static vector<double> IntegrateFunc_x;
static vector<double> IntegrateFunc_w;

///////////////////////////////////////////////////////////
//************** Integration of Functions ***************//
///////////////////////////////////////////////////////////


////////////////////////////////////////////////////////////////////////////////////
double IntegrateFunc(GenericFunc& func,double xmin,double xmax
       ,double perc,double dxinc,double dxmax,unsigned short glorder)
// --- Integration adaptative ---
//     Integrate[func[x], {x,xmin,xmax}]
// ..xmin,xmax are the integration limits
// ..dxinc is the searching increment x = xmin+i*dxinc
//         if <0  npt = int(|dxinc|)   (si<2 alors npt=100)
//                et dxinc = (xmax-xmin)/npt
// ..dxmax is the maximum possible increment (if dxmax<=0 no test)
// ---
// Partition de [xmin,xmax] en intervalles [x(i),x(i+1)]:
// on parcourt [xmin,xmax] par pas de "dxinc" : x(i) = xmin + i*dxinc
// on cree un intervalle [x(i),x(i+1)]
//     - si |f[x(i+1)] - f[x(i)]| > perc*|f[[x(i)]|
//     - si |x(i+1)-x(i)| >= dxmax  (si dxmax>0.)
// Dans un intervalle [x(i),x(i+1)] la fonction est integree
//   avec une methode de gauss-legendre d'ordre "glorder"
{
  double signe = 1.;
  if(xmin>xmax) {double tmp=xmax; xmax=xmin; xmin=tmp; signe=-1.;}

  if(glorder==0) glorder = 4;
  if(perc<=0.) perc=0.1;
  if(dxinc<=0.) {int n=int(-dxinc); if(n<2) n=100; dxinc=(xmax-xmin)/n;}
  if(glorder != IntegrateFunc_GlOrder) {
    IntegrateFunc_GlOrder = glorder;
    Compute_GaussLeg(glorder,IntegrateFunc_x,IntegrateFunc_w,0.,1.);
  }

  // Recherche des intervalles [x(i),x(i+1)]
  int_4 ninter = 0;
  double sum = 0., xbas=xmin, fbas=func(xbas), fabsfbas=fabs(fbas);
  for(double x=xmin+dxinc; x<xmax+dxinc/2.; x += dxinc) {
    double f = func(x);
    double dx = x-xbas;
    bool doit = false;
    if( x>xmax ) {doit = true; x=xmax;}
    else if( dxmax>0. && dx>dxmax ) doit = true;
    else if( fabs(f-fbas)>perc*fabsfbas ) doit = true;
    if( ! doit ) continue;
    double s = 0.;
    for(unsigned short i=0;i<IntegrateFunc_GlOrder;i++)
      s += IntegrateFunc_w[i]*func(xbas+IntegrateFunc_x[i]*dx);
    sum += s*dx;
    xbas = x; fbas = f; fabsfbas=fabs(fbas); ninter++;
  }
  //cout<<"Ninter="<<ninter<<endl;

  return sum*signe;
}

////////////////////////////////////////////////////////////////////////////////////
double IntegrateFuncLog(GenericFunc& func,double lxmin,double lxmax
         ,double perc,double dlxinc,double dlxmax,unsigned short  glorder)
// --- Integration adaptative ---
// Idem IntegrateFunc but integrate on logarithmic (base 10) intervals:
//    Int[ f(x) dx ] = Int[ x*f(x) dlog10(x) ] * log(10)
// ..lxmin,lxmax are the log10() of the integration limits
// ..dlxinc is the searching logarithmic (base 10) increment lx = lxmin+i*dlxinc
// ..dlxmax is the maximum possible logarithmic (base 10) increment (if dlxmax<=0 no test)
// Remarque: to be use if "x*f(x) versus log10(x)" looks like a polynomial
//            better than "f(x) versus x"
// ATTENTION: la fonction func que l'on passe en argument
//            est "func(x)" et non pas "func(log10(x))"
{
  double signe = 1.;
  if(lxmin>lxmax) {double tmp=lxmax; lxmax=lxmin; lxmin=tmp; signe=-1.;}

  if(glorder==0) glorder = 4;
  if(perc<=0.) perc=0.1;
  if(dlxinc<=0.) {int n=int(-dlxinc); if(n<2) n=100; dlxinc=(lxmax-lxmin)/n;}
  if(glorder != IntegrateFunc_GlOrder) {
    IntegrateFunc_GlOrder = glorder;
    Compute_GaussLeg(glorder,IntegrateFunc_x,IntegrateFunc_w,0.,1.);
  }

  // Recherche des intervalles [lx(i),lx(i+1)]
  int_4 ninter = 0;
  double sum = 0., lxbas=lxmin, fbas=func(pow(10.,lxbas)), fabsfbas=fabs(fbas);
  for(double lx=lxmin+dlxinc; lx<lxmax+dlxinc/2.; lx += dlxinc) {
    double f = func(pow(10.,lx));
    double dlx = lx-lxbas;
    bool doit = false;
    if( lx>lxmax ) {doit = true; lx=lxmax;}
    else if( dlxmax>0. && dlx>dlxmax ) doit = true;
    else if( fabs(f-fbas)>perc*fabsfbas ) doit = true;
    if( ! doit ) continue;
    double s = 0.;
    for(unsigned short i=0;i<IntegrateFunc_GlOrder;i++) {
      double y = pow(10.,lxbas+IntegrateFunc_x[i]*dlx);
      s += IntegrateFunc_w[i]*y*func(y);
    }
    sum += s*dlx;
    lxbas = lx; fbas = f; fabsfbas=fabs(fbas); ninter++;
  }
  //cout<<"Ninter="<<ninter<<endl;

  return M_LN10*sum*signe;
}

////////////////////////////////////////////////////////////////////////////////////
/*
Integration des fonctions a une dimension y=f(x) par la Methode de Gauss-Legendre.
  --> Calcul des coefficients du developpement pour [x1,x2]
| INPUT:
|  x1,x2 : bornes de l'intervalle (dans nbinteg.h -> x1=-0.5 x2=0.5)
|  glorder = degre n du developpement de Gauss-Legendre
| OUTPUT:
|  x[] = valeur des abscisses ou l'on calcule (dim=n)
|  w[] = valeur des coefficients associes
| REMARQUES:
|  - x et w seront dimensionnes a n.
|  - l'integration est rigoureuse si sur l'intervalle d'integration
|    la fonction f(x) peut etre approximee par un polynome
|    de degre 2*m (monome le + haut x**(2*n-1) )
*/
void Compute_GaussLeg(unsigned short glorder,vector<double>& x,vector<double>& w,double x1,double x2)
{
  if(glorder==0) return;
  int n = (int)glorder;
  x.resize(n,0.); w.resize(n,0.);

   int m,j,i;
   double z1,z,xm,xl,pp,p3,p2,p1;

   m=(n+1)/2;
   xm=0.5*(x2+x1);
   xl=0.5*(x2-x1);
   for (i=1;i<=m;i++)  {
      z=cos(M_PI*(i-0.25)/(n+0.5));
      do {
         p1=1.0;
         p2=0.0;
         for (j=1;j<=n;j++) {
            p3=p2;
            p2=p1;
            p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
         }
         pp=n*(z*p1-p2)/(z*z-1.0);
         z1=z;
         z=z1-p1/pp;
      } while (fabs(z-z1) > 3.0e-11);  // epsilon
      x[i-1]=xm-xl*z;
      x[n-i]=xm+xl*z;
      w[i-1]=2.0*xl/((1.0-z*z)*pp*pp);
      w[n-i]=w[i-1];
   }
}