c gauleg.for     P145 Numerical Recipies in Fortran
c compute x(i) and w(i)  i=1,n  Legendre ordinates and weights
c on interval -1.0 to 1.0 (length is 2.0)
c use ordinates and weights for Gauss Legendre integration
c
      subroutine gaulegf(x1, x2, x, w, n)
      integer n, i, j, m
      double precision x1, x2, x(100), w(100), eps
      double precision p1, p2, p3, pp, xl, xm, z, z1
      parameter (eps=3.d-14)
      parameter (debug=0)
      
      m = (n+1)/2
      xm = 0.5d0*(x2+x1)
      xl = 0.5d0*(x2-x1)
      do i=1,m
        z = cos(3.141592654d0*(i-0.25d0)/(n+0.5d0))
1       continue
        p1 = 1.0d0
        p2 = 0.0d0
        do j=1,n
          p3 = p2
          p2 = p1
          p1 = ((2.0d0*j-1.0d0)*z*p2-(j-1.0d0)*p3)/j
        end do
        pp = n*(z*p1-p2)/(z*z-1.0d0)
        z1 = z
        z = z1 - p1/pp
        if (abs(z-z1) .gt. eps) goto 1
        x(i) = xm - xl*z
        x(n+1-i) = xm + xl*z
        w(i) = (2.0d0*xl)/((1.0d0-z*z)*pp*pp)
        w(n+1-i) = w(i)
      end do
      return
      end

      program gauleg
      integer i, j
      double precision x(100), w(100), sum, a, b

      print *, 'test gauleg.for on interval -1.0 to 1.0 ordinates',
     c ', weights'
      do i=1,15
        call gaulegf(-1.0d0, 1.0d0, x, w, i)
        sum = 0.0d0
        do j=1,i
          print *, 'x(',j,')=', x(j), '  w(',j,')=', w(j)
          sum = sum + w(j)
        end do
        print *, '                           sum= ', sum
        print *, ' '
      end do

      a = 0.5d0
      b = 1.0d0
      print *, 'test gauleg on integral(sin(x),',a,'..',b,')'
      do i=2,10
        call gaulegf(a, b, x, w, i)
        sum = 0.0d0
        do j=1,i
          if(debug>0) then
            print *, 'x(',j,')=', x(j), '  w(',j,')=', w(j)
          end if
          sum = sum +w(j)*sin(x(j))
        end do
        print *, i, ' integral (0.5,1.0) sin(x) dx = ', sum
      end do
      print *, 'Maple says 0.3372802560'
      print *, ' '

      a = 0.5d0
      b = 5.0d0
      print *, 'test gauleg on integral(exp(x),',a,'..',b,')'
      do i=2,15
        call gaulegf(a, b, x, w, i)
        sum = 0.0d0
        do j=1,i
          if(debug>0) then
            print *, 'x(',j,')=', x(j), '  w(',j,')=', w(j)
          end if
          sum = sum + w(j)*exp(x(j))
        end do
        print *, i, ' integral (0.5,5.0) exp(x) dx = ', sum
      end do
      print *, 'Maple says 146.7644378'
      print *, ' '

      a = 0.5d0
      b = 5.0d0
      print *, 'test gauleg on integral(mess(x),',a,'..',b,')'
      do i=2,30
        call gaulegf(a, b, x, w, i)
        sum = 0.0d0
        do j=1,i
          sum = sum + w(j)*(((x(j)**x(j))**x(j))*(x(j)*
     c                (dlog(x(j))+1.0d0)+x(j)*dlog(x(j))))
        end do
        print *, i, ' integral (0.5,5.0) mess(x) dx = ', sum
      end do
      print *, 'mess(x) = ((x^x)^x)*(x*(log(x)+1)+x*log(x)))'
      print *, 'Maple says 2.980232239E+17'
      print *, ' '
      print *, 'Done.'
      end