cocomice
11/5/2014 - 4:42 PM

linear interpolation in C style

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;
}