/**************************************************************************
 *{@C
 *      Copyright:      1988-2024 Paul Obermeier (obermeier@poSoft.de)
 *
 *      Module:         Utilities
 *      Filename:       UT_MathHardy.c
 *
 *      Author:         Paul Obermeier
 *
 *      Description:    Routines to interpolate smoothly between irregularly
 *                      distributed points in two dimensions.
 *
 *                      The functions defined in this file are derived from
 *                      the Fortran subroutines found on pages 177 to 179 in
 *
 *                              Helmut Spaeth
 *                              "Zweidimensionale Spline-
 *                              Interpolations-Algorithmen"
 *                              R. Oldenbourg Verlag, Muenchen, Wien, 1991
 *
 *                      The original routines in Spaeth's book have one
 *                      additional parameter, "q", which has been changed into
 *                      a constant, 0.5, in this implementation. In the book,
 *                      the examples demonstrating the routines' behavior
 *                      suggest that 0.5 is a good value for "q". This allows
 *                      for some optimizations - calls to "pow" can be replaced
 *                      by calls to the faster "sqrt".
 *                      Another usable value for "q" is -0.5. Here "sqrt" can
 *                      be used instead of "pow", too. The effect of the
 *                      data points is more "local" than with "q" set to
 *                      0.5, but the interpolating function is not as smooth.
 *
 *      Additional documentation:
 *                      None.
 *
 *      Exported functions:
 *                      UT_MathHardyCoeff1
 *                      UT_MathHardyCoeff2
 *                      UT_MathHardy1
 *                      UT_MathHardy2
 *                      UT_MathHardyGrad1
 *                      UT_MathHardyGrad2
 *
 **************************************************************************/

#include <math.h>

#include "UT_Compat.h"
#include "UT_Const.h"
#include "UT_Macros.h"
#include "UT_Error.h"
#include "UT_Memory.h"
#include "UT_Math.h"

/* Set up a system of linear equations and solve it in order to find the
   coefficients, "R2" and "a", of the function interpolating between the "n"
   points given by "x", "y" and "z". Use a "q" value of 0.5. */

static UT_Bool find_R2_a_1
        (Int32 n, const Float32 *x, const Float32 *y,
         const Float64 *z, Float64 *R2, Float64 *a)
{
    Int32   i, j, *perm;
    Float64 xmin, ymin, xmax, ymax, R, tmp, xi, yi, *B, *lu;

    if (n < 1) {
        return UT_True;
    }
    if (n == 1) {
        a[0] = z[0];
        return UT_True;
    }
    if (!(B = UT_MemTempArray (n * n, Float64)) ||
        !(lu = UT_MemTempArray (n * n, Float64)) ||
        !(perm = UT_MemTempArray (n, Int32))) {
        return UT_False;
    }
    xmin = xmax = x[0];
    ymin = ymax = y[0];
    for (i = 1; i < n; ++i) {
        if (x[i] < xmin)
            xmin = x[i];
        else if (x[i] > xmax)
            xmax = x[i];
        if (y[i] < ymin)
            ymin = y[i];
        else if (y[i] > ymax)
            ymax = y[i];
    }
    *R2 = 0.1 * UT_MAX (xmax - xmin, ymax - ymin);
    R = sqrt (*R2);
    for (i = 0; i < n; ++i) {
        UT_MATRIX (B, n, i, i) = R;
        xi = x[i];
        yi = y[i];
        for (j = i + 1; j < n; ++j) {
            UT_MATRIX (B, n, i, j) = UT_MATRIX (B, n, j, i) =
                sqrt (UT_SQR (x[j] - xi) + UT_SQR (y[j] - yi) + *R2);
        }
    }
    return UT_MathGaussSolve (UT_False, n, B, lu, perm, z, a, &tmp);
}

/* Set up a system of linear equations and solve it in order to find the
   coefficients, "R2" and "a", of the function interpolating between the "n"
   points given by "x", "y" and "z". Use a "q" value of -0.5. */

