-- aquad.adb  package body
package body aquade is

-- aquad3  from book page 301, with modifications
-- integral :=  aquad3(f'access, 0.1, 2.0, 0.001);

  function aquad3(f : funct ; xmin, xmax, eps : Real;
                  i,j,n:Integer; xg:Real_Vector)
               return real is
  maxtop : Integer := 40;
  type stack_type is array (integer range <>, integer range <>) of real;
  stack : stack_type(1..maxtop,1..7);  -- a place to store/retrieve
  top : Integer;                       -- top of stack
  a, b, c, d, e : real;
  Fa, Fb, Fc, Fd, Fe : real;
  h1, h2, hs : real;
  Sab, Sac, Scb, S2ab : real;
  tol : real; -- eps
  val, value : real;
  ns : Integer := 9;

  procedure store(s1, s2, s3, s4, s5, s6, S7: real) is
  begin
    stack(top,1) := S1;
    stack(top,2) := S2;
    stack(top,3) := S3;
    stack(top,4) := S4;
    stack(top,5) := S5;
    stack(top,6) := S6;
    stack(top,7) := S7;
  end store;

  procedure retrieve(s1, s2, s3, s4, s5, s6, S7 : out real) is
  begin
    s1 := stack(top,1);
    s2 := stack(top,2);
    s3 := stack(top,3);
    s4 := stack(top,4);
    s5 := stack(top,5);
    s6 := stack(top,6);
    s7 := stack(top,7);
  end retrieve;

  function Sn(F0, F1, F2, h: real) return real is
  begin
    return h*(F0 + 4.0*F1 + F2)/3.0;
    -- error term 1/12 h^3 f^(2)(c)
  end Sn;

  function RS(F0, F1, F2, F3, F4, h : real) return real is
  begin
    return h*(14.0*F0 +64.0*F1 + 24.0*F2 + 64.0*F3 + 14.0*F4)/45.0;
    -- error term  8/945  h^7 f^(8)(c)
  end RS;

begin

  top := 1;
  value := 0.0;
  tol := eps;
  a := xmin;
  hs := (xmax-xmin)/Real(ns);
  b := a + hs;
  for k in 1..ns loop
    h1 := (b-a)/2.0;
    c := a + h1;
    Fa := f(a,i,j,n,xg);
    Fc := f(c,i,j,n,xg);
    Fb := f(b,i,j,n,xg);
    Sab := Sn(Fa, Fc, Fb, h1);
    store(a, Fa, Fc, Fb, h1, tol, Sab);
    top := top+1;
    a := b;
    b := a + Hs;
  end loop;
  while top > 1 loop
    top := top-1;
    retrieve(a, Fa, Fc, Fb, h1, tol, Sab);
    c := a + h1;
    b := a + 2.0*h1;
    h2 := h1/2.0;
    d := a + h2;
    e := a + 3.0*h2;
    Fd := f(d,i,j,n,xg);
    Fe := f(e,i,j,n,xg);
    Sac := Sn(Fa, Fd, Fc, h2);
    Scb := Sn(Fc, Fe, Fb, h2);
    S2ab := Sac + Scb;
    if abs(S2ab-Sab) < tol then
      val := RS(Fa, Fd, Fc, Fe, Fb, h2);
      value := value + Val;
    else
      h1 := h2;
      tol := tol/2.0;
      store(a, Fa, Fd, Fc, h1, tol, Sac);
      top := top+1;
      store(c, Fc, Fe, Fb, h1, tol, Scb);
      top := top+1;
    end if;
    if top >= maxtop then exit; end if;
  end loop;
  return value;
end aquad3;

end aquade; -- package body