```/*
FITRAT.I
Polynomial and rational function interpolators.

\$Id\$
*/
/*    Copyright (c) 1994.  The Regents of the University of California.

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.

*/
{
/* 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 */
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++) {
ho= dx(1:nx-m,..);
hp= dx(1+m:nx,..);
den= (c(2:nx-m+1,..)-d(1:nx-m,..))/(ho-hp);
d(1:nx-m,..)= hp*den;
c(1:nx-m,..)= ho*den;
yp+= y;
}

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.

*/
{
/* 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 */
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++) {
di= d(1:nx-m,);
ci= c(2:nx-m+1,);
t= dx(1:nx-m,)*di/dx(1+m:nx,);
dd= (ci-di)/(t-ci);
d(1:nx-m,)= ci*dd;
c(1:nx-m,)= t*dd;
yp+= y;
}

} else {
yp= [];
}

if (keep) {
extern yperr;
yperr= merge(yerr, y, exact);
}

return merge(yexact, yp, exact);
}

fitrat_tiny= 1.e-30;
```