128 Chapter 3.Interpolation and Extrapolation evaluation for the final column spline.(Recall that a spline construction is a process of order N,while a spline evaluation is only of order log N-and that is just to find the place in the table!) Here is a routine to precompute the auxiliary second-derivative table: void splie2(float xian],float x2an],float *ya,int m,int n,float **y2a) Given an m by n tabulated function ya[1..m][1..n],and tabulated independent variables x2a[1..n],this routine constructs one-dimensional natural cubic splines of the rows of ya and returns the second-derivatives in the array y2a[1..m][1..n].(The array x1a[1..m]is included in the argument list merely for consistency with routine splin2. 5常 void spline(float x[,float y,int n,float ypl,float ypn,float y20); 81 int j; for (j=1;j<=m;j++) sp1ine(x2a,ya[j],n,1.0e30,1.0e30,y2a[j]) Values 1x1030 signal a nat- ural spline. (If you want to interpolate on a sub-block of a bigger matrix,see $1.2.) RECIPES After the above routine has been executed once,any number of bicubic spline interpolations can be performed by successive calls of the following routine: 令 #include "nrutil.h" Press. void splin2(float xia[],float x2a[],float **ya,float **y2a,int m,int n, float x1,float x2,float *y) Given x1a,x2a,ya,m,n as described in splie2 and y2a as produced by that routine;and given a desired interpolating point x1,x2;this routine returns an interpolated function value y by bicubic spline interpolation. IENTIFIC void spline(float x[],float y[],int n,float yp1,float ypn,float y20); void splint(float xa[],float yan],float y2a,int n,float x,float ty); 6 int j; float *ytmp,*yytmp; ytmp=vector(1,m); yytmp=vector(1,m); Perform m evaluations of the row splines constructed by for (j=1;j<=m;j++) splie2,using the one-dimensional spline evaluator splint(x2a,ya[j],y2a[j],n,x2,&yytmp[j]); splint. spline(x1a,yytmp,m,1.0e30,1.0e30,ytmp); Construct the one-dimensional col- 10621 splint(xia,yytmp,ytmp,m,x1,y); umn spline and evaluate it. free_vector(yytmp,1,m); uction, Numerical Recipes 43108 free_vector(ytmp,1,m); (outside North Software. 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),825.2. Kinahan,B.F.,and Harm,R.1975,Astrophysica/Journal,vol.200.pp.330-335. Johnson,L.W.,and Riess,R.D.1982,Numerica/Analysis,2nd ed.(Reading,MA:Addison- Wesley).85.2.7. Dahlquist,G.,and Bjorck,A.1974,Numerica/Methods (Englewood Cliffs,NJ:Prentice-Hall). $7.7.128 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). evaluation for the final column spline. (Recall that a spline construction is a process of order N, while a spline evaluation is only of order log N — and that is just to find the place in the table!) Here is a routine to precompute the auxiliary second-derivative table: void splie2(float x1a[], float x2a[], float **ya, int m, int n, float **y2a) Given an m by n tabulated function ya[1..m][1..n], and tabulated independent variables x2a[1..n], this routine constructs one-dimensional natural cubic splines of the rows of ya and returns the second-derivatives in the array y2a[1..m][1..n]. (The array x1a[1..m] is included in the argument list merely for consistency with routine splin2.) { void spline(float x[], float y[], int n, float yp1, float ypn, float y2[]); int j; for (j=1;j<=m;j++) spline(x2a,ya[j],n,1.0e30,1.0e30,y2a[j]); Values 1×1030 signal a nat- } ural spline. (If you want to interpolate on a sub-block of a bigger matrix, see §1.2.) After the above routine has been executed once, any number of bicubic spline interpolations can be performed by successive calls of the following routine: #include "nrutil.h" void splin2(float x1a[], float x2a[], float **ya, float **y2a, int m, int n, float x1, float x2, float *y) Given x1a, x2a, ya, m, n as described in splie2 and y2a as produced by that routine; and given a desired interpolating point x1,x2; this routine returns an interpolated function value y by bicubic spline interpolation. { 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 j; float *ytmp,*yytmp; ytmp=vector(1,m); yytmp=vector(1,m); Perform m evaluations of the row splines constructed by splie2, using the one-dimensional spline evaluator splint. for (j=1;j<=m;j++) splint(x2a,ya[j],y2a[j],n,x2,&yytmp[j]); spline(x1a,yytmp,m,1.0e30,1.0e30,ytmp); Construct the one-dimensional col￾splint(x1a,yytmp,ytmp,m,x1,y); umn spline and evaluate it. free_vector(yytmp,1,m); free_vector(ytmp,1,m); } 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. Kinahan, B.F., and Harm, R. 1975, Astrophysical Journal, vol. 200, pp. 330–335. Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis, 2nd ed. (Reading, MA: Addison￾Wesley), §5.2.7. Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), §7.7