text: /*
   * PiecewiseLinearInterpolation.cpp
   *
   *  Created on: Sep 7, 2010
   *      Author: TF
   * \copyright
   * Copyright (c) 2020, OpenGeoSys Community (http://www.opengeosys.org)
   *            Distributed under a Modified BSD License.
   *              See accompanying file LICENSE.txt or
   *              http://www.opengeosys.org/project/license
   */
  
  #include "PiecewiseLinearInterpolation.h"
  #include "binarySearch.h"
  #include "swap.h"
  
  #include <iostream>
  
  namespace MathLib
  {
  PiecewiseLinearInterpolation::PiecewiseLinearInterpolation(
      const std::vector<double>& supporting_points,
      const std::vector<double>& values_at_supp_pnts)
      : _supporting_points(supporting_points),
        _values_at_supp_pnts(values_at_supp_pnts)
  {
  }
  
  PiecewiseLinearInterpolation::PiecewiseLinearInterpolation(
      const std::vector<double>& supporting_points,
      const std::vector<double>& values_at_supp_pnts,
      const std::vector<double>& points_to_interpolate,
      std::vector<double>& values_at_interpol_pnts)
      : _supporting_points(supporting_points),
        _values_at_supp_pnts(values_at_supp_pnts)
  {
      //	std::cout << "PiecewiseLinearInterpolation::PiecewiseLinearInterpolation
      // support_points, values_at_supp_pnts: "
      //<< "\n";
      //	for (size_t k(0); k<supporting_points.size(); k++) {
      //		std::cout << supporting_points[k] << " " << values_at_supp_pnts[k]
      //<<
      //"\n";
      //	}
      //	std::cout << "\n";
      values_at_interpol_pnts.clear();
      for (size_t k(0); k < points_to_interpolate.size(); k++)
          values_at_interpol_pnts.push_back(
              this->getValue(points_to_interpolate[k]));
  }
  
  PiecewiseLinearInterpolation::~PiecewiseLinearInterpolation() {}
  
  double PiecewiseLinearInterpolation::getValue(double pnt_to_interpolate)
  {
      // search interval that has the point inside
      size_t interval_idx(std::numeric_limits<size_t>::max());
      for (size_t k(1); k < _supporting_points.size() &&
                        interval_idx == std::numeric_limits<size_t>::max();
           k++)
          if (_supporting_points[k - 1] <= pnt_to_interpolate &&
              pnt_to_interpolate <= _supporting_points[k])
              interval_idx = k - 1;
  
      // compute linear interpolation polynom: y = m * x + n
      long double m((_values_at_supp_pnts[interval_idx + 1] -
                     _values_at_supp_pnts[interval_idx]) /
                    (_supporting_points[interval_idx + 1] -
                     _supporting_points[interval_idx]));
      //	double m ((_values_at_supp_pnts[interval_idx] -
      //_values_at_supp_pnts[interval_idx+1]) /
      //(_supporting_points[interval_idx] - _supporting_points[interval_idx+1]));
      //	double n (_values_at_supp_pnts[interval_idx+1] - m *
      //_supporting_points[interval_idx+1]);
      long double n(_values_at_supp_pnts[interval_idx] -
                    m * _supporting_points[interval_idx]);
  
      return m * pnt_to_interpolate + n;
  }
  }  // namespace MathLib
  