1 /*
2  * $Id: fitrat.i,v 1.2 2010-04-18 10:33:38 thiebaut Exp $
3  * Polynomial and rational function interpolators.
4  */
5 /* Copyright (c) 2005, The Regents of the University of California.
6  * All rights reserved.
7  * This file is part of yorick (http://yorick.sourceforge.net).
8  * Read the accompanying LICENSE file for details.
9  */
10 
11 func fitpol(y, x, xp, keep=)
12 /* DOCUMENT yp= fitpol(y, x, xp)
13          or yp= fitpol(y, x, xp, keep=1)
14      is an interpolation routine similar to interp, except that fitpol
15      returns the polynomial of degree numberof(X)-1 which passes through
16      the given points (X,Y), evaluated at the requested points XP.
17 
18      The X must either increase or decrease monotonically.
19 
20      If the KEEP keyword is present and non-zero, the external variable
21      yperr will contain a list of error estimates for the returned values
22      yp on exit.
23 
24      The algorithm is taken from Numerical Recipes (Press, et. al.,
25      Cambridge University Press, 1988); it is called Neville's algorithm.
26      The rational function interpolator fitrat is better for "typical"
27      functions.  The Yorick implementaion requires numberof(X)*numberof(XP)
28      temporary arrays, so the X and Y arrays should be reasonably small.
29 
30    SEE ALSO: fitrat, interp
31  */
32 {
33   /* find lower and upper indices in x of the interval containing xp */
34   nx= numberof(x);
35   u= digitize(xp, x);
36   l= max(u-1, 1);
37   u= min(u, nx);
38 
39   /* initialize yp to the y(x) closest to xp */
40   mask= (abs(x(l)-xp)<=abs(x(u)-xp));
41   ns= l*mask + u*(!mask);
42   yp= y(ns);
43   ns--;
44 
45   c= d= y= array(double(y), dimsof(xp));
46   dx= x - xp(-,);
47   if (dimsof(ns)(1)) {
48     is= ns;
49     is(*)= nx*indgen(0:numberof(xp)-1);
50   } else {
51     is= 0;
52   }
53   for (m=1 ; m<nx ; m++) {
54     ho= dx(1:nx-m,..);
55     hp= dx(1+m:nx,..);
56     den= (c(2:nx-m+1,..)-d(1:nx-m,..))/(ho-hp);
57     d(1:nx-m,..)= hp*den;
58     c(1:nx-m,..)= ho*den;
59     mask= (2*ns>=(nx-m));
60     y= d(max(ns,1)+is)*mask + c(ns+1+is)*(!mask); /* this is error estimate */
61     yp+= y;
62     ns-= mask;
63   }
64 
65   if (keep) {
66     extern yperr;
67     yperr= y;
68   }
69 
70   return yp;
71 }
72 
73 func fitrat(y, x, xp, keep=)
74 /* DOCUMENT yp= fitrat(y, x, xp)
75          or yp= fitrat(y, x, xp, keep=1)
76      is an interpolation routine similar to interp, except that fitpol
77      returns the diagonal rational function of degree numberof(X)-1 which
78      passes through the given points (X,Y), evaluated at the requested
79      points XP.  (The numerator and denominator polynomials have equal
80      degree, or the denominator has one larger degree.)
81 
82      The X must either increase or decrease monotonically.  Also, this
83      algorithm works by recursion, fitting successively to consecutive
84      pairs of points, then consecutive triples, and so forth.
85      If there is a pole in any of these fits to subsets, the algorithm
86      fails even though the rational function for the final fit is non-
87      singular.  In particular, if any of the Y values is zero, the
88      algorithm fails, and you should be very wary of lists for which
89      Y changes sign.  Note that if numberof(X) is even, the rational
90      function is Y-translation invariant, while numberof(X) odd generally
91      results in a non-translatable fit (because it decays to y=0).
92 
93      If the KEEP keyword is present and non-zero, the external variable
94      yperr will contain a list of error estimates for the returned values
95      yp on exit.
96 
97      The algorithm is taken from Numerical Recipes (Press, et. al.,
98      Cambridge University Press, 1988); it is called the Bulirsch-Stoer
99      algorithm.  The Yorick implementaion requires numberof(X)*numberof(XP)
100      temporary arrays, so the X and Y arrays should be reasonably small.
101 
102    SEE ALSO: fitpol, interp
103  */
104 {
105   /* find lower and upper indices in x of the interval containing xp */
106   nx= numberof(x);
107   u= digitize(xp, x);
108   l= max(u-1, 1);
109   u= min(u, nx);
110 
111   /* initialize yp to the y(x) closest to xp */
112   mask= (abs(x(l)-xp)<=abs(x(u)-xp));
113   ns= l*mask + u*(!mask);
114   yp= y(ns);
115 
116   exact= (x(ns)==xp);
117   list= where(exact);
118   if (numberof(list)) {
119     yexact= yp(list);
120     yerr= 0.0*yexact;
121   }
122 
123   list= where(!exact);
124   if (numberof(list)) {
125     yp= yp(list);
126     xp= xp(list);
127     ns= ns(list);
128 
129     c= d= y= array(double(y), numberof(xp));
130     d+= fitrat_tiny;  /* I question whether this ever rescues a
131                          failed fit, but I'll leave it here anyway.  */
132     dx= x - xp(-,);
133     ns--;
134     is= nx*indgen(0:numberof(xp)-1);
135     for (m=1 ; m<nx ; m++) {
136       di= d(1:nx-m,);
137       ci= c(2:nx-m+1,);
138       t= dx(1:nx-m,)*di/dx(1+m:nx,);
139       dd= (ci-di)/(t-ci);
140       d(1:nx-m,)= ci*dd;
141       c(1:nx-m,)= t*dd;
142       mask= (2*ns>=(nx-m));
143       y= d(max(ns,1)+is)*mask + c(ns+1+is)*(!mask);
144       yp+= y;
145       ns-= mask;
146     }
147 
148   } else {
149     yp= [];
150   }
151 
152   if (keep) {
153     extern yperr;
154     yperr= merge(yerr, y, exact);
155   }
156 
157   return merge(yexact, yp, exact);
158 }
159 
160 fitrat_tiny= 1.e-30;
161