get_legendre_grids_from_cheb Subroutine

  subroutine get_legendre_grids_from_cheb (x1, x2, zero, wgt)
    use constants, only: pi => dpi
    real, intent (in) :: x1, x2
    real, dimension (:), intent (out) :: zero, wgt
    integer :: i, nn, nh
    double precision :: xold, xnew, pold, pnew
    double precision, dimension (:), allocatable :: zz, ww

    nn = size(zero)
    nh = (nn+1)/2
    allocate (zz(nh))
    allocate (ww(nh))

    ! search zero from 1 to 0 using chebyshev grid points
    ! this is O(nn^2) operations
    xold = cos(pi / (2*(nn+1)))
    pold = legendre_p(nn,xold)
    do i=1, nh
       xnew = cos(pi*(2*i+1)/(2.0*(nn+1)))
       pnew = legendre_p(nn,xnew)
       call find_zero_bisect_newton (nn, xold, xnew, pold, pnew, zz(i))
       xold = xnew
       pold = pnew
    end do
    ! invert them to give zeros in (-1,0]
    zz(1:nn/2) = -zz(1:nn/2)

    ! weight from formula
!    ww = dble(2.0) / (dble(1.0) - zz**2) / legendre_pp(nn,zz,dble(0.0))**2
    ww = dble(2.0) / (dble(1.0) - zz**2) / legendre_pp(nn,zz)**2

    ! rescale (x2 may be smaller than x1)
    zero(1:nh) = real(zz * (x2-x1) + (x1+x2)) / 2.0
    zero(nh+1:nn) = real(zz(nn/2:1:-1) * (x1-x2) + (x1+x2)) / 2.0
    wgt(1:nh) = real(ww * abs(x2-x1) / 2.0)
    wgt(nh+1:nn) = wgt(nn/2:1:-1)

    deallocate (zz)
    deallocate (ww)

    ! roundoff correction
!!$    if (abs(sum(wgt)/abs(x2-x1)) - 1.0 > epsilon(wgt)) then
!!$       print *, 'roundoff correction occurred'
    if (weight_roundoff_correction) then
       if (mod(nn,2)==0) then
          wgt(nh) = abs(x2-x1) / 2.0 - sum(wgt(1:nh-1))
          wgt(nh+1) = wgt(nh)
       else
          wgt(nh) = abs(x2-x1) - sum(wgt(1:nh-1)) * 2.0
       end if
    end if

    call check_legendre_zero (x1, x2, zero)
    call check_legendre_weights (abs(x2-x1), wgt)

    if (debug) then
       print *, 'get_legendre_grids_from_cheb: sum of weights = ', sum(wgt)
       print *, 'get_legendre_grids_from_cheb: section length = ', abs(x2-x1)
    end if

  end subroutine get_legendre_grids_from_cheb