Interpolation.c

Go to the documentation of this file.

#include <stdio.h>
#include <conio.h>
#include <math.h>
#include "Marq.h"
#include "nrutil.h"
//#include "driver.h"
/* 
INTERPOLATION routines

  Danbo 12/02/04

*/
#define TINY 1.0e-25 // A small number.
#define FREERETURN {free_vector(d,1,n);free_vector(c,1,n);return;}
#define PI 3.1415926 

/*
Here is a routine for polynomial interpolation or extrapolation from N input
points. Note that the input arrays are assumed to be unit-offset. If you have
zero-offset arrays, remember to subtract 1 (see §1.2):
*/

/*Given arrays xa[1..n] and ya[1..n], and given a value x, this routine 
    returns a value y, and an error estimate dy.
    If P(x) is the polynomial of degree N - 1 such that P(xai) = yai, i =
    1, . . . , n, then the returned value y = P(x).
*/
void polint(float xa[], float ya[], int n, float x, float *y, float *dy)
{
    int i,m,ns=1;
    float den,dif,dift,ho,hp,w;
    float *c,*d;
    dif=fabs(x-xa[1]);
    c=vector(1,n);
    d=vector(1,n);
    for (i=1;i<=n;i++)
    {                       //Here we .nd the index ns of the closest table entry,
        if((dift=fabs(x-xa[i])) < dif)
        {
            ns=i;
            dif=dift;
        }
        c[i]=ya[i]; //and initialize the tableau of c’s and d’s.
        d[i]=ya[i];
    }
    *y=ya[ns--]; //This is the initial approximation to y.
    for (m=1;m<n;m++)
    {                   //For each column of the tableau,
        for(i=1;i<=n-m;i++)
        {           //we loop over the current c’s and d’s and update them.     
            ho=xa[i]-x;
            hp=xa[i+m]-x;
            w=c[i+1]-d[i];
            if((den=ho-hp) == 0.0) //This error can occur only if two input 
                                    //xa’s are (to within roundo.) identical.
                nrerror("Error in routine polint");
            den=w/den;
            d[i]=hp*den; //Here the c’s and d’s are updated.
            c[i]=ho*den;
        }
        *y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));
        /*
        After each column in the tableau is completed, we decide which correction,
        c or d, we want to add to our accumulating value of y,
        i.e., which path to take through the tableau forking up or down.
        We do this in such a way as to take the most  straight
        line route through the tableau to its apex, updating ns accordingly ù       to keep track of where we are. 
        This route keeps the partial approximations centered (insofar as possible)
        on the target x. The last dy added is thus the error indication.
        */
    }
    free_vector(d,1,n);
    free_vector(c,1,n);
}


/*Given arrays xa[1..n] and ya[1..n], and given a value of x, this routine returns a value of
y and an accuracy estimate dy. The value returned is that of the diagonal rational function,
evaluated at x, which passes through the n points (xai, yai), i = 1...n.
*/
void ratint(float xa[], float ya[], int n, float x, float *y, float *dy)

{
    int m,i,ns=1;
    float w,t,hh,h,dd,*c,*d;
    c=vector(1,n);
    d=vector(1,n);
    hh=fabs(x-xa[1]);
    for (i=1;i<=n;i++)
    {
        h=fabs(x-xa[i]);
        if (h == 0.0)
        {
            *y=ya[i];
            *dy=0.0;
            FREERETURN
        }
        else if(h < hh)
        {
            ns=i;
            hh=h;
        }
        c[i]=ya[i];
        d[i]=ya[i] + TINY; //The TINY part is needed to prevent a rare zero-over-zero condition.
    }
    *y=ya[ns--];
    for (m=1;m<n;m++)
    {
        for (i=1;i<=n-m;i++)
        {
            w=c[i+1]-d[i];
            h=xa[i+m]-x; //h will never be zero, since this was 
                         //tested in the initializing loop.
            t=(xa[i]-x)*d[i]/h;
            dd=t-c[i+1];
            if (dd == 0.0)
                nrerror("Error in routine ratint");
                //This error condition indicates that the interpolating function has a pole at the
                //requested value of x.
            dd=w/dd;
            d[i]=c[i+1]*dd;
            c[i]=t*dd;
        }
        *y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));
    }
    FREERETURN
}


