/* Array tests */

#include <math.h>
#include "test.h"

#define D2R(DEG) ((DEG)*(M_PI/180.0))

#define X_STEP	5
#define N_Y	73

static double x1;

static void
remove_func(vscalar *s)
{
    v_destroy(s);
}

static double
sin_cos(double x2)
{
    return sin(D2R(x1)) * cos(D2R(x2));
}

static int
interpolate_test(varray * a, float x, float y)
{
    float v, v1, v2, fraction;
    vscalar *s1, *s2;
    vfunc *f1, *f2;
    int pts;

    /* Interpolate in the array  */
    if ((pts = va_interp(a, x, &s1, &s2, &fraction)) < 0)
	return 0;

    f1 = vs_pget(s1);
    f2 = vs_pget(s2);

    v1 = vf_value(f1, y);
    if (pts > 1) {
	v2 = vf_value_linear(f1, y);
	v = v1 + fraction * (v2 - v1);
    } else {
	v = v1;
    }

    /* Find error in interpolation */
    x1 = x;
    v -= sin_cos(y);

    /* Make sure error is small enough */
    return fabs(v) < 0.1;
}

static int
extrapolate_test(varray * a, float x, float y)
{
    float v, v1, v2, fraction;
    vscalar *s1, *s2;
    vfunc *f1, *f2;
    int pts;

    /* Interpolate in the array  */
    if ((pts = va_extrapolate(a, x, &s1, &s2, &fraction)) < 0)
	return 0;

    f1 = vs_pget(s1);
    f2 = vs_pget(s2);

    v1 = vf_value(f1, y);
    if (pts > 1) {
	v2 = vf_value_extrapolate(f1, y);
	v = v1 + fraction * (v2 - v1);
    } else {
	v = v1;
    }

    /* Find error in interpolation */
    x1 = x;
    v -= sin_cos(y);

    /* Make sure error is small enough */
    return fabs(v) < 0.1;
}

static int
lookup_test(varray *a, float x, float y,
	    enum v_interptype method, float tolerance)
{
    vscalar *points[4];
    float fractions[4];
    float v;
    int i;

    /* Interpolate in the array */
    va_lookup(a, x, method, points, fractions);

    v = 0.0;
    for (i = 0; i < 4; i++) {
	float frac = fractions[i];
	if (frac != 0.0) {
	    vfunc *f = vs_pget(points[i]);
	    
	    /* Use same interpolation method on function as on array */
	    vf_interpolate(f, method);

	    /* Interpolate & add */
	    v += frac * vf_value(f,y);
	}
    }

    /* Find error in the interpolation */
    x1 = x;
    v -= sin_cos(y);

    return fabs(v) < tolerance;
}



static int
test_replace(varray * a, float x)
{
    vfunc *f, *same;

    x1 = x;
    f = vf_create();
    vf_add_func(f, sin_cos, 0.0, 360.0, N_Y * 2);
    va_pstore(a, x1, f);

    /* Make sure its inserted  */
    same = va_pget(a, x1);

    return (same == f);
}

int
main(void)
{
    float t1, t2, t3;
    varray *a1, *a2;
    vfunc *f;
    int i;

    TEST_START("array");

    TEST("creation", (a1 = va_create()) != NULL);

    /* Create a 2D table using functions & the array */
    for (i = 0; i <= 360; i += X_STEP) {
	x1 = (double) i;

	f = vf_create();
	vf_add_func(f, sin_cos, 0.0, 360.0, N_Y);

	va_pstore(a1, x1, f);
    }

    /* Test freeze/thaw */
    TEST_NOFAIL("freezing",
                v_freeze_file(a1, TEST_OUTPUT));

    TEST_NOFAIL("thawing",
                (a2 = v_thaw_file(TEST_OUTPUT)) != NULL);

    /* Test for interpolation */
    for (i = 0; i < 5; i++) {
	t1 = v_randreal(0.0, 360.0);
	t2 = v_randreal(0.0, 360.0);

	TEST("interpolation", interpolate_test(a1, t1, t2));
	TEST("extrapolation", extrapolate_test(a1, t1, t2));

	TEST("lookup-nearest",
             lookup_test(a1, t1, t2, V_INTERP_NEAREST, 0.05));

	TEST("lookup-linear",
             lookup_test(a1, t1, t2, V_INTERP_LINEAR, 0.01));

	TEST("lookup-extrapolate",
             lookup_test(a1, t1, t2, V_INTERP_EXTRAPOLATE, 0.01));

	TEST("lookup-lagrange",
             lookup_test(a1, t1, t2, V_INTERP_LAGRANGE, 0.001));
    }

    va_remove_func(a1, remove_func);
    TEST("replacement", test_replace(a1, 179.0));
    TEST("replacement", test_replace(a1, 180.0));

    TEST("find", (i = va_find(a1, 179.0)) >= 0);
    TEST("delete", va_delete(a1, i));

    TEST("destruction", (va_destroy(a1), 1));

    TEST_FINISH;
}