#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.;

  }
}