gauss_legendre_unit Subroutine

  subroutine gauss_legendre_unit(order, node, weight)
    integer(i32), intent(in) :: order
    real(dp), allocatable, intent(out) :: node(:), weight(:)
    integer :: i, j, midpoint
    real(dp) :: z, z_previous, polynomial, derivative, p0, p1, p2

    if (order < 2_i32) error stop 'panel oracle quadrature order must be >= 2.'
    allocate (node(order), weight(order))
    midpoint = (order + 1)/2
    do i = 1, midpoint
      z = cos(pi*(real(i, dp) - 0.25_dp)/(real(order, dp) + 0.5_dp))
      do
        p0 = 1.0_dp
        p1 = z
        do j = 2, order
          p2 = ((2.0_dp*real(j, dp) - 1.0_dp)*z*p1 - (real(j, dp) - 1.0_dp)*p0)/real(j, dp)
          p0 = p1
          p1 = p2
        end do
        polynomial = merge(1.0_dp, p1, order == 0_i32)
        derivative = real(order, dp)*(z*p1 - p0)/(z*z - 1.0_dp)
        z_previous = z
        z = z_previous - polynomial/derivative
        if (abs(z - z_previous) <= 8.0_dp*epsilon(1.0_dp)) exit
      end do
      node(i) = 0.5_dp*(1.0_dp - z)
      node(order + 1 - i) = 0.5_dp*(1.0_dp + z)
      weight(i) = 1.0_dp/((1.0_dp - z*z)*derivative*derivative)
      weight(order + 1 - i) = weight(i)
    end do
  end subroutine gauss_legendre_unit