program spherical
      implicit none
      integer*4 l,m
      real*8 t,f,re,im,fp,pi

      call InitFac
      pi=4.0d0*datan(1.0d0)
      fp=4.0d0*pi

      t=4.0d0
      f=-6.0d0
      call    sph(0, 0,t,f,re,im)
      write(6,900)0, 0,re,im,1.0d0/dsqrt(fp),0.0d0
900   format(2i3,4f12.6)
      call    sph(1, 0,t,f,re,im)
      write(6,900)1, 0,re,im,0.5d0*dsqrt(3.0d0/pi)*cos(t),0.0d0
      call    sph(1, 1,t,f,re,im)
      write(6,900)1, 1,re,im,-0.5d0*dsqrt(3.0d0/pi/2.0d0)*sin(t)*cos(f)
     1                      ,-0.5d0*dsqrt(3.0d0/pi/2.0d0)*sin(t)*sin(f)
      call    sph(1,-1,t,f,re,im)
      write(6,900)1,-1,re,im,+0.5d0*dsqrt(3.0d0/pi/2.0d0)*sin(t)*cos(f)
     1                      ,-0.5d0*dsqrt(3.0d0/pi/2.0d0)*sin(t)*sin(f)
      call    sph(2, 0,t,f,re,im)
      write(6,900)2, 0,re,im,+0.25d0*dsqrt(5.0d0/pi)*(3.0d0*cos(t)**2-
     1                        1.0d0),0.0d0
      call    sph(2, 1,t,f,re,im)
      write(6,900)2, 1,re,im,-0.5d0*dsqrt(15.0d0/pi/2.0d0)*
     1sin(t)*cos(t)*cos(f),  +0.5d0*dsqrt(15.0d0/pi/2.0d0)*
     1sin(t)*cos(t)*sin(f) 
      call    sph(2,-1,t,f,re,im)
      write(6,900)2,-1,re,im, 0.5d0*dsqrt(15.0d0/pi/2.0d0)*
     1sin(t)*cos(t)*cos(f),  -0.5d0*dsqrt(15.0d0/pi/2.0d0)*
     1sin(t)*cos(t)*sin(f) 
      call    sph(2, 2,t,f,re,im)
      write(6,900)2, 2,re,im,   0.25d0*dsqrt(15.0d0/pi/2.0d0)*
     1(sin(t))**2*cos(2.0d0*f), 0.25d0*dsqrt(15.0d0/pi/2.0d0)*
     1(sin(t))**2*sin(2.0d0*f) 
      call    sph(2,-2,t,f,re,im)
      write(6,900)2,-2,re,im,   0.25d0*dsqrt(15.0d0/pi/2.0d0)*
     1(sin(t))**2*cos(2.0d0*f),-0.25d0*dsqrt(15.0d0/pi/2.0d0)*
     1(sin(t))**2*sin(2.0d0*f) 
      call    sph(6, 5,t,f,re,im)
      write(6,900)6, 5,re,im,  
     1-3.0d0/16.0d0*dsqrt(77.0d0*13.0d0/pi/4.0d0)*
     1(sin(t))**5*cos(t)*cos(5.0d0*f),
     1-3.0d0/16.0d0*dsqrt(77.0d0*13.0d0/pi/4.0d0)*
     1(sin(t))**5*cos(t)*sin(5.0d0*f)
      call    sph(6,-5,t,f,re,im)
      write(6,900)6,-5,re,im,  
     1 3.0d0/16.0d0*dsqrt(77.0d0*13.0d0/pi/4.0d0)*
     1(sin(t))**5*cos(t)*cos(5.0d0*f),
     1-3.0d0/16.0d0*dsqrt(77.0d0*13.0d0/pi/4.0d0)*
     1(sin(t))**5*cos(t)*sin(5.0d0*f)

c     read(5,*)l,m,t,f
c     call sph(l,m,t,f,re,im)
c     write(6,900)l,m,re,im

      end

      subroutine sph(l,m,t,f,re,im)
