source: Sophya/trunk/SophyaLib/NTools/slininterp.h@ 3951

#ifndef SLININTERP_H_SEEN
#define SLININTERP_H_SEEN
/*    
   --- SOPHYA software - NTools module ---
   --- Class SLinInterp1D : Simple linear 1D interpolation class
   R. Ansari , C. Magneville   2000-2010 
   (C) UPS+LAL IN2P3/CNRS     (C) IRFU-SPP/CEA 
*/

#include "machdefs.h"
#include <iostream>
#include <vector>
#include <string>
#include <exception>
#include <stdexcept>

#include "genericfunc.h"

using namespace std;

//-------------------------------------------
// --- Class SLinInterp1D : 
//  Simple linear 1D interpolation class 
//-------------------------------------------

namespace SOPHYA {

class SLinInterp1D : public GenericFunc {
public :
  //! Default constructor - represent the function y=x
  SLinInterp1D(); 
  // Regularly spaced points in X with Y values defined by yreg         
  SLinInterp1D(double xmin, double xmax, vector<double>& yreg);
  //  Interpolate to a finer regularly spaced grid, from xmin to xmax with npt points if (npt>0)
  SLinInterp1D(vector<double>& xs, vector<double>& ys, double xmin=1., double xmax=-1., size_t npt=0); 

  virtual ~SLinInterp1D() { }
        
  double XMin() const { return xmin_; }
  double XMax() const { return xmax_; }
  double DeltaX()  { return dx_; }
  inline double X(long i) const {return xmin_ + i*dx_;}  // returns x value for index i

  // --------------------------------------------------------------
  //! Return the interpolated Y value as a function of X 
  double YInterp(double x) const ;
  // To compile without SOPHYA : Comment the following line 
  //! Return the interpolated Y value as a function of X 
  virtual inline double operator()(double x) {  return YInterp(x); }
  // --------------------------------------------------------------
        
  // Define the interpolation points through a set of regularly spaced points on X, between xmin and xmax 
  void DefinePoints(double xmin, double xmax, vector<double>& yreg);
  // Interpolate to a finer regularly spaced grid, from xmin to xmax with npt points 
  void DefinePoints(vector<double>& xs, vector<double>& ys, double xmin=1., double xmax=-1., size_t npt=0); 

  // Read  Y's (one/line) for regularly spaced X's from file and call DefinePoints(xmin, xmax, yreg)
  size_t ReadYFromFile(string const& filename, double xmin, double xmax, char clm='#');
  // Read pairs of X Y (one pair/line) from file and call DefinePoints(xs, ys ...)
  size_t ReadXYFromFile(string const& filename, double xmin=1., double xmax=-1., size_t npt=0, char clm='#');
  
  vector<double>& GetVX()  { return xs_; }
  vector<double>& GetVY()  { return ys_; }

  //! Print the object (interpolation points) on "ostream". lev=0,1
  ostream& Print(ostream& os,int lev=0) const ;
  //! Print the object (interpolation points) on "cout". lev=0,1
  inline ostream& Print(int lev=0) const { return Print(cout, lev); } 

protected:
  vector<double> yreg_, xs_, ys_;  // interpolated y value for regularly spaced x 
  double xmin_, xmax_, dx_;        // dx is spacing of finer grid of x's
  size_t ksmx_;                    // Maximum index value in xs_, ys_
  size_t npoints_;                 // Number of regularly spaced points, xmax not included 
};

/*! operator << overloading - Prints the interpolator object on the stream \b s */
inline ostream& operator << (ostream& s, SLinInterp1D const& a) 
{ a.Print(s,0);  return s; }

} // namespace SOPHYA

#endif   /* SLININTERP_H__SEEN */