/*Given arrays x[1..n] and y[1..n] containing a tabulated function, i.e., yi = f(xi), with
x1 < x2 < .. . < xN, and given values yp1 and ypn for the first derivative of the interpolating
function at points 1 and n, respectively, this routine returns an array y2[1..n] that contains
the second derivatives of the interpolating function at the tabulated points xi. If yp1 and/or
ypn are equal to 1 × 1030 or larger, the routine is signaled to set the corresponding boundary
condition for a natural spline, with zero second derivative on that boundary.
*/
void spline(float x[], float y[], int n, float yp1, float ypn, float y2[])
{
    int i,k;
    float p,qn,sig,un,*u;
    u=vector(1,n-1);
    if (yp1 > 0.99e30) //The lower boundary condition is set either to be “natural”
        y2[1]=u[1]=0.0;
    else
    {       //or else to have a specified first derivative.
        y2[1] = -0.5;
        u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
    }
    for (i=2;i<=n-1;i++)
    { //This is the decomposition loop of the tridiagonal algorithm.
      //    y2 and u are used for temporary storage of the decomposed factors.
        sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
        p=sig*y2[i-1]+2.0;
        y2[i]=(sig-1.0)/p;
        u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
        u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
    }
    if (ypn > 0.99e30) //The upper boundary condition is set either to be natural
        qn=un=0.0;
    else
    {   //or else to have a specified first derivative.
        qn=0.5;
        un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
    }
    y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
    for (k=n-1;k>=1;k--) //This is the backsubstitution loop of the tridiagonal algorithm.
        y2[k]=y2[k]*y2[k+1]+u[k];
    
    free_vector(u,1,n-1);
}
/*
It is important to understand that the program spline is called only once to
process an entire tabulated function in arrays xi and yi. Once this has been done,
values of the interpolated function for any value of x are obtained by calls (as many
as desired) to a separate routine splint (for “spline interpolation”):
*/

/*Given the arrays xa[1..n] and ya[1..n], which tabulate a function (with the xai’s in order),
and given the array y2a[1..n], which is the output from spline above, and given a value of
x, this routine returns a cubic-spline interpolated value y.
*/
void splint(float xa[], float ya[], float y2a[], int n, float x, float *y)
{
    void nrerror(char error_text[]);
    int klo,khi,k;
    float h,b,a;
    /*
    We will find the right place in the table by means of
    bisection. This is optimal if sequential calls to this
    routine are at random values of x. If sequential calls
    are in order, and closely spaced, one would do better
    to store previous values of klo and khi and test if
    they remain appropriate on the next call.
    */
    klo=1;
    khi=n;
    while (khi-klo > 1)
    {
        k=(khi+klo) >> 1;
        if (xa[k] > x)
            khi=k;
        else klo=k;
    }   //klo and khi now bracket the input value of x.
    h=xa[khi]-xa[klo];
    if (h == 0.0)
        nrerror("Bad xa input to routine splint"); //The xa’s must be distinct.
    a=(xa[khi]-x)/h;
    b=(x-xa[klo])/h; //Cubic spline polynomial is now evaluated.
    *y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
}