static UT_Bool find_R2_a_2
        (Int32 n, const Float32 *x, const Float32 *y,
         const Float64 *z, Float64 *R2, Float64 *a, Float32 r)
{
    Int32   i, j, *perm;
    Float64 xmin, ymin, xmax, ymax, R, tmp, xi, yi, *B, *lu;

    if (n < 1) {
        return UT_True;
    }
    if (n == 1) {
        a[0] = z[0];
        return UT_True;
    }
    if (!(B = UT_MemTempArray (n * n, Float64)) ||
        !(lu = UT_MemTempArray (n * n, Float64)) ||
        !(perm = UT_MemTempArray (n, Int32))) {
        return UT_False;
    }
    xmin = xmax = x[0];
    ymin = ymax = y[0];
    for (i = 1; i < n; ++i) {
        if (x[i] < xmin)
            xmin = x[i];
        else if (x[i] > xmax)
            xmax = x[i];
        if (y[i] < ymin)
            ymin = y[i];
        else if (y[i] > ymax)
            ymax = y[i];
    }
    *R2 = 0.1 * UT_MAX (xmax - xmin, ymax - ymin) * r;
    R = 1.0 / sqrt (*R2);
    for (i = 0; i < n; ++i) {
        UT_MATRIX (B, n, i, i) = R;
        xi = x[i];
        yi = y[i];
        for (j = i + 1; j < n; ++j) {
            UT_MATRIX (B, n, i, j) = UT_MATRIX (B, n, j, i) =
                1.0 / sqrt (UT_SQR (x[j] - xi) + UT_SQR (y[j] - yi) + *R2);
        }
    }
    return UT_MathGaussSolve (UT_False, n, B, lu, perm, z, a, &tmp);
}

/***************************************************************************
 *[@e
 *      Name:           UT_MathHardyCoeff1
 *
 *      Usage:          Compute the coefficients for a function which
 *                      interpolates smoothly between irregularly distributed
 *                      data points in two dimensions, using Hardy's method.
 *
 *      Synopsis:       UT_Bool UT_MathHardyCoeff1
 *                              (Int32 n,
 *                              const Float32 *x, const Float32 *y,
 *                              const Float64 *z,
 *                              Float64 *R2, Float64 *a)
 *
 *      Description:    "N" is the number of data points at which the value of
 *                      the interpolating function to be found is known. "x",
 *                      "y", "z" and "a" must be vectors with "n" elements.
 *                      "N" must not be less than one.
 *                      Vectors "x" and "y" contain the coordinates of the "n"
 *                      points for which the value of the function is given.
 *                      Vector "z" contains these values.
 *                      "Make_UT_MathHardy1" and "UT_MathHardyCoeff2 compute the
 *                      coefficients of the interpolating function and store
 *                      them in in "R2" and in vector "a".
 *
 *                      "Make_UT_MathHardy1" uses a "q" value of 0.5.
 *                      The resulting function is very smooth, but changing one
 *                      data point changes the interpolating function globally,
 *                      not only near the changed data point.
 *
 *                      Note: Due to numerical instability "UT_MathHardyCoeff1"
 *                      can become unreliable if the number of data points is
 *                      too high. It is recommended to keep the number of
 *                      data points below 100.
 *
 *      Return value:   UT_True if successful, UT_False if the interpolating
 *                      function's coefficients cannot be computed.
 *
 *      See also:       UT_MathHardyCoeff2
 *                      UT_MathHardy1
 *                      UT_MathHardyGrad1
 *
 ***************************************************************************/

UT_Bool UT_MathHardyCoeff1
        (Int32 n,
         const Float32 *x, const Float32 *y,
         const Float64 *z,
         Float64 *R2, Float64 *a)
{
    UT_MemState memstate;
    UT_Bool     success;

    memstate = UT_MemRemember ();
    success = find_R2_a_1 (n, x, y, z, R2, a);
    UT_MemRestore (memstate);
    return success;
}

