-- nuderiv_test_exp.adb  non uniformly spaced derivative coefficients

with Ada.Text_IO;
use  Ada.Text_IO;
with Ada.Numerics;
use  Ada.Numerics;
with Ada.Numerics.Long_Elementary_Functions;
use  Ada.Numerics.Long_Elementary_Functions;
with real_Arrays;
use  real_Arrays;
with nuderiv; -- procedure under test
with rderiv;  -- for comparison

procedure nuderiv_test_exp is
  npoint : integer;
  n : Integer := 12;
  X : real_vector(1..n);
  Y : real_vector(1..n);
  Yp : real_vector(1..n);
  Ypp : real_vector(1..n);
  Yppp : real_vector(1..n);
  Ypppp : real_vector(1..n);
  Yppppp : real_vector(1..n);
  C, Cu : real_vector(1..n);
  nuyp, nyp : real; -- any order
  Xmin, Xmax, hx : real; -- for comparison test
  Err1, Err2, Err3, Err4 , Err5 : Real;
  Max_Err1, Max_Err2, Max_Err3, Max_Err4, Max_Err5 : Real := 0.0;
  nx : integer;
  h : Real;

  package Fltio is new Float_Io(Long_Float);
  use Fltio;
  package Intio is new integer_Io(integer);
  use Intio;

  -- test function and check functions
  function yf5(x : real) return real is -- test function 5
  begin
    return exp(x);
  end yf5;

  function ypf5(x : real) return real is -- derivative of test function 5
  begin
    return exp(x);
  end ypf5;

  function yppf5(x : real) return real is -- 2nd derivative of test function 5
  begin
    return exp(x);
  end yppf5;

  function ypppf5(x : real) return real is -- 3rd derivative of test function5
  begin
    return exp(x);
  end ypppf5;

  function yppppf5(x : real) return real is -- 4th derivative of test function5
  begin
    return exp(x);
  end yppppf5;

  function ypppppf5(x : real) return real is -- 5th derivative of test function5
  begin
    return exp(x);
  end ypppppf5;

