# test_gaulegf.jl     from 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
#
using Printf

function gaulegf(x1, x2, x, w, n)
  # n = size(x1,2) = size(x2,1) = size(x2,2) = size(x,1) = size(w,1)
  # integer :: i, j, m
  # double precision :: p1, p2, p3, pp, xl, xm, z, z1
  pp = 1.0
  eps=3.0e-14
      
  m = (n+1)/2
  xm = 0.5*(x2+x1)
  xl = 0.5*(x2-x1)
  for i in 1:m
    z = cos(3.141592654*(i-0.25)/(n+0.5))
    z1 = 0.0
    while abs(z-z1)>eps
      p1 = 1.0
      p2 = 0.0
      for j in 1:n
        p3 = p2
        p2 = p1
        p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
      end # j
      pp = n*(z*p1-p2)/(z*z-1.0)
      z1 = z
      z = z1 - p1/pp
    end # end while
    x[Int(i)] = xm - xl*z
    x[Int(n+1-i)] = xm + xl*z
    w[Int(i)] = (2.0*xl)/((1.0-z*z)*pp*pp)
    w[Int(n+1-i)] = w[Int(i)]
  end # end i
  return nothing
end # end gaulegf

function vecprt(name,V) # works for row and column vector
  n = size(V,1)
  if n>1
    for i in 1:n
      @printf "%s[%d]=%f \n" name i  V[i]
    end
    println(" ")
  else
    n = size(V,2)
    for i in 1:n
      @printf "%s[1,%d]=%f \n" name i  V[1,i]
    end
    println(" ")
  end # end if else
  return nothing
end # end vecprt

function f(x)
  return x*x + 2.0*x + 3.0
end # end f(x)

println("test_gaulegf running")
n = 5
x1 = 0.0
x2 = 1.0
@printf "n=%d  x1=%f  x2=%f \n" n x1 x2
x = zeros(Float64, n)
w = zeros(Float64, n)

gaulegf(x1, x2, x, w, n)
vecprt("x",x)
vecprt("w",w)
println("integral=sum i in 1:n w[i]*f(x[i])")
println(" ")

intf = 0.0
for i in 1:n
  global intf
  intf = intf + w[i]*f(x[i])
end 
@printf "intf=%f \n" intf

n = 6
x1 = -1.0
x2 = 3.0
@printf "n=%d  x1=%f  x2=%f \n" n x1 x2
x = zeros(Float64, n)
w = zeros(Float64, n)
gaulegf(x1, x2, x, w, n)
vecprt("x",x)
vecprt("w",w)

intf = 0.0
for i in 1:n
  global intf
  intf = intf + w[i]*f(x[i])
end 
@printf "intf=%f \n" intf

println("test_gaulegf finished")
# end test_gaulegf.jl