# test_gaulegf.rb  compute Gauss Legendre integration weights and ordinates
# integrate from x=a to x=b using n evaluations of the function f(x)
# usage: n = 20
#        a = -1.0
#        b = 1.0
#        x = Array.new(n+1) computed by gaulegf
#        w = Array.new(n+1)
#        gaulegf(a, b, x, w, n)
#        area = 0.0;
#        for i in 1..n         #   yes, 1..n 
#          area += w[i]*f(x[i])

def gaulegf(a, b, x, w, n)
  puts "in gaulegf, a=#{a}, b=#{b}, n=#{n}"
  eps = 3.0E-14;
  
  m = (n+1)/2;
  xm = 0.5*(b+a);
  xl = 0.5*(b-a);
  for i in 1..m
    z = Math.cos(3.141592654*(i-0.25)/(n+0.5))
    while 1==1
      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
      pp = n*(z*p1-p2)/(z*z-1.0)
      z1 = z
      z = z1 - p1/pp
      if (z-z1).abs <= eps
        break
      end
    end
    x[i] = xm - xl*z
    x[n+1-i] = xm + xl*z
    w[i] = 2.0*xl/((1.0-z*z)*pp*pp)
    w[n+1-i] = w[i]
  end
end # gaulegf

def f(p)       # test function to integrate
  return p*p;
end

puts "test_gauleg.rb running"
a = 1.0
b = 2.0
n = 4
area = 0.0
x = Array.new(n+1)
w = Array.new(n+1)

puts "calling gaulegf(1.0, 2.0, x, w, n)"
gaulegf(a, b, x, w, n)

puts "x[1]=#{x[1]}, x[2]=#{x[2]}"
puts "w[1]=#{w[1]}, w[2]=#{w[2]}"

for i in 1..n
  puts "x[#{i}]=#{x[i]}, w[#{i}]=#{w[i]}"
end
area = 0.0;
for i in 1..n 
  area += w[i]*f(x[i]);
end
puts "area=#{area}"
puts "error=#{area-(7.0/3.0)}"

puts "test_gaulegf.rb finished"