/* Driver for polint routine */
int testpolintshort(void)
{

    int i,n,nfunc;
    float dy,f,x,y,*xa,*ya;

    printf("generation of interpolation tables\n");
    printf(" ... sin(x) 0<x<PI\n");
    printf(" ... exp(x) 0<x<1 \n");
    printf("how many entries go in these tables?\n");
    if (scanf("%d",&n) == EOF) return 1;

    for (nfunc=1;nfunc<=2;nfunc++)
    {
        xa=vector(1,n);
        ya=vector(1,n) ;

        
        if (nfunc == 1)
        {
            printf("\nsine function from 0 to PI\n");
            printf("\n   x     y   \n");
            for (i=1;i<=n;i++)
            {
                xa[i] =i*PI/n;
                ya[i] =sin(xa[i]);
                printf("%12.6f %12.6f \n", xa[i],ya[i]);

            }

        }
        else if( nfunc == 2)
        {
            printf("\nexponential function from 0 to 1\n");
            printf("\n   x     y   \n");
            for (i=1;i<=n;i++)
            {
                xa[i] =i*1.0/n;
                ya[i]=exp(xa[i]) ;
                printf("%12.6f %12.6f \n", xa[i],ya[i]);

            }
        }

        else
        {
            break;
        }
        
        printf("\n%9s %13s %16s %13s\n", "x", "f(x)", "interpolated", "error");
//      for (i=1;i<=10;i++)
//      {
            if (nfunc == 1)
            {
                //x=(-0.05+i/10.0)*PI;
                x=(0.5)*PI;
                f=sin(x);
            }
            else if (nfunc == 2)
            {
                //x=(-0.05+i/10.0);
                x=(0.5);
                f=exp(x);
            }
/*          else if (nfunc == 3)
            {
                x=(-0.05+i/10.0);
                f=x*x ;
            }
*/          
            //polint(xa,ya,n,x,&y,&dy);
            //polint(&xa[4],&ya[4],3,x,&y,&dy);
            //printf("%12.6f %12.6f %12.6f %4s %11f\n", x,f,y,"    ",dy);
//      }
            polint(&xa[n/2],&ya[n/2],5,x,&y,&dy);
            printf("%12.6f %12.6f %12.6f %4s %11f\n", x,f,y,"    ",dy);

    free_vector(ya,1,n);
    free_vector(xa,1,n);


    printf("\n***********************************\n");
    printf("press RETURN\n");
    (void) getchar();

        

    }


return 0;

}

int testpolintfull(void)
{

    int i,n,nfunc;
    float dy,f,x,y,*xa,*ya;

    printf("generation of interpolation tables\n");
    printf(" ... sin(x) 0<x<PI\n");
    printf(" ... exp(x) 0<x<1 \n");
    printf("how many entries go in these tables?\n");
    if (scanf("%d",&n) == EOF) return 1;
    xa=vector(1,n);
    ya=vector(1,n) ;
    for (nfunc=1;nfunc<=2;nfunc++)
    {
        if (nfunc == 1)
        {
            printf("\nsine function from 0 to PI\n");
            for (i=1;i<=n;i++)
            {
                xa[i] =i*PI/n;
                ya[i] =sin(xa[i]);

            }

        }
        else if( nfunc == 2)
        {
            printf("\nexponential function from 0 to 1\n");
            for (i=1;i<=n;i++)
            {
                xa[i] =i*1.0/n;
                ya[i]=exp(xa[i]) ;
            }
        }
        else
        {
            break;
        }
        
        printf("\n%9s %13s %16s %13s\n", "x", "f(x)", "interpolated", "error");
        for (i=1;i<=10;i++)
        {
            if (nfunc == 1)
            {
                x=(-0.05+i/10.0)*PI;
                f=sin(x);
            }
            else if (nfunc == 2)
            {
                x=(-0.05+i/10.0);
                f=exp(x);
            }
            
            polint(xa,ya,n,x,&y,&dy);
            printf("%12.6f %12.6f %12.6f %4s %11f\n", x,f,y,"    ",dy);
        }

    printf("\n***********************************\n");
    printf("press RETURN\n");
    (void) getchar();

        

    }
free_vector(ya,1,n);
free_vector(xa,1,n);
return 0;

}

