//
// $Id: Interpolate.cc,v 1.3 2005/04/19 22:46:55 thea Exp $
//
#include <vector>
#include <functional>
#include <algorithm>
#include "Interpolate.hh"
ClassImp(Interpolate)
Bool_t row_less( vector<Double_t> x, vector<Double_t> y) {
return (x[0] < y[0]);
}
//______________________________________________________________________________
Interpolate::Interpolate(const string &fn, size_t nval): fFilename(fn) {
//
// Contructor
//
NumbersFileParser parser(fFilename, nval+1);
size_t nRows=parser.GetCol(0).size();
for(size_t i=0; i<nRows; i++) {
vector<Double_t> row = parser.GetRow(i);
fValues.push_back(row);
}
sort(fValues.begin(), fValues.end(), row_less);
fXmin = fValues[0][0];
fXmax = fValues[nRows-1][0];
for (size_t i(0); i<nval+1; i++)
fUnits.push_back(1.);
}
//______________________________________________________________________________
Interpolate::~Interpolate() {
//
// Destructor
//
}
//______________________________________________________________________________
Double_t Interpolate::GetValue(Double_t x, size_t col) {
//
// Return interpolation in x
//
// if x is out of range return last value
if(x > fXmax) {
MsgForm(EsafMsg::Warning,"Parameter out of range [%.3f,%.f] (file=%s , x = %f)",
fXmin, fXmax, fFilename.c_str(), x);
return fValues[fValues.size()-1][col+1];
} else if(x < fXmin) {
MsgForm(EsafMsg::Warning,"Parameter out of range [%.3f,%.3f] (file=%s , x = %f)",
fXmin, fXmax, fFilename.c_str(), x);
return fValues[0][col+1];
}
// x in range
size_t i;
for(i=1; i<fValues.size(); i++)
if(fValues[i][0]>=x) break;
i--;
Double_t x1=fValues[i][0];
Double_t y1=fValues[i][col+1];
Double_t x2=fValues[i+1][0];
Double_t y2=fValues[i+1][col+1];
// linear interpolation
Double_t value= (y2-y1)/(x2-x1)*(x-x1)+y1;
return value;
}
//______________________________________________________________________________
void Interpolate::SetXUnit(Double_t unit) {
//
// Set the unit of measure for x variable
//
Double_t ratio = unit/fUnits[0];
for (size_t r(0); r < fValues.size(); r++)
fValues[r][0] = fValues[r][0]*ratio;
fUnits[0] = unit;
fXmin = fValues[0][0];
fXmax = fValues[fValues.size()-1][0];
}
//______________________________________________________________________________
void Interpolate::SetUnit(Double_t unit, size_t nval) {
//
// Set the unit of measure for n-th variable
//
if ( nval >= fValues[0].size() )
Msg(EsafMsg::Panic) << "SetUnit: variable number out of range "<< MsgDispatch;
// data index goes from 1 to fValues[0].size()
nval++;
Double_t ratio = unit/fUnits[nval];
for (size_t r(0); r < fValues.size(); r++)
fValues[r][nval] = fValues[r][nval]*ratio;
fUnits[nval] = unit;
}