function legendre,n,x
;+
; ROUTINE:  legendre
;
; PURPOSE:  compute legendre polynomial
;
; USEAGE:   result=legendre(n,x)
;
; INPUT:    
;  n        order of polynomial
;  x        argument of polynomial
;
;
; DISCUSSION:
;           use recursion relation to compute legendre polynomial
;           for orders greater than 2:
;
;           n*P(n,x)=(2n-1)*x*P(n-1,x)-(n-1)*P(n-2,x)
;
;           P(0,x)=1.  P(1,x)=x
;
;  
; EXAMPLE:  
;
;    x=findrng(-1,1,10000)
;    f=legendre(20,x)
;    plot,x,f
;
;; find zeroes of legendre polynomial
;
;    r=roots(f,0.) & xxx=interpolate(x,r) & fff=interpolate(f,r)
;    oplot,xxx,fff,psym=2
;    print,xxx
;
;; print zeroes of rescaled legendre polynomial n=10,20,2
;
;    for n=10,20,2 do begin & f=legendre(n,x) & r=roots(f,0.) & print,n,min(abs(cos(!dtor*39.)-(interpolate(x,r)+1.)/2.)) & endfor
;
; AUTHOR:   Paul Ricchiazzi                        06 Apr 98
;           Institute for Computational Earth System Science
;           University of California, Santa Barbara
;           paul@icess.ucsb.edu
;
; REVISIONS:
;
;-
;
valmm=1.
valm=x

if n eq 0 then return,replicate(valmm,n_elements(x))
if n eq 1 then return,valm

for i=2,n do begin
  ri=float(i)
  val=(2.*ri-1.)*x*valm-(ri-1.)*valmm
  val=val/ri
  valmm=valm
  valm=val
endfor

return,val
end