int testpolintlambda(void)
{

    int i,n, l = 1, linit;
    float dy,f,x,y;
    FILE *nf;
    float *xa,*ya;
    xa=dvector(1,137);
    ya=dvector(1,137) ;


    if((nf = fopen("In.dat", "r")) == NULL)
        nrerror("In.dat not found\n");
    while(!feof(nf))
    {
        fscanf(nf, "%f%f", &xa[l], &ya[l]);
        printf("%d %12.6f %12.6f \n", l, xa[l],ya[l]);

        l++;
    }
    fclose(nf);
    linit = xa[1];

    nf = fopen("Is.dat", "w");

    for (i=1; i<l-2;i++)
    {
        polint(&xa[i], &ya[i],2,(float) i+linit,&y,&dy);
        printf("\n%d %12.6f %12.6f ", i, xa[i],ya[i]);
        printf("\n%d %d %12.6f %12.6f", i, i+linit,y,dy);
        printf("\n**************");

        fprintf(nf, "%12.6lf %12.6lf %d  %12.6lf \n",xa[i],ya[i], i+linit, y);

    }

    fclose(nf);

}

int testsplint(void)
{

    int i,n,nfunc;
    float f,x,y,yp1,ypn,*xa,*ya,*y2;

    printf("generation of interpolation tables\n");
    printf(" ... sin(x) 0<x<PI\n");
    printf(" ... exp(x) 0<x<1 \n");
    printf("how many entries go in these tables?\n");
    if (scanf("%d",&n) == EOF) return 1;
    xa=vector(1,n);
    ya=vector(1,n) ;
    y2=vector(1,n) ;
    for (nfunc=1;nfunc<=2;nfunc++)
    {
        if (nfunc == 1)
        {
            printf("\nsine function from 0 to PI\n");
            for (i=1;i<=n;i++)
            {
                xa[i] =i*PI/n;
                ya[i] =sin(xa[i]);

            }
            yp1 = cos(xa[1]);
            ypn = cos(xa[n]);
        }
        else if( nfunc == 2)
        {
            printf("\nexponential function from 0 to 1\n");
            for (i=1;i<=n;i++)
            {
                xa[i] =i*1.0/n;
                ya[i]=exp(xa[i]) ;
            }
            yp1 = exp(xa[1]);
            ypn = exp(xa[n]);

        }
        else
        {
            break;
        }
        
        spline(xa,ya,n,yp1,ypn,y2);

        printf("\n%9s %13s %16s %13s\n", "x", "f(x)", "interpolated", "error");
        for (i=1;i<=10;i++)
        {
            if (nfunc == 1)
            {
                x=(-0.05+i/10.0)*PI;
                f=sin(x);
            }
            else if (nfunc == 2)
            {
                x=(-0.05+i/10.0);
                f=exp(x);
            }
            
            splint(xa,ya,y2,n,x,&y);
            printf("%12.6f %12.6f %12.6f \n", x,f,y);
        }

    printf("\n***********************************\n");
    printf("press RETURN\n");
    (void) getchar();

        

    }
free_vector(ya,1,n);
free_vector(xa,1,n);
return 0;

}

