NB. Associated Legendre polynomials P_l^m(x), x=theta=th
NB. These are renormalized so multiplication by e^(i m \phi) gives 
NB. Y_lm(\theta,\phi), See "Numerical Recipes" 3rd ed p 292-295.
NB. Usage: (l m)&plegendre th (with the &, th can be an array) 
Plegendre=: 4 : 0
'l m'=. x
xx=. cos y
pmm=. 1
   if. m>0 do.
      omx2=. (1-xx)*(1+xx)
      fact=. k=. 1 
      whilst. k<:m do.
	 pmm=. pmm*omx2*fact%(fact+1)
	 fact=. fact+2
	 k=. k+1
      end. 
   end.   
      pmm=. %: (1+2*m)*pmm%4p1
      if. 2&|m do. pmm=. -pmm end.
      if. l=m do. pmm return. end.
      pmmp1=. pmm*xx*%:(3+2*m)
      if. l=m+1 do. pmmp1 return. end. 
      oldfact=. %:(3+2*m)
      ll=. m+2
       whilst. ll<:l do.
	 l2=. ll*ll 
         fact=. %: (_1+4*l2)%(l2-m*m)
         pll=. fact*(xx*pmmp1)-(pmm%oldfact)
         oldfact=. fact
         pmm=. pmmp1
         pmmp1=. pll 
         ll=. ll+1 
       end.
       pll 
)

NB. Example:  int01=: i.@>: % ]
NB. xx=. _1+2*int01 2000
NB. z=. (15 0)Plegendre xx
NB. load'plot'
NB. 'pensize 4' plot xx;z

eim=: [: ^0j1*]  NB. eim z = ^0j1*z
Ythp=: 4 : 0
'l m'=. x
'th phi'=. y
(x&Plegendre th)*eim m*phi
)

NB. This provides P_l^m(x) without renomalization.
NB. Beware that the results become huge for large m.
plgndr=: 4 : 0
'l m'=. x
prod=. 1.0
j=. l+1-m
while. j<: l+m do.
   prod=. prod * j
   j=. j+1
end.
(4p1*prod%1+2*l)*(x&Plegendre y)
)