/***************************************************************************
 *[@e
 *      Name:           UT_MathHardyCoeff2
 *
 *      Usage:          Compute the coefficients for a function which
 *                      interpolates smoothly between irregularly distributed
 *                      data points in two dimensions, using Hardy's method.
 *
 *      Synopsis:       UT_Bool UT_MathHardyCoeff2
 *                              (Int32 n,
 *                              const Float32 *x, const Float32 *y,
 *                              const Float64 *z,
 *                              Float64 *R2, Float64 *a,
 *                              Float32 r)
 *
 *      Description:    "N" is the number of data points at which the value of
 *                      the interpolating function to be found is known. "x",
 *                      "y", "z" and "a" must be vectors with "n" elements.
 *                      "N" must not be less than one.
 *                      Vectors "x" and "y" contain the coordinates of the "n"
 *                      points for which the value of the function is given.
 *                      Vector "z" contains these values.
 *                      "Make_UT_MathHardy1" and "UT_MathHardyCoeff2 compute the
 *                      coefficients of the interpolating function and store
 *                      them in in "R2" and in vector "a".
 *
 *                      "Make_UT_MathHardy2" uses a "q" value of -0.5.
 *                      The resulting function is much less smooth than one
 *                      computed with "UT_MathHardyCoeff1", but with 
 *                      "UT_MathHardyCoeff2" the effect of changing one data
 *                      point is much more limited to the area around the
 *                      changed data point. 
 *                      The interpolating function's smoothness can be tuned
 *                      somewhat by changing the "r" parameter. Usually "r"
 *                      should be set to 1.0.
 *                      The interpolating function becomes smoother for greater
 *                      "r" values. "R" must always be greater than 0.0.
 *
 *      Return value:   UT_True if successful, UT_False if the interpolating
 *                      function's coefficients cannot be computed.
 *
 *      See also:       UT_MathHardyCoeff1
 *                      UT_MathHardy2
 *                      UT_MathHardyGrad2
 *
 ***************************************************************************/

UT_Bool UT_MathHardyCoeff2
        (Int32 n,
         const Float32 *x, const Float32 *y,
         const Float64 *z,
         Float64 *R2, Float64 *a,
         Float32 r)
{
    UT_MemState memstate;
    UT_Bool     success;

    memstate = UT_MemRemember ();
    success = find_R2_a_2 (n, x, y, z, R2, a, r);
    UT_MemRestore (memstate);
    return success;
}

/***************************************************************************
 *[@e
 *      Name:           UT_MathHardy1
 *
 *      Usage:          Evaluate an interpolation function whose
 *                      coefficients have been computed using
 *                      "UT_MathHardyCoeff1".
 *
 *      Synopsis:       Float32 UT_MathHardy1
 *                              (Int32 n,
 *                              const Float32 *x, const Float32 *y,
 *                              Float64 R2, const Float64 *a,
 *                              Float32 u, Float32 v)
 *
 *      Description:    "n", "x" and "y" have the same meaning as in routines
 *                      "UT_MathHardyCoeff1" and "UT_MathHardyCoeff2".
 *                      "A" and "R2" are the coefficients computed by
 *                      "UT_MathHardyCoeff1" or "UT_MathHardyCoeff2".
 *                      "U" and "v" define the point where the function will
 *                      be evaluated.
 *
 *      Return value:   The value of the function described by "n", "R2",
 *                      "x", "y" and "a" at point (u, v).
 *
 *      See also:       UT_MathHardy2
 *                      UT_MathHardyCoeff1
 *                      UT_MathHardyGrad1
 *
 ***************************************************************************/

Float32 UT_MathHardy1 (Int32 n, const Float32 *x, const Float32 *y,
               Float64 R2, const Float64 *a,
               Float32 u, Float32 v)
{
    Int32   i;
    Float32 result;

    if (n == 1) {
        return (Float32) a[0];
    }
    result = 0.0;
    for (i = 0; i < n; ++i) {
        result += (Float32) (a[i] * 
                             sqrt (UT_SQR(u-x[i]) + UT_SQR(v-y[i]) + R2));
    }
    return result;
}

/***************************************************************************
 *[@e
 *      Name:           UT_MathHardy2
 *
 *      Usage:          Evaluate an interpolation function whose
 *                      coefficients have been computed using
 *                      "UT_MathHardyCoeff1" or "UT_MathHardyCoeff2".
 *
 *      Synopsis:       Float32 UT_MathHardy2
 *                              (Int32 n,
 *                              const Float32 *x, const Float32 *y,
 *                              Float64 R2, const Float64 *a,
 *                              Float32 u, Float32 v)
 *
 *      Description:    "N", "x" and "y" have the same meaning as in routines
 *                      "UT_MathHardyCoeff1" and "UT_MathHardyCoeff2".
 *                      "A" and "R2" are the coefficients computed by
 *                      "UT_MathHardyCoeff1" or "UT_MathHardyCoeff2".
 *                      "U" and "v" define the point where the function will
 *                      be evaluated.
 *
 *      Return value:   The value of the function described by "n", "R2",
 *                      "x", "y" and "a" at point (u, v).
 *
 *      See also:       UT_MathHardy1
 *                      UT_MathHardyCoeff2
 *                      UT_MathHardyGrad2
 *
 ***************************************************************************/

