linearInterpolation.hpp

/*
  Copyright 2009, 2010, 2013 SINTEF ICT, Applied Mathematics.
  Copyright 2009, 2010 Statoil ASA.

  This file is part of the Open Porous Media project (OPM).

  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef OPM_LINEARINTERPOLATION_HEADER_INCLUDED
#define OPM_LINEARINTERPOLATION_HEADER_INCLUDED


#include <algorithm>
#include <tuple>
#include <vector>

namespace Opm
{

inline int tableIndex(const std::vector<double>& table, double x)
{
    if (table.size() < 2)
        return 0;

    const auto lower = std::lower_bound(table.begin(), table.end(), x);

    if (lower == table.end())
        return table.size()-2;
    else if (lower == table.begin())
        return 0;
    else
      return std::distance(table.begin(), lower)-1;
}

inline std::pair<double, int>
linearInterpolationSlope(const std::vector<double>& xv,
                         const std::vector<double>& yv,
                         const double x)
{
    const auto i = Opm::tableIndex(xv, x);
    return { (yv[i + 1] - yv[i]) / (xv[i + 1] - xv[i]), i };
}


inline double linearInterpolationDerivative(const std::vector<double>& xv,
                                            const std::vector<double>& yv, double x)
{
    // Extrapolates if x is outside xv
    return linearInterpolationSlope(xv, yv, x).first;
}

inline double linearInterpolation(const std::vector<double>& xv,
                                  const std::vector<double>& yv, double x)
{
    // Extrapolates if x is outside xv
    const auto& [t, i] = linearInterpolationSlope(xv, yv, x);
    return t * (x - xv[i]) + yv[i];
}

inline double linearInterpolationNoExtrapolation(const std::vector<double>& xv,
                                                 const std::vector<double>& yv, double x)
{
    // Return end values if x is outside xv
    if (x < xv.front()) {
        return yv.front();
    }
    if (x > xv.back()) {
        return yv.back();
    }

    return linearInterpolation(xv, yv, x);
}

inline double linearInterpolation(const std::vector<double>& xv,
                                  const std::vector<double>& yv,
                                  double x, int& ix1)
{
    // Extrapolates if x is outside xv
    double t;
    std::tie(t, ix1) = linearInterpolationSlope(xv, yv, x);
    return t * (x - xv[ix1]) + yv[ix1];
}

} // namespace Opm

#endif // OPM_LINEARINTERPOLATION_HEADER_INCLUDED