source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/tools/src/Interpolate.cc @ 114

//
// $Id: Interpolate.cc 2804 2008-10-09 12:10:06Z biktem $
//
#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;
}