Index: /trunk/SophyaLib/NTools/functab.cc
===================================================================
--- /trunk/SophyaLib/NTools/functab.cc	(revision 2618)
+++ /trunk/SophyaLib/NTools/functab.cc	(revision 2618)
@@ -0,0 +1,125 @@
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include <iostream>
+#include <math.h>
+#include <string.h>
+#include <string>
+
+#include "pexceptions.h"
+#include "functab.h"
+
+/*!
+  \class SOPHYA::FuncTab
+  \ingroup NTools
+  Function tabulation
+*/
+
+/* --Methode-- */
+//! Constructor
+/*!
+  \param func : function to be tabulated
+  \param ntab : number of tabulated points
+  \param xmin xmax : tabulate from xmin to max
+  \verbatim
+  - Tabul[1] contient func(xmin), Tabul[NTab] contient func(xmax)
+  - ATTENTION: Tabul[i]   contient func(xmin + (i-1)*Dx)
+               Tabul[i+1] contient func(xmin + i*Dx)
+  - Tabul[0] et Tabul[NTab+1] servent pour interpoler rapidement
+    Par defaut Tabul[0] contient Tabul[1] (=func(xmin))
+               Tabul[NTab+1] contient Tabul[NTab] (=func(xmax))
+    Ces valeurs peuvent etre changees avec les methodes SetLimVal
+    Si la fonction existe pour x<xmin et x>xmax alors le mieux est
+       d'utiliser la methode SetLimVal(void):
+            Tabul[0]      = func(xmin-dx)
+            Tabul[NTab+1] = func(xmax+dx)
+    Si la fonction n'existe pas pour x<xmin et x>xmax utiliser
+            la methode SetLimVal(val0,valn):
+            Tabul[0]      = val0
+            Tabul[NTab+1] = valn
+  \endverbatim
+
+ */
+FuncTab::FuncTab(double (*func)(double const),uint_4 ntab,r_8 xmin,r_8 xmax)
+{
+  cout<<"func="<<func<<endl;
+  if(ntab<=1) ntab = 10;
+  if(xmin>=xmax) {xmin=0.; xmax=1.;}
+  NTab = ntab;
+  XMin = xmin;
+  XMax = xmax;
+  Func = func;
+  Dx = (r_8)(XMax-XMin)/(r_8)(NTab-1);
+  Tabul = new r_8[NTab+2];
+  for(int_4 i=0;i<NTab+2;i++) Tabul[i] = 0.;
+  PInterpL = PInterpP[0] = PInterpP[1] = NULL;
+  if(!Tabul) return;
+  for(int_4 i=0;i<NTab;i++) Tabul[i+1] = func(XMin + i*Dx);
+  Tabul[0] = Tabul[1];
+  Tabul[NTab+1] = Tabul[NTab];
+}
+
+/* --Methode-- */
+//! Destructor
+FuncTab::~FuncTab(void)
+{
+  if(Tabul) {delete [] Tabul; Tabul=NULL;}
+  if(PInterpL) {delete [] PInterpL; PInterpL=NULL;}
+}
+
+/* --Methode-- */
+//! Set automatically values for the index [0] and [NTab+1]
+/*! Set automatically values for the index [0] and [NTab+1]
+  That is needed for fast interpolation
+  Value are set for xmin-dx and xmax+dx
+ */
+void FuncTab::SetLimVal(void)
+{
+  Tabul[0] = Func(XMin-Dx);
+  Tabul[NTab+1] = Func(XMax+Dx);
+}
+
+
+/* --Methode-- */
+//! Give values for the index [0] and [NTab+1]
+/*! Set automatically values for the index [0] and [NTab+1]
+  That is needed for fast interpolation
+ */
+void FuncTab::SetLimVal(r_8 val0,r_8 valn)
+{
+  Tabul[0] = val0;
+  Tabul[NTab+1] = valn;
+}
+
+/* --Methode-- */
+//! Construc table for linear and parabolic interpolation
+/*!
+  \param typ : 0=nearest tabulated value, 1=linear, 2=parabolic
+  \warning For speed consideration, no protection is done
+           in the corresponding Value?() inline routine.
+ */
+void FuncTab::SetInterp(uint_2 typ)
+{
+  if(typ==1) {  // Linear Interpolation between two consecutive values
+
+    if(PInterpL) delete [] PInterpL;
+    PInterpL = new r_8[NTab+2];
+    for(int_4 i=0;i<NTab+1;i++)
+      PInterpL[i] = (Tabul[i+1]-Tabul[i])/Dx;
+    PInterpL[NTab+1] = 0.;
+
+  } else if(typ==2) {
+
+    if(PInterpP[0]) delete [] PInterpP[0];
+    if(PInterpP[1]) delete [] PInterpP[1];
+    PInterpP[0] = new r_8[NTab+2];
+    PInterpP[1] = new r_8[NTab+2];
+    for(int_4 i=1;i<=NTab;i++) {
+      PInterpP[0][i] = (Tabul[i+1]+Tabul[i-1]-2.*Tabul[i])/(2.*Dx*Dx);
+      PInterpP[1][i] = (Tabul[i+1]-Tabul[i-1])/(2.*Dx);
+    }
+    PInterpP[0][0]=PInterpP[0][NTab+1]=PInterpP[1][0]=PInterpP[1][NTab+1]=0.;
+
+  }
+}
Index: /trunk/SophyaLib/NTools/functab.h
===================================================================
--- /trunk/SophyaLib/NTools/functab.h	(revision 2618)
+++ /trunk/SophyaLib/NTools/functab.h	(revision 2618)
@@ -0,0 +1,56 @@
+#ifndef FUNCTAB_SEEN
+#define VFUNCTAB_SEEN
+
+#include "machdefs.h"
+#include "tvector.h"
+#include "tarray.h"
+
+namespace SOPHYA {
+
+class FuncTab {
+public:
+  FuncTab(double (*func)(double const),uint_4 nbin,r_8 xmin,r_8 xmax);
+  virtual ~FuncTab(void);
+
+  //! Give values for the index [0] and [NTab+1]
+  void SetLimVal(void);
+  void SetLimVal(r_8 val0,r_8 valn);
+
+  //! Interpolation type
+  void SetInterp(uint_2 typ=0);
+
+  //! L'index de la valeur tabulee la plus proche de x (entre 1 et NTab)
+  inline uint_4 IndexNear(r_8 x)
+    {uint_4 i=(uint_4)((x-XMin)/Dx+0.5) +1; return (i>NTab)? NTab: i;}
+  //! L'index de la valeur tabulee juste inferieure a x (entre 1 et NTab)
+  inline uint_4 IndexLow(r_8 x)
+    {uint_4 i=(uint_4)((x-XMin)/Dx) +1; return (i>NTab)? NTab: i;}
+
+  //! L'abscisse pour la valeur i de l'index
+  inline r_8 X(uint_4 i) {return XMin+(i-1)*Dx;}
+
+  //! Le valeur tabulee la plus proche de x
+  inline r_8 Value(r_8 x) {return Tabul[IndexNear(x)];}
+
+  //! L'interpolation lineaire dans le bin (attention no-protection for time saving)
+  inline r_8 ValueL(r_8 x)
+    {uint_4 i=IndexLow(x);
+     return Tabul[i]+PInterpL[i]*(x-X(i));}
+
+  //! L'interpolation parabolique (attention no-protection for time saving)
+  inline r_8 ValueP(r_8 x)
+  {uint_4 i=IndexLow(x); x -= X(i);
+     return Tabul[i]+PInterpP[0][i]*x*x+PInterpP[1][i]*x;}
+
+protected:
+  double (*Func)(double const);
+  uint_4 NTab;
+  r_8 XMin,XMax,Dx;
+  r_8 *Tabul;
+  r_8 *PInterpL;
+  r_8 *PInterpP[2];
+};
+
+} // Fin du namespace
+
+#endif
