Interpolation.c File Reference

#include <stdio.h>
#include <conio.h>
#include <math.h>
#include "Marq.h"
#include "nrutil.h"

Include dependency graph for Interpolation.c:

Go to the source code of this file.

Macros
#define	FREERETURN   {free_vector(d,1,n);free_vector(c,1,n);return;}
#define	PI   3.1415926
#define	TINY   1.0e-25

Functions
void	polint (float xa[], float ya[], int n, float x, float *y, float *dy)
void	ratint (float xa[], float ya[], int n, float x, float *y, float *dy)
void	spline (float x[], float y[], int n, float yp1, float ypn, float y2[])
void	splint (float xa[], float ya[], float y2a[], int n, float x, float *y)
int	testpolintfull (void)
int	testpolintlambda (void)
int	testpolintshort (void)
int	testsplint (void)
int	testsplintrc (void)

Macro Definition Documentation

FREERETURN

#define FREERETURN   {free_vector(d,1,n);free_vector(c,1,n);return;}

Definition at line 14 of file Interpolation.c.

Referenced by ratint().

PI

#define PI   3.1415926

Definition at line 15 of file Interpolation.c.

Referenced by testpolintfull(), testpolintshort(), and testsplint().

TINY

#define TINY   1.0e-25

Definition at line 13 of file Interpolation.c.

Referenced by ratint().

Function Documentation

polint()

void polint	(	float	xa[],
		float	ya[],
		int	n,
		float	x,
		float *	y,
		float *	dy
	)

Definition at line 28 of file Interpolation.c.

References free_vector(), nrerror(), and vector().

Referenced by testpolintfull(), testpolintlambda(), and testpolintshort().

{
    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);
}

        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);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ratint()

void ratint	(	float	xa[],
		float	ya[],
		int	n,
		float	x,
		float *	y,
		float *	dy
	)

Definition at line 82 of file Interpolation.c.

References FREERETURN, nrerror(), TINY, and vector().

{
    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
}

Here is the call graph for this function:

spline()

void spline	(	float	x[],
		float	y[],
		int	n,
		float	yp1,
		float	ypn,
		float	y2[]
	)

Definition at line 138 of file Interpolation.c.

References free_vector(), p, and vector().

Referenced by testsplint(), and testsplintrc().

{
    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);
}

Here is the call graph for this function:

Here is the caller graph for this function:

splint()

void splint	(	float	xa[],
		float	ya[],
		float	y2a[],
		int	n,
		float	x,
		float *	y
	)

Definition at line 183 of file Interpolation.c.

References nrerror().

Referenced by testsplint(), and testsplintrc().

{
    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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

testpolintfull()

int testpolintfull

(

void

)

Definition at line 308 of file Interpolation.c.

References free_vector(), PI, polint(), and vector().

{

    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;

}

Here is the call graph for this function:

testpolintlambda()

int testpolintlambda

(

void

)

Definition at line 379 of file Interpolation.c.

References dvector(), nrerror(), and polint().

{

    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);

}

Here is the call graph for this function:

testpolintshort()

int testpolintshort

(

void

)

Definition at line 214 of file Interpolation.c.

References free_vector(), PI, polint(), and vector().

{

    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;

}

Here is the call graph for this function:

testsplint()

int testsplint

(

void

)

Definition at line 419 of file Interpolation.c.

References free_vector(), PI, spline(), splint(), and vector().

{

    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;

}

Here is the call graph for this function:

testsplintrc()

int testsplintrc

(

void

)

Definition at line 497 of file Interpolation.c.

References free_vector(), nrerror(), spline(), splint(), and vector().

{

    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;

}

Here is the call graph for this function: