/*
 File :       interp.c
 Brief :      Interpolation functions.
 Author :     CLS
 */
#include <assert.h>

#include "interp.h"
#include "fes_int.h"

/*
 _linearWeighting

 Linear Weighting.

 x    Abscissa where interpolation will occur
 x1   Abscissa corresponding to the first value
 x2   Abscissa corresponding to the second value
 w1   Linear weight w1
 w2   Linear weight w2
 */
static void _linearWeighting(
        const double x,
        const double x_1,
        const double x_2,
        double* w_1,
        double* w_2)
{
    assert(x >= x_1 - EPSILON);
    assert(x <= x_2 + EPSILON);

    if (EQUALS(x_1, x_2) || EQUALS(x, x_1))
    {
        *w_1 = 1.0;
        *w_2 = 0.0;
    }
    else if (EQUALS(x, x_2))
    {
        *w_1 = 0.0;
        *w_2 = 1.0;
    }
    else
    {
        *w_1 = (x_2 - x) / (x_2 - x_1);
        *w_2 = (x - x_1) / (x_2 - x_1);
    }
}

/*
 _sumWeighting

 Sum Weighting.

 x_12
 w_12
 s
 w
 */
static void _sumWeighting(
        const double x_12,
        const double w_12,
        double* s,
        double* w)
{
    if (!EQUALS(x_12, DV))
    {
        *s += w_12 * x_12;
        *w += w_12;
    }
}

/*
 */
double bilinearInterp(
        const double x_1,
        const double x_2,
        const double y_1,
        const double y_2,
        const double value_11,
        const double value_21,
        const double value_12,
        const double value_22,
        const double x,
        const double y)
{
    double w = 0.0;
    double s = 0.0;
    double w_x1;
    double w_x2;
    double w_y1;
    double w_y2;

    _linearWeighting(x, x_1, x_2, &w_x1, &w_x2);
    _linearWeighting(y, y_1, y_2, &w_y1, &w_y2);

    _sumWeighting(value_11, w_x1 * w_y1, &s, &w);
    _sumWeighting(value_12, w_x1 * w_y2, &s, &w);
    _sumWeighting(value_21, w_x2 * w_y1, &s, &w);
    _sumWeighting(value_22, w_x2 * w_y2, &s, &w);

    return w == 0.0 ? DV : s / w;
}