110 Chapter 3.Interpolation and Extrapolation float *c,*d; dif=fabs(x-xa[1]); c=vector(1,n); d=vector(1,n); for(1=1;1<=n;1++){ Here we find the index ns of the closest table entry, if (dift=fabs(x-xa[i]))<dif){ ns=1; 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, 83g for(i=1;i<=n-m;1+)[ we loop over the current c's and d's and update ho=xa[i]-x; them. hp=xa[i+m]-x; w=c[i+1]-d[1]: if (den=ho-hp)==0.0)nrerror("Error in routine polint"); This error can occur only if two input xa's are (to within roundoff)identical. 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, America Press. 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) 9 Programs on the target x.The last dy added is thus the error indication. free_vector(d,1,n); SCIENTIFIC free_vector(c,1,n); to dir Quite often you will want to call polint with the dummy arguments xa and ya replaced by actual arrays with offsets.For example,the construction 192 COMPUTING (ISBN polint (&xx[14],&yy [14],4,x,y,dy)performs 4-point interpolation on the tab- ulated values xx[15..18],yy [15..18].For more on this,see the end of $3.4. Fuurggoglrion Numerical Recipes 10621 43106 CITED REFERENCES AND FURTHER READING: (outside Abramowitz,M.,and Stegun,I.A.1964,Handbook of Mathematical Functions,Applied Mathe- matics Series,Volume 55 (Washington:National Bureau of Standards;reprinted 1968 by North Software. Dover Publications,New York),$25.2. Stoer.J.,and Bulirsch,R.1980,Introduction to Numerical Analysis(New York:Springer-Verlag). 62.1. Gear,C.W.1971,Numerical Initial Value Problems in Ordinary Differential Equations(Englewood Cliffs,NJ:Prentice-Hall),$6.1.110 Chapter 3. Interpolation and Extrapolation Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). float *c,*d; dif=fabs(x-xa[1]); c=vector(1,n); d=vector(1,n); for (i=1;i<=n;i++) { Here we find 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 ho=xa[i]-x; them. hp=xa[i+m]-x; w=c[i+1]-d[i]; if ( (den=ho-hp) == 0.0) nrerror("Error in routine polint"); This error can occur only if two input xa’s are (to within roundoff) identical. 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); } Quite often you will want to call polint with the dummy arguments xa and ya replaced by actual arrays with offsets. For example, the construction polint(&xx[14],&yy[14],4,x,y,dy) performs 4-point interpolation on the tab￾ulated values xx[15..18], yy[15..18]. For more on this, see the end of §3.4. CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathe￾matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 by Dover Publications, New York), §25.2. Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag), §2.1. Gear, C.W. 1971, Numerical Initial Value Problems in Ordinary Differential Equations (Englewood Cliffs, NJ: Prentice-Hall), §6.1