CbmBilinearSplineFunction.h

Go to the documentation of this file.

/*
 *====================================================================
 *
 *  Class for 2D spline function used by CbmModels package
 *  
 *  Authors: V.Vovchenko
 *
 *  e-mail : 
 *
 *====================================================================
 *
 *  2D spline class
 *
 *====================================================================
 */
 
#ifndef BILINEARSPLINEFUNCTION_H
#define BILINEARSPLINEFUNCTION_H
#include "CbmSplineFunction.h"
 
class BilinearSplineFunction {
  std::vector<double> xs;
  std::vector<SplineFunction> xspls;
 
public:
  BilinearSplineFunction(void) : xs(), xspls() {
    xs.resize(0);
    xspls.resize(0);
  }
  BilinearSplineFunction(const std::vector<double>& x,
                         const std::vector<double>& y,
                         const std::vector<double>& vals)
    : xs(), xspls() {
    setData(x, y, vals);
  }
  void setData(const std::vector<double>& x,
               const std::vector<double>& y,
               const std::vector<double>& vals) {
    if (x.size() > 0) {
      xs.resize(0);
      xspls.resize(0);
      double cx = -1e50;
      for (unsigned int i = 0; i < x.size(); ++i) {
        if (fabs(x[i] - cx) > 1e-6) {
          xspls.push_back(SplineFunction());
          xs.push_back(x[i]);
          xspls[xspls.size() - 1].add_val(y[i], vals[i]);
          cx = x[i];
        } else {
          xspls[xspls.size() - 1].add_val(y[i], vals[i]);
        }
      }
    }
  }
  double Eval(double x, double y) const {
    if (xs.size() < 2) return -1.;
    unsigned int indx = 0;
    std::vector<double>::const_iterator it =
      lower_bound(xs.begin(), xs.end(), x);
    indx     = distance(xs.begin(), it);
    int ind1 = 0, ind2 = 0;
    if (indx == 0) {
      ind1 = 0;
      ind2 = 1;
    } else if (indx == xs.size()) {
      ind1 = indx - 2;
      ind2 = indx - 1;
    } else {
      ind1 = indx - 1;
      ind2 = indx;
    }
    double f1v = xspls[ind1].f(y);
    double f2v = xspls[ind2].f(y);
    return f1v + (x - xs[ind1]) * (f2v - f1v) / (xs[ind2] - xs[ind1]);
  }
  ~BilinearSplineFunction(void) {}
};
 
#endif