;+ ; NAME: ; LEGCHEB ; ; 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: ; Compute Legendre polynomial coefficents from Chebyshev coefficients ; ; MAJOR TOPICS: ; Curve and Surface Fitting, Special Functions ; ; CALLING SEQUENCE: ; b = LEGCHEB(a) ; ; DESCRIPTION: ; ; This routine computes the coefficients of a Legendre polynomial ; expansion when the Chebyshev expansion is known. ; ; Users can determine the Chebyshev expansion coefficients using a ; routine like CHEBFIT, CHEBCOEF or CHEBGRID. Then, if the Legendre ; expansion is needed instead, this conversion routine should be ; used. Evaluation of the Legendre series can be performed using ; the POLYLEG function in the IDL Astronomy Library. ; ; Internally, the computational precision is double precision. ; This routine relies upon the algorithm of Piessens (1974). ; ; INPUTS: ; ; A - a vector, the coefficients of the Chebyshev series of the ; desired function. ; ; RETURNS: ; ; The vector B, which contains the coefficients of the Legendre ; polynomial expansion. Both A and B will have the same number of ; elements and data type. ; ; KEYWORD PARAMETERS: ; ; NONE ; ; EXAMPLE: ; ; ;; Compute the Chebyshev series coefficients of 1/(2-X) on [-1,1] ; A = CHEBCOEF('1d/(2d - X)', /expr) ; ; ;; Convert to Legendre series coefficients ; B = LEGCHEB(A) ; ; REFERENCES: ; ; Abramowitz, M. & Stegun, I., 1965, *Handbook of Mathematical ; Functions*, 1965, U.S. Government Printing Office, Washington, ; D.C. (Applied Mathematical Series 55) ; Piessens, R. 1974, Comm. ACM, v. 17, p. 25 (TOMS 473) ; ; MODIFICATION HISTORY: ; Written and documented, CM, 25 Sep 2002 ; ; $Id: legcheb.pro,v 1.1 2002/09/25 21:12:35 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. ;- function legcheb, a, reset=reset common legcheb_common, ink n1 = n_elements(a) n = n1-1 ;; Internal routine: reset the common block if keyword_set(reset) then begin ink = 0 & dummy = temporary(ink) return, -1d endif ;; A common block is used to store the matrix which converts from ;; one representation to another. This matrix needs to be expanded ;; if the input vector (size N1) is longer than the size of the ;; matrix (NINK x NINK). nink = sqrt(n_elements(ink)) if nink LT n1 then begin ;; If we've never been called before, initialize the 2x2 array ;; with hard coded numbers if n_elements(ink) EQ 0 then begin ink = reform([[2d,0d],[0d, 2d/3d]], 2,2) nink = 2 endif ;; Insert the new array into the old array inknew = dblarr(n1,n1) inknew(0:nink-1,0:nink-1) = ink ink = reform(inknew,n1,n1) ;; Compute the diagonal components, based on recurrence relation ;; in Piessens (1974) for nn = nink, n1-1 do begin ;; Recurrence relation for diagonal components ;; ORIG: ink(nn,nn) = ink(nn-1,nn-1)*4d*nn^2/(2d*nn+1)/(2d*nn) ink(nn,nn) = ink(nn-1,nn-1)*2d*nn/(2d*nn+1) endfor ;; Special case: first row, because of a 0/0 condition kk = dindgen(n1/2)*2 ink(kk,0) = -2d/(kk*kk - 1) ;; Fill remaining columns of the INK array, using the recurrence ;; of eqn 6 in Piessens for nn = 1, n1-1 do begin ;; KSTART is the starting index, which could be larger than ;; NINK because the elements left of the diagonal are always ;; zero. The NINK-1 vs NINK-2 logic is because the cells ;; alternate nonzero quantities. kstart = nink-1 if ink(kstart,nn) EQ 0 then kstart = nink-2 kstart = kstart>nn ;; Recurrence used here. for kk = kstart, n1-3, 2 do if ink(kk,nn) NE 0 then begin ink(kk+2,nn) = (ink(kk,nn) $ * (double((kk-1)*kk - nn*(nn+1)) $ / ((kk+3)*(kk+2) - nn*(nn+1))) $ * (double(kk+2)/kk)) end endfor endif ;; Extract relevant portion of INK matrix nn = dindgen(n1) mat = ink(0:n1-1,0:n1-1) ;; Keep same data type for A and B b = reform(a)*0. ;; Compute B by matrix multiplication b(*) = (mat ## a(*)) ;; Apply final multiplicative factor, the normalization of the ;; Legendre polynomials. return, b*(0.5d + nn) end