c     spherical function Yl,m(t,f) ... r, im, real and imaginaty parts
      implicit none
      integer*4 l,m
      real*8 t,f,re,im,n,norm,pol,p
      norm=n(l,m)
      pol=p(m,l,t)
      re=norm*pol*cos(dble(m)*f)
      im=norm*pol*sin(dble(m)*f)
      return
      end

      subroutine ct(l,m,t,f,rect,imct)
c     tensor operators
      implicit none
      integer*4 l,m
      real*8 pi,t,f,rect,imct,f,re,im
      pi=4.0d0*datan(1.0d0)
      call sph(l,m,t,f,re,im)
      f=dsqrt(4.0d0*pi/dble(2*l+1))
      rect=f*re
      imct=f*im
      return
      end

      function p(m,l,t)
c     associated Legendre polynomial              
      implicit none
      integer*4 nfac0
      parameter (nfac0=100)
      integer*4 l,m,ll,i,mp
      real*8 t,z,p,sgn,pmm,pmmp,po
      real*8 dm,d0,df
      common/dpar/dm,d0,df(nfac0)
      real*8 f0,ff
      common/npar/f0,ff(nfac0)
      real*8,allocatable::pp(:)
      if(l.eq.0.and.m.eq.0)then
       p=1.0d0
       return
      endif
      z=cos(t)
c     pp is p for positive m:
      mp=abs(m)
c     Pm,m
      pmm=df(2*mp-1)*(1-z**2)**(dble(mp)/2.0d0)
      if(mp.eq.l)then
       po=pmm
      else
c      Pm,m+1
       pmmp=z*dble(2*mp+1)*pmm
       if(mp+1.eq.l)then
        po=pmmp
       else
        allocate(pp(l-mp+1))
        pp(1)=pmm
        pp(2)=pmmp
        do 1 ll=mp+2,l
        i=ll-mp+1
1       pp(i)=
     1  (z*dble(2*ll-1)*pp(i-1)-dble(ll+mp-1)*pp(i-2))/dble(ll-mp)
        po=pp(l-mp+1)
       endif
      endif
      if(m.gt.0)then
       p=po
      else
       if(mod(mp,2).eq.0)then
        sgn= 1.0d0
       else
        sgn=-1.0d0
       endif
       p=sgn*ff(l-mp)/ff(l+mp)*po
      endif
      return
      end

      function n(l,m)
      implicit none
      real*8 n,fp
      integer*4 l,m
      integer*4 nfac0
      parameter (nfac0=100)
      real*8 f0,factorial
      common/npar/f0,factorial(nfac0)
      fp=4.0d0*4.0d0*datan(1.0d0)
c     n=dsqrt(dble(2*l+1)*factorial(l-abs(m))/fp/factorial(l+abs(m)))
      n=dsqrt(dble(2*l+1)*factorial(l-m)/fp/factorial(l+m))
      return
      end

      Function Fac(n)
      implicit none
      integer*4 n,ii
      real*8 Fac,bb
      bb=1.0d0
      do 1 ii=2 , n
1     bb=bb*dble(ii)
      fac=bb
      return
      end

      Function DFac(n)
c     double factorial = n!!
      implicit none
      integer*4 n,ii
      real*8 DFac,bb
      if(n.eq.-1.or.n.eq.0)then
       bb=1.0d0
      else
       bb=dble(n)
       if(mod(n,2).ne.0)then
        do 1 ii=n-2,1,-2
1       bb=bb*dble(ii)
       else
        do 2 ii=n-2,2,-2
2       bb=bb*dble(ii)
       endif
      endif
      dfac=bb
      return
      end

      subroutine Initfac
      implicit none
      integer*4 nfac0,i
      parameter (nfac0=100)
      real*8 f0,factorial,fac
      common/npar/f0,factorial(nfac0)
      real*8 dm,d0,df,dfac
      common/dpar/dm,d0,df(nfac0)
      do  1 i=0 , nfac0
1     factorial(i)=fac(i)
      do  2 i=-1, nfac0
2     df(i)=dfac(i)
      return
      end