;+ ; NAME: ; CHEBCOEF ; ; 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: ; Estimate Chebyshev polynomial coefficients of a function on an interval ; ; MAJOR TOPICS: ; Curve and Surface Fitting ; ; CALLING SEQUENCE: ; p = CHEBCOEF(FUNC, PRIVATE, FUNCTARGS=functargs, /DOUBLE, /EXPRESSION, $ ; PRECISION=prec, ERROR=err, NMAX=nmax, INTERVAL=interval, $ ; REDUCE_ALGORITHM=, STATUS=) ; ; DESCRIPTION: ; ; CHEBCOEF estimates the coefficients for a finite sum of Chebyshev ; polynomials approximating the function FUNC(x) over an interval. ; The user can choose the desired precision and maximum number of ; chebyshev coefficients. ; ; This routine is intended for functions which can be evaluated to ; full machine precision at arbitrary abcissae, and which are smooth ; enough to ensure that the coefficients are a decreasing sequence. ; For already-tabulated or potentially noisy data, the routines ; CHEBGRID or CHEBFIT should be used instead. ; ; The function to be approximated may either be the name of an IDL ; function (the default behavior), or an IDL expression (using the ; /EXPRESSION keyword). ; ; The procedure uses a modified form of the classic algorithm for ; determining the coefficients, which relies the orthogonality ; relation for Chebyshev polynomials. The interval [a,b] is ; subdivided successively into sets of subintervals of length ; 2^(-k)*(b-a),(k = 0,1,2...). After each subdivision the ; orthogonality properties of the Chebyshev polynomials with respect ; to summation over equally-spaced points are used to compute two ; sets of approximate values of the coefficients cj, one set ; computed using the end-points of the subintervals, and one set ; using the mid-points. Certain convergence requirements must be ; met before terminating. If the routine fails to converge with 64 ; coefficents, then the current best-fitting coefficients are ; returned, along with an error estimate in the ERROR keyword. ; CHEBCOEF never returns more than 64 coefficients. ; ; The coefficients may be further refined. If the keyword ; REDUCE_ALGORITHM is set to a value of 1, then any high order ; coefficients below a certain threshold are discarded. If ; REDUCE_ALGORITHM is set to 2 (the default), then all coefficients ; below the threshold are discarded rather than just the high order ; ones. The threshold is determined by the PRECISION keyword. ; ; INPUTS: ; ; FUNC - a scalar string, the name of the function to be ; approximated, or an IDL string containing an expression to ; be approximated (if /EXPRESSION is set). ; ; PRIVATE - any optional variable to be passed on to the function to ; be integrated. For functions, PRIVATE is passed as the ; second positional parameter; for expressions, PRIVATE can ; be referenced by the variable 'P'. CHEBCOEF does not ; examine or alter PRIVATE. ; ; 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: ; ; DOUBLE - if set, then computations are done in double precision ; rather than single precision. ; ; ERROR - upon return, this keyword contains an estimate of the ; maximum absolute error in the approximation. ; ; EXPRESSION - if set, then FUNC is an IDL expression to be ; approximated, rather than the name of a function. ; ; FUNCTARGS - A structure which contains the parameters to be passed ; to the user-supplied function specified by FUNCT via ; the _EXTRA mechanism. This is the way you can pass ; additional data to your user-supplied function without ; using common blocks. By default, no extra parameters ; are passed to the user-supplied function. ; ; 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 coefficients to be ; estimated. This number may not exceed 64. ; Default: 64 ; ; PRECISION - a scalar, the requested precision in the ; approximation. Any terms which do not contribute ; significantly, as defined by this threshold, are ; discarded. If the function to be estimated 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 ; ; STATUS - upon return, this keyword contains information about the ; status of the approximation. A value of -1 indicates bad ; input values; a value of 0 indicates the required ; accuracy was not obtained; a value of 1 indicates ; success. ; ; 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 E406 ; 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 ; Changed docs slightly, CM, 25 Mar 2002 ; ; $Id: chebcoef.pro,v 1.6 2002/05/03 18:40:27 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. ;- ;; Evaluate a user-supplied expression function chebcoef_eval, x, p, expression=expr, _EXTRA=extra y = 0 cmd = 'Y = '+expr dummy = execute(cmd) return, y end function chebcoef, f0, priv, functargs=fa, double=double, error=err, $ nmax=nmax, interval=interval, precision=prec0, $ expression=expr, reduce_algorithm=redalg0, $ status=status, indices=igood if n_params() EQ 0 then begin message, 'USAGE:', /info message, 'P = CHEBCOEF(FUNCT, [PRIV,] INTERVAL=[a,b], NMAX=...)', /info return, !values.d_nan endif sz = size(f0) err = -1 if sz(sz(0)+1) NE 7 OR n_elements(f0) NE 1 then begin NO_FUNCT: message, 'ERROR: FUNCT must be a scalar string', /info return, 0 endif ;; Check for empty string f = strtrim(f0(0),2) if f EQ '' then goto, NO_FUNCT ;; Prepare for EXPRESSION if requested if keyword_set(expr) then begin f = 'CHEBCOEF_EVAL' fa = {expression: strtrim(f0(0),2)} endif else begin f = strtrim(f0(0),2) endelse ;; 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 approximating '+f, /info message, !err_string, /info errmsg = 0 if NOT keyword_set(expr) then begin f1 = byte(strupcase(strtrim(f0(0),2))) ca = (byte('A'))(0) cz = (byte('Z'))(0) c0 = (byte('0'))(0) c9 = (byte('9'))(0) c_ = (byte('_'))(0) wh = where((f1 GE ca AND f1 LE cz) EQ 0 AND f1 NE c_ $ AND (f1 GE c0 AND f1 LE c9) EQ 0, ct) if ct GT 0 OR (f1(0) GE c0 AND f1(0) LE c9) then begin message, ('FUNCT appears to be an expression. Did you '+$ 'intend to pass the /EXPRESSION keyword?'), /info errmsg = 1 endif endif if errmsg EQ 0 then $ message, ('Please verify that function works and conforms to '+$ 'the documentation'), /info ier = -1L return, 0L endif endif if n_elements(prec0) EQ 0 then prec = 1e-7 else prec = prec0(0) zero = prec*0. if keyword_set(double) then zero = 0D if n_elements(interval) LT 2 then interval = zero + [-1., 1.] if n_elements(redalg0) EQ 0 then redalg = 2 else redalg = floor(redalg0(0)) status = -1 a = interval(0) b = interval(1) hf = zero + 0.5 eps = prec z1 = zero + 1 z2 = zero + 2 sz = size(zero) if sz(sz(0)+1) EQ 5 then pi = !dpi else pi = !pi x0 = [a, b] if n_elements(priv) GT 0 then begin if n_elements(fa) GT 0 then fv = call_function(f, x0, priv, _EXTRA=fa) $ else fv = call_function(f, x0, priv) endif else begin if n_elements(fa) GT 0 then fv = call_function(f, x0, _EXTRA=fa) $ else fv = call_function(f, x0) endelse ALFA=HF*(B-A) BETA=HF*(B+A) C1=fv(0) C2=fv(1) AC = [C2+C1, C2-C1] BC = AC*0 for i = 1, 7 do begin I1=2^(I-1) I2=I1-1 I3=2*I1 C1=Z2/I1 C2=PI/I1 jj = dindgen(i2+1) x = alfa*cos((jj+hf)*c2)+beta if n_elements(priv) GT 0 then begin if n_elements(fa) GT 0 then fv = call_function(f, x, priv, $ _EXTRA=fa) $ else fv = call_function(f, x, priv) endif else begin if n_elements(fa) GT 0 then fv = call_function(f, x, _EXTRA=fa) $ else fv = call_function(f, x) endelse c = fv ;; Compute B-coefficients for j = 0L, i2 do begin F1=J*C2 F2=-HF*F1 C3=2*COS(F1) A2=zero A1=zero A0=C(I2) for K = I2-1,0L,-1 do begin A2=A1 A1=A0 A0=C(K)+C3*A1-A2 endfor BC(J)=C1*(A0*COS(F1+F2)-A1*COS(F2)) BC(I1)=zero endfor c = hf*[ac(0:i1-1)+bc(0:i1-1), rotate(ac(0:i1)-bc(0:i1),2)] cc = abs(c) cmx = max(cc) if (CMX GT 0) THEN begin CMX=1/CMX CC(I3)=HF*CC(I3) A0=CC(I2)*CMX A1=CC(I1)*CMX for J = I1+2,I3 do begin A2=CC(J)*CMX IF(A0 LE EPS AND A1 LE EPS AND A2 LE EPS) THEN $ goto, CHEB9 A0=A1 A1=A2 endfor ENDIF ;; DOUBLE THE NUMBER OF COEFFICIENTS. if i LT 7 then begin ac = c(0:i3) bc = ac*0 endif endfor ;; REQUIRED ACCURACY NOT OBTAINED NC=64 DELTA=total(abs(c(60:nc))) message, 'WARNING: Required accuracy not obtained', /info status = 0 goto, CLEANUP CHEB9: ;; REQUIRED ACCURACY OBTAINED ;; SUM NEGLECTED TERMS IN EXPANSION status = 1 DELTA=total(cc(j:i3)) ;; CHECK IF FURTHER REDUCTION OF COEFFICIENTS IS POSSIBLE. NC=J-1 REST=EPS-DELTA IF (REST GT 0) AND redalg GT 0 THEN begin while (CC(NC) LT REST) do begin DELTA=DELTA+CC(NC) REST=REST-CC(NC) NC=NC-1 endwhile ENDIF CLEANUP: C(0)=HF*C(0) p = c(0:nc) rest = eps - delta if redalg EQ 2 then begin wh = where(cc(0:nc) LT prec, ct) i = ct-1 while (i GE 0) AND ((rest GT 0) OR (status EQ 1)) do begin delta = delta + cc(wh(i)) rest = rest - cc(wh(i)) p(wh(i)) = 0 i = i - 1 endwhile endif DONE: igood = where(p NE 0) err = delta RETURN, p end