c===========================================================
c     Demonstrates sample use of SLATEC routine dxlegf
c     to computs normalized Legendre polynomials and 
c     associated Legendre functions.
c
c     % maple
c     > Digits := 16:
c     > evalf();
c
c     SUBROUTINE DXLEGF (DNU1, NUDIFF, MU1, MU2, THETA, ID, 
c    &                   PQA, IPQA, IERROR)
c===========================================================
      program     tdxlegf

      implicit    none

      integer     iargc
      integer     i4arg
      logical     ltrace

      integer     maxrange 
      parameter ( maxrange = 100 )

      real*8      pqa(maxrange)
      integer     ipqa(maxrange)
      integer     nudiff, mu1, mu2, id
      real*8      dnu1
      integer     ierror
   
      real*8      theta, dtheta
      integer     i, ntheta

c-----------------------------------------------------------
c     Argument parsing.
c-----------------------------------------------------------
      if( iargc() .lt. 1 ) go to 900
      dnu1 = i4arg(1,-99)
      if( dnu1 .eq. -99 ) go to 900
      ntheta = i4arg(2,101)
      ltrace = iargc() .gt. 2

      if( ltrace ) then
         write(0,*) 'tdxlegf: nu = ', dnu1
         write(0,*) 'tdxlegf: ntheta = ', ntheta
      end if
c-----------------------------------------------------------
c     Compute and output Legendre polynomial, case 4a with 
c     mu1 = mu2 = 0, order of polynomial is nu = dnu1.
c-----------------------------------------------------------
      nudiff  = 0 
      mu1     = 0 
      mu2     = 0 
      id      = 4

      ntheta  = 101
      dtheta  = 2.0d0 * atan(1.0d0)  / (ntheta - 1)
      
      theta   = dtheta
      do i = 1 , ntheta - 1
         call DXLEGF (DNU1, NUDIFF, MU1, MU2, THETA, ID,
     &                   PQA, IPQA, IERROR) 
         if( ltrace ) then
            write(0,*) 'tdxlegf: ierror  = ', ierror
            write(0,*) 'tdxlegf: ipqa(1) = ', ipqa(1)
            write(0,*) 
         end if
         if( ierror .eq. 0  .and. ipqa(1) .eq. 0 ) then
            theta = theta + dtheta
            write(*,*) cos(theta), pqa(1)
         end if
      end do

      stop

 900  continue
         write(0,*) 'usage: tdxlegf <nu> [<ntheta> <trace>]'
      stop

      end
