All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Pages
linearInterpolation.hpp
1 /*
2  Copyright 2009, 2010, 2013 SINTEF ICT, Applied Mathematics.
3  Copyright 2009, 2010 Statoil ASA.
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef OPM_LINEARINTERPOLATION_HEADER_INCLUDED
22 #define OPM_LINEARINTERPOLATION_HEADER_INCLUDED
23 
24 
25 #include <vector>
26 #include <algorithm>
27 
28 namespace Opm
29 {
30 
31  inline int tableIndex(const std::vector<double>& table, double x)
32  {
33  // Returns an index in an ordered table such that x is between
34  // table[j] and table[j+1]. If x is out of range, first or last
35  // interval is returned; Binary search.
36  int n = table.size() - 1;
37  if (n < 2) {
38  return 0;
39  }
40  int jl = 0;
41  int ju = n;
42  bool ascend = (table[n] > table[0]);
43  while (ju - jl > 1) {
44  int jm = (ju + jl)/2; // Compute a midpoint
45  if ( (x >= table[jm]) == ascend) {
46  jl = jm; // Replace lower limit
47  } else {
48  ju = jm; // Replace upper limit
49  }
50  }
51  return jl;
52  }
53 
54 
55  inline double linearInterpolationDerivative(const std::vector<double>& xv,
56  const std::vector<double>& yv, double x)
57  {
58  // Extrapolates if x is outside xv
59  int ix1 = tableIndex(xv, x);
60  int ix2 = ix1 + 1;
61  return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1]);
62  }
63 
64  inline double linearInterpolation(const std::vector<double>& xv,
65  const std::vector<double>& yv, double x)
66  {
67  // Extrapolates if x is outside xv
68  int ix1 = tableIndex(xv, x);
69  int ix2 = ix1 + 1;
70  return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1];
71  }
72 
73  inline double linearInterpolationNoExtrapolation(const std::vector<double>& xv,
74  const std::vector<double>& yv, double x)
75  {
76  // Return end values if x is outside xv
77  if (x < xv.front()) {
78  return yv.front();
79  }
80  if (x > xv.back()) {
81  return yv.back();
82  }
83 
84  int ix1 = tableIndex(xv, x);
85  int ix2 = ix1 + 1;
86  return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1];
87  }
88 
89  inline double linearInterpolation(const std::vector<double>& xv,
90  const std::vector<double>& yv,
91  double x, int& ix1)
92  {
93  // Extrapolates if x is outside xv
94  ix1 = tableIndex(xv, x);
95  int ix2 = ix1 + 1;
96  return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1];
97  }
98 
99 
100 
101 } // namespace Opm
102 
103 
104 
105 
106 
107 #endif // OPM_LINEARINTERPOLATION_HEADER_INCLUDED