;+
; 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