! adapted from www.cs.umbbc.edu/~squire/download/gauleg.f90
! gauleg.f90     P145 Numerical Recipes in Fortran
! compute x(i) and w(i)  i=1,n  Legendre ordinates and weights
! on interval -1.0 to 1.0 (length is 2.0)
! use ordinates and weights for Gauss Legendre integration
! MODIFIED April 15, 2020 by Riccardo Rende , Maria Peressi
!
subroutine gaulegf(x1, x2, x, w, n)
  implicit none
  integer, intent(in) :: n
  integer, parameter :: dp=kind(1.d0)
  real(dp), parameter :: pigr=acos(-1.d0)
  real(dp), intent(in) :: x1, x2
  real(dp), dimension(n), intent(out) :: x, w
  integer :: i, j, m, k
  real(dp) :: p1, p2, p3, pp, xl, xm, z, z1
  real(dp), parameter :: eps=3.d-14

  m = (n+1)/2
  xm = 0.5_dp*(x2+x1)
  xl = 0.5_dp*(x2-x1)
  do i=1,m
   z = cos(pigr*(i-0.25_dp)/(n+0.5_dp))
!  z = cos(3.14159265358979d0*(i-0.25d0)/(n+0.5d0))
!  z = cos(3.1415926535897931d0*(i-0.25d0)/(n+0.5d0))
!  z = cos(3.141592653589793238d0*(i-0.25d0)/(n+0.5d0))  ! no ???
  z1 = 0.0_dp
!    do while(abs(z-z1) > eps)  ! commento questo
   do k = 1, 200 ! aggiunto questo
      p1 = 1.0_dp
      p2 = 0.0_dp
      do j=1,n
        p3 = p2
        p2 = p1
        p1 = ((2.0_dp*j-1.0_dp)*z*p2-(j-1.0_dp)*p3)/j
      end do
      pp = n*(z*p1-p2)/(z*z-1.0_dp)
      z1 = z
      z = z1 - p1/pp
            if (abs(z-z1) < eps) exit   ! aggiunto questo
            if (k == 200) then                                     !aggiunto
            print*,' convergence not reached, stopping....'; STOP  !aggiunto
            end if                                                 !aggiunto
    end do
    x(i) = xm - xl*z
    x(n+1-i) = xm + xl*z
    w(i) = (2.0_dp*xl)/((1.0_dp-z*z)*pp*pp)
    w(n+1-i) = w(i)
  end do
end subroutine gaulegf

program gauleg
  implicit none
  integer :: i, j
  integer, parameter :: dp=kind(1.d0)
  real(dp), dimension(100) :: x, w
  real(dp) :: sum, a, b
  integer, parameter :: debug=0

  print *, 'test gauleg.f90 on interval -1.0 to 1.0 ordinates, weights'
  do i=1,15
    call gaulegf(-1.0_dp, 1.0_dp, 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 *, '             integrate(1.0, from -1.0 to 1.0)= ', sum
    print *, ' '
  end do

end program gauleg