VERB_code_2.2  2
 All Classes Namespaces Files Functions Variables Typedefs Macros Pages
polilinear.cpp
Go to the documentation of this file.
1 /**
2  * Polilinear interpolation sources.
3  * \file polilinear.cpp
4  * \author Y.Shprits
5  */
6 
7 #include "../../Exceptions/error.h"
8 
9 /** Do polinomial interpolation between 3 points for all data exepts the boundaries.
10  * Linear interpolation on the boundaries.
11  * Boundary values for extrapolation.
12  */
13 double polilinear( double *xa, double *ya, int n, double x, double lb, double ub)
14 {
15  int i = 0;
16 
17  while (i <= n-1 && x > xa[i]) {
18  i = i + 1;
19  }
20 
21 // if ((x <= xa[0])) {
22  if (i <= 0) {
23  return lb;
24  } else if (i <= 1) {
25  // linear
26  double a = (x - xa[i-1])/(xa[i] - xa[i-1]);
27  return ya[i-1] + a*(ya[i] - ya[i-1]);
28  } else if (i >= n-1) {
29  return ub;
30 // return -21;
31  } else if (i >= n-2) {
32  // linear
33  double a = (x - xa[i-1])/(xa[i] - xa[i-1]);
34  return ya[i-1] + a*(ya[i] - ya[i-1]);
35 // double a = (x - xa[i])/(xa[i-1] - xa[i]);
36 // return ya[i] + a*(ya[i-1] - ya[i]);
37 /* } else if (i == n-1) {
38  //polinom
39  double A = (x - xa[i-1])*(x - xa[i]) / (xa[i-2] - xa[i-1]) / (xa[i-2] - xa[i]);
40  double B = (xa[i-2] - x)*(x - xa[i]) / (xa[i-2] - xa[i-1]) / (xa[i-1] - xa[i]);
41  double C = (xa[i-2] - x)*(xa[i-1] - x) / (xa[i-1] - xa[i]) / (xa[i-2] - xa[i]);
42  return A*ya[i-2] + B*ya[i-1] + C*ya[i];*/
43  } else if (i >= 1 && i <= n-2) {
44  //polinom
45  double A1 = (x - xa[i])*(x - xa[i+1]) / (xa[i-1] - xa[i]) / (xa[i-1] - xa[i+1]);
46  double B1 = (xa[i-1] - x)*(x - xa[i+1]) / (xa[i-1] - xa[i]) / (xa[i] - xa[i+1]);
47  double C1 = (xa[i-1] - x)*(xa[i] - x) / (xa[i] - xa[i+1]) / (xa[i-1] - xa[i+1]);
48  double res = A1*ya[i-1] + B1*ya[i] + C1*ya[i+1];
49 // if (res > 99 || (!(res >= 0) && !(res <= 0))) {
50 // cout << "111";
51 // }
52  return res;
53  } else {
54  //cout << "interpolation error";
55  //exit(0);
56  throw error_msg("POLINOM_ERROR", "Polinomial interpolation error");
57  }
58 }