-- check_test_4d.adb  check every derivative function up to 4th order
--                    from package test_4d, 69 tests
-- this computes all 15 cases of possible derivatives for 4 independent
-- variables up to total 4th derivative, given Xorder, Yorder, Zorder, Torder:
-- Four cases: One partial Uxxxx, Uxxx, Uxx, Ux,  also for y, z, and t      16
-- Six cases:  Two partials xy, xz, xt, yz, yt, zt to order 1,1 1,2 1,3
--                                                          2,1 2,2 3,1     36
-- Four cases: Three partials xyz, xyt, xzt, yzt to order 1,1,1  2,1,1
--                                                        1,2,1  1,1,2      16
-- One case:   Four partials  Uxyzt                                          1
-- see pde44_eq.adb for eval_derivative that also handles boundary values
with Ada.Text_IO;
use  Ada.Text_IO;
with real_arrays;
use  real_arrays;
with test_4d;
use  test_4d;
with deriv;
with rderiv;

procedure check_test_4d is

  xg, yg, zg, tg : real_vector(1..10);
  hx, hy, hz, ht : real;
  nx, ny, nz, nt : integer;
  maxerror, gmaxerror : real;
  xmin, xmax, ymin, ymax, zmin, zmax, tmin, tmax : real;
  type funct is access function (x : real; y : real; z : real; t :real)
                                return real;

  procedure check(name : string;
                  F : funct;
                  points : integer; xord : integer;
                  yord : integer; zord : integer; tord : integer) is
    cx, cy, cz, ct : real_vector(1..10);
    p1, p2, p3 : real_vector(1..10);
    analytic, discrete : real;
    error : real;
    x, y, z, t : real;

  begin
    maxerror := 0.0;
    for i in 1..points loop
      rderiv(xord, points, i, hx, cx);
      x := xg(i);
      for ii in 1..points loop
        rderiv(yord, points, ii, hy, cy);
        y := yg(ii);
        for iii in 1..points loop
          rderiv(zord, points, iii, hz, cz);
          z := zg(iii);
          for iiii in 1..points loop
            rderiv(tord, points, iiii, ht, ct);
            t := tg(iiii);
            analytic := f(x, y, z, t);
            discrete := 0.0;

            -- four cases of one partial x, y, z, t to any order
            if xord/=0  and  yord=0  and  zord=0  and  tord=0 then
              for j in 1..points loop
                discrete := discrete + cx(j)*u(xg(j),y,z,t);
              end loop; -- j

            elsif xord=0  and  yord/=0  and  zord=0  and  tord=0 then
              for jj in 1..points loop
                discrete := discrete + cy(jj)*u(x,yg(jj),z,t);
              end loop; -- jj

            elsif xord=0  and  yord=0  and  zord/=0  and  tord=0 then
              for jjj in 1..points loop
                discrete := discrete + cz(jjj)*u(x,y,zg(jjj),t);
              end loop; -- jjj

            elsif xord=0  and  yord=0  and  zord=0  and  tord/=0 then
              for jjjj in 1..points loop
                discrete := discrete + ct(jjjj)*u(x,y,z,tg(jjjj));
              end loop; -- jjjj

            -- six cases of two partials xy, xz, xt, yz, yt, zt to any order
            elsif xord/=0  and  yord/=0  and  zord=0  and  tord=0 then
              for j in 1..points loop
                p1(j) := 0.0;
                for jj in 1..points loop
                   p1(j) := p1(j) + cy(jj)*u(xg(j),yg(jj),z,t);
                end loop; -- jj
                discrete := discrete + cx(j)*p1(j);
              end loop; -- j

            elsif xord/=0  and  yord=0  and  zord/=0  and  tord=0 then
              for j in 1..points loop
                p1(j) := 0.0;
                for jjj in 1..points loop
                  p1(j) := p1(j) + cz(jjj)*u(xg(j),y,zg(jjj),t);
                end loop; -- jjj
                discrete := discrete + cx(j)*p1(j);
              end loop; -- j

            elsif xord/=0  and  yord=0  and  zord=0  and  tord/=0 then
            for j in 1..points loop
              p1(j) := 0.0;
              for jjjj in 1..points loop
                 p1(j) := p1(j) + ct(jjjj)*u(xg(j),y,z,tg(jjjj));
              end loop; -- jjjj
              discrete := discrete + cx(j)*p1(j);
            end loop; -- j

            elsif xord=0  and  yord/=0  and  zord/=0  and  tord=0 then
              for jj in 1..points loop
                p1(jj) := 0.0;
                for jjj in 1..points loop
                  p1(jj) := p1(jj) + cz(jjj)*u(x,yg(jj),zg(jjj),t);
                end loop; -- jjj
                discrete := discrete + cy(jj)*p1(jj);
              end loop; -- jj

            elsif xord=0  and  yord/=0  and  zord=0  and  tord/=0 then
              for jj in 1..points loop
                p1(jj) := 0.0;
                for jjjj in 1..points loop
                   p1(jj) := p1(jj) + ct(jjjj)*u(x,yg(jj),z,tg(jjjj));
                end loop; -- jjjj
                discrete := discrete + cy(jj)*p1(jj);
              end loop; -- jj

            elsif xord=0  and  yord=0  and  zord/=0  and  tord/=0 then
              for jjj in 1..points loop
                p1(jjj) := 0.0;
                for jjjj in 1..points loop
                   p1(jjj) := p1(jjj) + ct(jjjj)*u(x,y,zg(jjj),tg(jjjj));
                end loop; -- jjjj
                discrete := discrete + cz(jjj)*p1(jjj);
              end loop; -- jjj

            -- four cases of three partials xyz, xyt, xzt, yzt to any order
            elsif xord/=0  and  yord/=0  and  zord/=0  and  tord=0 then
              for j in 1..points loop
                p1(j) := 0.0;
                for jj in 1..points loop
                  p2(jj) := 0.0;
                  for jjj in 1..points loop
                    p2(jj) := p2(jj) + cz(jjj)*u(xg(j),yg(jj),zg(jjj),t);
                  end loop; -- jjj
                  p1(j) := p1(j) + cy(jj)*p2(jj);
                end loop; -- jj
                discrete := discrete + cx(j)*p1(j);
              end loop; -- j

            elsif xord/=0  and  yord/=0  and  zord=0  and  tord/=0 then
              for j in 1..points loop
                p1(j) := 0.0;
                for jj in 1..points loop
                  p2(jj) := 0.0;
                  for jjjj in 1..points loop
                    p2(jj) := p2(jj) + ct(jjjj)*u(xg(j),yg(jj),z,tg(jjjj));
                  end loop; -- jjjj
                  p1(j) := p1(j) + cy(jj)*p2(jj);
                end loop; -- jj
                discrete := discrete + cx(j)*p1(j);
              end loop; -- j

            elsif xord/=0  and  yord=0  and  zord/=0  and  tord/=0 then
              for j in 1..points loop
                p1(j) := 0.0;
                for jjj in 1..points loop
                  p2(jjj) := 0.0;
                  for jjjj in 1..points loop
                    p2(jjj) := p2(jjj) + ct(jjjj)*u(xg(j),y,zg(jjj),tg(jjjj));
                  end loop; -- jjjj
                  p1(j) := p1(j) + cz(jjj)*p2(jjj);
                end loop; -- jjj
                discrete := Discrete + cx(j)*p1(j);
              end loop; -- j

            elsif xord=0  and  yord/=0  and  zord/=0  and  tord/=0 then
              for jj in 1..points loop
                p1(jj) := 0.0;
                for jjj in 1..points loop
                  p2(jjj) := 0.0;
                  for jjjj in 1..points loop
                    p2(jjj) := p2(jjj) + ct(jjjj)*u(x,yg(jj),zg(jjj),tg(jjjj));
                  end loop; -- jjjj
                  p1(jj) := p1(jj) + cz(jjj)*p2(jjj);
                end loop; -- jjj
                discrete := discrete + cy(jj)*p1(jj);
              end loop; -- jj

            -- one case of four partials
            else
              for j in 1..points loop
                p1(j) := 0.0;
                for jj in 1..points loop
                  p2(jj) := 0.0;
                  for jjj in 1..points loop
                    p3(jjj) := 0.0;
                    for jjjj in 1..points loop
                      p3(jjj) := p3(jjj) +ct(jjjj)*u(xg(j),yg(jj),zg(jjj),tg(jjjj));
                    end loop; -- jjjj
                    p2(jj) := p2(jj) + cz(jjj)*p3(jjj);
                    end loop; -- jjj
                  p1(j) := p1(j) + cy(jj)*p2(jj);
                  end loop; -- jj
                discrete := Discrete + cx(j)*p1(j);
                end loop; -- j
              end if;
              error := analytic-discrete;
              if abs(Error)>maxerror then
                Maxerror := abs(error);
              end if;
              if maxerror>0.001 and i=0 and ii=0 and iii=0 and iiii=0 then
                put_line("analytic="&real'image(Analytic)&
                         ", discrete="&real'image(Discrete)&
                         ", error="&real'image(error));
              end if;
            end loop; -- iiii for t
          end loop; -- iii for z
        end loop; -- ii for y
      end loop; -- i for x
      put_line("                                 "&name&" maxerror="&
               real'image(maxerror));
      if maxerror>gmaxerror then
        gmaxerror := maxerror;
      end if;
  end Check;

begin

  put_line("check_test_4d.adb running");
  -- a test case, all do not need to be same
  xmin := -1.0;
  xmax := 1.0;
  nx := 5;
  hx := (xmax-xmin)/real(nx-1);
  for i in 1..nx loop
    xg(i) := xmin + real(i-1)*hx;
  end loop;
  ymin := -1.0;
  ymax := 1.0;
  ny := 5;
  hy := (ymax-ymin)/real(ny-1);
  for ii in 1..ny loop
    yg(ii) := ymin + real(ii-1)*hy;
  end loop;
  zmin := -1.0;
  zmax := 1.0;
  nz := 5;
  hz := (zmax-zmin)/real(nz-1);
  for iii in 1..nz loop
    zg(iii) := zmin + real(iii-1)*hz;
  end loop;
  tmin := -1.0;
  tmax := 1.0;
  nt := 5;
  ht := (tmax-tmin)/real(nt-1);
  for iiii in 1..nt loop
    tg(iiii) := tmin + real(iiii-1)*ht;
  end loop;
  gmaxerror := 0.0;
  -- check all 69 cases of 4 independent variables for total order 1 through 4
  check("Ux   ", Ux'access,    5, 1, 0, 0, 0);
  check("Uy   ", Uy'access,    5, 0, 1, 0, 0);
  check("Uz   ", Uz'access,    5, 0, 0, 1, 0);
  check("Ut   ", Ut'access,    5, 0, 0, 0, 1);
  check("Uxx  ", Uxx'access,   5, 2, 0, 0, 0);
  check("Uxy  ", Uxy'access,   5, 1, 1, 0, 0);
  check("Uxz  ", Uxz'access,   5, 1, 0, 1, 0);
  check("Uxt  ", Uxt'access,   5, 1, 0, 0, 1);
  check("Uyy  ", Uyy'access,   5, 0, 2, 0, 0);
  check("Uyz  ", Uyz'access,   5, 0, 1, 1, 0);
  check("Uyt  ", Uyt'access,   5, 0, 1, 0, 1);
  check("Uzz  ", Uzz'access,   5, 0, 0, 2, 0);
  check("Uzt  ", Uzt'access,   5, 0, 0, 1, 1);
  check("Utt  ", Utt'access,   5, 0, 0, 0, 2);
  check("Uxxx ", Uxxx'access,  5, 3, 0, 0, 0);
  check("Uxxy ", Uxxy'access,  5, 2, 1, 0, 0);
  check("Uxxz ", Uxxz'access,  5, 2, 0, 1, 0);
  check("Uxxt ", Uxxt'access,  5, 2, 0, 0, 1);
  check("Uxyy ", Uxyy'access,  5, 1, 2, 0, 0);
  check("Uxyz ", Uxyz'access,  5, 1, 1, 1, 0);
  check("Uxyt ", Uxyt'access,  5, 1, 1, 0, 1);
  check("Uxzz ", Uxzz'access,  5, 1, 0, 2, 0);
  check("Uxzt ", Uxzt'access,  5, 1, 0, 1, 1);
  check("Uxtt ", Uxtt'access,  5, 1, 0, 0, 2);
  check("Uyyy ", Uyyy'access,  5, 0, 3, 0, 0);
  check("Uyyz ", Uyyz'access,  5, 0, 2, 1, 0);
  check("Uyyt ", Uyyt'access,  5, 0, 2, 0, 1);
  check("Uyzz ", Uyzz'access,  5, 0, 1, 2, 0);
  check("Uyzt ", Uyzt'access,  5, 0, 1, 1, 1);
  check("Uytt ", Uytt'access,  5, 0, 1, 0, 2);
  check("Uzzz ", Uzzz'access,  5, 0, 0, 3, 0);
  check("Uzzt ", Uzzt'access,  5, 0, 0, 2, 1);
  check("Uztt ", Uztt'access,  5, 0, 0, 1, 2);
  check("Uttt ", Uttt'access,  5, 0, 0, 0, 3);

  check("Uxxxx", Uxxxx'access, 5, 4, 0, 0, 0);
  check("Uxxxy", Uxxxy'access, 5, 3, 1, 0, 0);
  check("Uxxxz", Uxxxz'access, 5, 3, 0, 1, 0);
  check("Uxxxt", Uxxxt'access, 5, 3, 0, 0, 1);
  check("Uxxyy", Uxxyy'access, 5, 2, 2, 0, 0);
  check("Uxxyz", Uxxyz'access, 5, 2, 1, 1, 0);
  check("Uxxyt", Uxxyt'access, 5, 2, 1, 0, 1);
  check("Uxxzz", Uxxzz'access, 5, 2, 0, 2, 0);
  check("Uxxzt", Uxxzt'access, 5, 2, 0, 1, 1);
  check("Uxxtt", Uxxtt'access, 5, 2, 0, 0, 2);
  check("Uxyyy", Uxyyy'access, 5, 1, 3, 0, 0);
  check("Uxyyz", Uxyyz'access, 5, 1, 2, 1, 0);
  check("Uxyyt", Uxyyt'access, 5, 1, 2, 0, 1);
  check("Uxyzz", Uxyzz'access, 5, 1, 1, 2, 0);
  check("Uxyzt", Uxyzt'access, 5, 1, 1, 1, 1);
  check("Uxytt", Uxytt'access, 5, 1, 1, 0, 2);
  check("Uxzzz", Uxzzz'access, 5, 1, 0, 3, 0);
  check("Uxzzt", Uxzzt'access, 5, 1, 0, 2, 1);
  check("Uxztt", Uxztt'access, 5, 1, 0, 1, 2);
  check("Uxttt", Uxttt'access, 5, 1, 0, 0, 3);
  check("Uyyyy", Uyyyy'access, 5, 0, 4, 0, 0);
  check("Uyyyz", Uyyyz'access, 5, 0, 3, 1, 0);
  check("Uyyyt", Uyyyt'access, 5, 0, 3, 0, 1);
  check("Uyyzz", Uyyzz'access, 5, 0, 2, 2, 0);
  check("Uyyzt", Uyyzt'access, 5, 0, 2, 1, 1);
  check("Uyytt", Uyytt'access, 5, 0, 2, 0, 2);
  check("Uyzzz", Uyzzz'access, 5, 0, 1, 3, 0);
  check("Uyzzt", Uyzzt'access, 5, 0, 1, 2, 1);
  check("Uyztt", Uyztt'access, 5, 0, 1, 1, 2);
  check("Uyttt", Uyttt'access, 5, 0, 1, 0, 3);
  check("Uzzzz", Uzzzz'access, 5, 0, 0, 4, 0);
  check("Uzzzt", Uzzzt'access, 5, 0, 0, 3, 1);
  check("Uzztt", Uzztt'access, 5, 0, 0, 2, 2);
  check("Uzttt", Uzttt'access, 5, 0, 0, 1, 3);
  check("Utttt", Utttt'access, 5, 0, 0, 0, 4);
  put_line("check_test_4d.adb finished, max maxerror="&
           real'image(gmaxerror));

end check_test_4d;