C **********************************************************************
C *                                                                    *
C *  LEGENDRE and SHERICAL-HARMONICS FUNCTIONS (legendre.f)            *
C *                                                                    *
C * Library module containing the following functions and subroutines: *
C *                                                                    *
C *   plgndr:  computes the associated Legendre polynomial P_l^m (x).  *
C *            Here m and l are integers satisfying 0<=m<=l, while x   *
C *            lies in the range -1<=x<=1   (AJ, 8/96)                 *
C *   dlegendre: computes P_n(cos(theta)) where P_n(x) is the Legendre *
C *              polynomial of order n.  (MB, 8/96)                    *
C *   sphharm: Computes the spherical harmonics  Y_l^m (theta,phi).    *
C *            Here m and l are integers satisfying 0<=l and -l<=m<=+l *
C *            general spherical harmonics function with               *
C *            Condon-Shortley phase. (MCH 97)                         *
C *   algndr:  computes the associated Legendre polynomial P_l^m (x).  *
C *            Similar to plgndr but computes an arry containing all   *
C *            P_l^m (x) for m<=l<=lmax. This subroutine computes      *
C *            the L^2-mormalized functions, if switch=.true. is set.  *
C *                                                                    *
C **********************************************************************


C-----------------------------------------------------------------------
C                FUNCTION   PLGNDR
C
C PURPOSE:
C   Computes the associated Legendre polynomial P_l^m (x). Here m and l
c   are integers satisfying 0<=m<=l, while x lies in the range
C   -1<=x<=1
C INPUT PARAMETERS:
C   l: integer (see above)
C   m: integer (see above)
C   x: real*8  (see above)
C FUNCTION VALUE:
C     value of the "plgndr" function as defined above
C
C PHASE CONVENTION:
C   The phase convention differs by a factor (-1)**m from most textbooks
C   as e.g. Messiah (Quantum Mechanics) or Zare (Angular Momentum).
C-----------------------------------------------------------------------

      function plgndr(l,m,x)

      implicit none 

      integer l,m
      real*8  plgndr,x

      integer i,ll
      real*8  fact,pll,pmm,pmmp1,somx2

      if (m.lt.0.or.m.gt.l.or.abs(x).gt.1) then
         write(6,*) 'ERROR: bad arguments in plgndr ',m,l,x
         write(2,*) 'ERROR: bad arguments in plgndr ',m,l,x
         stop
      endif

C --- Compute P_m^m ---

      pmm=1.
      if (m.gt.0) then
         somx2=sqrt((1.d0-x)*(1.d0+x))
         fact=1.
         do i=1,m
            pmm=-pmm*fact*somx2
            fact=fact+2.d0
         enddo
      endif
      if (l.eq.m) then
         plgndr=pmm
      else

C --- Compute P_m+1^m ---

         pmmp1=x*(2*m+1)*pmm
         if (l.eq.m+1) then
            plgndr=pmmp1
         else 

C --- Compute P_l^m, l>m+1 ---

            pll=0.0
            do ll=m+2,l
               pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m)
               pmm=pmmp1
               pmmp1=pll
            enddo
            plgndr=pll
         endif
      endif

      return
      end

C ----------------------------------------------------------------------
C                            FUNCTION DLEGENDRE
C
C PURPOSE:
C     Computes the value P_n(cos(theta)) of the Legendre polynomial 
C                 P_n(x) := 1/(2^n n!) d^n/dx^n (x^2-1)^n
C     using the recursion relation
C               nP_n(x) = (2n-1)xP_(n-1)(x) - (n-1)P_(n-2)(x),
C                        P_1(x) = x,  P_0(x) = 1.
C     The definition of the Legendre polynomials implies that
C           Integral_(-1)^(+1) P_n(x) P_m(x) dx = 2/(2n+1) delta_mn.
C     All calculations are performed with double precision. Note that
C     for |x| > 1 there is no guaranty that the above recurrence
C     relation is stable; the special form of the argument, cos(theta),
C     ensures that the routine works correctly for any parameter theta.
C INPUT PARAMETERS:
C     n,theta: parameters as stated above
C FUNCTION VALUE:
C     value of the "dlegendre" function as defined above
C ----------------------------------------------------------------------

      real*8 function dlegendre (n,theta)

      implicit none

      real*8  theta,x,old,veryold
      integer n,i

