Changeset 3196 in Sophya for trunk/Cosmo/SimLSS/geneutils.cc

Timestamp:

Apr 3, 2007, 4:06:50 PM (18 years ago)

Author:

cmv

Message:

les AGN selon C.Jackson, une premiere approche simplifiee, recodage from Jim Rich. cmv 03/04/2007

File:

• trunk/Cosmo/SimLSS/geneutils.cc (modified) (2 diffs)

trunk/Cosmo/SimLSS/geneutils.cc

@@ -150 +150 @@
 }
 
+//----------------------------------------------------
+double InterpTab(double x0,vector<double>& X,vector<double>& Y,unsigned short typint)
+// Interpole in x0 the table Y = f(X)
+//           X doit etre ordonne par ordre croissant (strictement)
+// typint = 0 : nearest value
+//          1 : linear interpolation
+//          2 : parabolique interpolation
+{
+ long n = X.size();
+ if(Y.size()<n || n<2)
+   throw ParmError("InterpTab_Error :  incompatible size between X and Y tables!");
+
+ if(x0<X[0] || x0>X[n-1]) return 0.;
+ if(typint>2) typint = 0;
+
+ // Recherche des indices encadrants par dichotomie
+ long k, klo=0, khi=n-1;
+ while (khi-klo > 1) {
+   k = (khi+klo) >> 1;
+   if (X[k] > x0) khi=k; else klo=k;
+ }
+
+ // Quel est le plus proche?
+ k = (x0-X[klo]<X[khi]-x0) ? klo: khi;
+
+ // On retourne le plus proche
+ if(typint==0) return Y[k];
+
+ // On retourne l'extrapolation lineaire
+ if(typint==1 || n<3)
+   return Y[klo] + (Y[khi]-Y[klo])/(X[khi]-X[klo])*(x0-X[klo]);
+
+ // On retourne l'extrapolation parabolique
+ if(k==0) k++; else if(k==n-1) k--;
+ klo = k-1; khi = k+1;
+ double x1 = X[klo]-X[k], x2 = X[khi]-X[k];
+ double y1 = Y[klo]-Y[k], y2 = Y[khi]-Y[k];
+ double den = x1*x2*(x1-x2);
+ double a = (y1*x2-y2*x1)/den;
+ double b = (y2*x1*x1-y1*x2*x2)/den;
+ x0 -= X[k];
+ return Y[k] + (a*x0+b)*x0;;
+
+}
+
 //-------------------------------------------------------------------
 int FuncToHisto(GenericFunc& func,Histo& h,bool logaxex)
@@ -252 +297 @@
   return n;
 }
+
+//-------------------------------------------------------------------
+
+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];
+   }
+}