gaussquad_rosetta Subroutine

    subroutine gaussquad_rosetta(n, r1, r2)

        ! This code was originally found at the following website:
        !  http://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature#Fortran

        integer,  intent(in)    :: n
        real(wp), intent(out)   :: r1(n), r2(n)


        integer                 :: k
        real(wp), parameter     :: pi = 4._wp*atan(1._wp)
        real(wp)                :: x, f, df, dx
        integer                 :: i,  iter
        real(wp), allocatable   :: p0(:), p1(:), tmp(:)

        p0 = [1._wp]
        p1 = [1._wp, 0._wp]

        do k = 2, n
            tmp = ((2*k-1)*[p1,0._wp]-(k-1)*[0._wp, 0._wp,p0])/k
            p0 = p1; p1 = tmp
        enddo

        outer: do i = 1, n
            x = cos(pi*(i-0.25_wp)/(n+0.5_wp))
            inner: do iter = 1, 10
                f = p1(1); df = 0._wp
                deep: do k = 2, size(p1)
                    df = f + x*df
                    f  = p1(k) + x * f
                enddo deep
                dx =  f / df
                x = x - dx
                if (abs(dx)<10*epsilon(dx)) exit
            enddo inner
            r1(i) = x
            r2(i) = 2/((1-x**2)*df**2)
        enddo outer
        return
    end subroutine gaussquad_rosetta