int testsplintrc(void)
{

    int i,n,nfunc, linit, l = 1;
    float f,x,y,yp1,ypn;
    float *xa,*ya,*y2;
    FILE *nf;
    double a,b;

    if((nf = fopen("In.dat", "r")) == NULL)
        nrerror("In.dat not found\n");
    while(!feof(nf))
    {
        fscanf(nf, "%lf%lf", &a, &b);
        l++;
    }
    l--;
    fclose(nf);

    xa=vector(1,l);
    ya=vector(1,l) ;
    y2=vector(1,l) ;


    if((nf = fopen("In.dat", "r")) == NULL)
        nrerror("In.dat not found\n");


    for (i=1; i<l; i++)
    {
        fscanf(nf, "%f%f", &xa[i], &ya[i]);
        printf("%d %12.6f %12.6f \n", i, xa[i],ya[i]);

    }
    printf("**********************************************\n");

    fclose(nf);

/*
    l=1;
    while(!feof(nf))
    {
        fscanf(nf, "%f%f", &xa[l], &ya[l]);
        printf("%d %12.6f %12.6f \n", l, xa[l],ya[l]);
        l++;

    }
    printf("**********************************************\n");
    l--;
    fclose(nf);
*/
    linit = xa[1];


    yp1 = 1e30;
    ypn = 1e30;
    
        
        spline(xa,ya,l,yp1,ypn,y2);


    nf = fopen("Is.dat", "w");

    for (i=1; i<l;i++)
    {

        splint(&xa[i], &ya[i],&y2[i],3,(float) i+linit,&y);

        //polint(&xa[i], &ya[i],3,(float) i+linit,&y,&dy);
        printf("\n%d %12.6f %12.6f %d %d %12.6f ", i, xa[i],ya[i],i, i+linit,y);
        //printf("\n%d %d %12.6f ", i, i+linit,y);
        //printf("\n**************");

        fprintf(nf, "%12.6lf %12.6lf %d  %12.6lf \n",xa[i],ya[i], i+linit, y);

    }

    fclose(nf);

free_vector(ya,1,n);
free_vector(xa,1,n);
return 0;

}

        to keep track of where we are. 
        This route keeps the partial approximations centered (insofar as possible)
        on the target x. The last dy added is thus the error indication.
        */
    }
    free_vector(d,1,n);
    free_vector(c,1,n);
}


/*Given arrays xa[1..n] and ya[1..n], and given a value of x, this routine returns a value of
y and an accuracy estimate dy. The value returned is that of the diagonal rational function,
evaluated at x, which passes through the n points (xai, yai), i = 1...n.
*/
void ratint(float xa[], float ya[], int n, float x, float *y, float *dy)

{
    int m,i,ns=1;
    float w,t,hh,h,dd,*c,*d;
    c=vector(1,n);
    d=vector(1,n);
    hh=fabs(x-xa[1]);
    for (i=1;i<=n;i++)
    {
        h=fabs(x-xa[i]);
        if (h == 0.0)
        {
            *y=ya[i];
            *dy=0.0;
            FREERETURN
        }
        else if(h < hh)
        {
            ns=i;
            hh=h;
        }
        c[i]=ya[i];
        d[i]=ya[i] + TINY; //The TINY part is needed to prevent a rare zero-over-zero condition.
    }
    *y=ya[ns--];
    for (m=1;m<n;m++)
    {
        for (i=1;i<=n-m;i++)
        {
            w=c[i+1]-d[i];
            h=xa[i+m]-x; //h will never be zero, since this was 
                         //tested in the initializing loop.
            t=(xa[i]-x)*d[i]/h;
            dd=t-c[i+1];
            if (dd == 0.0)
                nrerror("Error in routine ratint");
                //This error condition indicates that the interpolating function has a pole at the
                //requested value of x.
            dd=w/dd;
            d[i]=c[i+1]*dd;
            c[i]=t*dd;
        }
        *y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));
    }
    FREERETURN
}


