-- test_lsfit3.adb  using stuff from  pde_nl33.adb  for test

with Ada.Text_IO;
use  Ada.Text_IO;
with Ada.Numerics.Long_Elementary_Functions;
use  Ada.Numerics.Long_Elementary_Functions;
with Real_Arrays; -- types real, real_vector, real_matrix
use  Real_Arrays;
with Array3d;
use  Array3d;
with Integer_Arrays; -- types integer_vector
use  Integer_Arrays;
with Lsfit;
use  Lsfit;

procedure Test_Lsfit3 is
  Nx : Integer := 5;      -- number of data in each dimension
  Ny : Integer := 5;
  Nz : Integer := 5;
  Xg : Real_Vector(1..Nx);    -- x grid, does not need to be uniform
  Yg : Real_Vector(1..Ny);    -- y grid
  Zg : Real_Vector(1..Nz);    -- z grid
  Ug : Real_Matrix3d(1..Nx,1..Ny,1..Nz); -- computed solution at xyz grid,
  Ua : Real_matrix3d(1..nx,1..ny,1..nz); -- known includes boundary
  Hx, Hy, Hz : real;
  xmin : Real := -1.0;
  xmax : Real :=  1.0;
  ymin : Real := -1.0;
  ymax : Real :=  1.0;
  zmin : Real := -1.0;
  zmax : Real :=  1.0;
  X, Y, Z, U, Ue : Real := 0.0;
  Last : Boolean := False;
  Err, Maxerr : Real := 0.0;

  function ub(x:Real; y:Real; z:Real) return real is -- boundary
  begin
    return x*x*x*x+2.0*y*y*y*y+3.0*z*z*z*z+4.0*x*y*z+5.0*x*y+z+
           6.0*x + 7.0*y + 8.0;
  end ub;

  function U1(X:Real; Y:Real; Z:real) return Real is -- data
  begin
    return 1.0 + 2.0*X + 3.0*Y + 4.0*Z;
  end U1;

  package Real_IO is new Float_IO(Real);
  use Real_IO;

begin
  Put_Line("test_lsfit3.adb");

  hx := (xmax-xmin)/real(nx-1);
  Put_Line("nx="&Integer'Image(nx)&", step="&Real'Image(hx));
  for i in 1..nx loop
    xg(i) := xmin+Real(i-1)*hx;
  end loop;
  hy := (ymax-ymin)/real(ny-1);
  Put_Line("ny="&Integer'Image(ny)&", step="&Real'Image(hy));
  for ii in 1..ny loop
    yg(ii) := ymin+Real(ii-1)*hy;
  end loop;
  hz := (zmax-zmin)/real(nz-1);
  Put_Line("nz="&Integer'Image(nz)&", step="& Real'Image(hz));
  for iii in 1..nz loop
    zg(iii) := zmin+Real(iii-1)*hz;
  end loop;

  Fit_Init(1, 3); -- first power, three independent variables

  Last := False;
  for i in 1..nx loop
    X := Xg(i);
    for ii in 1..ny loop
      Y := Yg(Ii);
      for iii in 1..nz loop
        Z := Zg(iii);
        U := U1(X,Y,Z);
        if i=nx and ii=ny and iii=nz then
           Last := True;
        end if;
        Fit3_Data(X, Y, Z, U, Last);
      end loop;
    end loop;
  end loop;
  -- built, now check
  for i in 1..nx loop
    X := Xg(i);
    for ii in 1..ny loop
      Y := Yg(ii);
      for iii in 1..nz loop
        Z := Zg(iii);
        U := U1(X,Y,Z);
        Ue := Eval3_Poly(X, Y, Z);
        Err := abs(U-Ue);
        if Err>Maxerr then Maxerr := Err; end if;
      end loop;
    end loop;
  end loop;
  Put_Line("power 1 test, maxerr="&Real'Image(Maxerr));

  for i in 1..nx loop
    for ii in 1..ny loop
      for iii in 1..nz loop
        Ua(i,ii,iii) := ub(xg(i),yg(ii),zg(iii));
        Put("xg=");
        Put(xg(i),3,5,0);
        Put(", yg=");
        Put(yg(ii),3,5,0);
        Put(", zg=");
        Put(zg(iii),3,5,0);
        Put(", Ua("&Integer'Image(i)&","&Integer'Image(ii)&","&
            ","&Integer'Image(iii)&")=");
        Put(Ua(i,ii,iii),3,5,0);
        New_Line;
        if i=1 or i=nx or ii=1 or ii=ny or iii=1 or iii=nz then
          ug(i,ii,iii) := ub(xg(i),yg(ii),zg(iii));
        else
          ug(i,ii,iii) := 0.0;
        end if;
      end loop;
    end loop;
  end loop;
  New_Line;

  Fit_Init(4,3); -- max power, dimension

  Fit3_Bpoly(Nx, Ny, Nz, Xg, Yg, Zg, Ug);

  -- test fit
  for i in 2..nx-1 loop
     X := Xg(i);
     for ii in 2..ny-1 loop
      Y := Yg(ii);
      for iii in 2..nz-1 loop
        Z := Zg(iii);
        U := Eval3_Poly(X, Y, Z);
        Put("expected=");
        Put(Ua(i,ii,iii),8,3,0);
        Put("  fit=");
        Put(U,8,3,0);
        Put("  error="&Real'Image(Ua(i,ii,iii)-U));
        New_Line;
      end loop;
    end loop;
  end loop;
  New_Line;
  Put_Line("test_lsfit3.adb finished");
end Test_Lsfit3;