source: Sophya/trunk/Cosmo/SimLSS/geneutils.h@ 3248

#ifndef GENEUTILS_SEEN
#define GENEUTILS_SEEN

#include "machdefs.h"
#include <math.h>
#include "genericfunc.h"
#include "histos.h"
#include "tvector.h"
#include "cspline.h"

#include <vector>

namespace SOPHYA {

//----------------------------------------------------
class InterpFunc {
public:
  InterpFunc(double xmin,double xmax,vector<double>& y);
  virtual ~InterpFunc(void) { }

  double XMin(void) {return _xmin;}
  double XMax(void) {return _xmax;}
  inline double X(long i) {return _xmin + i*_dx;}

  //! Retourne l'element le plus proche de f donnant y=f(x)
  inline double operator()(double x)
         {
         x -= _xmin;
         long i = long(x/_dx+0.5);  // On prend le "i" le plus proche
         if(i<0) i=0; else if(i>=_nm1) i=_nm1-1;
         return _y[i];
         }

  // idem operator(double) et retourne
  // ok==0 si valeur trouvee, 1 si x<xmin, 2 si x>xmax
  inline double operator()(double x,unsigned short& ok)
    {ok=0; if(x<_xmin) ok=1; else if(x>_xmax) ok=2; return (*this)(x);}

  //! Retourne l'interpolation lineaire de f donnant y=f(x)
  // ok==0 si valeur trouvee, 1 si x<xmin, 2 si x>xmax
  double Linear(double x,unsigned short& ok);

  //! Retourne l'interpolation parabolique de f donnant y=f(x)
  // ok==0 si valeur trouvee, 1 si x<xmin, 2 si x>xmax
  double Parab(double x,unsigned short& ok);

protected:
  double _xmin,_xmax,_dx;
  long _nm1;  // n-1
  vector<double>& _y;
};

//----------------------------------------------------
class InverseFunc {
public:
  InverseFunc(vector<double>& x,vector<double>& y);
  virtual ~InverseFunc(void);
  int ComputeLinear(long n,vector<double>& xfcty);
  int ComputeParab(long n,vector<double>& xfcty);
  double YMin(void) {return _ymin;}
  double YMax(void) {return _ymax;}
protected:
  inline void find_in_y(double x,long& klo,long& khi)
    {
      long k;
      klo=0, khi=_y.size()-1;
      while (khi-klo > 1) {
        k = (khi+klo) >> 1;
        if (_y[k] > x) khi=k; else klo=k;
      }
    }

  vector<double>& _x;
  vector<double>& _y;
  double _ymin,_ymax,_dy;
};

}

//----------------------------------------------------
double InterpTab(double x0,vector<double>& X,vector<double>& Y,unsigned short typint=0);

//----------------------------------------------------
int FuncToHisto(GenericFunc& func,Histo& h,bool logaxex=false);
int FuncToVec(GenericFunc& func,TVector<r_8>& h,double xmin,double xmax,bool logaxex=false);

//----------------------------------------------------
double AngSol(double dtheta,double dphi,double theta0=M_PI/2.);
double AngSol(double dtheta);

//----------------------------------------------------
unsigned long PoissRandLimit(double mu,double mumax=10.);

//----------------------------------------------------
double IntegrateFunc(GenericFunc& func,double xmin,double xmax
         ,double perc=0.1,double dxinc=-1.,double dxmax=-1.,unsigned short glorder=4);

double IntegrateFuncLog(GenericFunc& func,double lxmin,double lxmax
         ,double perc=0.1,double dlxinc=-1.,double dlxmax=-1.,unsigned short glorder=4);

void Compute_GaussLeg(unsigned short glorder,vector<double>& x,vector<double>& w,double x1=0.,double x2=1.);

#endif