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