GNU-Darwin Web
i/fitrat
 Home Manual Packages Global Index Keywords Quick Reference ``` /* FITRAT.I Polynomial and rational function interpolators. \$Id\$ */ /* Copyright (c) 1994. The Regents of the University of California. All rights reserved. */ func fitpol (y, x, xp, keep=) /* DOCUMENT yp= fitpol(y, x, xp) -or- yp= fitpol(y, x, xp, keep=1) is an interpolation routine similar to interp, except that fitpol returns the polynomial of degree numberof(X)-1 which passes through the given points (X,Y), evaluated at the requested points XP. The X must either increase or decrease monotonically. If the KEEP keyword is present and non-zero, the external variable yperr will contain a list of error estimates for the returned values yp on exit. The algorithm is taken from Numerical Recipes (Press, et. al., Cambridge University Press, 1988); it is called Neville's algorithm. The rational function interpolator fitrat is better for "typical" functions. The Yorick implementaion requires numberof(X)*numberof(XP) temporary arrays, so the X and Y arrays should be reasonably small. SEE ALSO: fitrat, interp */ { /* find lower and upper indices in x of the interval containing xp */ nx= numberof(x); u= digitize(xp, x); l= max(u-1, 1); u= min(u, nx); /* initialize yp to the y(x) closest to xp */ mask= (abs(x(l)-xp)<=abs(x(u)-xp)); ns= l*mask + u*(!mask); yp= y(ns); ns--; c= d= y= array(double(y), dimsof(xp)); dx= x - xp(-,); if (dimsof(ns)(1)) { is= ns; is(*)= nx*indgen(0:numberof(xp)-1); } else { is= 0; } for (m=1 ; m=(nx-m)); y= d(max(ns,1)+is)*mask + c(ns+1+is)*(!mask); /* this is error estimate */ yp+= y; ns-= mask; } if (keep) { extern yperr; yperr= y; } return yp; } func fitrat (y, x, xp, keep=) /* DOCUMENT yp= fitrat(y, x, xp) -or- yp= fitrat(y, x, xp, keep=1) is an interpolation routine similar to interp, except that fitpol returns the diagonal rational function of degree numberof(X)-1 which passes through the given points (X,Y), evaluated at the requested points XP. (The numerator and denominator polynomials have equal degree, or the denominator has one larger degree.) The X must either increase or decrease monotonically. Also, this algorithm works by recursion, fitting successively to consecutive pairs of points, then consecutive triples, and so forth. If there is a pole in any of these fits to subsets, the algorithm fails even though the rational function for the final fit is non- singular. In particular, if any of the Y values is zero, the algorithm fails, and you should be very wary of lists for which Y changes sign. Note that if numberof(X) is even, the rational function is Y-translation invariant, while numberof(X) odd generally results in a non-translatable fit (because it decays to y=0). If the KEEP keyword is present and non-zero, the external variable yperr will contain a list of error estimates for the returned values yp on exit. The algorithm is taken from Numerical Recipes (Press, et. al., Cambridge University Press, 1988); it is called the Bulirsch-Stoer algorithm. The Yorick implementaion requires numberof(X)*numberof(XP) temporary arrays, so the X and Y arrays should be reasonably small. SEE ALSO: fitpol, interp */ { /* find lower and upper indices in x of the interval containing xp */ nx= numberof(x); u= digitize(xp, x); l= max(u-1, 1); u= min(u, nx); /* initialize yp to the y(x) closest to xp */ mask= (abs(x(l)-xp)<=abs(x(u)-xp)); ns= l*mask + u*(!mask); yp= y(ns); exact= (x(ns)==xp); list= where(exact); if (numberof(list)) { yexact= yp(list); yerr= 0.0*yexact; } list= where(!exact); if (numberof(list)) { yp= yp(list); xp= xp(list); ns= ns(list); c= d= y= array(double(y), numberof(xp)); d+= fitrat_tiny; /* I question whether this ever rescues a failed fit, but I'll leave it here anyway. */ dx= x - xp(-,); ns--; is= nx*indgen(0:numberof(xp)-1); for (m=1 ; m=(nx-m)); y= d(max(ns,1)+is)*mask + c(ns+1+is)*(!mask); yp+= y; ns-= mask; } } else { yp= []; } if (keep) { extern yperr; yperr= merge(yerr, y, exact); } return merge(yexact, yp, exact); } fitrat_tiny= 1.e-30; ```