-- nuderiv_test.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 is
  npoint : integer;
  X : real_vector(1..5) := (1.0, 2.0, 4.0, 7.0, 8.5);
  Y : real_vector(1..5);
  Yp : real_vector(1..5);
  Ypp : real_vector(1..5);
  Yppp : real_vector(1..5);
  Ypppp : real_vector(1..5);
  C, Cu : real_vector(1..5);
  nuyp, nyp : real; -- any order
  Xmin, Xmax, hx : real; -- for comparison test
  nx : integer;

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

  -- test functions and check functions
  function yf4(x : real) return real is -- test function 4
  begin
    return 3.0 -0.5* x + 0.2*x*x - 0.1*x*x*x;
  end yf4;

  function ypf4(x : real) return real is -- derivative of test function 4
  begin
    return -0.5 + 0.4*x - 0.3*x*x;
  end ypf4;

  function yppf4(x : real) return real is -- 2nd derivative of test function 4
  begin
    return 0.4 - 0.6*x;
  end yppf4;

  function ypppf4(x : real) return real is -- 3rd derivative of test function4
  begin
    return -0.6;
  end ypppf4;

  function yf5(x : real) return real is -- test function 5
  begin
    return 3.0 -0.5* x + 0.1*x*x - 0.05*x*x*x + 0.002*x*x*x*x;
  end yf5;

  function ypf5(x : real) return real is -- derivative of test function 5
  begin
    return -0.5 + 0.2*x - 0.15*x*x + 0.008*x*x*x;
  end ypf5;

  function yppf5(x : real) return real is -- 2nd derivative of test function 5
  begin
    return 0.2 - 0.3*x + 0.024*x*x;
  end yppf5;

  function ypppf5(x : real) return real is -- 3rd derivative of test function5
  begin
    return -0.3 + 0.048*x;
  end ypppf5;

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

begin
  put_line("nuderiv_test.adb running ");
  npoint := 4;
  put_line(" testing  npoint="&integer'Image(npoint));
  for i in 1..npoint loop
    Y(i) := yf4(X(i));
    Yp(i) := ypf4(X(i));
    Ypp(i) := yppf4(X(i));
    Yppp(i) := ypppf4(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);
    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(", err=");
        put(yp(point)-nuyp,3,3);
        new_line;
      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(", err=");
        put(ypp(point)-nuyp,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(", err=");
        put(yppp(point)-nuyp,3,3);
        new_line;
      end if;
    end loop;
  end loop;

  npoint := 5;
  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));
    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);
    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(", err=");
        put(yp(point)-nuyp,3,3);
        new_line;
      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(", err=");
        put(ypp(point)-nuyp,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(", err=");
        put(yppp(point)-nuyp,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(", err=");
        put(ypppp(point)-nuyp,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
end nuderiv_test;