00001
00006 #include "bpm_nr.h"
00007
00008 double nr_quadinterpol( double x,
00009 double x1, double x2, double x3,
00010 double y1, double y2, double y3 ) {
00011
00012 double d1, d2, d3, d4;
00013 double a, b, c;
00014
00015
00016 if( ( fabs( y1 - y2 ) < 1.0e-15 ) &&
00017 ( fabs( y2 - y3 ) < 1.0e-15 ) ) return y2;
00018
00019
00020
00021
00022 d1 = x1*y2 + x2*y3 + x3*y1 - x3*y2 - y3*x1 - x2*y1;
00023 d2 = x1*x1*y2 + x2*x2*y3 + x3*x3*y1 - x3*x3*y2 - y3*x1*x1 - x2*x2*y1;
00024 d3 = x1*x1*x2 + x2*x2*x3 + x3*x3*x1 - x3*x3*x2 - x3*x1*x1 - x2*x2*x1;
00025 d4 = x1*x1*x2*y3 + x2*x2*x3*y1 + x1*y2*x3*x3
00026 - x3*x3*x2*y1 - x1*x1*x3*y2 - x2*x2*x1*y3;
00027
00028 if( fabs( d3 ) < 1.e-15 ) {
00029
00030 return y2;
00031 }
00032
00033
00034
00035 a = - d1 / d3;
00036 b = d2 / d3;
00037 c = d4 / d3;
00038
00039
00040 return a*x*x + b*x + c;
00041 }