/usr/share/gnudatalanguage/astrolib/polyleg.pro is in gdl-astrolib 2018.02.16+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 | 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
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 size(x,/TNAME) EQ 'DOUBLE' then begin
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
|