/*Given arrays x[1..n] and y[1..n] containing a tabulated function, i.e., yi = f(xi), with
x1 < x2 < .. . < xN, and given values yp1 and ypn for the first derivative of the interpolating
function at points 1 and n, respectively, this routine returns an array y2[1..n] that contains
the second derivatives of the interpolating function at the tabulated points xi. If yp1 and/or
ypn are equal to 1 × 1030 or larger, the routine is signaled to set the corresponding boundary
condition for a natural spline, with zero second derivative on that boundary.
*/
void spline(float x[], float y[], int n, float yp1, float ypn, float y2[])
{
    int i,k;
    float p,qn,sig,un,*u;
    u=vector(1,n-1);
    if (yp1 > 0.99e30) //The lower boundary condition is set either to be “natural”
        y2[1]=u[1]=0.0;
    else
    {       //or else to have a specified first derivative.
        y2[1] = -0.5;
        u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
    }
    for (i=2;i<=n-1;i++)
    { //This is the decomposition loop of the tridiagonal algorithm.
      //    y2 and u are used for temporary storage of the decomposed factors.
        sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
        p=sig*y2[i-1]+2.0;
        y2[i]=(sig-1.0)/p;
        u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
        u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
    }
    if (ypn > 0.99e30) //The upper boundary condition is set either to be natural
        qn=un=0.0;
    else
    {   //or else to have a specified first derivative.
        qn=0.5;
        un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
    }
    y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
    for (k=n-1;k>=1;k--) //This is the backsubstitution loop of the tridiagonal algorithm.
        y2[k]=y2[k]*y2[k+1]+u[k];
    
    free_vector(u,1,n-1);
}
/*
It is important to understand that the program spline is called only once to
process an entire tabulated function in arrays xi and yi. Once this has been done,
values of the interpolated function for any value of x are obtained by calls (as many
as desired) to a separate routine splint (for “spline interpolation”):
*/

/*Given the arrays xa[1..n] and ya[1..n], which tabulate a function (with the xai’s in order),
and given the array y2a[1..n], which is the output from spline above, and given a value of
x, this routine returns a cubic-spline interpolated value y.
*/
void splint(float xa[], float ya[], float y2a[], int n, float x, float *y)
{
    void nrerror(char error_text[]);
    int klo,khi,k;
    float h,b,a;
    /*
    We will find the right place in the table by means of
    bisection. This is optimal if sequential calls to this
    routine are at random values of x. If sequential calls
    are in order, and closely spaced, one would do better
    to store previous values of klo and khi and test if
    they remain appropriate on the next call.
    */
    klo=1;
    khi=n;
    while (khi-klo > 1)
    {
        k=(khi+klo) >> 1;
        if (xa[k] > x)
            khi=k;
        else klo=k;
    }   //klo and khi now bracket the input value of x.
    h=xa[khi]-xa[klo];
    if (h == 0.0)
        nrerror("Bad xa input to routine splint"); //The xa’s must be distinct.
    a=(xa[khi]-x)/h;
    b=(x-xa[klo])/h; //Cubic spline polynomial is now evaluated.
    *y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
}

/* Driver for polint routine */
int testpolintshort(void)
{

    int i,n,nfunc;
    float dy,f,x,y,*xa,*ya;

    printf("generation of interpolation tables\n");
    printf(" ... sin(x) 0<x<PI\n");
    printf(" ... exp(x) 0<x<1 \n");
    printf("how many entries go in these tables?\n");
    if (scanf("%d",&n) == EOF) return 1;

    for (nfunc=1;nfunc<=2;nfunc++)
    {
        xa=vector(1,n);
        ya=vector(1,n) ;

        
        if (nfunc == 1)
        {
            printf("\nsine function from 0 to PI\n");
            printf("\n   x     y   \n");
            for (i=1;i<=n;i++)
            {
                xa[i] =i*PI/n;
                ya[i] =sin(xa[i]);
                printf("%12.6f %12.6f \n", xa[i],ya[i]);

            }

        }
        else if( nfunc == 2)
        {
            printf("\nexponential function from 0 to 1\n");
            printf("\n   x     y   \n");
            for (i=1;i<=n;i++)
            {
                xa[i] =i*1.0/n;
                ya[i]=exp(xa[i]) ;
                printf("%12.6f %12.6f \n", xa[i],ya[i]);

            }
        }

        else
        {
            break;
        }
        
        printf("\n%9s %13s %16s %13s\n", "x", "f(x)", "interpolated", "error");
//      for (i=1;i<=10;i++)
//      {
            if (nfunc == 1)
            {
                //x=(-0.05+i/10.0)*PI;
                x=(0.5)*PI;
                f=sin(x);
            }
            else if (nfunc == 2)
            {
                //x=(-0.05+i/10.0);
                x=(0.5);
                f=exp(x);
            }
/*          else if (nfunc == 3)
            {
                x=(-0.05+i/10.0);
                f=x*x ;
            }
*/          
            //polint(xa,ya,n,x,&y,&dy);
            //polint(&xa[4],&ya[4],3,x,&y,&dy);
            //printf("%12.6f %12.6f %12.6f %4s %11f\n", x,f,y,"    ",dy);
//      }
            polint(&xa[n/2],&ya[n/2],5,x,&y,&dy);
            printf("%12.6f %12.6f %12.6f %4s %11f\n", x,f,y,"    ",dy);

    free_vector(ya,1,n);
    free_vector(xa,1,n);


    printf("\n***********************************\n");
    printf("press RETURN\n");
    (void) getchar();

        

    }


return 0;

}

