1 /**
2 @file
3 @ingroup nr
4 */
5
6 #include "bpm_nr.h"
7
nr_quadinterpol(double x,double x1,double x2,double x3,double y1,double y2,double y3)8 double nr_quadinterpol( double x,
9 double x1, double x2, double x3,
10 double y1, double y2, double y3 ) {
11
12 double d1, d2, d3, d4;
13 double a, b, c;
14
15 // if the 3 points are equal, set the interpolated point equal as well
16 if( ( fabs( y1 - y2 ) < 1.0e-15 ) &&
17 ( fabs( y2 - y3 ) < 1.0e-15 ) ) return y2;
18
19 // construct the equation of the parabola from the determinant
20 // d1->4 are the 4 3x3 minor determinants
21 // d1 * x^2 - d2 * x + d3 * y - d4 = 0
22 d1 = x1*y2 + x2*y3 + x3*y1 - x3*y2 - y3*x1 - x2*y1;
23 d2 = x1*x1*y2 + x2*x2*y3 + x3*x3*y1 - x3*x3*y2 - y3*x1*x1 - x2*x2*y1;
24 d3 = x1*x1*x2 + x2*x2*x3 + x3*x3*x1 - x3*x3*x2 - x3*x1*x1 - x2*x2*x1;
25 d4 = x1*x1*x2*y3 + x2*x2*x3*y1 + x1*y2*x3*x3
26 - x3*x3*x2*y1 - x1*x1*x3*y2 - x2*x2*x1*y3;
27
28 if( fabs( d3 ) < 1.e-15 ) {
29 //bpm_warning( "Collinearity issue in nr_quadinterpol()", __FILE__, __LINE__ );
30 return y2;
31 }
32
33 // now calculate a,b,c, don't forget the - signs !
34 // y = a * x^2 + b*x + c
35 a = - d1 / d3;
36 b = d2 / d3;
37 c = d4 / d3;
38
39 // return interpolated value at the point x
40 return a*x*x + b*x + c;
41 }
42