-- test_deriv.adb  compute formulas for numerical derivatives
--             order is order of derivative, 1 = first derivative, 2 = second
--             points is number of points where value of function is known
--             f(x1), f(x2), f(x2) ... f(x points)
--             term is the point where derivative is computer
--             f'(x1) = (a(1)*f(x1) + a(2)*f(x2)
--                       + ... + a(points)*f(x points)) /h^order

with Ada.text_io;
use  Ada.text_io;
with Ada.Numerics;
use  Ada.Numerics;
with Real_Arrays;
use  Real_Arrays;
with Deriv;

procedure test_deriv is
  package realio is new Float_IO(Real);
  use realio;

  norder : integer := 6;
  npoints : integer := 12;
  a : real_vector(1..50);         -- coefficients of f(x0), f(x1), ...

  function pow(x : real; pwr : integer) return real is
    val : real := 1.0;
  begin
    for i in 0..pwr-1 loop
      val := val * x;
    end loop;
    return val;
  end pow;

  procedure check(order : integer;
                  npoints : integer;
                  point : integer;
                  a : real_vector) is
    pwr : integer; -- simple test using x^pwr
    isum : Real;
    sum, x, term, err, xpoint, aderiv : real;
    h : real;
  begin
    for degree in npoints-2..npoints+1 loop
      h := 0.125;
      if order>3 then h := h/2.0; end if; -- cumulative
      if order>4 then h := h/2.0; end if;
      if order>5 then h := h/2.0; end if;
      if order>6 then h := h/2.0; end if;
      pwr := degree;
      sum := 0.0;
      x := 1.0;
      for i in 1..npoints loop
        term := a(i)*pow(x, pwr);
        sum := sum + term;
        x := x + h;
      end loop;
      sum := sum / pow(h, order);

      xpoint := 1.0+real(point-1)*h;
      aderiv := 1.0; -- exact derivative at point
      for i in 0..order-1 loop
        aderiv := aderiv * real(pwr);
        pwr := pwr -1;
      end loop;
      aderiv := aderiv * pow(xpoint, pwr);
      err := aderiv - sum;
      put("          err=" & real'image(err) & ", h=");
      put(h,1,4,0);
      put(", aderiv=");
      put(aderiv,4,1,0);
      put(", x=");
      put(xpoint,1,4,0);
      Put(", pwr=" & integer'image(degree));
      new_line;
    end loop;
    isum := 0.0;
    for i in 1..npoints loop
      isum := isum + a(i);
    end loop;
    if(abs(isum)>1.0e-12) then
      put_Line("ERROR  sum a(i) /= 0.0, is "&real'image(isum));
    end if;
  end check;

begin -- test_deriv

  put_line("Formulas for numerical computation of derivatives ");
  put_line("  f^(3) means third derivative ");
  put_line("  f^(2)(x(1)) means second derivative computed at point x(1) ");
  put_line("  f(x(2)) means function evaluated at point x(2) ");
  put_line("  Points are x(0), x(1)=x(0)+h, x(2)=x0+2h, ... ");
  put_line("  The error term gives the power of h and ");
  put_line("  derivative order with z chosen for maximum value.");
  put_line("  The error terms, err=, are for test polynomial sum(x^pwr) ");
  put_line("  with h = 0.125 or smaller, and order = 'pwr' ");

  for order in 1..norder loop
    npoints := npoints -2; -- fewer for larger derivatives
    if npoints<=order+1 then npoints := order + 2; end if;
    for points in order+1..npoints loop
      for term in 1..points loop
        put_line("computing order=" & integer'image(order) &
                 ", npoints=" & integer'image(points) &
                 ", at term=" & integer'image(term));
        deriv(order, points, term, a);

        for jterm in 1..points loop
          if jterm=1 then
            put_line("f^(" & integer'image(order) &
                     ")(x(" & integer'image(term) &
                     ")= " & real'image(a(jterm)) &
                     " * f(x(" & integer'image(jterm) & ")) + ");
          elsif jterm=points then
            put_line("              " & real'image(a(jterm)) &
                     " * f(x(" & integer'image(jterm) &
                     ") ) + O(h^" & integer'image(points-order) &
                     ")f^(" & integer'image(points) & ")(z) ");
          else
            put_line("              " & real'image(a(jterm)) &
                     " * f(x(" & integer'image(jterm) & ") + ");
          end if;
        end loop; -- jterm
        check(order, points, term, a);
        new_line;
      end loop; -- term
      new_line;
    end loop;  -- points
    new_line;
  end loop; -- order
  -- extra checks
  put_line("numerical accuracy fails at 4,?  multiple precision required");
  new_line;
  put_line("order=4, points=9, term=1 ");
  deriv(4, 9, 1, a);
  check(4, 9, 1, a);
  put_line("order=4, points=9, term=2 ");
  deriv(4, 9, 2, a);
  check(4, 9, 2, a);
  put_line("order=4, points=10, term=1 ");
  deriv(4, 10, 1, a);
  check(4, 10, 1, a);
  put_line("order=4, points=10, term=2 ");
  deriv(4, 10, 2, a);
  check(4, 10, 2, a);
  new_line;
  put_line("order=4, points=11, term=1 ");
  deriv(4, 11, 1, a);
  check(4, 11, 1, a);
  put_line("order=4, points=11, term=2 ");
  deriv(4, 11, 2, a);
  check(4, 11, 2, a);
  put_line("order=4, points=12, term=1 ");
  deriv(4, 12, 1, a);
  check(4, 12, 1, a);
  put_line("order=4, points=12, term=2 ");
  deriv(4, 12, 2, a);
  check(4, 12, 2, a);
  put_line("order=4, points=13, term=1 ");
  deriv(4, 13, 1, a);
  check(4, 13, 1, a);
  put_line("order=4, points=13, term=2 ");
  deriv(4, 13, 2, a);
  check(4, 13, 2, a);
  put_line("order=4, points=14, term=1 ");
  deriv(4, 14, 1, a);
  check(4, 14, 1, a);
  put_line("order=4, points=14, term=2 ");
  deriv(4, 14, 2, a);
  check(4, 14, 2, a);
  put_line("many formulas checked against Abramowitz and Stegun,");
  put_line("Handbook of Mathematical Functions Table 25.2");
end test_deriv;