get_laguerre_grids Subroutine

  subroutine get_laguerre_grids (zero, wgt)
    use mp, only: mp_abort
    use file_utils, only: error_unit
    implicit none
    real, dimension (:), intent (out) :: zero
    real, dimension (:), intent (out) :: wgt
    logical :: error=.false.
    integer :: i, j, n, nzero
    double precision :: x, delx, pold, pnew
    double precision, dimension (:), allocatable :: zz

    n = size (zero)
    allocate (zz(n))
    zz = 0.0

    if (n > 180) then
       write (error_unit(),*) 'ERROR: can''t get so many laguerre grid points'
       write (error_unit(),*) 'ERROR: size(zero)= ', n
       error = .true.
    end if

    ! search zero from 0 to xmax using evenly spaced grid points
    if (n==1) then

       zero(1) = 1.0
       wgt(1) = 1.0
       return

    else

       nzero = 0
       pold = laguerre_l(n,dble(0.0))
       delx = 0.001
       do i=1, 1000      ! up to x=1
          x = delx*i
!          print*, x
          pnew = laguerre_l(n,x)
!          if (pold*pnew < epsilon(0.0)) then
          if (pold*pnew < 0.0) then
             nzero = nzero+1
             call find_zero (n, x-delx, x, pold, pnew, zz(nzero))
          end if
          pold=pnew
       end do

       do j=0, 3
          delx = delx * 10.
          do i=1, 900
             x = delx*(i+100)
!             print*, x
             pnew = laguerre_l(n,x)
             if (pold*pnew < 0.0) then
                nzero = nzero+1
                call find_zero (n, x-delx, x, pold, pnew, zz(nzero))
             end if
             if (nzero == n) exit
             pold=pnew
          end do
          if (nzero == n) exit
       end do

    end if

    zero = zz
    wgt = real(zz / (n+1)**2 / laguerre_l(n+1,zz)**2)

    deallocate (zz)

    ! roundoff correction
    if (weight_roundoff_correction) then
       i = sum(maxloc(wgt))
       wgt(i) = 1.0 - sum(wgt(1:i-1)) - sum(wgt(i+1:n))
    end if

    ! check number of found zeros
    if (nzero < n) then
       write (error_unit(),*) 'ERROR in laguerre: didn''t find all zeros'
       do i=1, n
          write (error_unit(),*) i, zero(i)
       end do
       call mp_abort("ERROR in laguerre: didn't find all zeros")
    end if

    if (error) call mp_abort ('STOP in get_laguerre_grids')

  end subroutine get_laguerre_grids