j.davis - october 1992 - subroutine legnds
modified 1993 January - J.F.Chandler
derived from: f amuchastegui - september 1969 - subroutine legndr
evaluation of scaled legendre polynomials and scaled legendre functions
using recursion formulas

a(n,h) = sqrt{((2*n + 1)*(n-h))/((2*n - 1)*(n+h))}
b(n,h) = sqrt{((2*n + 1)*(n+h-1))/((2*n - 1)*(n+h))}     h>1
b(n,1) = sqrt{((2*n + 1)*(n)*2)/(n + 1)}
d(n,h) = sqrt{(n + h)*(n - h + 1)}                       h>1
d(n,1) = sqrt{(n + 1)*(n)*2*(2*n + 1)}

legnds  evaluates the legendre polynomials, its
derivatives and also the associated legendre functions
the legendre polynomials are defined as:
pn(z) = (1 / (n!*2**n) ) * (nth derivat of (z*z-1)**n)
and the associated legendre functions as
p(n,h) = ((sqrt(1-z*z))**h) * ph(n) ,
where ph(n) is the  hth derivative of  pn(z)
we can see that for h=0  p(n,0)=pn(z)
The Legendre functions are scaled by the factor
sqrt({2*(2*n+1)*(n-h)!}/{(n+h)!}.  The Legendre polynomials
are left unscaled in this routine.

the formulae are from  "tables of integrals, series
and products" - gradsteyn & ryzhik,pags 1004-27
$$$$$$$ the factor (-1)**h  for the legendre function used in this
reference has been set = 1 to use the convention of
the smithsonian astr. observatory journal
the formulae for  p'(z)  and  p''(z)  are from "differential
equations with applications", ritger & rose, page 223

evaluation of  p(z)
the first polynomial in the array  leg(n)  is the
second order polynomial, since the first order is not
used in PEP

evaluation of p(n,h)                 -- scaled from LEGNDR routine
the order of the legendre functions in the array gleg,
as well as the order of the partials in gleg1, gleg2 is
p(2,1), p(2,2), p(3,1), p(3,2), p(3,3), ....... , p(n,n)
The following recursive formulas were used for p(n,h) based on the
modification that p(n,h) is scaled for h>0 and unscaled for h=0
p(n,h) = a(n,h)*z*p(n-1,h) + zz*b(n,h)*p(n-1,h-1)
p(n,n) = zz*b(n,n)*p(n-1,n-1)

evaluation of  p'(n,h)                 -- scaled from LEGNDR routine
The following recursive formulas were used for p'(n,h) based on the
modification that p'(n,h) is scaled for h>0 and unscaled for h=0
p'(n,h) = h*z/zz^2 * p(n,h) - d(n,h)/zz * p(n, h-1)

evaluation of  p''(n,h)                 -- scaled from LEGNDR routine
The following recursive formulas were used for p''(n,h) based on the
modification that p''(n,h) is scaled for h>0 and unscaled for h=0
p''(n,h) = h*(1+z^2)/zz^4 * p(n,h) + h*z/zz^2 * p'(n,h)
- d(n,h)/zz^3 * {z * p(n,h-1) + zz^2 * p'(n, h-1)}