-- gaulegf.adb     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

with Ada.Text_IO;
use Ada.Text_IO;
with Ada.Numerics.Long_Elementary_Functions;
use Ada.Numerics.Long_Elementary_Functions;
with Real_Arrays;
use Real_Arrays;

procedure gaulegf(x1 : real; x2 : real;
                  x : in out Real_Vector; w : in out Real_Vector;
                  n : integer) is
  m : integer;
  eps : real := 3.0E-14;
  p1, p2, p3, pp, xl, xm, z, z1 : real;
begin
  m := (n+1)/2;
  xm := 0.5*(x2+x1);
  xl := 0.5*(x2-x1);
  for i in 1..m loop
    z := cos(3.141592654*(real(i)-0.25)/(real(n)+0.5));
    loop
      p1 := 1.0;
      p2 := 0.0;
      for j in 1..n loop
        p3 := p2;
        p2 := p1;
        p1 := ((2.0*real(j)-1.0)*z*p2-(real(j)-1.0)*p3)/real(j);
      end loop;
      pp := real(n)*(z*p1-p2)/(z*z-1.0);
      z1 := z;
      z := z1 - p1/pp;
      exit when abs(z-z1) <= eps;
    end loop;
    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 loop;
end gaulegf;