begin
  put_line("nuderiv_test_exp.adb running ");
  npoint := n;
  H := 0.1;
    for I in 1..N loop
    X(I) := H;
    H := H*1.33; -- increasing spacing
  end loop;

  put_line(" testing  npoint="&integer'Image(npoint));
  for i in 1..npoint loop
    Y(i) := yf5(x(i));
    Yp(i) := ypf5(x(i));
    Ypp(i) := yppf5(x(i));
    Yppp(i) := ypppf5(x(i));
    Ypppp(i) := yppppf5(x(i));
    Yppppp(i) := ypppppf5(x(i));
    put("x(");
    put(I,1);
    put(")=");
    put(X(i),3,3);
    new_line;
    Put("     y=");
    Put(Y(i),3,3);
    Put(", y'=");
    Put(Yp(i),3,3);
    Put(", y''=");
    Put(Ypp(i),3,3);
    Put(", y'''=");
    Put(Yppp(i),3,3);
    Put(", y''''=");
    Put(Ypppp(i),3,3);
    Put(", y'''''=");
    Put(Yppppp(i),3,3);
    new_line;
  end loop;
  for order in 1..npoint-1 loop
    for point in 1..npoint loop
      Put_line("testing order="&integer'image(order)&
               ", npoint="&integer'image(npoint)&
               ", point="&integer'image(point));
      nuderiv(order, npoint, point, X, C);
      for i in 1..npoint loop
        put("nu c(");
        put(i,1);
        put(")=");
        put(C(i),3,3);
        new_line;
      end loop;
      nuyp := 0.0;
      for i in 1..npoint  loop
        nuyp := nuyp + C(i)*Y(i);
      end loop;
      if order=1 then
        put_line("order="&integer'image(order)&
                 ", npoint="&integer'Image(npoint)&
                 ", point="&Integer'Image(point));
        put("    yp=");
        put(yp(point),3,3);
        put(", nuyp=");
        put(nuyp,3,3);
        put(", err1=");
        err1 := yp(point)-nuyp;
        put(err1,3,3);
        new_line;
        if abs(err1) > Max_Err1 then Max_Err1 := abs(err1); end if;
      elsif order=2 then
        put_line("order="&integer'image(order)&
                 ", npoint="&integer'Image(npoint)&
                 ", point="&Integer'Image(point));
        put("    ypp=");
        put(ypp(point),3,3);
        put(", nuyp=");
        put(nuyp,3,3);
        put(", err2=");
        err2 := ypp(point)-nuyp;
        put(err2,3,3);
        new_line;
        if abs(err2) > Max_Err2 then Max_Err2 := abs(err2); end if;
        put(err2,3,3);
        new_line;
      elsif order=3 then
        put_line("order="&integer'image(order)&
                 ", npoint="&integer'Image(npoint)&
                 ", point="&Integer'Image(point));
        put("    yppp=");
        put(yppp(point),3,3);
        put(", nuyp=");
        put(nuyp,3,3);
        put(", err3=");
        err3 := yppp(point)-nuyp;
        put(err3,3,3);
        new_line;
        if abs(err3) > Max_Err3 then Max_Err3 := abs(err3); end if;
        put(err3,3,3);
        new_line;
      elsif order=4 then
        put_line("order="&integer'image(order)&
                 ", npoint="&integer'Image(npoint)&
                 ", point="&Integer'Image(point));
        put("    ypppp=");
        put(ypppp(point),3,3);
        put(", nuyp=");
        put(nuyp,3,3);
        put(", err4=");
        err4 := ypppp(point)-nuyp;
        put(err4,3,3);
        new_line;
        if abs(err4) > Max_Err4 then Max_Err4 := abs(err4); end if;
        put(err4,3,3);
        new_line;
      elsif order=5 then
        put_line("order="&integer'image(order)&
                 ", npoint="&integer'Image(npoint)&
                 ", point="&Integer'Image(point));
        put("    yppppp=");
        put(yppppp(point),3,3);
        put(", nuyp=");
        put(nuyp,3,3);
        put(", err5=");
        err5 := yppppp(point)-nuyp;
        put(err5,3,3);
        new_line;
        if abs(err5) > Max_Err5 then Max_Err5 := abs(err5); end if;
        put(err5,3,3);
        new_line;
      end if;
    end loop;
  end loop;
  New_Line;
  Put_Line("Compare accuracy of nuderiv and nderiv on uniform spaced");
  Xmin := 0.0;
  Xmax := 1.0;
  nx := 5;
  Hx := (Xmax-Xmin)/real(nx-1);
  for i in 1..nx loop
    X(i) := Xmin + Hx*Real(i-1);
    Y(i) := X(i)**4 + 2.0*X(I)**3 + 3.0*X(I)**2 + 4.0*X(I) + 5.0;
    Put_Line("x="&Real'Image(X(i))&", y="&Real'Image(Y(i)));
  end loop;
  for i in 1..nx loop -- point
    nuderiv(4, nx, i, X,  C);
    rderiv( 4, nx, i, Hx, Cu);
    nuyp := 0.0;
    nyp  := 0.0;
    Put_Line("testing order 4, 5 points, at point i="&Integer'Image(i));
    for j in 1..nx loop
      nuyp := nuyp + Cu(j)*Y(j);
      nyp  := nyp  + C(j)*Y(j);
    end loop; -- j
    Put("nu=");
    Put(nuyp,2,3,0);
    Put(", nd=");
    Put(nyp,2,3,0);
    Put(", actual=");
    Put(24.0,2,3,0);
    Put_Line(", diff="&Real'Image(Nuyp-Nyp));
  end loop; -- i
  Put_Line("first  order max_err1="&Real'Image(Max_Err1));
  Put_Line("second order max_err2="&Real'Image(Max_Err2));
  Put_Line("third  order max_err3="&Real'Image(Max_Err3));
  Put_Line("fourth order max_err4="&Real'Image(Max_Err4));
  Put_Line("fifth  order max_err5="&Real'Image(Max_Err5));
end nuderiv_test_exp;