```/*
* cheby.i
* Chebyshev polynomial approximation routines
* after Numerical Recipes (Press et. al.) section 5.6
*/

func cheby_fit (f, x, n)
/* DOCUMENT fit = cheby_fit(f, interval, n)
*       or fit = cheby_fit(f, x, n)
*   returns the Chebyshev fit (for use in cheby_eval) of degree N
*   to the function F on the INTERVAL (a 2 element array [a,b]).
*   In the second form, F and X are arrays; the function to be
*   fit is the piecewise linear function of xp interp(f,x,xp), and
*   the interval of the fit is [min(x),max(x)].
*
*   The return value is the array [a,b, c0,c1,c2,...cN] where [a,b]
*   is the interval over which the fit applies, and the ci are the
*   Chebyshev coefficients.  It may be useful to use a relatively
*   large value of N in the call to cheby_fit, then to truncate the
*   resulting fit to fit(1:3+m) before calling cheby_eval.
*
*/
{
a = double(min(x));
b = max(x);
++n;
p = (pi/n) * span(0.5,n-0.5,n);
c = cos(p*indgen(0:n-1)(-,));
p = a + 0.5*(b-a)*(c(,2)+1.);
if (is_array(f)) p = interp(f,x, p);
else for (i=1 ; i<=n ; ++i) p(i) = f(p(i));
return grow([a,b], (2./n) * (p(+)*c(+,)));
}

func cheby_eval (fit, x)
/* DOCUMENT cheby_eval(fit, x)
*   evaluates the Chebyshev fit (from cheby_fit) at points X.
*   the return values have the same dimensions as X.
*
*/
{
x = interp([-2.,2.],fit(1:2), x);
a = b = 0.;
for (i=numberof(fit) ; i>2 ; --i) {
c = b;
b = a;
a = x*b - c + fit(i);
}
return 0.5*(a-c);
}

func cheby_integ (fit, x0)
/* DOCUMENT cheby_integ(fit)
*       or cheby_integ(fit, x0)
*   returns Chebyshev fit to the integral of the function of the
*   input Chebyshev FIT.  If X0 is given, the returned integral will
*   be zero at X0 (which should be inside the fit interval fit(1:2)),
*   otherwise the integral will be zero at x=fit(1).
*
*/
{
if (is_void(x0)) x0 = fit(1);
f = fit;
c = 0.25*(fit(2)-fit(1));
n = numberof(fit) - 2;
if (n>2) f(4:n+1) = c * (fit(3:n)-fit(5:n+2))/indgen(n-2);
f(0) = c * fit(n+1)/(n-1);
f(3) = 0.;
f(3) = -2.*cheby_eval(f, x0);
return f;
}

func cheby_deriv (fit)
/* DOCUMENT cheby_deriv(fit)
*   returns Chebyshev fit to the derivative of the function of the
*   input Chebyshev FIT.
*
*/
{
n = numberof(fit) - 2;
if (n<2) return [fit(1),fit(2),0.];
f = fit(1:-1);
f(0) = 2.*(n-1)*fit(0);
if (n>2) f(-1) = 2.*(n-2)*fit(-1);
for (i=-2 ; i>1-n ; --i) f(i) = f(i+2) + 2.*(i+n-1)*fit(i);
return (2./(fit(2)-fit(1))) * f;
}
```