interptest.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
#include <sp/spLib.h>
 
int main(int argc, char **argv)
{
    long k;
    spDVector x, y;
    spDVector xi, yi, yi2, yi3;
 
    /* See:
     *   http://oshiete1.goo.ne.jp/kotaeru.php3?q=1319182
     *   http://blog.goo.ne.jp/supurainhikaku/e/76266ba89fdf8e890f3d6be48172d8a3
     */
    x = xdvalloc(4);
    x->data[0] = -1.0;
    x->data[1] = 0.0;
    x->data[2] = 1.0;
    x->data[3] = 2.0;
    y = xdvzeros(x->length);
    y->data[1] = 1.0;
    xi = xdvinit(0.0, 0.125, 1.0);
    yi = xdvspline(x, y, xi, SP_TRUE);
    yi2 = xdvspline(x, y, xi, SP_FALSE);
#if 0
    yi3 = xdvsplinecl(x, y, xi, 0.0, 0.0);
#else
    yi3 = xdvinterp1l(x, y, xi);
#endif
    
#if 1
    for (k = 0; k < xi->length; k++) {
        printf("%f %f %f %f\n", xi->data[k], yi->data[k], yi2->data[k], yi3->data[k]);
    }
    printf("\n");
#endif
    
    xdvfree(x);
    xdvfree(y);
    xdvfree(xi);
    xdvfree(yi);
    xdvfree(yi2);
    xdvfree(yi3);
 
    /* from MATLAB 'doc interp1' */
    x = xdvinit(1900.0, 10.0, 1990.0);
    y = xdvalloc(10);
    y->data[0] = 75.995; y->data[1] = 91.972;
    y->data[2] = 105.711; y->data[3] = 123.203;
    y->data[4] = 131.669; y->data[5] = 150.697;
    y->data[6] = 179.323; y->data[7] = 203.212;
    y->data[8] = 226.505; y->data[9] = 249.633;
 
    xi = xdvinit(1900.0, 2.5, 2000.0);
    yi = xdvinterp1s(x, y, xi);
    yi2 = xdvinterp1l(x, y, xi);
 
#if 1
    for (k = 0; k < xi->length; k++) {
        printf("%f %f %f\n", xi->data[k], yi->data[k], yi2->data[k]);
    }
#endif
 
    xdvfree(x);
    xdvfree(y);
    xdvfree(xi);
    xdvfree(yi);
    xdvfree(yi2);
    
    return 0;
}