;+ ; NAME: ; CHEBGRID ; ; 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 a grid ; ; MAJOR TOPICS: ; Curve and Surface Fitting ; ; CALLING SEQUENCE: ; p = CHEBGRID(T, X, [ DXDT, NPOINTS=, NPOLY=, NGRANULE= , $ ; RMS=, DRMS=, RESIDUALS=, DRESIDUALS= , $ ; XMATRIX=, DXMATRIX=, RESET=, ; DERIV_WEIGHT= ] ) ; ; DESCRIPTION: ; ; CHEBGRID estimates the coefficients for a finite sum of Chebyshev ; polynomials approximating a continuous tabulated function over an ; interval. The function (and optionally its derivative) must be ; tabulated on a regularly sampled grid. The implementation of this ; function is taken from a method described by X. X. Newhall, used ; in estimating coefficients for ephemerides in the solar system. ; ; The tabulated function is assumed to be continuous over the entire ; interval. A Chebyshev series is fitted to the function over small ; segments, called granules. The size of each granule, the number ; of points in each granule, and the number of Chebyshev polynomials ; are all configurable. ; ; Users may specify either the function alone, or the function and ; its first derivative. By also giving the tabulated derivative, a ; more accurate Chebyshev polynomial can be developed. Aside from ; the constraints mentioned in the next paragraph, the polynomial ; that is returned is the best-fit polynomial in a least-squares ; sense. ; ; Here is a definition of terms: ; ; GRANULE - a single continuous fitted segment. The length of the ; granule, NGRANULE, is specified in units of the tabulated ; grid size. Because of the continuity requirements developed ; below, granules will always overlap at their endpoints. ; Thus, then length of a granule should be a factor of ; N_ELEMENTS(X)-1. For simple functions over short intervals, ; the granule size can be equal to N_ELEMENTS(X)-1 ; ; NUMBER OF POINTS the number of points, NPOINTS, within a ; granule to be fitted to the polynomial, not necessarily ; equal to the granule size. The greater the number of ; points, the more computation time and storage is required. ; This number *must* be a factor of NGRANULE. Typically ; NPOINTS is a number between 8 and 12. Because of the ; single-point overlap between granules (see below), the ; actual number of points per fit is NPOINTS+1. ; ; NUMBER OF POLYNOMIALS the number of Chebyshev polynomial terms, ; NPOLYNOMIAL, to be fitted per granule. The greater the ; number of polynomial terms, the more computation time and ; storage is required, but also the greater the approximating ; precision of the fit. ; ; The particular set of Chebyshev polynomial coefficients developed ; by this function have some special properties. If both the ; function and its derivative are specified, then the value and ; derivative of the interpolating polynomial at the granule ; endpoints will be exactly equal to the tabulated endpoint values. ; This feature allows many approximations to be strung together ; piecewise, and the function value and first derivative will be ; continuous across granule boundaries. ; ; If only the function value is specified, then only the function ; value will be continuous at the granule endpoints, and not the ; derivative. ; ; An extensive set of statistics are computed to assess the quality ; of the Chebyshev polynomial fit. The keywords RESIDUALS and ; DRESIDUALS return the residuals of the fit after subtracting the ; interpolation. The RMS and DRMS keywords return the root mean ; squared deviations between data and model. ; ; If the user does not know how many granules, points, or polynomial ; coefficients to use, then he or she should try several ; combinations and see which minimizes the r.m.s. value with the ; fewest number of coefficients. ; ; If the XMATRIX and DXMATRIX keywords are passed, then CHEBGRID ; attempts to avoid recomputing several of the matrices it uses in ; estimating the coefficients. If multiple calls to CHEBGRID are to ; be made, some compution time savings can be made. In the first ; call CHEBGRID the required matrices are computed and returned. In ; subsequent calls, CHEBGRID detects the XMATRIX and DXMATRIX ; keyword values and uses those values if it can. ; ; The user can also estimate their own coefficients. The matrices ; returned are (NPOINTS+1)x(NPOLYNOMIAL). The coefficients from a ; NPOINTS+1 tabulation, X, are found by: ; ; PCHEB = XMATRIX ## X + DXMATRIX ## DXDT ; ; if derivative information is known, or ; ; PCHEB = XMATRIX ## X ; ; if no derivative information is known. [ Note: the matrices are ; different, depending on whether derivative information is known or ; not. ] ; ; ; INPUTS: ; ; T - array of regularly sampled *independent* variables. The number ; of elements in T should be a multiple of NGRANULE, plus one. ; ; X - array of regularly sampled *dependent* variables. The number ; of elements in X should be equal to the number of elements in ; T. ; ; DXDT - optionally, a tabulated array of first derivatives of X ; with respect to T, at the same grid points. ; ; KEYWORD PARAMETERS: ; ; NGRANULE - size of a "granule", in grid intervals. NGRANULE must ; be at least 2, and a factor of N_ELEMENTS(T)-1. ; Default: 8 ; ; NPOINTS - number of points per granule that are fitted. NPOINTS ; must be at least 2, and a factor of NGRANULE. ; Default: NGRANULE ; ; NPOLYNOMIAL - number of Chebyshev polynomial terms per fit. ; NPOLYNOMIAL must be at least 2 and less than ; 2*(NPOINTS+1), when derivative information is ; specified; or less than NPOINTS+1, when no ; derivative information is specified. ; Default: 7 ; ; RESIDUALS - upon return, an array of size N_ELEMENTS(T), with ; residuals of the tabulated function minus the ; interpolated function. ; ; DRESIDUALS - same as RESIDUALS, but for the function's first ; derivative. ; ; RMS - upon return, the root mean square of the function value ; residuals. ; ; DRMS - same as RMS, but for the function's first derivative. ; ; XMATRIX - upon return, the matrix used to compute Chebyshev ; polynomial coefficients from the function value. ; ; Upon input, CHEBGRID determines if XMATRIX will apply to ; the data, and if so, XMATRIX is reused rather than ; computed. If XMATRIX cannot be reused, then it is ; computed afresh, and the new value is returned in the ; XMATRIX keyword. ; ; The user should not modify the contents of this array. ; ; DXMATRIX - same as XMATRIX, but for the function's first ; derivative. ; ; RESET - if set, force a recomputation of XMATRIX and/or DXMATRIX. ; ; DERIV_WEIGHT - amount of weight to give to function derivative, ; relative to the function value. ; Default: 0.16d ; ; ; RETURNS: ; ; An array of coefficient values. The dimensions of the array are ; NPOLYNOMIALxNSEGS, where NSEGS is the number of granules in the ; entire interval. ; ; ; EXAMPLE: ; ; ;; Estimate Chebyshev coefficients for the function SIN(X), on the ; ;; interval [-1,+1]. ; xx = dindgen(9)/4d - 1d ;; Regular grid from -1 to 1 (9 points) ; yy = sin(xx) ;; Function values, sin(x), ... ; dy = cos(xx) ;; ... and derivatives ; ; ;; Estimate coefficients using CHEBGRID (single granule of 8 intervals) ; p = chebgrid(xx, yy, dy, npoints=8, ngranule=8, npoly=10) ; ; xxx = dindgen(1001)/500 - 1d ;; New grid for testing ; res = sin(xxx) - chebeval(xxx, p) ; plot, xxx, res ; ; ;; Same as example above, except extended range to [-1, +15], ; using eight granules. ; xx2 = dindgen(65)/4d - 1 ; yy2 = sin(xx2) ; dy2 = cos(xx2) ; p = chebgrid(xx2, yy2, dy2, ngranule=8, npoint=8, npoly=10) ; help, p ; P DOUBLE = Array[10, 8] ; ;; (i.e., 10 polynomial coefficients over 8 granules) ; ; ; REFERENCES: ; ; Abramowitz, M. & Stegun, I., 1965, *Handbook of Mathematical ; Functions*, 1965, U.S. Government Printing Office, Washington, ; D.C. (Applied Mathematical Series 55) ; Newhall, X. X. 1989, Celestial Mechanics, 45, p. 305-310 ; ; MODIFICATION HISTORY: ; Written, CM, Feb 2002 ; Documented, CM, 24 Mar 2002 ; Corrected documentation, CM, 28 Apr 2002 ; Typo correction, CM, 10 Oct 2002 ; ; $Id: chebgrid.pro,v 1.5 2002/11/07 00:15:23 craigm Exp $ ; ;- ; Copyright (C) 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. ;- ;; Utility function: compute XMATRIX and DXMATRIX using Newhall approach pro chebpcmat, npts, npoly, xmat, vmat, dweight=weight0 ;; n0 is the number of intervals in Cheb approx. n0 = npts - 1 if n_elements(weight0) EQ 0 then $ weight = 0.16d $ else $ weight = weight0(0) tmat = dblarr(npoly, npts) tdot = tmat cj = dblarr(npoly) xj = 1d - 2d*dindgen(npts)/n0 for i = 0, npoly-1 do begin cj(*) = 0 & cj(i) = 1 tmat(i,*) = chebeval(xj, cj, deriv=v) tdot(i,*) = v endfor ;; Form matrix T*W tw = dblarr(2,npts,npoly) tw(0,*,*) = transpose(tmat) tw(1,*,*) = transpose(tdot) * weight ;; Form matrix T*WT twt = reform(tw(0,*,*),npts,npoly) ## tmat + $ reform(tw(1,*,*),npts,npoly) ## tdot tw = reform(tw, 2*npts, npoly, /overwrite) twt = reform(twt, npoly, npoly, /overwrite) ;; Augment matrix T*W to get matrix C2 c2 = dblarr(2*npts,npoly+4) c2(*,0:npoly-1) = tw c2(0,npoly) = 1 & c2(1,npoly+1) = 1 c2(2*npts-2,npoly+2) = 1 & c2(2*npts-1,npoly+3) = 1 ;; Augment matrix T*WT to get the matrix C1 c1 = dblarr(npoly+4,npoly+4) c1(0:npoly-1,0:npoly-1) = twt c1(0:npoly-1,npoly+0) = tmat(*,0) c1(0:npoly-1,npoly+1) = tdot(*,0) c1(0:npoly-1,npoly+2) = tmat(*,npts-1) c1(0:npoly-1,npoly+3) = tdot(*,npts-1) c1(npoly:*,0:npoly-1) = transpose(c1(0:npoly-1,npoly:*)) ;; Compute matrix C1^(-1) c1inv = invert(c1) ;; Compute matrix C1^(-1) C2 c1c2 = c1inv ## c2 c1c2 = reform(c1c2, 2,npts,npoly+4) c1c2 = reverse(c1c2,2) c1c2 = reform(c1c2, 2*npts,npoly+4) ii = lindgen(npts)*2 xmat = c1c2(ii,0:npoly-1) ;; Split into terms multiplying Y and VY vmat = c1c2(ii+1,0:npoly-1) return end ;; Utility function: compute XMATRIX only, using only the constraint ;; on the function values at the endpoints. pro chebpcmat_xonly, npts, npoly, xmat ;; n0 is the number of points in Cheb approx. n0 = npts - 1 tmat = dblarr(npoly, npts) cj = dblarr(npoly) xj = 1d - 2d*dindgen(npts)/n0 for i = 0, npoly-1 do begin cj(*) = 0 & cj(i) = 1 tmat(i,*) = chebeval(xj, cj, deriv=v) endfor ;; Augment matrix T to get matrix C2 c2 = dblarr(npts,npoly+2) c2(*,0:npoly-1) = transpose(tmat) c2(0,npoly) = 1 c2(npts-1,npoly+1) = 1 ;; Augment matrix T*WT to get the matrix C1 c1 = dblarr(npoly+2,npoly+2) c1(0:npoly-1,0:npoly-1) = transpose(tmat) ## tmat c1(0:npoly-1,npoly+0) = tmat(*,0) c1(0:npoly-1,npoly+1) = tmat(*,npts-1) c1(npoly:*,0:npoly-1) = transpose(c1(0:npoly-1,npoly:*)) ;; Compute matrix C1^(-1) c1inv = invert(c1) ;; Compute matrix C1^(-1) C2 c1c2 = c1inv ## c2 c1c2 = reform(c1c2, npts,npoly+2) c1c2 = reverse(c1c2,1) xmat = c1c2(*,0:npoly-1) return end function chebgrid, t, x, dxdt, ngranule=ngran0, npoints=npts0, $ npolynomial=npoly0, deriv_weight=dweight0, $ rms=rms, drms=drms, residuals=resid, dresiduals=dresid, $ xmatrix=xmatrix, dxmatrix=dxmatrix, reset=reset ;; Default processing if n_elements(ngran0) EQ 0 then ngran = 8 $ else ngran = round(ngran0(0)) > 2 if n_elements(npts0) EQ 0 then npts = ngran $ else npts = round(npts0(0)) > 2 if n_elements(npoly0) EQ 0 then npoly = 7 $ else npoly = round(npoly0(0)) > 2 ;; Error checking if ngran LT npts then begin message, 'ERROR: Granule size ('+strtrim(ngran,2)+') is too '+ $ 'small for number of samples ('+strtrim(npts,2)+')' return, !values.d_nan endif ;; Be sure NGRAN is a multiple of NPTS - or not. Instead, a warning ;; message is printed in the loop. ; if abs(double(ngran)/npts - round(ngran/npts)) GT 1d-5 then begin ; message, 'ERROR: NPOINTS must be a multiple of NGRANULE' ; return, !values.d_nan ; endif ;; Be sure we are solving a least-squares problem. If the number of ;; polynomials is too great then it becomes underconstrained, not ;; overconstrained. if n_elements(dxdt) GT 0 then begin if npoly GE 2*(npts+1) then $ message, 'ERROR: NPOLYNOMIAL must be less than 2*(NPOINTS+1)' endif else begin if npoly GE npts+1 then $ message, 'ERROR: NPOLYNOMIAL must be less than NPOINTS+1' endelse ;; Begin size checking of input matrices - we may be able to use the ;; previously computed version. szx = size(xmatrix) szv = size(dxmatrix) ;; Cases: recompute because existing X matrix is wrong size; ;; recompute because existing V matrix is wrong size; ;; recompute because a V matrix was passed, but no DXDT was redo_x = (szx(0) NE 2 OR szx(1) NE npts+1 OR szx(2) NE npoly) redo_v = (n_elements(dxdt) GT 0 AND $ (szv(0) NE 2 OR szv(1) NE npts+1 OR szv(2) NE npoly)) no_v = (n_elements(dxdt) EQ 0 AND n_elements(dxmatrix) GT 0) ;; Actual recomputation of matrices if redo_x OR redo_v OR no_v OR keyword_set(reset) then begin COMPUTE_CHEBMAT: xmatrix = 0 & dummy = temporary(xmatrix) dxmatrix = 0 & dummy = temporary(dxmatrix) if n_elements(dxdt) GT 0 then $ chebpcmat, npts+1, npoly, xmatrix, dxmatrix, dweight=dweight0 $ else $ chebpcmat_xonly, npts+1, npoly, xmatrix endif rms = 0.*x(0) drms = rms chebm = dblarr(npoly, (n_elements(x)-1)/ngran) resid = x*0. dresid = resid ispan = lindgen(npts+1)*(ngran/npts) imax = max(ispan) ng = 0L for ibase = 0, n_elements(x)-1, ngran do begin if ibase EQ n_elements(x)-1 then goto, DONE if n_elements(x)-ibase LT ngran+1 then begin nlost = n_elements(x)-ibase message, 'WARNING: last '+strtrim(nlost,2)+' elements of X '+$ 'were discarded because they formed only a fractional granule.', $ /info goto, DONE endif tspan = [t(ibase), t(ibase+imax)] tgran = t(ibase:ibase+imax-1)-t(ibase) dt = tspan(1) - tspan(0) tspan = tspan - tspan(0) ;; Compute the X portion of the coefficients xgran = x(ibase+ispan) chebi = xmatrix ## xgran ;; Compute the DXDT portion if it is available if n_elements(dxdt) GT 0 then begin dxgran = dxdt(ibase+ispan) * dt/2. chebi = chebi + dxmatrix ## dxgran ;; Statistics - V first, then X comes later xmod = chebeval(tgran, chebi, interval=tspan, derivative=dxmod) ;; DXDT portion of statistics dresid(ibase:ibase+imax-1) = dxdt(ibase:ibase+imax-1) - dxmod diff_dx = (dxdt(ibase:ibase+imax-1) - dxmod)^2 drms = drms + total(diff_dx) endif else begin ;; Statistics - X only xmod = chebeval(tgran, chebi, interval=tspan) endelse ;; Finish statistics with X portion resid(ibase:ibase+imax-1) = x(ibase:ibase+imax-1) - xmod diff_x = ( x(ibase:ibase+imax-1) - xmod)^2 rms = rms + total(diff_x) ;; Append to existing coefficient list chebm(*,ng) = chebi(*) ng = ng + 1L endfor DONE: ;; Final adjustments to statistics rms = sqrt( rms / ngran) if n_elements(dxdt) GT 0 then drms = sqrt(drms / ngran) return, chebm end