function polyleg,x,coeff
;+
; NAME:
;       POLYLEG
;
; PURPOSE:
;       Evaluate a Legendre polynomial with specified coefficients.
; EXPLANATION:
;       Meant to be used analogously to the POLY function in the IDL User's
;       Library distribution.
;
; CALLING SEQUENCE:
;       Result = POLYLEG( X, C )        
;
; INPUTS:
;       X - input variable, scalar or vector    
;       C - vector of Legendre polynomial coefficients. 
; OUTPUTS:
;       POLYLEG returns a result equal to:
;               C(0) + C(1)*P_1(x) + C(2)*P_2(x) + ...
;
;       where P_j(x) is the jth Legendre polynomial.   The output will have
;       the same dimensions as the input X variable.
;
; EXAMPLE:
;       If x = [0.5, 1.0] and C = [2.4, 1.3, 2.5] then
;       print, polyleg(x, c)    ====> [2.7375, 6.20]
;
;       The result can be checked using the first 3 Legendre polynomial terms
;       C(0) + C(1)*x + C(2)*(0.5*(3*x^2-1))
; METHOD:
;       Uses the recurrence relation of Legendre polynomials
;               (n+1)*P_n+1(x) = (2n+1)*x*P_n(x) - n*P_n-1(x)
;       evaluated with the Clenshaw recurrence formula, see Numerical Recipes
;       by Press et al. (1992), Section 5.5
;
; REVISION HISTORY:
;       Written W. Landsman   Hughes STX Co.        April, 1995    
;       Fixed for double precision  W. Landsman     May, 1997
;       Converted to IDL V5.0   W. Landsman   September 1997
;-
 On_error,2
 
 if N_params() LT 2 then begin
        print,'Syntax - result = POLYLEG( X, Coeff)'
        return, -1
 endif

 sz_x = size(x)
 N= N_elements(coeff) -1
 M = N_elements(x)

 case N of 
 0: return, replicate( coeff, M)
 1: return, x* coeff[1] + coeff[0]
 else:
 endcase

; If X is double then compute in double; otherwise compute in real

 if sz_x[sz_x[0] + 1] EQ 5 then begin        ;Fixed May 1997
        y = dblarr( M, N+2)
        jj = dindgen(N) + 2.0d
 endif else begin
        y = fltarr( M, N+2 )
        jj = findgen(N) + 2.
 endelse

 beta1 =  -jj / (jj+1)
 for j = N,1,-1 do begin

        alpha = (2*j + 1.)*x/float(j + 1.) 
        y[0,j-1] = alpha*y[*,j] + beta1[j-1]*y[*,j+1] + coeff[j]
 endfor

 return, -0.5*y[*,1] + x*y[*,0] + coeff[0]
 end