/*
 * 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.
 *
 * SEE ALSO: cheby_eval, cheby_integ, cheby_deriv
 */
{
  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.
 *
 * SEE ALSO: cheby_fit
 */
{
  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).
 *
 * SEE ALSO: cheby_fit, cheby_deriv
 */
{
  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.
 *
 * SEE ALSO: cheby_fit, cheby_integ
 */
{
  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;
}