;+
; NAME:
;   CHEBFIT
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Fit Chebyshev polynomial coefficients to a tabulated function
;
; MAJOR TOPICS:
;   Curve and Surface Fitting
;
; CALLING SEQUENCE:
;   p = CHEBFIT(X, Y, ERR, INTERVAL=interval, NMAX=nmax,
;               PRECISION=prec, /EVEN, /ODD, REDUCE_ALGORITHM=)
;
; DESCRIPTION:
;
;   CHEBFIT fits a series of Chebyshev polynomials to a set of
;   tabulated and possibly noisy data points.  The functions MPFIT and
;   CHEBEVAL, available from the above web page, must also be in your
;   IDL path for this function to work properly.  The user can choose
;   the desired precision and maximum number of chebyshev
;   coefficients.
;
;   This function is intended for use on already-tabulated data which
;   are potentially noisy.  The user should never expect more than
;   NPOINTS terms, where NPOINTS is the number of (x,y) pairs.  For
;   functions which can be evaluated to full machine precision at
;   arbitrary abcissae, the routine CHEBCOEF should be used instead.
;   For exact data tabulated on a regular grid, the routine CHEBGRID
;   should be tried.
;
;   The user can also specify that the function is even or odd, using
;   the keywords EVEN or ODD.  This saves computation time because
;   certain terms in the expansion can be ignored.  For the purposes
;   of this function even and odd refer to the symmetry about the
;   center of the interval.
;
;   The algorithm is employed in three steps.  In the first step, the
;   coefficients are estimated at a crude level.  In the second step,
;   it is determined whether certain coefficients are deemed
;   "ignoreable", i.e., they do not contribute significantly to the
;   function and are discarded.  The operation of this step is
;   determined by the REDUCE_ALGORITHM keyword.  Finally, the
;   remaining "good" coefficients are re-fitted to achieve the best
;   fit.
;
; INPUTS:
;
;   X, Y - the x- and y- tabulated values to be fitted.
;
;   ERR - (optional) the y-error bar associated with each (x,y) pair.
;         Default: 1
;
; RETURNS:
;
;   An array of Chebyshev coefficients which can be passed to
;   CHEBEVAL.  NOTE: the convention employed here is such that the
;   constant term in the expansion is P(0)*T0(x) (i.e., the convention
;   of Luke), and not P(0)/2 * T0(x).
;
; KEYWORD PARAMETERS:
;
;   EVEN, ODD - if set, then the fitting routine assumes the function
;               is even or odd, about the center of the interval.
;
;   INTERVAL - a 2-element vector describing the interval over which
;              the polynomial is to be evaluated.
;              Default: [-1, 1]
;
;   NMAX - a scalar, the maximum number of polynomial terms to be
;          fitted at one time.
;          Default: 16
;
;   PRECISION - a scalar, the requested precision in the fit.  Any
;               terms which do not contribute significantly, as
;               defined by this threshold, are discarded.  If the
;               function to be fitted is not well-behaved, then the
;               precision is not guaranteed to reach the desired
;               level.
;               Default: 1E-7
;
;   REDUCE_ALGORITHM - a scalar integer, describes how insignificant
;               terms are removed from the fit.  If 0, then all terms
;               are kept, and none are dicarded.  If 1, then only
;               trailing terms less than PRECISION are discarded.  If
;               2, then both trailing and intermediate terms less than
;               PRECISION are discarded.
;               Default: 2
;
; EXAMPLE:
;
;   x = dindgen(1000)/100     ; Range of 0 to 10
;   y = cos(x) + randomn(seed,1000)*0.01  ; Function with some noise
;   p = chebfit(x, y, interval=[0d,10d])
;   plot, x, y - chebeval(x,p, interval=[0d,10d])
;
; REFERENCES:
;
;   Abramowitz, M. & Stegun, I., 1965, *Handbook of Mathematical
;     Functions*, 1965, U.S. Government Printing Office, Washington,
;     D.C. (Applied Mathematical Series 55)
;   CERN, 1995, CERN Program Library, Function E407
;   Luke, Y. L., *The Special Functions and Their Approximations*,
;     1969, Academic Press, New York
;
; MODIFICATION HISTORY:
;   Written and documented, CM, June 2001
;   Copyright license terms changed, CM, 30 Dec 2001
;   Added usage message, CM, 20 Mar 2002
;   Slight docs change, CM, 25 Mar 2002
;
;  $Id: chebfit.pro,v 1.7 2003/07/20 05:53:44 craigm Exp $
;
;-
; Copyright (C) 2001, 2002, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-

