-- gauleg.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;
procedure gauleg is
  sum, a, b : long_float;
  type Vector is array (integer range <>) of long_float;
  x, w : Vector(1..100);

  procedure gaulegf(x1 : long_float; x2 : long_float;
                    x : in out Vector; w : in out Vector;  n : integer) is
    m : integer;
    eps : long_float := 3.0E-14;
    p1, p2, p3, pp, xl, xm, z, z1 : long_float;
  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*(long_float(i)-0.25)/(long_float(n)+0.5));
      loop
        p1 := 1.0;
        p2 := 0.0;
        for j in 1..n loop
          p3 := p2;
          p2 := p1;
          p1 := ((2.0*long_float(j)-1.0)*z*p2-(long_float(j)-1.0)*p3)/
                long_float(j);
        end loop;
        pp := long_float(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;
begin
  put_line("test gauleg.adb on interval -1.0 to 1.0  ordinates, weights");
  for i in 1..15 loop
    gaulegf(-1.0, 1.0, x, w, i);
    sum := 0.0;
    for j in 1..i loop
      put_line("x(" & integer'image(j) & ")=" & long_float'image(x(j)) &
               "w(" & integer'image(j) & ")=" & long_float'image(w(j)));
      sum := sum + w(j);
    end loop;
    put_line("              integral(1.0, -1.0..1.0)=" & long_float'image(sum));
    new_line;
  end loop;

  a := 0.5;
  b := 1.0;
  for i in 2..10 loop
    gaulegf(a, b, x, w, i);
    sum := 0.0;
    for j in 1..i loop
      sum := sum + w(j)*sin(x(j));
    end loop;
    put_line("integral (0.5,1.0) sin(x) dx = " & long_float'image(sum));
  end loop;
  put_line("-cos(1.0)+cos(0.5) = " & long_float'image(-cos(1.0)+cos(0.5)));
  put_line("Maple says 3.372802560E-01");
  new_line;

  a := 0.5;
  b := 5.0;
  for i in 2..10 loop
    gaulegf(a, b, x, w, i);
    sum := 0.0;
    for j in 1..i loop
      sum := sum + w(j)*exp(x(j));
    end loop;
    put_line("integral (0.5,5.0) exp(x) dx = " & long_float'image(sum));
  end loop;
  put_line("exp(5.0)-exp(0.5) = " & long_float'image(exp(5.0)-exp(0.5)));
  put_line("Maple says 1.467644378E+02");
  new_line;

  a := 0.5;
  b := 5.0;
  for i in 2..30 loop
    gaulegf(0.5, 5.0, x, w, i);
    sum := 0.0;
    for j in 1..i loop
      sum := sum + w(j)*(((x(j)**x(j))**x(j))*(x(j)*(log(x(j))+1.0)+x(j)*log(x(j))));
    end loop;
    put_line("integral (0.5,5.0) mess(x) dx = " & long_float'image(sum));
  end loop;
  put_line("((5.0**5.0)**5.0)-(0.5**0.5)**0.5 = " &
           long_float'image((5.0**5.0)**5.0-(0.5**0.5)**0.5));
  put_line("Maple says 2.980232239E+17");
  new_line;
  put_line("Done.");
end gauleg;