Two-dimensional linear interpolation function

The function LinearInterp2D() implements linear interpolation on a two-dimensional surface.

LinearInterp2D()'s matrix argument aadZ is passed as an array of pointers to one-dimensional arrays of doubles. Each of these arrays of doubles represents one column of data. LinearInterp2D() returns 1 if it succeeds, 0 if it fails.

Add the code to the end of Tutorial1.cpp.

static int FindInterpPoints(int cPoints, double* adValues, double dInterp, int* cc);

// Interpolate in a 2 dimensional array
// Returns 1 for success, 0 for failure
int LinearInterp2D(int cX, double* adX, int cY, double *adY, double** aadZ, 
                   double dXInterp, double dYInterp, double* pdZInterp)
{
    int ccX[2], ccY[2], i;
    double adZInterpX[2];

    // Find the right poionts to use for interpolation
    if (!FindInterpPoints(cX, adX, dXInterp, ccX))
        return 0;
    if (!FindInterpPoints(cY, adY, dYInterp, ccY))
        return 0;

    // Interpolate in the X dimension
    for (i = 0; i < 2; i++) {
        if (ccX[0] == ccX[1])
            adZInterpX[i] = aadZ[ccY[i]][ccX[0]];
        else 
            adZInterpX[i] = aadZ[ccY[i]][ccX[0]]
                          + ( ( aadZ[ccY[i]][ccX[1]] - aadZ[ccY[i]][ccX[0]] )
                            * ( dXInterp - adX[ccX[0]] )
                            / ( adX[ccX[1]] - adX[ccX[0]] ) );
    }

    // Interpolate in the Y dimension
    if (ccY[0] == ccY[1])
        *pdZInterp = adZInterpX[0];
    else
        *pdZInterp = adZInterpX[0]
                   + ( ( adZInterpX[1] - adZInterpX[0] )
                     * ( dYInterp - adY[ccY[0]] )
                     / ( adY[ccY[1]] - adY[ccY[0]] ) );
    return 1;
}

// Locate the interpolation points in an ordered unique array
// Puts the coordinates of the two points in cc[] 
// Returns 1 for success, 0 for failure
static int FindInterpPoints(int cPoints, double* adValues, double dInterp, int* cc)
{
    int i;

    // Too few points
    if (cPoints < 1)
        return 0;
    // Only one point
    if (cPoints == 1) {
        cc[0] = cc[1] = 0;
        return 1;
    }
    // Only two points
    if (cPoints == 2) {
        cc[0] = 0;
        cc[1] = 1;
        return 1;
    }
    // Find position of dInterp in adValues[]
    for (i = 0; i < cPoints; i++) {
        // Check that adValues[] is strongly ordered
        if (i > 0 && adValues[i] <= adValues[i - 1])
            return 0;
        // Look for an exact match
        if (dInterp == adValues[i]) {
            cc[0] = cc[1] = i;
            return 1;
        }
        // Find upper bound
        if (dInterp < adValues[i]) {
            // Extrapolate left
            if (i == 0) {
                cc[0] = 0;
                cc[1] = 1;
            }
            // Interpolate
            else {
                cc[0] = i - 1;
                cc[1] = i;
            }
            return 1;
        }
    }
    // Extrapolate right
    cc[0] = cPoints - 2;
    cc[1] = cPoints - 1;
    return 1;
}

Next: Creating a matrix argument >>