;; Compute residuals for MPFIT
function chebfit_eval, p, interval=interval, nterms=nterms, igood=igood, $
                       _EXTRA=extra

  common chebfit_common, x, y, err

  if n_elements(igood) EQ 0 then begin
      p1 = p
  endif else begin
      p1 = replicate(p(0)*0, nterms)
      p1(igood) = p
  endelse

  ;; Compute the Chebyshev polynomial
  f = chebeval(x, p1, interval=interval)

  ;; Compute the deviates, applying either errors or weights
  if n_elements(err) GT 0 then begin
      result = (y-f)/err
  endif else if n_elements(wts) GT 0 then begin
      result = (y-f)*wts
  endif else begin
      result = (y-f)
  endelse
  
  ;; Make sure the returned result is one-dimensional.
  result = reform(result, n_elements(result), /overwrite)
  return, result
end

function chebfit, x, y, err, nmax=nterms0, interval=interval, $
                  precision=prec, even=even, odd=odd, quiet=quiet, $
                  initialize=init, reduce_algorithm=redalg0, $
                  indices=igood, nocatch=nocatch, $
                  yfit=yfit, perror=perror, bestnorm=bestnorm, dof=dof

  if n_params() EQ 0 then begin
      message, 'USAGE:', /info
      message, 'P = CHEBFIT(X, Y, ERR, INTERVAL=[a,b], NMAX=, ...)', /info
      return, !values.d_nan
  endif

  if n_elements(nterms0) EQ 0 then nterms = 16L $
  else                             nterms = floor(nterms0(0)) > 2L
  nterms = nterms < n_elements(x)
  if n_elements(interval) LT 2 then interval = [-1., 1.]
  if n_elements(prec) EQ 0 then prec = 1.e-7
  if n_elements(redalg0) EQ 0 then redalg = 2 else redalg = floor(redalg0(0))
  if n_elements(quiet) EQ 0 then quiet = 1

  ;; Handle error conditions gracefully
  if NOT keyword_set(nocatch) then begin
      catch, catcherror
      if catcherror NE 0 then begin
          catch, /cancel
          message, 'Error detected while fitting', /info
          message, !err_string, /info
          ier = -1L
          return, 0L
      endif
  endif

  if n_elements(p) LT nterms OR keyword_set(init) then begin
      p = replicate(x(0)*0 + 1, nterms) / (findgen(nterms)+1)^2
      p(0) = total(y)/n_elements(y)
      ;; If mean is *exactly* zero, then shift it off slightly
      if p(0) EQ 0 then p(0) = sqrt(total(y^2))/n_elements(y)/10
  endif
  p0 = p
  igood = lindgen(nterms)

  if keyword_set(even) OR keyword_set(odd) then $
    igood = lindgen(n_elements(p)/2)*2 + keyword_set(odd)
  nt = min([nterms, max(igood)+1])

  ;; Cancel out old common entries
  common chebfit_common, xc, yc, errc
  xc   = 0 & dummy = temporary(xc)
  yc   = 0 & dummy = temporary(yc)
  errc = 0 & dummy = temporary(errc)
  
  xc = x
  yc = y
  if n_elements(err) GT 0 then begin
      errc = err
  endif

  fa = {interval: interval, igood: igood, nterms: nt}
  p1 = mpfit('CHEBFIT_EVAL', p0(igood), functargs=fa, maxiter=5, quiet=quiet)
  p0(igood) = p1

  ;; Look for and remove the insignificant terms from the fit
  if redalg GT 0 then begin
      wh = where(abs(p1) GT prec(0), ct)
      if ct EQ 0 then begin
          ALL_ZERO:
          message, 'WARNING: no significant Chebyshev terms were detected', $
            /info
          p = p0*0
          return, 0L
      endif
      if max(wh) LT n_elements(igood)-1 then begin
          imax = max(wh)
          igood = igood(0:imax)
          p1 = p1(0:imax)
      endif

      if redalg EQ 2 then begin
          wh = where(abs(p1) GT 0.1*prec, ct)
          if ct EQ 0 then goto, ALL_ZERO
          igood = igood(wh)
          p1 = p1(wh)
      endif
  endif

  nt = min([nterms, max(igood)+1])
  fa = {interval: interval, igood: igood, nterms: nt}
  p2 = mpfit('CHEBFIT_EVAL', p1, functargs=fa, maxiter=10, quiet=quiet, $
             perror=dp2, bestnorm=bestnorm, dof=dof)

  xc = 0 & yc = 0 & errc = 0
  p = p0*0 
  perror = p
  p(igood) = p2
  perror(igood) = dp2

  if arg_present(yfit) then $
    yfit = chebeval(x, p, interval=interval)

  return, p
end