source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/optics/src/Mspline.cc @ 114

/* -------------------------------------------------------------------------
 *   spline.c
 *
 *   --- cubic spline interpolation for non-periodic functions
 *        "The latest algorithm collections in C" (1991) (in Japanese)
 *        Haruhiko Okumura, Professor (Computer Science), Ph.D.
 *        Matsusaka University
 *
 *        modified by N.Sakaki
 *
 * $Id: spline.c,v 1.1 2003/03/27 07:38:45 sakaki Exp sakaki $
 * -------------------------------------------------------------------------
 */
/*
 */
#include <stdlib.h>
#include "Mspline.hh"

Int_t Mcsp_maketable(Double_t x[], Double_t y[], Double_t z[], size_t size)
{
    Int_t i;
    Double_t t;
    Double_t *h, *d;

    h=(Double_t *)malloc(sizeof(Double_t)*size);
    d=(Double_t *)malloc(sizeof(Double_t)*size);
    if(h==NULL || d==NULL)
      return -1;

    z[0] = 0;  z[size-1] = 0;  /* y"(x) / 6 at both end points */
    for(i=0; i<(Int_t)(size-1); i++)
      {
        h[i  ] =  x[i+1] - x[i];
        d[i+1] = (y[i+1] - y[i])/h[i];
      }

    z[1] = d[2] - d[1] - h[0]*z[0];
    d[1] = 2.0*(x[2]-x[0]);

    for (i=1; i<(Int_t)(size-2); i++)
      {
        t = h[i]/d[i];
        z[i+1] = d[i+2] - d[i+1] - z[i]*t;
        d[i+1] = 2.0*(x[i+2]-x[i]) - h[i]*t;
    }

    z[size-2] -= h[size-2]*z[size-1];
    for (i=size-2; i>0; i--)
        z[i] = (z[i] - h[i]*z[i+1]) /d[i];

    free(h);
    free(d);

    return 0;
}

Double_t Mcspline(Double_t t, Double_t x[], Double_t y[], Double_t z[], size_t size)
{
    Int_t i, j, k;
    Double_t d, h;

    i = 0;  j = size - 1;
    while (i<j)
      {
        k = (i+j)/2;
        if (x[k]<t)
      i = k+1;
    else
      j = k;
      }

    if (i>0)
      i--;
    h = x[i+1] - x[i];
    d = t - x[i];

    return (( (z[i+1] - z[i])*d/h + z[i]*3.0) * d
        + ( (y[i+1] - y[i])/h - (z[i]*2. + z[i+1])* h) )*d + y[i];
}

#if 0
#include <stdio.h>

#define  N  5
Double_t x[N] = {   0,      10,     20,     30,     40   },
       y[N] = { 610.66, 1227.4, 2338.1, 4244.9, 7381.2 },
       z[N];
Int_t main()
{
    Int_t i;

    csp_maketable(x, y, z, N);
    for (i = 10; i <= 30; i++)
        printf("%3d %6.1f\n", i, cspline(i, x, y, z, N));
    return EXIT_SUCCESS;
}
#endif