int testpolintfull(void)
{

    int i,n,nfunc;
    float dy,f,x,y,*xa,*ya;

    printf("generation of interpolation tables\n");
    printf(" ... sin(x) 0<x<PI\n");
    printf(" ... exp(x) 0<x<1 \n");
    printf("how many entries go in these tables?\n");
    if (scanf("%d",&n) == EOF) return 1;
    xa=vector(1,n);
    ya=vector(1,n) ;
    for (nfunc=1;nfunc<=2;nfunc++)
    {
        if (nfunc == 1)
        {
            printf("\nsine function from 0 to PI\n");
            for (i=1;i<=n;i++)
            {
                xa[i] =i*PI/n;
                ya[i] =sin(xa[i]);

            }

        }
        else if( nfunc == 2)
        {
            printf("\nexponential function from 0 to 1\n");
            for (i=1;i<=n;i++)
            {
                xa[i] =i*1.0/n;
                ya[i]=exp(xa[i]) ;
            }
        }
        else
        {
            break;
        }
        
        printf("\n%9s %13s %16s %13s\n", "x", "f(x)", "interpolated", "error");
        for (i=1;i<=10;i++)
        {
            if (nfunc == 1)
            {
                x=(-0.05+i/10.0)*PI;
                f=sin(x);
            }
            else if (nfunc == 2)
            {
                x=(-0.05+i/10.0);
                f=exp(x);
            }
            
            polint(xa,ya,n,x,&y,&dy);
            printf("%12.6f %12.6f %12.6f %4s %11f\n", x,f,y,"    ",dy);
        }

    printf("\n***********************************\n");
    printf("press RETURN\n");
    (void) getchar();

        

    }
free_vector(ya,1,n);
free_vector(xa,1,n);
return 0;

}

int testpolintlambda(void)
{

    int i,n, l = 1, linit;
    float dy,f,x,y;
    FILE *nf;
    float *xa,*ya;
    xa=dvector(1,137);
    ya=dvector(1,137) ;


    if((nf = fopen("In.dat", "r")) == NULL)
        nrerror("In.dat not found\n");
    while(!feof(nf))
    {
        fscanf(nf, "%f%f", &xa[l], &ya[l]);
        printf("%d %12.6f %12.6f \n", l, xa[l],ya[l]);

        l++;
    }
    fclose(nf);
    linit = xa[1];

    nf = fopen("Is.dat", "w");

    for (i=1; i<l-2;i++)
    {
        polint(&xa[i], &ya[i],2,(float) i+linit,&y,&dy);
        printf("\n%d %12.6f %12.6f ", i, xa[i],ya[i]);
        printf("\n%d %d %12.6f %12.6f", i, i+linit,y,dy);
        printf("\n**************");

        fprintf(nf, "%12.6lf %12.6lf %d  %12.6lf \n",xa[i],ya[i], i+linit, y);

    }

    fclose(nf);

}

