;+
; NAME: 
;       ALEG
; 
; PURPOSE:
;       This procedure calculates the associated Legendre functions P(l,n)
;       at a given point x, up to a maximum order P(nmax,nmax). It is based
;       on the varying order recurrence relation given by Abramovitz & Stegun,
;       Handbook of Mathematical Functions, 1970 (eq. 8.5.3, p. 333). The
;       Schmidt normalization is used here (see Akasofu & Chapman, Solar
;       Terrestrial Physics, 1972, p. 77), in accordance with usual practice
;       in geomagnetic studies. The procedure also gives d P(n,m) d theta, 
;       where theta = arccos(x), by taking the derivative of the same
;       recurrence relation. Finally, the procedure gives the value of
;	P(n,m)/sin(theta).
;
; CATEGORY:
;       Mathematical procedures.
;
; CALLING SEQUENCE:
;       ALEG, x, nmax, p, q, s
;
; INPUTS:
;          X: the point at which P(n,m) is to be evaluated.
;
;       NMAX: maximum order of P(n,m).
;
; KEYWORDS:
;       None.
;
; OUTPUTS:
;       P: an array of size (lmax + 1) x (lmax + 1), containing the function values
;          (P(n,m) = 0 for n < m).
;       Q: an array of size (lmax + 1) x (lmax + 1), containing the values of 
;          d P(n,m) / d theta.
;	S: an array of size lmax + 1) x (lmax + 1), containing the values of
;	   P(n,m)/sin(theta).
;
; COMMON BLOCKS:
;       None.
;
; SIDE EFFECTS:
;       None.
;
; RESTRICTIONS:
;       None.
;
; MODIFICATION HISTORY:
;       September 6, 1993. Written by A. Louro.
;-

pro aleg, x, nmax, p, q, s

p = fltarr(nmax + 1, nmax + 1)
q = fltarr(nmax + 1, nmax + 1)
s = fltarr(nmax + 1, nmax + 1)
aux = sqrt(1. - x * x)
sr2 = sqrt(2.)

p(0,0) = sr2 				; Set initial values for
p(1,0) = sr2 * x			; recurrence procedure
q(0,0) = 0.
q(1,0) = - sr2 * aux

for n = 1, nmax do begin 		; Calculate diagonal values P(n,n)
  c1 = sqrt(1. - 1. / (2. * n))
  p(n,n) = c1 * aux * p(n - 1, n - 1)
  q(n,n) = c1 * (x * p(n - 1, n - 1) + aux * q(n -1, n - 1))
  s(n,n) = c1 * p(n - 1, n - 1)
endfor

for m = 0, nmax - 1 do begin		; Fill in rest of array
  mm = max([m, 1])
  for n = mm, nmax - 1 do begin
    c2 = 2. * n + 1.
    c3 = sqrt((n + m) * (n - m))
    c4 = sqrt((n - m + 1) * (n + m + 1))
    p(n + 1, m) = (c2 * x * p(n,m) - c3 * p(n - 1, m)) / c4
    q(n + 1, m) = (c2 * (x * q(n,m) - aux * p(n,m)) - $
		   c3 * q(n - 1, m)) / c4
    s(n + 1, m) = (c2 * x * s(n,m) - c3 * s(n - 1, m)) / c4
  endfor
endfor

for n = 0, nmax do begin		; "Correct" (n,0) terms by
  p(n, 0) = p(n, 0) / sr2		; dividing by sqrt(2)
  q(n, 0) = q(n, 0) / sr2
endfor
return
end