linear interpolation in C style
double interp_lin(int array_size, double X[], double Y[], double x){
int dim = array_size-1;
double y;
// extreme cases (x<X(0) or x>X(end): extrapolation
if(x <= X[0]){
y = ( Y[1] - Y[0] ) / ( X[1] - X[0] )*( x - X[0] ) + Y[0] ;
return y;
}
if(x >= X[dim]){
y = Y[dim] + ( Y[dim] - Y[dim-1] ) / ( X[dim] - X[dim-1] ) * ( x - X[dim] );
return y;
}
// otherwise
// [ x - X(A) ] / [ X(B) - x ] = [ y - Y(A) ] / [ Y(B) - y ]
// y = [ Y(B)*x - X(A)*Y(B) + X(B)*Y(A) - x*Y(A) ] / [ X(B) - X(A) ]
double delta = 0.0;
double min_d = INFINITY;
int j = -99, k;
signed int i ;
for( i=0; i<dim+1; i++ ){
if(X[i] == x){
y = Y[i];
return y;
}
delta = fabs( X[i] - x ) ;
if(delta < min_d){
min_d = delta ;
j = i;
}
}
if(X[j] < x){
k = j;
}else{
k = j-1;
}
double a = (Y[k+1] - Y[k]) / (X[k+1] - X[k]) ;
double b = Y[k] - a*X[k];
y = a*x + b;
return y;
}