int testsplint(void)
{

    int i,n,nfunc;
    float f,x,y,yp1,ypn,*xa,*ya,*y2;

    printf("generation of interpolation tables\n");
    printf(" ... sin(x) 0<x<PI\n");
    printf(" ... exp(x) 0<x<1 \n");
    printf("how many entries go in these tables?\n");
    if (scanf("%d",&n) == EOF) return 1;
    xa=vector(1,n);
    ya=vector(1,n) ;
    y2=vector(1,n) ;
    for (nfunc=1;nfunc<=2;nfunc++)
    {
        if (nfunc == 1)
        {
            printf("\nsine function from 0 to PI\n");
            for (i=1;i<=n;i++)
            {
                xa[i] =i*PI/n;
                ya[i] =sin(xa[i]);

            }
            yp1 = cos(xa[1]);
            ypn = cos(xa[n]);
        }
        else if( nfunc == 2)
        {
            printf("\nexponential function from 0 to 1\n");
            for (i=1;i<=n;i++)
            {
                xa[i] =i*1.0/n;
                ya[i]=exp(xa[i]) ;
            }
            yp1 = exp(xa[1]);
            ypn = exp(xa[n]);

        }
        else
        {
            break;
        }
        
        spline(xa,ya,n,yp1,ypn,y2);

        printf("\n%9s %13s %16s %13s\n", "x", "f(x)", "interpolated", "error");
        for (i=1;i<=10;i++)
        {
            if (nfunc == 1)
            {
                x=(-0.05+i/10.0)*PI;
                f=sin(x);
            }
            else if (nfunc == 2)
            {
                x=(-0.05+i/10.0);
                f=exp(x);
            }
            
            splint(xa,ya,y2,n,x,&y);
            printf("%12.6f %12.6f %12.6f \n", x,f,y);
        }

    printf("\n***********************************\n");
    printf("press RETURN\n");
    (void) getchar();

        

    }
free_vector(ya,1,n);
free_vector(xa,1,n);
return 0;

}

int testsplintrc(void)
{

    int i,n,nfunc, linit, l = 1;
    float f,x,y,yp1,ypn;
    float *xa,*ya,*y2;
    FILE *nf;
    double a,b;

    if((nf = fopen("In.dat", "r")) == NULL)
        nrerror("In.dat not found\n");
    while(!feof(nf))
    {
        fscanf(nf, "%lf%lf", &a, &b);
        l++;
    }
    l--;
    fclose(nf);

    xa=vector(1,l);
    ya=vector(1,l) ;
    y2=vector(1,l) ;


    if((nf = fopen("In.dat", "r")) == NULL)
        nrerror("In.dat not found\n");


    for (i=1; i<l; i++)
    {
        fscanf(nf, "%f%f", &xa[i], &ya[i]);
        printf("%d %12.6f %12.6f \n", i, xa[i],ya[i]);

    }
    printf("**********************************************\n");

    fclose(nf);

/*
    l=1;
    while(!feof(nf))
    {
        fscanf(nf, "%f%f", &xa[l], &ya[l]);
        printf("%d %12.6f %12.6f \n", l, xa[l],ya[l]);
        l++;

    }
    printf("**********************************************\n");
    l--;
    fclose(nf);
*/
    linit = xa[1];


    yp1 = 1e30;
    ypn = 1e30;
    
        
        spline(xa,ya,l,yp1,ypn,y2);


    nf = fopen("Is.dat", "w");

    for (i=1; i<l;i++)
    {

        splint(&xa[i], &ya[i],&y2[i],3,(float) i+linit,&y);

        //polint(&xa[i], &ya[i],3,(float) i+linit,&y,&dy);
        printf("\n%d %12.6f %12.6f %d %d %12.6f ", i, xa[i],ya[i],i, i+linit,y);
        //printf("\n%d %d %12.6f ", i, i+linit,y);
        //printf("\n**************");

        fprintf(nf, "%12.6lf %12.6lf %d  %12.6lf \n",xa[i],ya[i], i+linit, y);

    }

    fclose(nf);

free_vector(ya,1,n);
free_vector(xa,1,n);
return 0;

}