;+ ; NAME: ; CHEBEVAL ; ; 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: ; Evaluate a Chebyshev polynomial on an interval, given the coefficients ; ; MAJOR TOPICS: ; Curve and Surface Fitting ; ; CALLING SEQUENCE: ; y = CHEBEVAL(X, P, INTERVAL=interval, DERIVATIVE=deriv) ; ; DESCRIPTION: ; ; CHEBEVAL computes the values of a Chebyshev polynomial function at ; specified abcissae, over the interval [a,b]. The user must supply ; the abcissae and the polynomial coefficients. The function is of ; the form: ; ; N ; y(x) = Sum p_n T_n(x*) x in [a,b] ; i=0 ; ; Where T_n(x*) are the orthogonal Chebyshev polynomials of the ; first kind, defined on the interval [-1,1] and p_n are the ; coefficients. The scaled variable x* is defined on the [-1,1] ; interval such that (x*) = (2*x - a - b)/(b - a), and x is defined ; on the [a,b] interval. ; ; The derivative of the function may be computed simultaneously ; using the DERIVATIVE keyword. ; ; The is some ambiguity about the definition of the first ; coefficient, p_0, namely, the use of p_0 vs. the use of p_0/2. ; The p_0 definition of Luke is used in this function. ; ; INPUTS: ; ; X - a numerical scalar or vector, the abcissae at which to ; evaluate the polynomial. If INTERVAL is specified, then all ; values of X must lie within the interval. ; ; P - a vector, the Chebyshev polynomial coefficients, as returned ; by CHEBFIT or CHEBCOEF. ; ; RETURNS: ; ; An array of function values, evaluated at the abcissae. The ; numeric precision is the greater of X or P. ; ; KEYWORD PARAMETERS: ; ; DERIVATIVE - upon return, a vector containing the derivative of ; the function at each abcissa is returned in this ; keyword. ; ; INTERVAL - a 2-element vector describing the interval over which ; the polynomial is to be evaluated. ; Default: [-1, 1] ; ; EXAMPLE: ; ; x = dindgen(1000)/100 ; Range of 0 to 10 ; p = chebcoef('COS(x)', /expr, interval=[0d, 10d]) ;; Compute coefs ; y = chebeval(x, p, interval=[0d,10d]) ;; Eval Cheby poly ; plot, x, y - cos(x) ; Plot residuals ; ; 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 ; Return a vector even when P has one element, CM, 22 Nov 2004 ; Fix bug in evaluation of derivatives, CM, 22 Nov 2004 ; ; $Id: chebeval.pro,v 1.6 2004/11/22 07:08:00 craigm Exp $ ; ;- ; Copyright (C) 2001, 2002, 2004, 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. ;- function chebeval, x0, p, derivative=v, interval=interval0 if n_params() EQ 0 then begin message, 'USAGE:', /info message, 'F = CHEBEVAL(X, P, INTERVAL=[a,b], DERIVATIVE=DFDX)', /info return, !values.d_nan endif zero = x0(0)*0 + 0. & one = zero + 1 y = zero v = zero if n_elements(interval0) LT 2 then begin t = x0 a = -one & b = +one endif else begin a = interval0(0)+0. & b = interval0(1)+0. if a EQ b then return, t*0 t = (2*x0 - a - b)/(b-a) endelse p0 = one + t*0 p1 = t v0 = zero + t*0 v1 = one + t*0 v2 = 4.*t t = 2.*t n = n_elements(p) if arg_present(v) then begin for i = 0L, n-1, 2 do begin if i EQ n-1 then begin p1 = zero v1 = zero endif j = (i+1)<(n-1) y = temporary(y) + p(i)*p0 + p(j)*p1 v = temporary(v) + p(i)*v0 + p(j)*v1 ;; Advance to the next set of Chebyshev polynomials. For ;; velocity we need to keep the next orders around ;; momentarily. p2 = t*p1 - p0 p3 = t*p2 - p1 v2 = t*v1 - v0 + 2*p1 v3 = t*v2 - v1 + 2*p2 p0 = temporary(p2) & p1 = temporary(p3) v0 = temporary(v2) & v1 = temporary(v3) endfor endif else begin for i = 0L, n-1, 2 do begin if i EQ n-1 then p1 = zero j = (i+1)<(n-1) y = temporary(y) + p(i)*p0 + p(j)*p1 ;; Advance to the next set of Chebyshev polynomials. For no ;; derivative, we can re-use old variables. p0 = t*p1 - temporary(p0) p1 = t*p0 - temporary(p1) endfor endelse v = temporary(v)*2/(b-a) return, y end