Float32 UT_MathHardy2 (Int32 n, const Float32 *x, const Float32 *y,
               Float64 R2, const Float64 *a,
               Float32 u, Float32 v)
{
    Int32   i;
    Float32 result;

    if (n == 1) {
        return (Float32) a[0];
    }
    result = 0.0;
    for (i = 0; i < n; ++i) {
        result += (Float32) (a[i] /
                             sqrt (UT_SQR(u-x[i]) + UT_SQR(v-y[i]) + R2));
    }
    return result;
}

/***************************************************************************
 *[@e
 *      Name:           UT_MathHardyGrad1
 *
 *      Usage:          Find the gradient of an interpolation function whose
 *                      coefficients have been computed using 
 *                      "UT_MathHardyCoeff1" or "UT_MathHardyCoeff2".
 *
 *      Synopsis:       void UT_MathHardyGrad1
 *                              (Int32 n,
 *                              const Float32 *x, const Float32 *y,
 *                              Float64 R2, const Float64 *a,
 *                              Float32 u, Float32 v,
 *                              Float32 *dx, Float32 *dy)
 *
 *      Description:    "N", "x" and "y" have the same meaning as in routines
 *                      "UT_MathHardyCoeff1" and "UT_MathHardyCoeff2".
 *                      "A" and "R2" are the coefficients computed by
 *                      "UT_MathHardyCoeff1" or "UT_MathHardyCoeff2".
 *                      "U" and "v" define the point where the gradient is to
 *                      be computed.
 *
 *                      The x and y components of the gradient at point (u, v)
 *                      are returned in "dx" and "dy".
 *
 *      Return value:   None.
 *
 *      See also:       UT_MathHardyGrad2
 *                      UT_MathHardyCoeff1
 *                      UT_MathHardy1
 *
 ***************************************************************************/

void UT_MathHardyGrad1
        (Int32 n, const Float32 *x, const Float32 *y,
         Float64 R2, const Float64 *a,
         Float32 u, Float32 v, Float32 *dx, Float32 *dy)
{
    Int32   i;
    Float64 tmp;

    *dx = 0.0;
    *dy = 0.0;
    if (n == 1) {
        return;
    }
    for (i = 0; i < n; ++i) {
        tmp = a[i] / sqrt (UT_SQR (u - x[i]) + UT_SQR (v - y[i]) + R2);
        *dx += (Float32) (tmp * (u - x[i]));
        *dy += (Float32) (tmp * (v - y[i]));
    }
    return;
}

/***************************************************************************
 *[@e
 *      Name:           UT_MathHardyGrad2
 *
 *      Usage:          Find the gradient of an interpolation function whose
 *                      coefficients have been computed using
 *                      "UT_MathHardyCoeff1" or "UT_MathHardyCoeff2".
 *
 *      Synopsis:       void UT_MathHardyGrad2
 *                              (Int32 n,
 *                              const Float32 *x, const Float32 *y,
 *                              Float64 R2, const Float64 *a,
 *                              Float32 u, Float32 v,
 *                              Float32 *dx, Float32 *dy)
 *
 *      Description:    "N", "x" and "y" have the same meaning as in routines
 *                      "UT_MathHardyCoeff1" and "UT_MathHardyCoeff2".
 *                      "A" and "R2" are the coefficients computed by
 *                      "UT_MathHardyCoeff1" or "UT_MathHardyCoeff2".
 *                      "U" and "v" define the point where the gradient is to
 *                      be computed.
 *
 *                      The x and y components of the gradient at point (u, v)
 *                      are returned in "dx" and "dy".
 *
 *      Return value:   None.
 *
 *      See also:       UT_MathHardyGrad1
 *                      UT_MathHardyCoeff2
 *                      UT_MathHardy2
 *
 ***************************************************************************/

void UT_MathHardyGrad2
        (Int32 n, const Float32 *x, const Float32 *y,
         Float64 R2, const Float64 *a,
         Float32 u, Float32 v, Float32 *dx, Float32 *dy)
{
    Int32   i;
    Float64 tmp;

    *dx = 0.0;
    *dy = 0.0;
    if (n == 1) {
        return;
    }
    for (i = 0; i < n; ++i) {
        tmp = sqrt (UT_SQR (u - x[i]) + UT_SQR (v - y[i]) + R2);
        tmp = -a[i] / (tmp * tmp * tmp);
        *dx += (Float32) (tmp * (u - x[i]));
        *dy += (Float32) (tmp * (v - y[i]));
    }
    return;
}