C --- COMPUTE FUNCTION VALUE ---

      x = cos(theta)
      dlegendre = 0.0
      if (n .eq. 0) then
         dlegendre = 1.0d0
      elseif (n .eq. 1) then
         dlegendre = x
      elseif (n .ge. 2) then
         old = 1.0d0
         dlegendre = x
         do i = 2,n
            veryold = old
            old = dlegendre
            dlegendre = (dble(2*i-1)*x*old-dble(i-1)*veryold)/dble(i)
         enddo
      else
         print*,'ERROR: Invalid Legendre order in dlegendre :',n
         stop
      endif

      return
      end

C-----------------------------------------------------------------------
C                FUNCTION  SPHHARM
C
C PURPOSE:
C   Computes the spherical harmonics  Y_l^m (theta,phi). Here m and l
C   are integers satisfying 0<=l and -l<=m<=+l
C - general spherical harmonics function with Condon-Shortley phase -
C INPUT PARAMETERS:
C   l: integer (see above)
C   m: integer (see above)
C   theta,phi : real*8 
C FUNCTION VALUE:
C     value of the "sphharm" function as defined above
C CALL TO:
C function plgndr (l,m,x=cos(theta)) 
c-----------------------------------------------------------------------
      complex*16 function sphharm(l,m,theta,phi)

      implicit none
      real*8     theta,phi,pi,ass,plgndr,fac,x
      integer    i,im,m,l

      pi  = 4.d0*atan(1.d0)
      im  = iabs(m)
      x=cos(theta)

      fac=1.d0
      do i=l-im+1,l+im
         fac=fac*i
      enddo
      ass=plgndr(l,im,x)
      if (m.lt.0 .and. mod(im,2).eq.1 ) ass=-ass

      sphharm=sqrt((2*l+1)/(fac*4.d0*pi))*ass*cdexp(dcmplx(0.d0,m*phi))

      return
      end


C-----------------------------------------------------------------------
C                SUBROUTINE   ALGNDR
C
C PURPOSE:
C   Similar to PLGNDR, but stores all P_l^m (x) values in array plm.
C   l runs from m to lmax. 
C   Computes the associated Legendre polynomial P_l^m (x). Here m and l
c   are integers satisfying 0<=m<=l, while x lies in the range
C   -1<=x<=1
C INPUT PARAMETERS:
C   lmax: integer (see above)
C   m: integer (see above),  0<=m<=lmax
C   x: real*8  (see above),  x = cos(theta)
C   switch : logical , if .true. then compute L^2-mormalized functions.
C OUTPUT PARAMETERS:
C     Array plm(k), k=1,..,lmax+1-m.  l = m-1+k.
C
C PHASE CONVENTION:
C   The phase differs by a factor (-1)**m from most textbooks as
C   e.g. Messiah (Quantum Mechanics) or Zare (Angular Momentum).
C-----------------------------------------------------------------------

      subroutine algndr(lmax,m,x,plm,switch)

      implicit none 

      logical switch
      integer lmax,m,i,ll
      real*8  x,fact,pll,pmm,pmmp1,somx2
      real*8  plm(*)

      if (m.lt.0 .or. m.gt.lmax .or. abs(x).gt.1) then
         write(6,*) 'bad arguments in algndr  - stop -'
         write(2,*) 'bad arguments in algndr  - stop -'
         stop
      end if

C --- Compute P_m^m ---

      pmm=1.
      fact=0
      if (m.gt.0) then
         somx2=sqrt((1.-x)*(1.+x))
         fact=1.
         do i=1,m
            pmm=-pmm*fact*somx2
            fact=fact+2.
         enddo
      endif
      plm(1)=pmm
      if (switch)  then
         fact = 1.d0
         do i = 2, 2*m
            fact = fact*i
         end do
         plm(1)= plm(1)*sqrt((m+0.5d0)/fact)
      end if
      if (lmax.eq.m) return

C --- Compute P_m+1^m ---

      pmmp1=x*(2*m+1)*pmm
      plm(2)=pmmp1
      if (switch) then
         fact = fact*(2*m+1)
         plm(2)= plm(2)*sqrt((m+1.5d0)/fact)
      end if
      if (lmax.eq.m+1) return

C --- Compute P_l^m, l>m+1 ---

      do ll=m+2,lmax
         pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m)
         pmm=pmmp1
         pmmp1=pll
         plm(ll-m+1)=pll
      enddo

      if (switch) then
         do ll = m+2,lmax
            fact = fact*(ll+m)/dble(ll-m)
            plm(ll-m+1)=plm(ll-m+1)*sqrt((ll+0.5d0)/fact)
         end do
      end if
      return
      end

C ----------------------------------------------------------------------









