/* Function tests */

#include "test.h"

#define NPOINTS 20
#define VALUE   M_PI / 3
#define E_VALUE 2 * M_PI + (M_PI / 100)
#define SMALL   1.0e-6

#define MAX_PTS 100
#define STEP	(M_PI / MAX_PTS)
#define TOL	(MAX_PTS * 0.0002)

#define IN_CYCLE(v) if (v < 0.0 || v > 360.0) ok = 0

#define NCYCLE  6
#define REDUCE  20.0

static int test_integral(void);
static int test_differential(double cycle);

int
main(void)
{
    vfunc *f, *integral, *reduced;
    double x, fx, f0, f1, f2, f3;
    int i, ok;
    vlist *l;
    FILE *fp;

    TEST_START("func");

    TEST("creation 1",
         (f = vf_func(0.0, 0.0,
                      1.0, 3.0,
                      2.0, 6.0,
                      VF_END)) != NULL &&
         vf_point_count(f) == 3);

    TEST("destruction",
         (vf_destroy(f), 1));

    f = vf_create();
    for (i = 0; i < NPOINTS; i++) {
        x = 2 * M_PI * i / (NPOINTS - 1);
        vf_add_point(f, x, sin(x));
    }

    TEST("creation 2", vf_point_count(f) == NPOINTS);

#if 0
    printf("Linear: %g\n", vf_value_linear(f, VALUE));
    printf("Lagrange: %g\n", vf_value_lagrange(f, VALUE));
    printf("Spline: %g\n", vf_value_spline(f, VALUE));
    printf("Nearest: %g\n", vf_value_nearest(f, VALUE));
#endif

    TEST("linear 1", V_ABS(vf_value_linear(f, VALUE) - 0.859206) < SMALL);
    TEST("linear 2", V_ABS(vf_value_linear(f, E_VALUE) - 0.0) < SMALL);

    TEST("lagrange", V_ABS(vf_value_lagrange(f, VALUE) - 0.865892) < SMALL);

    TEST("spline", V_ABS(vf_value_spline(f, VALUE) - 0.866015) < SMALL);

    TEST("nearest", V_ABS(vf_value_nearest(f, VALUE) - 0.837166) < SMALL);

    TEST("extrapolate 1", V_ABS(vf_value_extrapolate(f, VALUE) - 0.859206) < SMALL);
    TEST("extrapolate 2", V_ABS(vf_value_extrapolate(f, E_VALUE) - 0.030846) < SMALL);

    TEST("defined 1", vf_defined(f, 2 * M_PI - 0.1));

    TEST("defined 2", !vf_defined(f, 2 * M_PI + 0.1));

    TEST("maximum", V_ABS(vf_maximum(f) - 2 * M_PI) < SMALL);

    TEST("integral", test_integral());
    TEST("differential (normal)", test_differential(-1.0));
    TEST("differential (cyclic)", test_differential(+0.5));

    vf_destroy(f);

    f = vf_create();
    for (x = 0; x < 5; x++)
        vf_add_point(f, x, x);

    TEST("X point list",
         (l = vf_points_x(f)) != NULL &&
         V_ABS(vl_dget(l, 1) - 1) < SMALL);

    TEST("FX point list",
         (l = vf_points_fx(f)) != NULL &&
         V_ABS(vl_dget(l, 4) - 4) < SMALL);

    TEST("deletion",
         vf_delete_point(f, 1) &&
         vf_point_count(f) == 4 &&
         vf_get_point(f, 3, &x, &fx) &&
         V_ABS(fx - 4) < SMALL);

    vf_empty(f);
    TEST("emptying", vf_point_count(f) == 0);

    for (i = 0, x = 0.0; i < 24; i++, x += 15.0)
        vf_add_point(f, x, fmod(400.0 * sin(x * M_PI / 180.0), 360.0));

    vf_set_cycle(f, 360.0);

    for (ok = 1, x = 0.0; x <= 360.0; x += 5.0) {
        f0 = vf_value_nearest(f, x);
	f1 = vf_value_linear(f, x);
	f2 = vf_value_lagrange(f, x);
	f3 = vf_value_spline(f, x);

#if 0
        printf("%8g -> %12g\t%12g\t%12g\t%12g\n", x, f0, f1, f2, f3);
#endif
	IN_CYCLE(f0);
	IN_CYCLE(f1);
	IN_CYCLE(f2);
	IN_CYCLE(f3);
    }

    TEST("cyclic interpolation", ok);

    f = vf_create();

    vf_add_point(f, 0.0, 0.0);
    vf_add_point(f, 1.0, 1.0);
    vf_add_point(f, 2.0, 2.0);
    vf_add_point(f, 3.0, 2.0);
    vf_add_point(f, 4.0, 2.0);
    vf_add_point(f, 5.0, 2.0);
    vf_add_point(f, 6.0, 1.0);
    vf_add_point(f, 7.0, 0.0);

    reduced = vf_reduce(f, REDUCE);

    TEST("reduction", reduced != NULL && vf_point_count(reduced) == 4);

    TEST("reduction 2",
         vf_reduce(reduced, REDUCE) == NULL);

    TEST_FINISH;
}

/* Test the integrate function */
static int
test_integral(void)
{
    int i;
    double x, fx, error;
    vfunc *sin_func, *cos_func;

    /* Create cosine function */
    cos_func = vf_create();
    for (i = 0; i < MAX_PTS; i++) {
	x = STEP * i;
	vf_add_point(cos_func, x, cos(x));
    }

    /* Integrate it */
    sin_func = vf_integral(cos_func);

    /* Test accuracy */
    for (error = 0.0, i = 0; i < MAX_PTS; i++) {
	x = STEP * i;
	fx = vf_value(sin_func, x);
	error += fabs(fx - sin(x));
    }

    /* Recycle der memory */
    vf_destroy(sin_func);
    vf_destroy(cos_func);

    /* Did it work */
    return error < TOL;
}

/* Test the differentiate function */
static int
test_differential(double cycle)
{
    int i;
    double x, fx, error;
    vfunc *sin_func, *cos_func;

    /* Create sine function */
    sin_func = vf_create();
    for (i = 0; i < MAX_PTS; i++) {
	x = STEP * i;
	vf_add_point(sin_func, x, sin(x));
    }

    /* Diffentiate it */
    vf_set_cycle(sin_func, cycle);
    cos_func = vf_differential(sin_func);

    /* Test accuracy */
    for (error = 0.0, i = 0; i < MAX_PTS; i++) {
	x = STEP * i;
	fx = vf_value(cos_func, x);
	error += fabs(fx - cos(x));
    }

    /* Recycle der memory */
    vf_destroy(sin_func);
    vf_destroy(cos_func);

    /* Did it work */
    return error < TOL;
}